Revisiting nonfactorizable contributions to the factorization-forbidden decays of $B$ meson to charmonium

Motivated by the observed large rates of $B\rightarrow (\chi_{c0}, \chi_{c2}, h_c)K$ decays by $BABAR$ and Belle collaborations, we investigate the nonfactorizable contributions to these factorization-forbidden decays, which can occur via a gluon exchange between the $c\bar c$ system and the spectator quark. Our numerical results demonstrate that the spectator contributions are able to produce a large branching ratio consistent with experiments. As a by-product, we also study the Cabibbo-suppressed decays, such as $B\rightarrow (\chi_{c0}, \chi_{c2}, h_c)\pi$ and the U-spin related $B_s$ decay, which receive less attention from both theory and experiment. The calculated branching ratios reach the order of $10^{-6}$, which are in the scope of Belle-II and LHCb experiments. Besides, the $CP$-asymmetry parameters are also calculated for these decays. The obtained results are compared with available experimental data and numbers from other predictions. We also investigate the sources of theoretical uncertainties in our calculation.


I. INTRODUCTION
In the Standard Model (SM) of particle physics, the charmonium decays of B meson arise from the quark-level process b → qcc with q = d, s involving tree and penguin amplitudes. If the charmonium states are produced from a cc system with the orbital angular momentum L = 0, such as η c and J/ψ. These decays belong to the color-suppressed category and receive large nonfactorizable corrections. For the orbital excitation of the cc assignments with L = 1, the spin and orbit interaction between the charm-anticharm quarks pair can create four P -wave charmonium states, namely, χ c0 , χ c1 , χ c2 , and h c . Except for the χ c1 modes, which are allowed under the factorization hypothesis, others are prohibited due to the V − A structure of the weak vertex [1]. However, these processes could occur via an gluon exchange between the charmonium system and the quark in other mesons, which induces the so-called nonfactorizable contributions. Therefore, the factorization-forbidden decays of B meson to a P -wave charmonium state can provide valuable insights into the nonfactorizable mechanism.
Experimentally, the first observation of the decay B + → χ c0 K + was reported by the Belle Collaboration [2] using χ c0 decays to pion or kaon pair, and later confirmed by the BABAR Collaboration [3]. Subsequently, both the BABAR [4][5][6][7] and Belle Collaborations [8,9] have performed the improved measurements. The current world averages of the absolute branching ratio B → χ c0 K can reach the order of 10 −4 [10], which is of the same order in magnitude as that of the factorization-allowed χ c1 mode. The corresponding vector K * mode have been searched and observed by BABAR Collaboration [11,12] with a similar rate. Another factorization-inhibited decay B + → χ c2 K + has been measured by the Belle [13,14] and BABAR [15] with an average branching ratio of (1.1 ± 0.4) × 10 −5 [10], which is an order of magnitude smaller than that of χ c0 mode. Several Collaborations have also searched for the process B → h c K [16][17][18]. The current branching ratio upper limit on its branching ratio is 3.8 × 10 −5 at the 90% confidence level [10], obtained in the h c search by Belle [16]. Very recently, the Belle Collaboration performed a search for the decays B + → h c K + and B 0 → h c K 0 S [19]. They found the evidence for the former with a 4.8σ significance, and set the upper limit on the latter. These results are comparable to B → χ c2 K but below the measured rate for B → χ c0 K. Most recently, both the Belle [20] and BABAR [21] present the measurement of the absolute branching fractions of B + → X cc K + , where X cc refers to one charmonium state.
The measured branching ratios in these factorization-inhibited decays are surprisingly larger than the expectations from factorization, which suggests that nonfactorizable contributions in B decays to charmonium can be sizable. Many attempts have been made to resolve above puzzles. Melić analyze the soft nonfactorizable effects to B → (η c , J/ψ, χ c0,c1 )K decays [1] by using the light-cone sum rules (LCSR) approach. The calculated branching ratio of the B → χ c0 K, however, are too small to accommodate the data due to the larger cancellation between the twist-3 and FIG. 1: The nonfactorizable spectator diagrams for the decays B → χcP , where P stands for a final state hadron pion or kaon.
twist-4 pieces in the nonfactorizable contributions. A similar conclusion is also drawn in [22] with the same method. The exclusive B decays to P -wave charmonium states have been studied in the framework of QCD factorization (QCDF) [23][24][25][26][27][28]. They found the soft contributions may be large since there exist infrared divergences arising from the vertex corrections as well as endpoint singularities due to the leading twist spectator corrections. However, M. Beneke pointed out that the concerned decays may occur via the color octet mechanism [29], and the endpoint divergence in hard spectator-scattering factorizes can be absorbed into color octet operator matrix elements [30]. The B → χ c0 K and B → h c K decays [31,32] have been investigated in the perturbative QCD (PQCD) approach based on the k T factorization theorem [33][34][35][36], in which the vertex corrections are ignored and the endpoint singularity is cured by including the parton's transverse momentum. Without endpoint singularity, the PQCD approach so far, has been successfully applied to the studies of various factorization-allowed charmonium decays of B meson [37][38][39][40][41][42][43][44][45][46][47][48][49]. Most predictions are consistent with the current experiments. In another approach [50,51], by considering intermediate charmed meson rescattering effects, the authors claimed that the rescattering effects could provide the large part of the B → χ c0 K, h c K amplitudes. Generally, the nonfactorizable dynamics include both the vertex corrections and spectator amplitudes. We agree on the comment made in [31]: the PQCD formalism for the vertex corrections require the transverse momenta dependent charmonium meson wave functions as the necessary nonperturbative inputs, which are not yet available completely. The authors in [26] confirm that the vertex corrections are numerically small, and the spectator corrections are large and dominant in the QCDF. Furthermore, the previous PQCD calculations [31,32] demonstrate the spectator contributions are sufficient to account for the measurements. Motivated by the same idea, we have reason to believe that it is appropriate to analyze other factorization-inhibited decays, such as B → χ c2 K, in the PQCD approach. Moreover, we update the P -wave charmonium distribution amplitudes (DAs) according to our previous work [52], where the new universal nonperturbative objects are successful in describing the P -wave charmonium decays of B/B c meson [52][53][54].
The layout of the paper is as follows: In Sec. II, we give our theoretical formulae based on the PQCD framework. The input parameters together with the numerical results and discussions are presented in Sec. III. Finally, we give our summary in Sec. IV. The relevant meson distribution amplitudes are relegated to the Appendix.

II. ANALYTIC FORMULAS
For the factorization-inhibited B decays in question, only the nonfactorization spectator diagrams contribute, which are displayed in Fig 1. In the rest frame of the B meson, we define P i and k i with i = 1, 2, 3 to be the four-momenta of the mesons and quarks in the initial and final states, as indicated in Fig 1. In the light-cone coordinates, they are parametrized as [54] where r = m/M represents the ratio of charmonium mass m to B meson mass M . We neglect the light pseudoscalar meson mass for simplicity. The x i and k iT represent the parton longitudinal momentum fractions and the parton transverse momentua, respectively. Similar to the vector charmonium, the longitudinal polarization vectors ǫ L of an axial-vector charmonium can be defined as which satisfy the normalization ǫ 2 L = −1 and the orthogonality ǫ L · P 2 = 0. For a tensor charmonium, the polarization tensor ǫ µν (λ) with helicity λ can be constructed via the polarization vector ǫ µ [55][56][57], whose detail expressions can be found in Refs [52,53].
The Hamiltonian referred to this work in the SM is written as [58]: where G F is the Fermi coupling constant, ξ c(t) = V * c(t)b V c(t)q is the production of Cabibbo-Kobayashi-Maskawa (CKM) matrix element. O i (µ) and C i (µ) are the effective local four quark operators and their QCD corrected Wilson coefficients at the renormalization scale µ, respectively. Their explicit expressions can be found in Ref. [58].
In the PQCD approach, the decay amplitudes are expressed as the convolution of the hard kernel H with the relevant meson light-cone wave functions Φ i where "Tr" denotes the trace over all dirac structures and color indices. t is the energy scale in hard function H. The meson wave functions Φ absorb the nonperturbative dynamics in the hadronization processes. The definitions of the Pwave charmonium wave functions refer to our previous work [52], while those for the B and light pseudoscalar mesons, one can consult Ref. [59]. We list the relevant meson distribution amplitudes and the corresponding parameters in the Appendix. As mentioned before, since only the nonfactorizable spectator diagrams contribute to the hard kernel H, the momentum convolution integral in Eq. (4) would involve all the initial and final states altogether. The perturbative calculations can be performed without endpoint singularity. In the following, we start to compute the decay amplitudes of the concerned decays. We mark subscript S, A, and T to denote a scalar, axial-vector, and tensor charmonium in the final states, respectively. The decay amplitudes for (V − A) ⊗ (V − A) operators read as with the color factor C f = 4/3, the mass ratios r c = m c /M , r p = m 0 /M , m c being the charm quark mass, while m 0 being the chiral scale parameter. For simplicity, in the above formulas, we have combined the two decay amplitudes in Fig. 1(a) and 1(b) because their hard kernel are symmetric under x 2 ↔ 1 − x 2 , and dropped those power-suppressed terms higher than r 2 orders. α and β are the virtuality of the internal gluon and quark, respectively, whose expressions are given as The hard scale t is chosen as the largest scale of the virtualities of the internal particles in the hard amplitudes, The functions h and the Sudakov factors S cd (t) are referred to Appendix A of Ref. [44] for details. For the decay amplitudes from the (V − A) ⊗ (V + A) operators, we need to do Fierz transformation to get the right color structure for factorization to work, they satisfy, By combining the contributions from different diagrams with corresponding Wilson coefficients, one obtains the total decay amplitudes as and the CP -averaged branching ratio is then expressed as whereĀ denote the corresponding charge conjugate decay amplitude, which can be obtained by conjugating the CKM elements in A.

III. NUMERICAL RESULTS
To estimate the contributions from the decay amplitudes, we need to specify various parameters entering our numerical analysis throughout this paper. We employ the meson lifetime τ B + = 1.638 ps, τ Bs = 1.51 ps, and τ B 0 = 1.51 ps [10]. The parameters relevant for the Wolfenstein parameters are λ = 0.22506, A = 0.811,ρ = 0.124, η = 0.356 [10]. We take M B = 5.28 GeV, M Bs = 5.37 GeV, m χc0 = 3.415 GeV, m χc2 = 3.556 GeV, m hc = 3.525 GeV for the meson masses [10] andm b (m b ) = 4.18 GeV,m c (m c ) = 1.275 GeV for the b and c quark "running" masses in the M S scheme. The chiral masses relates the pseudoscalar meson mass to the quark mass are set as m 0 = 1.6 ± 0.2 GeV [60]. The decay constants (GeV) can be extracted from other decay rates or evaluated from the QCD sum rules, which are summarized below [10,59,61] Using above input parameters and employing the analytic formulas in Sec. II, we derive the CP -averaged branching ratios with error bars, where the second equal-sign in each row denote the central value with all uncertainties added in quadrature. The theoretical errors correspond to the uncertainties due to (1) the shape parameters: ω b = 0.40 ± 0.04 GeV for the Decay amplitudes twist-2 twist-3 Total  ) and (b) for B + → (χc0, χc2, hc)K + decays. The results are given in units of 10 −3 GeV 3 .

Decay amplitudes
Aa  7), which we vary from 0.75t to 1.25t, and Λ (5) QCD = 0.25 ± 0.05 GeV. In general, The uncertainties induced by these parameter are comparable. The uncertainties stemming from the decay constants of charmonium states are not shown explicitly, whose effects on the branching ratios via the relation of B ∝ f 2 χc . We have checked the sensitivity of the branching ratios to the charm quark velocity v inside the charmonium DAs in Eq. (A4). The variation of v 2 in the range 0.25 to 0.35 leads to the difference of our results does not exceed 10%, which suggests that the relativistic corrections may be not important for these decays.
Moreover, because isospin is conserved in the heavy quark limit, we can obtain the branching ratios of neutral counterpart by multiplying the charged ones by the lifetime ratio τ B 0 /τ B + , From Eqs. (5)- (7), it can be seen that each decay amplitude receives contributions from both twist-2 and twist-3 DAs of the charmonium state, whose results are displayed separately in Table I, with all the input parameters taken at their central values. We simply symbolize them as "twist-2" and "twist-3", respectively, while the label "total" corresponds to the total contributions. The numerical results indicate that the dominant contribution comes from the twist-2 DA. This can be understood from the formulas in Eqs. (5)-(7), the contribution of the twist-3 DA is power-suppressed compared with the twist-2 one. Since we used the same asymptotic model of the twist-2 DA for the three charmonium states [see Eq. (A3)], these decay amplitudes are governed by their different decay constants. The relation among the decay constants f χc0 > f χc2 ∼ f hc in Eq. (13) roughly implies the hierarchy pattern as shown in Table I. For the twist-3 piece, the different asymptotic behaviors and tensor decay constants contribute to different values.
As discussed in Refs [62,63], for the non-leptonic two-body decay of B meson, if the emitted particle from the weak vertex is a pseudoscalar or vector meson, there is a destructive interference between the two nonfactorizable spectator diagrams in Fig.1 due to the symmetric twist-2 DAs of the emitted meson. However, for the case of this work, all the twist-2 DAs of the emitted charmonium state are antisymmetric under the exchange of the momentum fraction x 2 of the c quark and 1 − x 2 of thec quark [see Eq. (A3)], which reverses the constructive or destructive interference situation. To be more explicit, we present the values of two spectator amplitudes in Table II, where the A a and A b denote the decay amplitudes from Fig. 1(a) and Fig. 1(b), respectively. It can be seen that the constructive interference between A a and A b can enhance the total decay amplitudes. This is the main reason that the large decay rates for these factorization-forbidden modes are comparable to those naively factorizable decays.

Modes
This the infrared divergences arising from the vertex diagrams and endpoint divergences due to spectator amplitudes at the leading-twist order. The reason for their different numerical results mainly due to the treatment of these divergences.
In [26], the infrared and endpoint divergences are regularized by a non-zero gluon mass and the off-shellness of quarks, respectively. They find the B → χ c0 K is dominated by the spectator contribution, while the contributions arising from the vertex are numerically small. Subsequently, the same scheme is applying to the B → h c (χ c2 )K decays [27] except for neglecting the vertex corrections. In Ref. [30], the authors use space-time dimension as infrared regulator, while the endpoint divergence is absorbed into color octet operator matrix elements. The scheme in [25] is that a non-zero binding energy is introduced to regularize the infrared divergence, and the the spectator contributions are parameterized in a model-dependent way. In general, our results are more consistent with the QCDF predictions from [30] within the low charm quark mass region. Compared with previous PQCD calculations [31,32], we update the charmonium distribution amplitudes and some input parameters in this work. Our predictions on B(B → χ 0 K) are typically smaller than those of [31] and more close to the current experimental average [10]. For the h c mode, both the current PQCD result and the previous calculations in [27,32] are well consistent with the present data. Two earlier papers [50,51] studied the decays of B + → χ c0 K + and B + → h c K + by including the rescattering effects. The evaluated branching ratios are lie in the range of (1.1 ∼ 3.5) × 10 −4 and (2 ∼ 12) × 10 −4 , respectively. Although the rescattering effects can enhance the branching ratio of B + → χ c0 K + to match the data, the value of B(B + → h c K + ) may be overestimated due to uncontrollable theoretical uncertainties. From Table III, one finds that the B + → χ c2 K + channel is predicted to have 3 times larger branching ratio than the data, which is dominated by the Belle experiment [14]. BABAR gave the upper bounds (the values) [15], where the uncertainties are statistical and systematic, respectively. It seems that their central values are somewhat different, and suffer from sizeable statistical uncertainties. Our branching ratio for the neutral mode is to be comparable with the upper limit of BABAR [15]. It is worth noting that the predictions from QCDF on this mode have a relative big spread. For example, the QCDF prediction from [27] gave B(B + → χ c2 K + ) = 3.0 × 10 −6 , which is too small compared to the measured value. However, another QCDF prediction, of order 10 −4 [25], is too large. As pointed out in [25], the number can be adjusted to the right magnitude with an appropriate choice for the parameters of the spectator hard scattering contributions. The large numerical difference between two QCDF calculations is mainly due to the large twist-3 spectator contribution, which can be traced to the infrared behavior of the spectator interactions.
As mentioned above, the charged and neutral decay modes differ in the lifetimes of B 0 and B + in our formalism, the predicted branching ratios are almost in the same magnitude. However, the data of B(B 0 → χ c0 K 0 ) in Table III is two order of magnitude smaller than that of the charged one. This number was obtained by LHCb from the Dalitz plot (DP) analysis of the B 0 → K 0 S π + π − decays [64]. It shows a tension with the similar amplitude analysis performed by BABAR [7,65,66]. One can see that all of the model calculations in Table III are substantially larger than the LHCb data in magnitude. Improved measurements are certainly needed for this decay mode.
The CP asymmetry arises from the interference between tree and penguin amplitudes. Using the same definition in Ref. [54], we study direct CP violation A dir CP and mixing-induced CP asymmetry S for the decays in question. The two CP -violating parameters can be expressed as with λ f = η f e −2iβ (s)Ā A , and η f = 1(−1) for CP -even (odd) final states. β (s) is one of the angles of the unitarity triangle in the CKM quark mixing matrix [10]. Our numerical results are tabulated in Table IV, where the errors induced by the same sources as in Eq. (14). Since the hadronic parameter dependence canceled out in Eq. (18), these CP asymmetries are more sensitive to the heavy quark masses and the hard scale. It is obvious that direct CP violations for decays involving K are very small, of order 10 −4 , due to an almost null weak phase from the CKM matrix element V ts . For other channels, the weak phase from V td will enhance A dir CP substantially to the level of 10 −2 . Because the penguin contribution is small compared to the tree one, the mixing-induced CP asymmetry S in Eq. (18) is approximately given by −η f S ≈ sin(2β). It is clear from Table IV that the predicted S are not much deviation from the current world average values sin(2β) = 0.695 ± 0.019 and 2β s = (2.55 ± 1.15) × 10 −2 [10]. With more experimental data in the future, these modes can serve as an alternative places to extract CKM phases β (s) .
Experimentally, the world average of A CP (B + → χ c0 K + ) = −0.20 ± 0.18 due to the measurements −0.14 ± 0.15 +0.03 −0.06 [5], −0.96 ± 0.37 ± 0.04 [6] from BABAR and −0.065 ± 0.20 +0.035 −0.024 [9] from Belle. Note that the number from [6] has a large and non-Gaussian uncertainty and its difference from zero is not statistically significant. All measured direct CP violation are consistent with zero. Since LHCb has measured A CP to the accuracy of 10 −3 , it is conceivable that an observation of CP violation in other decays will be feasible in the near future. The mixing-induced CP asymmetry for B 0 → χ c0 K 0 S decay are measured by BABAR [65] with two solutions: where the last uncertainty represents the DP signal model dependence. We can find that the first experimental solution is more favored by our calculation.

IV. CONCLUSION
The factorization-forbidden decays of B meson to charmonium have been revisited in the PQCD formalism, which is free from the endpoint singularities. The charmonium distribution amplitudes are updated according to our previous study. We find the dominant contribution comes from leading-twist charmonium distribution amplitudes. The constructive interference between the two nonfactorizable spectator diagrams enhance the decay rate, which can be compatible with those factorization-allowed decays. The obtained branching ratios of the B + → χ c0 K + and B + → h c K + decays are essentially in agreement with the current data, whereas estimates of the χ c2 one is found to be larger, typically by a factor 3. For those not yet measured decays involving π orK in the final state, The calculated branching ratios will be further tested by the LHCb and Belle-II experiments in the near future. We further estimate the CP -violating parameters. As expected, the direct CP asymmetries in these channels are very small due to the suppressed penguin contributions. The mixing-induced CP asymmetry are close to sin(β (s) ), which suggests that these channels can give a cross-check on the measurement of the CKM phases β and β s .
We also have collected other theoretical results whenever available in Table III, and made a detailed comparison. The predicted branching ratios for most decay processes have similar magnitudes while they can differ by several factors for specific decay modes. In general, our predictions are more consistent with data than these earlier analyses.
We discussed the theoretical uncertainties arising from the hadronic parameters, such as ω, f B , and m 0 , in the meson wave function, the heavy quark masses and the hard scale, which can give sizable effects on the branching ratios, whereas the CP asymmetries are found to be relatively stable with respect to variations of hadronic parameters. The obtain results with reasonably accurate will be tested at the ongoing and forthcoming hadron colliders. where with v the charm quark velocity, which denote the relativistic corrections to the Coulomb wave functions [67]. We take v 2 = 0.3 for charmonium. The normalization constants N L,T,S can be determined by the corresponding normalization conditions [52]. The kaon and pion meson distribution amplitudes up to twist-3 are given by light-cone QCD sum rules [60,68] with the Gegenbauer polynomials