Probing new physics in class-I $B$-meson decays into heavy-light final states

With updated experimental data and improved theoretical calculations, several significant deviations are being observed between the Standard Model predictions and the experimental measurements of the branching ratios of $\bar{B}_{(s)}^0\to D_{(s)}^{(*)+} L^-$ decays, where $L$ is a light meson from the set $\{\pi,\rho,K^{(\ast)}\}$. Especially for the two channels $\bar{B}^0\to D^{+}K^-$ and $\bar{B}_{s}^0\to D_{s}^{+}\pi^-$, both of which are free of the weak annihilation contribution, the deviations observed can even reach 4-5$\sigma$. Here we exploit possible new-physics effects in these class-I non-leptonic $B$-meson decays within the framework of QCD factorization. Firstly, we perform a model-independent analysis of the effects from twenty linearly independent four-quark operators that can contribute, either directly or through operator mixing, to the quark-level $b\to c\bar{u} d(s)$ transitions. It is found that, under the combined constraints from the current experimental data, the deviations observed could be well explained at the $1\sigma$ level by the new-physics four-quark operators with $\gamma^{\mu}(1-\gamma_5)\otimes\gamma_{\mu} (1-\gamma_5)$ structure, and also at the $2\sigma$ level by the operators with $(1+\gamma_5)\otimes(1-\gamma_5)$ and $(1+\gamma_5)\otimes(1+\gamma_5)$ structures. However, the new-physics four-quark operators with other Dirac structures fail to provide a consistent interpretation, even at the $2\sigma$ level. Then, as two specific examples of model-dependent considerations, we discuss the case where the new-physics four-quark operators are generated by either a colorless charged gauge boson or a colorless charged scalar, with their masses fixed both at the $1$~TeV. Constraints on the effective coefficients describing the couplings of these mediators to the relevant quarks are obtained by fitting to the current experimental data.

Abstract: With updated experimental data and improved theoretical calculations, several significant deviations are being observed between the Standard Model predictions and the experimental measurements of the branching ratios ofB 0 (s) → D ( * )+ (s) L − decays, where L is a light meson from the set {π, ρ, K ( * ) }. Especially for the two channelsB 0 → D + K − and B 0 s → D + s π − , both of which are free of the weak annihilation contribution, the deviations observed can even reach 4-5σ. Here we exploit possible new-physics effects in these class-I non-leptonic B-meson decays within the framework of QCD factorization. Firstly, we perform a model-independent analysis of the effects from twenty linearly independent fourquark operators that can contribute, either directly or through operator mixing, to the quark-level b → cūd(s) transitions. It is found that, under the combined constraints from the current experimental data, the deviations observed could be well explained at the 1σ level by the new-physics four-quark operators with γ µ (1 − γ 5 ) ⊗ γ µ (1 − γ 5 ) structure, and also at the 2σ level by the operators with (1 + γ 5 ) ⊗ (1 − γ 5 ) and (1 + γ 5 ) ⊗ (1 + γ 5 ) structures. However, the new-physics four-quark operators with other Dirac structures fail to provide a consistent interpretation, even at the 2σ level. Then, as two specific examples of model-dependent considerations, we discuss the case where the new-physics four-quark operators are generated by either a colorless charged gauge boson or a colorless charged scalar, with their masses fixed both at the 1 TeV. Constraints on the effective coefficients describing the couplings of these mediators to the relevant quarks are obtained by fitting to the current experimental data.

Introduction
Flavor physics plays always an important role in testing the Standard Model (SM) of particle physics and probing new physics (NP) beyond it [1,2]. Here, the non-leptonic weak decays of B mesons are of particular interest, since they provide direct access to the fundamental parameters of the Cabibbo-Kobayashi-Maskawa (CKM) matrix [3,4] and further insight into the strong-interaction dynamics involved in these decays. Aiming at such a goal, the BaBar and Belle collaborations [5], as well as the LHCb experiment [6] have already performed many high-precision measurements of these kinds of decays [7,8]. In addition, new frontiers of precision are expected in the era of Belle II [9] and upgraded LHCb [10].
Confronted with the plethora of high precision measurements made by these dedicated experiments, we are forced to improve as much as possible the accuracy of theoretical predictions about these non-leptonic weak decays. Here the main challenge we are now facing is how to calculate reliably the hadronic matrix elements of four-quark operators in the effective weak Hamiltonian (see subsection 2.1). For a long time, the naive factorization (NF) assumption [11] and modifications thereof (see, e.g., refs. [12][13][14][15] and references therein) were used to estimate the non-leptonic B-decay amplitudes. Several more promising strategies built upon either the SU (3) flavor symmetry of strong interactions [16][17][18] or the factorization theorem, such as the QCD factorization (QCDF) [19][20][21] and its field theoretical formulation, the soft-collinear effective theory [22][23][24][25][26], as well as the perturbative QCD [27][28][29] approach, have been developed to study the same problem. Certain combinations of these approaches have also been adopted in, e.g., refs. [30][31][32].
In this paper, we shall consider the exclusive two-bodyB 0 (s) → D ( * )+ (s) L − decays, where L ∈ {π, ρ, K ( * ) }, within the QCDF framework. For these class-I non-leptonic decays, the spectator antiquark and other light degrees of freedom of the initialB 0 (s) mesons need to rearrange themselves only slightly to form the heavy D , and the light-cone distribution amplitude (LCDA), Φ L (u), of the light meson encode all the long-distance strong-interaction effects, both of which can be extracted from experimental data or calculated by using nonperturbative methods like QCD sum rules and/or lattice QCD. The hard kernels T ij (u) receive, on the other hand, contributions only from scales of O(m b ) in the heavy-quark limit and are therefore calculable perturbatively. At leading order (LO) in the strong coupling α s , eq. (1.1) reproduces the NF result, and both the next-to-leading-order (NLO) [20,34] and next-to-next-to-leading-order (NNLO) [35,36] corrections to T ij (u) are now known.
As all the four quark flavors in b → cūd (s) transitions are different from each other, these tree-level decays receive contributions from neither the penguin operators nor the penguin topology. There is also no color-suppressed tree topology in these class-I decays. At leading power in Λ QCD /m b , these decays are dominated by the color-allowed tree topology that receives only vertex corrections, while interactions with the spectator antiquark and the weak annihilation topology are both power-suppressed [20]. In fact, noting that the weak annihilation topology contributes only toB 0 → D ( * )+ π − andB s → D ( * )+ s K − , but not tō B 0 → D ( * )+ K − andB s → D ( * )+ s π − , one can use the ratios between the branching fractions of these two kinds of decays to probe the topology. Remarkably, the current experimental data shows already that the impact due to such a topology is negligible [37]. Other sources of power corrections, such as the higher-twist corrections to the light-meson LCDAs and the exchange of a single soft gluon between the B (s) → D ( * ) (s) transitions and the light meson, are also estimated to be quite small [20,38]. Therefore, these class-I non-leptonic decays are theoretically clean and the QCDF approach is expected to work well for them. However, with the updated input parameters, the SM predictions [36,38,39] are found to be generically higher than the current experimental measurements [7,8] of the branching ratios ofB 0 (s) → D predictions for the branching ratios of these class-I decays and their ratios with respect to the semi-leptonicB 0 (s) → D ( * )+ (s) −ν decay rates evaluated at q 2 = m 2 L , R ( * ) (s)L , and then discuss the NP effects both in a model-independent setup and in the case where the NP operators are generated by either a colorless charged gauge boson or a colorless charged scalar. Our conclusions are finally made in section 4. For convenience, the ranges for the NP Wilson coefficients C i (m b ) allowed by the ratios R ( * ) (s)L are given in the appendix.
2 Theoretical framework

Effective weak Hamiltonian
The class-IB 0 (s) → D ( * )+ (s) L − decays are mediated by the quark-level b → cūd(s) transitions. Once the top quark, the gauge bosons W ± and Z 0 , the Higgs boson, as well as other heavy degrees of freedom present in any extension of the SM are integrated out, the corresponding QCD amplitudes of the decays are computed most conveniently in the framework of effective weak Hamiltonian [46,47], which for the problem at hand reads 3 where G F is the Fermi constant, and V cb V * uq (q = d, s) the product of the CKM matrix elements. Q i (i = 1, 2) are the two SM four-quark current-current operators given in the Buchalla-Buras-Lautenbacher (BBL) basis [46], while the remaining ones in eq. (2.1) denote the full set of twenty linearly independent four-quark operators that can contribute, either directly or through operator mixing, to the quark-level b → cūd(s) transitions [48][49][50].
The NP four-quark operators can be further split into eight separate sectors, among which there is no mixing [50,51]. Firstly, the operators belonging to the two sectors V LL and V LR, which are relevant for tree-level contributions mediated by heavy charged gauge bosons in any extension of the SM, can be written, respectively, as [50,51] where α, β are the color indices, and Q V LL i are identical to the SM operators Q i given in the BBL basis [46]. Secondly, the operators belonging to the two sectors SLL and SLR, which are relevant for tree-level contributions generated by new heavy charged scalars, are given, respectively, by [50,51] Figure 1. Feynman rules for the couplings of a colorless charged gauge boson A + (upper) and a colorless charged scalar H + (lower) to an up-(i α ) and a down-type (j β ) quark, with the strengths normalized to that of the SM tree-level W + exchange, where g 2 is the SU (2) L gauge coupling and P L(R) = 1 2 (1 ∓ γ 5 ) denote the left-and right-handed chirality projectors.
Finally, the operators belonging to the four remaining chiralityflipped sectors V RR, V RL, SRR and SRL are obtained, respectively, from eqs. (2.2)-(2.5) by making the interchanges (1 ∓ γ 5 ) ↔ (1 ± γ 5 ). It should be noted that, due to parity invariance of strong interactions, the QCD ADMs of the chirality-flipped sectors are identical to that of the original sectors, simplifying therefore the renormalization group (RG) analysis of these operators [50].
The short-distance Wilson coefficients C i (µ) and C i (µ) in eq. (2.1) can be reliably calculated by using the RG-improved perturbation theory [46,47]. Explicit expressions up to the NNLO in α s for the SM part can be found, e.g., in ref. [52], and will be used throughout this paper. For the NP part, on the other hand, one can easily obtain the NLO results of C i (µ b ) evaluated at the typical scale µ b m b that is appropriate for the non-leptonic B-meson decays, by solving the RG equations satisfied by these short-distance Wilson coefficients, based on the one-and two-loop QCD ADMs of the NP four-quark operators [48][49][50], as well as the O(α s ) corrections to the matching conditions for C i (µ 0 ) evaluated at the NP scale µ 0 [51]. Here, for later convenience, we show in Fig. 1 the Feynman rules describing the couplings of both a colorless charged gauge boson A + and a colorless charged scalar H + to an up-(i α ) and a down-type (j β ) quark, the strengths of which have been normalized to that of the tree-level W + exchange in the SM. For further details about the matching and evolution procedures in the case of these tree-level mediators, the readers are referred to ref. [51] and references therein. Throughout this paper, we shall assume that  Figure 3. "Non-factorizable" vertex corrections to the hard kernels T ij (u) at the NLO in α s , where the other captions are the same as in Fig. 2.
the NP Wilson coefficients C i (µ) as well as the effective couplings ∆ L,R ij (A) and ∆ L,R ij (H) are all real, and take the same values for both the b → cūd and b → cūs transitions.

Calculation of one-loop vertex corrections
To obtain the non-leptonic B-decay amplitudes, we must also calculate the hadronic matrix elements of the local four-quark operators present in eq. (2.1). To this end, we shall adopt the QCDF approach [19][20][21], within the framework of which the hadronic matrix element of a four-quark operator satisfies the factorization formula given by eq. (1.1). For the SM contribution, the hard kernels T ij (u) have been calculated up to the NNLO in α s [35,36], and will be used throughout this paper. For the NP contribution, on the other hand, we shall calculate the one-loop vertex corrections to the hard kernels T ij (u), completing therefore a full NLO analysis of the class-I non-leptonicB 0 (s) → D ( * )+ (s) L − decays in the case where the short-distance Wilson coefficients of the four-quark operators are also known at the same order [51]. Such an NLO analysis in the NP sector is helpful for reducing the dependence of the final decay amplitudes on certain unphysical scale and renormalization scheme [51], as will be detailed in subsection 3.5.
As mentioned already in the last section, at leading power in Λ QCD /m b , these class-I non-leptonic decays are dominated by the color-allowed tree topology with the lowest-order Feynman diagram shown in Fig. 2, and the hard kernels T ij (u) receive only the "nonfactorizable" vertex corrections [20], with the corresponding one-loop Feynman diagrams shown in Fig. 3. It should be noted that, as the light quark-antiquark pair (ūq) has to be in a color-singlet configuration to produce an energetic light meson L in the leading Fock-state approximation, the hard kernels T ij (u) only receive non-vanishing contributions from the color-singlet operators starting at the zeroth order in α s and from the color-octet operators starting at the first order in α s , respectively. This implies that T ij (u) ∝ 1 + O(α 2 s ) + · · · for the color-singlet and T ij (u) ∝ O(α s ) + · · · for the color-octet operators, respectively. It is also observed that, although there exist both collinear and infrared divergences in each of the four vertex diagrams shown in Fig. 3, these divergences cancel when one sums over all the four diagrams, yielding therefore a finite and perturbatively calculable O(α s ) correction to the hard kernels T ij (u) [20,21]. Explicit evaluations of these one-loop vertex diagrams with insertions of the SM current-current operators in the Chetyrkin-Misiak-Münz (CMM) basis [53,54] can be found, e.g., in refs. [20,36,55]. Our results for the one-loop vertex corrections to the hard kernels T ij (u), which arise from insertions of the NP four-quark operators present in eq. (2.1) into the Feynman diagrams shown in Fig. 3, will be presented below. This, together with the NLO results of the NP Wilson coefficients C i (µ b ) [48][49][50][51], completes our analysis at the NLO in α s .
where q = p−p , 4 and the upper (lower) sign applies when L is a pseudoscalar (vector) meson. f L and Φ L denote respectively the decay constant and the leading-twist LCDA of the light meson L, while the reduced matrix elements D where C F = (N 2 c − 1)/(2N c ), with N c = 3 being the number of colors, and Here z = m c /m b ,ū = 1 − u, and the dilogarithm is defined by (2.10) In the limit z → 0, our results coincide with that for the charmless B-meson decays presented in refs. [21,56], where the four-quark operators are also defined in the BBL basis. In addition, we have also checked explicitly that, by using the relations among the short-distance Wilson coefficients [52,54,57] corresponding to the four-quark operators defined in the BBL and CMM bases respectively, our results for the hadronic matrix elements D given in the BBL basis, agree up to the NLO in α s with that presented in refs. [36,55,58], where the calculations are however performed with the four-quark operators defined in the CMM basis. To reproduce the results presented in ref. [20], on the other hand, one should keep in mind that the LO relations among the short-distance Wilson coefficients are used when transforming from one operator basis to another.
where the one-loop hard kernel T V LR (u, z) is now given by 14) It has been checked that, in the limit z → 0, the above results are also reduced to that for the charmless B-meson decays given in refs. [21,56].
where the parameters µ m are defined, respectively, as µ p (µ) = m 2 L /[m u (µ) + m q (µ)] for a pseudoscalar and µ v (µ) = m L f ⊥ L (µ)/f L for a vector meson, with m u,q (µ) being the running quark masses in the MS scheme and f ⊥ L (µ) the scale-dependent transverse decay constant of a vector meson. When all three-particle contributions are neglected, the twist-3 two-particle LCDA Φ p (u) is determined completely by the equations of motion, with its asymptotic form given exactly by Φ p (u) = 1 [56,59], while Φ v (u) is related to the twist-2 LCDA Φ ⊥ (u) of a transversely polarized vector meson by [56,60] Φ where the second equation is obtained by inserting the Gegenbauer expansion of Φ ⊥ (u), with α L n,⊥ (µ) being the Gegenbauer moments with α L 0,⊥ = 1, and P n (x) the Legendre polynomials. For further details about these hadronic parameters, the readers are referred to ref. [56] and references therein.
The reduced matrix elements of the scalar and pseudoscalar currents are related, respectively, to that of the vector and axial-vector currents by The one-loop hard kernel T SLL (u, z) reads where the one-loop hard kernel T T LL (u, z) is given by where the upper (lower) sign applies when L is a pseudoscalar (vector) meson, and the one-loop hard kernel T SLR (u, z) reads and It is noted that, in the limit z → 0, our results are consistent with that for the charmless B-meson decays presented in refs. [21,56,61].
The one-loop vertex corrections to the hard kernels T ij (u) with insertions of the chiralityflipped four-quark operators can be easily obtained from the results given above by changing, if necessary, the overall signs of the reduced matrix elements D + (s) |c · · · b|B 0 (s) and D * + (s) |c · · · b|B 0 (s) . It should be noted that our calculations of the hadronic matrix elements of these four-quark operators are performed in the naive dimensional regularization scheme with anti-commuting γ 5 in D = 4 − 2 dimensions, which matches exactly the one used for evaluations of the short-distance Wilson coefficients C i (µ) [50,51]. This ensures, therefore, the renormalization scheme and scale independence of the non-leptonic decay amplitudes up to the NLO in

Estimate of weak annihilation contribution
We now proceed to discuss the power-suppressed weak annihilation contribution to the class-IB 0 (s) → D ( * )+ (s) L − decays, with the corresponding Feynman diagrams shown in Fig. 4. It must be emphasized that, due to the presence of endpoint singularities, the weak annihilation topology cannot be computed self-consistently within the QCDF framework [20,21]. Nevertheless, we shall still follow the conventions used in refs. [21,56,62] to make an estimate of the weak annihilation effect in these class-I decays. Instead of considering all the four-quark operators present in eq. (2.1), we shall focus only on the SM current-current operators. Our purpose is to demonstrate that, even with the weak annihilation contribution taken into account, the deviations observed in the branching ratios of these class-I decays could not be explained in the SM, as will be shown numerically in subsection 3.2.
Following the same conventions as used in refs. [21,56,62], one can write the weak annihilation contribution to the decay amplitude of a class-I non-leptonic decay as with the upper sign applied when both final-state mesons are pseudoscalar or longitudinally polarized vector mesons, and the lower when one of them is a vector meson. The Wilson coefficient C 2 (µ h ) and the building blocks . Since the treatment of weak annihilation topology within the QCDF framework is model-dependent anyway, we shall assume that the building blocks A i 1 (µ h ) take the same expressions as for the charmless B-meson decays [21,56,62], although the asymptotic forms of the D ( * ) (s) LCDAs are quite different from that of a light charged charmless meson [20]. Explicitly, we have [21,56,62] when both final-state mesons are pseudoscalar, when one of them is a pseudoscalar and the other a vector meson, whereas , which means that we have assigned a 200% uncertainty to the default value obtained with A = 0. The ratios r M χ are defined as for a pseudoscalar (P ) and a vector (V ) meson respectively, where m 1,2 are the current quark masses of the two valence constituents of the meson considered. In view of the large uncertainties brought by the phenomenological parameters A and ϕ A , it is generally expected that such a model-dependent treatment should give the correct order of magnitude of the weak annihilation effect in B-meson decays into both the charmless [21,56,62] and the heavy-light final states [20,43]. In order to separate the weak annihilation contribution from that of the dominant color-allowed tree topology, we introduce the effective coefficients b 1 (D is the momentum of the two final-state mesons in the parent rest frame. Then, including also the LO contributions from the four-quark operators present in eq. (2.1), we can write the total decay amplitude of a given channel as [20,36]: 5 Note that, due to angular momentum conservation, the polarization vectors µ of D * + (s) and η µ of V in the final states take only the longitudinal part in eqs. (2.39) and (2.40). The decay amplitudes ofB 0 (s) → D * + (s) V − modes are more complicated and, to leading power in Λ QCD /m b , dominated also by the longitudinal polarization, with the transverse parts being suppressed by O(m V /m B (s) ); their explicit expressions could be found, e.g., in ref. [20]. The effective coefficients a 1 (D * + (s) L − ) can be expressed in terms of the shortdistance Wilson coefficients C i (µ) as well as the perturbatively calculable hard kernels T ij (u) convoluted with the light-meson LCDAs Φ L,m (u). For the SM contributions, both the NLO [20,34] and the NNLO [35,36] corrections to a 1 (D * + (s) L − ) are known. Combining our calculations of the one-loop vertex corrections to T ij (u) as well as the O(α s ) corrections to the matching conditions for the short-distance Wilson coefficients [51], the effective coefficients a 1 (D * + (s) L − ) associated with the complete set of NP four-quark operators present in eq. (2.1) are now known up to the NLO in α s .

Input parameters
To update the SM predictions ofB 0 (s) → D ( * )+ (s) L − decays presented in ref. [36], we should firstly update the theoretical input parameters, which include the strong coupling constant α s , the quark masses, the CKM matrix elements, the lifetimes of B 0 (s) mesons, as well as the hadronic parameters like the B (s) → D ( * ) (s) transition form factors and the decay constants and Gegenbauer moments of light mesons. We use the two-loop relation between pole and MS mass [74], to convert the top-quark pole mass m pole t to the scale-invariant mass m t (m t ). The three-loop running for α s is used throughout this paper. For convenience, we collect in Table 1 all the input parameters used throughout this paper. To obtain the theoretical uncertainty of an observable, we vary each input parameter within its 1σ range and then add each individual uncertainty in quadrature. In addition, we have included the uncertainty due to the variation of the renormalization scale For the B → D ( * ) transition form factors, we take the "L w≥1 +SR" fit results obtained in ref. [75], in which both O(Λ QCD /m b,c ) and O(α s ) contributions as well as the uncertainties in the predictions of the form-factor ratios at O(Λ QCD /m b,c ) are consistently included within the framework of heavy quark effective theory (HQET). 6 For the B s → D  Table 4. 6 For other similar analyses including the missing higher-order pieces in the relations among the HQET form factors, the readers are referred to refs. [76][77][78][79][80][81][82] and references therein.
decay rates are taken from the LHCb collaboration [85,86]. 7 As a comparison, we list in Table 2 our results for some of these transition form factors evaluated at q 2 = m 2 K − or q 2 = m 2 π − , together with the ones used in refs. [36,38]. It can be seen that our results for these transition form factors are all consistent within errors with that presented in ref. [38], while being much more precise than the ones used in ref. [36]. This justifies our choices of the form factors given in refs. [75,83,84], rather than adopting the more complete analysis performed in the heavy-quark-expansion framework, where the form-factor uncertainties including O(Λ 2 QCD /m 2 c ) corrections, the strong unitarity bounds, and a consistent treatment of the flavor symmetry (breaking) are all taken into account in the global fit [76,77].

Updated predictions for branching ratios
Our updated SM predictions for the branching ratios ofB 0 (s) → D ( * )+ (s) L − decays at different orders in α s are given in Table 3, together with the results obtained in refs. [36,38] as a comparison. The experimental data is taken from the Particle Data Group [7] and/or the Heavy Flavor Averaging Group [8]. For the decay modesB 0 s → D + s π − andB 0 s → D + s K − , we have also replaced the previous LHCb results by the latest updates [88] when performing the averages given in ref. [8]. It can be seen that our updated results are generally higher than the current experimental data, even at the LO, and the higher-order perturbative corrections always add constructively to the LO results. 8 Especially for the decay modes level and, after taking into account the theoretical and experimental uncertainties, the deviation can even reach about 4-5σ. It must be pointed out that such a large deviation has been observed for the first time in ref. [38], where the values of B (s) → D ( * ) (s) transition form factors were taken from ref. [76]. Compared with the results presented in ref. [36], our updated central values of the branching ratios ofB 0 decays are increased by about 16% for D + and 20% for D * + final states, respectively. This is mainly due to the following 7 It should be noted that the unitarity bounds applied in the lattice [83,84] and experimental [85,86] determinations of the Bs → D ( * ) s transition form factors are much looser than imposed in refs. [75][76][77]. 8 Due to the sizable decay width difference, ys ≡ ∆Γs 2Γs = 0.062 ± 0.004 [8], the measured branching ratios ofB 0 s decays should be multiplied by a factor 1 − y 2 s , when compared with the theoretical predictions [89]. 1.63 ± 0.50 Table 3.  [36,38] as a comparison. The column marked by NNLO # represents our results obtained with the weak annihilation contribution included. For the channelB 0 → D * + ρ − , only the longitudinal polarization amplitude is considered. The experimental data is taken from refs. [7,8], with the longitudinal polarization fraction ofB 0 → D * + ρ − decay taken from ref. [87]. For the decay modesB 0 s → D + s π − andB 0 s → D + s K − , the previous LHCb results have been superseded by the latest updates [88] when performing the averages given in ref. [8]. The measured branching ratios ofB 0 s decays should also be multiplied by a factor 1 − y 2 s , with y s = 0.062 ± 0.004 [8], when compared with the corresponding theoretical predictions [89].
two reasons: Firstly, our input of the CKM matrix element |V cb | is about 7.4% larger than the one used in ref. [36], where the value of |V cb | extracted from exclusive semi-leptonic Bmeson decays as of 2016 was used instead. Secondly, our inputs for the B → D and B → D * transition form factors, whose theoretical information available since the analysis of ref. [36] has been systematically taken into account [75], are now about 4.7% and 6.5% larger than the corresponding ones used in ref. [36], when evaluated at q 2 = 0. It is also observed that the theoretical uncertainties of the branching ratios ofB 0 K − decays are significantly reduced with respect to that obtained in ref. [36]. This is mainly due to the improved lattice determinations of the B s → D ( * ) s transition form factors [83,84]. Especially for the two channelsB 0 s → D * + s π − andB 0 s → D * + s K − , our updated results are about 70% larger than that given in ref. [36], due mainly to the increase of the transition  Table 4. Our estimates of the weak annihilation coefficients |b 1 (D  Table 2). On the other hand, our central values of the branching ratios are slightly larger for theB 0 (s) → D + (s) L − but smaller for thē B 0 (s) → D * + (s) L − decays than the corresponding ones presented in ref. [38], while being in perfect agreement within errors. This is also attributed mainly to the different inputs of the CKM matrix element |V cb | and the B (s) → D ( * ) (s) transition form factors. As can be seen from Table 3, our estimate of the weak annihilation contribution, although being plagued by large uncertainties due to the model parameters A and ϕ A , always contributes constructively to the dominant color-allowed tree amplitude, with its effect being less than 5% on the final branching ratios. Thus, the weak annihilation effect does not help to reconcile the deviations observed in the class-I non-leptonic decays [36,38]. To see the relative size of the weak annihilation contribution in these decays, we show in Table 4 our estimates of the effective coefficients |b 1 (D The latter can also be extracted from the measured ratios of branching fractions [8], [37]. The values obtained in such a way are shown in the last column of Table 4 as a comparison. It can be seen that both methods give similar magnitudes of the weak annihilation contributions, with our estimated results being positive while the extracted values negative, and no sign of an enhanced weak annihilation topology is shown in these class-I non-leptonic decays, as is generally expected both within the QCDF framework [20,43] and based on the SU (3)-flavor symmetry of strong interactions [37]. As a consequence, from now on, we shall no longer consider the weak annihilation contribution. with only a ∼ 11% precision [84]. Furthermore, different choices of |V cb | also affect the final results of the absolute branching ratios [36,38]. In order to minimize the impacts of these input parameters, one can consider the ratios of the non-leptonicB 0 (s) → D rates evaluated at q 2 = m 2 L , where refers to either an electron or a muon, and q 2 is the four-momentum squared transferred to the lepton pair. In this way, one obtains [12,20,90] which by construction are free of the uncertainty related to |V cb |. Neglecting the masses of light leptons, we have exactly X L = X * L = 1 for a vector meson L, valid both for the sum of and separately for the longitudinal and the transverse polarization of the D * + (s) mesons in the final state. This is due to the kinematic equivalence between the production of the lepton pair via the SM weak current with γ µ (1 − γ 5 ) structure in semi-leptonic decays and that of a vector meson with four-momentum q µ in non-leptonic decays [12,20]. For a light pseudoscalar meson L, on the other hand, X  Table 5, together with the available results presented in refs. [36,38]. In addition, we give in Table 6 the values of the ratios R ( * ) (s)L extracted from the current experimental data as well as our updated theoretical predictions at different orders in α s , which will be used later to analyze the NP effects in these class-I non-leptonic decays.
From Table 5, one can see that our results for the effective coefficients |a 1 (D ( * )+ (s) L − )| at the NNLO in α s are well consistent with that obtained in ref. [36], up to slight variations induced by the updated input parameters from α s (m Z ), the Gegenbauer moments, as well as the quark masses. 10 As emphasized already in refs. [20,36], an essentially universal value of |a 1 (D  Ref. [36] Ref. [38] Exp.   framework, which is however consistently higher than the central values fitted from the current experimental data. As shown in the last column of Table 6, the deviations observed inB 0 (s) → D ( * )+ (s) π − andB 0 (s) → D ( * )+ (s) K − decay modes are particularly remarkable, with some of them reaching even up to 4-5σ. This is attributed to the increased theoretical predictions [36] and, at the same time, the decreased experimental measurements [7,8] of the absolute branching ratios, together with their reduced uncertainties, as compared to the previous analysis performed at the NLO in α s within the same framework [20].
As pointed out already in refs. [36,38], it is quite difficult to understand the large deviations observed in these class-I non-leptonic B-meson decays in the SM, by simply considering the higher-order power and perturbative corrections to the decay amplitudes based on the QCDF approach [20,92]. Thus, as an alternative, we shall in the next subsections resort to possible NP explanations of these deviations, firstly in a model-independent setup by considering the NP effects from twenty linearly independent four-quark operators present in eq. (2.1), and then within two model-dependent scenarios where the NP four-quark operators are mediated by either a colorless charged gauge boson or a colorless charged scalar. See also refs. [40][41][42][43][44][45] for recent discussions along this line.  0.42 ± 0.14 2.5 Table 6. Theoretical and experimental values of the ratios R

Model-independent analysis
With our prescription for the effective weak Hamiltonian given by eq. (2.1), possible NP effects would be signaled by the non-vanishing NP Wilson coefficients C i that accompany the corresponding NP four-quark operators. As a model-independent analysis, we shall use the ten ratios R ( * ) (s)L collected in Table 6 to constrain these NP Wilson coefficients C i , both at the characteristic scale µ b = m b (low-scale scenario) and at the electroweak scale µ W = m W (high-scale scenario).

Low-scale scenario
Firstly, let us consider the case where only a single NP four-quark operator is present in eq. (2.1). The resulting constraints on the corresponding NP Wilson coefficient C i (m b ) from the ratios R  Table 7 given in the appendix. 11 In this case, the following observations can be made: (m b ) (right) from the ratios R π (red), R * π (blue), R ρ (yellow), R K (green), R * K (pink), R K * (orange), R sπ (purple), R sK (black), R * sπ (cyan), as well as R * sK (brown), respectively. The horizontal bounds represent the experimental ranges within 1σ (dark gray) and 2σ (light gray) error bars.
• As can be seen from Fig. 5, all the deviations observed inB 0 (s) → D ( * )+ (s) L − decays could be explained simultaneously by the two NP four-quark operators with γ µ (1 − γ 5 ) ⊗ γ µ (1 − γ 5 ) structure. As these two operators appear already in the SM, this means that we can account for the observed deviations collected in Table 6 by a shift to the SM Wilson coefficients C 1 and/or C 2 . The final allowed ranges for the NP Wilson coefficients C V LL One can see that the constraint on This is due to the fact that C V LL 2 (m b ) gives the leading contribution to the effective coefficients a 1 (D is suppressed by 1/N c at the LO and further by C F /4π at the NLO in α s , within the QCDF framework [20,36]. • From Figs. 6 and 7, one can see that the NP four-quark operators with either (1 + γ 5 ) ⊗ (1 − γ 5 ) or (1 + γ 5 ) ⊗ (1 + γ 5 ) structure could also be used to account for the observed deviations but now at the 2σ level, with the corresponding allowed ranges for the NP Wilson coefficients given, respectively, by The much weaker constraints on C are also due to the fact that the latter always provide the leading contributions to the hard kernels T ij (u). For the decay modes where L is a light pseudoscalar meson, the hadronic matrix elements of these (pseudo-)scalar four-quark operators, although being formally power-suppressed, would be chirally-enhanced by the factors 2µ p (µ)/(m b (µ) ∓ m c (µ)) and hence be not much suppressed numerically for realistic bottom-and charm-quark masses [21,56]. This explains the important role played by these (pseudo-)scalar four-quark operators in non-leptonic B-meson decays both  within the SM [21,56] 12 and in various NP models [42,61,[93][94][95][96][97].
• As can be seen from Table 7, the remaining NP four-quark operators with other Dirac structures present in eq. (2.1) are already ruled out by the combined constraints from the ten ratios R ( * ) (s)L collected in Table 6, even at the 2σ level. As the matrix elements L − |q(1±γ 5 )u|0 ≡ 0 for a light charged vector meson L − , there is no LO contribution to the decay amplitudes ofB 0 (s) → D + (s) ρ − andB 0 (s) → D + (s) K * − decays from the NP four-quark operators with (1 ± γ 5 ) ⊗ (1 ± γ 5 ) and (1 ± γ 5 ) ⊗ (1 ∓ γ 5 ) structures. This explains why the two ratios R ρ and R K * receive insignificant contributions from these operators (see also the third and the sixth plot in Figs. 6 and 7). For the NP fourquark operators with σ µν (1 ± γ 5 ) ⊗ σ µν (1 ± γ 5 ) structures, on the other hand, the ratios R π , R * π , R K and R * K receive only negligible contributions from the NP Wilson coefficients C SLL depend crucially on whether the final-state heavy mesons are D + (s) or D * + (s) , as shown in Figs. 8 and 9. As a consequence, the tensor four-quark operators also fail to provide a simultaneous explanation of the ten ratios R ( * ) (s)L collected in Table 6, even at the 2σ level.
• Due to the relatively larger experimental uncertainties of the three ratios R ρ , R * K , and R K * , their constraints on the NP Wilson coefficients are much weaker. More precise measurements of these decay modes are, therefore, expected from the LHCb [10] and Belle II [9] experiments, which will be helpful to further discriminate the NP contributions from We now consider the case where two NP four-quark operators with the same Dirac but different color structures are present in eq. (2.1), and allow the corresponding two NP Wilson coefficients to vary simultaneously. To obtain the allowed regions for the NP Wilson coefficients, we follow the strategies used in refs. [98,99]: each point in the NP parameter space corresponds to a theoretical range constructed for the ratios R ( * ) (s)L in the point, with the corresponding theoretical uncertainty taken also into account. If this range has overlap with the 2σ range of the experimental data on R ( * ) (s)L , this point is then assumed to be allowed. Here the theoretical uncertainty at each point in the NP parameter space is obtained in the same way as in the SM, i.e., by varying each input parameter within its respective range and then adding the individual uncertainty in quadrature. Such a treatment is motivated by the observation that, while the experimental data yields approximately a Gaussian distribution for the branching ratios ofB 0 (s) → D ( * )+ (s) L − decays, a theoretical calculation does not. As the latter depends on a set of hadronic input parameters like the heavy-toheavy transition form factors as well as the decay constants and Gegenbauer moments of the light mesons, for which no probability distribution is known, it is more suitable to assume that these theory parameters have no particular distribution but are only constrained in certain allowed ranges with an equal weighting, irrespective of how close they are from the edges of the allowed ranges [64,100].
In the case where two NP Wilson coefficients are present simultaneously, we show in Fig. 10 the allowed regions in the ( planes, under the combined constraints from the ten ratios R ( * ) (s)L varied within 2σ error bars. It is readily to see that, due to the partial cancellation between contributions from the two NP Wilson coefficients, the allowed regions for the NP parameter space become potentially larger than in the case where only one NP Wilson coefficient is present. In the presence of two NP four-quark operators with other Dirac structures, on the other hand, there exist no allowed regions for the corresponding NP Wilson coefficients that can provide a simultaneous explanation of the ratios R

High-scale scenario
From the point of view of constructing specific NP models and correlating the low-energy constraints with the direct searches performed at high-energy frontiers, it is also interesting to provide constraints on the NP Wilson coefficients C i (µ W ) that are given at the electroweak scale µ W = m W . To this end, we must take into account the RG evolution of these short-distance Wilson coefficients from µ W down to the low-energy scale µ b = m b , at which the hadronic matrix elements of the NP four-quark operators are evaluated. The most generic formulae for the RG equations satisfied by the NP Wilson coefficients C i (µ) can be written as where γ ij are the QCD ADMs of the NP four-quark operators, with their one-and two-loop results given already in refs. [48][49][50]. By solving eq. (3.4), one can then obtain the evolution matricesÛ (µ b , µ W ), which connect the Wilson coefficients at different scales [46,47]: where, once specific to our case with the effective weak Hamiltonian given by eq. (2.1), C is a two-dimensional column vector andÛ (µ b , µ W ) a 2 × 2 matrix for each V LL (V RR), V LR (V RL), SLR (SRL) sector, while C is a four-dimensional column vector andÛ (µ b , µ W ) a 4 × 4 matrix in the SLL (SRR) sector [50].
Here, instead of re-performing a detailed analysis of the NP effects at the electroweak scale, we focus only on the case where only a single NP four-quark operator is present in eq. (2.1), and investigate how the three solutions obtained in the low-scale scenario change when looked at the electroweak scale. Following the same way as in the low-scale scenario, we show in Figs. 11-13 the allowed ranges for the NP Wilson coefficients C i (m W ), under the constraints from the ratios R π , R * π , R ρ , R K , R * K , R K * , R sπ , R sK , R * sπ , as well as R * sK . It is found that, due to the RG evolution, the resulting allowed range for C V LL 12, 4.86], and should be therefore discarded as it is much larger than the SM case with C 1 (m W ) = 0 and C 2 (m W ) = 1 at the LO in α s [47], while the allowed range for C V LL 2 (m W ) remains almost the same as in the low-scale scenario (see eq. (3.2)), with under the combined constraints from the ten ratios R ( * ) (s)L at the 1σ level. On the other hand, the NP four-quark operators with either (1 + γ 5 ) ⊗ (1 − γ 5 ) or (1 + γ 5 ) ⊗ (1 + γ 5 ) structure, could still provide a reasonable explanation of the deviations observed inB 0 (s) → D ( * )+ (s) L − decays at the 2σ level, with the resulting allowed ranges for the NP Wilson coefficients given, respectively, by which, compared with the results obtained in the low-scale scenario (see eq. (3.3)), indicate a large RG evolution effect in these (pseudo-)scalar four-quark operators [50].

Model-dependent analysis
As found in the last subsection, the deviations observed inB 0 (s) → D  Table 6, we can then obtain constraints on the effective coefficients describing the couplings of these mediators to the relevant quarks (see Fig. 1).

Colorless charged gauge boson
Starting with the Feynman rules given in Fig. 1 and after integrating out the heavy colorless charged gauge boson A + , we can obtain the effective weak Hamiltonian describing the quark-level b → cūd(s) transitions mediated by A + [51]: where m A is the mass of the colorless charged gauge boson A + , and ∆ L,R i,j (A) represent the reduced couplings of A + to an up-and a down-type quark. The short-distance Wilson coefficients C i (µ b ) at the low-energy scale µ b = m b can be obtained through a two-step RG evolution [49,101] where the evolution matricesÛ (µ b , µ W ) andÛ (µ W , µ 0 ) are evaluated in an effective theory with f = 5 and f = 6 quark flavors, respectively. Analytic expressions for these evolution matrices can be found in ref. [101]. The matching conditions for the short-distance Wilson coefficients C i (µ 0 ), including the O(α s ) corrections, at the initial scale µ 0 = m A have been calculated in ref. [51]. Together with the one-loop vertex corrections to the hard kernels T ij (u) calculated in subsection 2.2, this enables us to perform a full NLO analysis of the NP effects in the class-I non-leptonicB 0 (s) → D where the one-loop hard kernels T V LL (u, z) and T V LR (u, z) are given, respectively, by eqs. (2.7) and (2.13). It can be seen that the scale dependence is reduced for the real part, but not for the imaginary part, since the latter vanishes at the LO in α s .  Table 6. The other captions are the same as in Fig. 5.
Specific to the case where the NP four-quark operators are mediated by a heavy colorless charged gauge boson A + , with its mass m A fixed at 1 TeV, we have generally four nonzero effective couplings, λ LL (A), λ LR (A), λ RR (A), and λ RL (A), which might be independent of each other. In order to simplify our analysis and reduce the number of free NP parameters, we shall consider the following three different scenarios: • In scenario I, we consider the case where only one effective coefficient is nonzero in eq. (3.8). Under the individual and combined constraints from the ratios R π , R * π , R ρ , R K , R * K , R K * , R sπ , R sK , R * sπ , as well as R * sK collected in Table 6, we can obtain the allowed ranges for this nonzero effective coefficient, which are shown in Figs at the 1σ level. Such a conclusion is also consistent with the recent observation made in ref. [44], which claims that part of the deviations can be reduced by a left-handed W model through a −10% shift in the b → cūd(s) decay amplitudes. All the other three cases are, however, ruled out already by the combined constraints from the ratios R ( * ) (s)L , even at the 2σ level.
• In scenario II, we consider the case where all the four effective coefficients are nonzero,  Table 6, in the scenarios where the left-and right-handed reduced couplings are symmetric (scenario II, left) and asymmetric (scenario III, right), defined respectively by eqs. (3.13) and (3.14). The other captions are the same as in Fig. 5. but with the additional left-right symmetric assumption on the reduced couplings [102]: for thē B 0 (s) → D + (s) π − andB 0 (s) → D + (s) K − decays. This explains why the ratios R (s)π and R (s)K are insensitive to the NP contributions, as shown in the left panel of Fig. 17.
• In scenario III, we consider instead the case where the left-and right-handed reduced couplings are asymmetric [102]: which implies that, while all the four effective couplings are still nonzero, they satisfy now the relation λ LL (A) = λ RR (A) = −λ LR (A) = −λ RL (A). As shown in the right panel of Fig. 17, such a scenario also fails to provide a simultaneous account for the ten ratios R ( * ) (s)L collected in Table 6, even at the 2σ level. Note that in this case the ratios R (s)π and R (s)K receive no contributions from the NP four-quark operators, which is now due to λ LL (A) = λ RR (A) and λ LR (A) = λ RL (A), resulting in therefore an exact cancellation between the hadronic matrix elements of Q for the decay modes involved.

Colorless charged scalar
Let us now proceed to discuss the case where the NP four-quark operators are generated by a heavy colorless charged scalar H + , with its mass m H fixed also at 1 TeV. The resulting effective weak Hamiltonian for the quark-level b → cūd(s) transitions mediated by such a charged scalar is now given by [51] H scalar and ∆ L,R i,j (H) represent the reduced couplings of H + to an up-and a down-type quark, as defined in Fig. 1. It should be noted that, at the matching scale µ 0 = m H , only the Wilson coefficients C SLL 2 (µ 0 ), C SLR 2 (µ 0 ), C SRR 2 (µ 0 ), and C SRL 2 (µ 0 ) are nonzero at the LO, while all the remaining ones appear firstly at the NLO in α s , with their explicit expressions given already in ref. [51]. To get their values at the low-energy scale µ b = m b , we should also perform a two-step RG evolution as in eq. (3.10), where the analytic formulae for the evolution matricesÛ (µ b , µ W ) andÛ (µ W , µ 0 ) can be found in ref. [101]. This, together with the O(α s ) vertex corrections to the hard kernels T ij (u) presented in subsection 2.2, makes it possible to investigate the NP effects on the class-I non-leptonicB 0 (s) → D  the effective coefficients r K χ a SLL 1 (D + K − ) (normalized by −λ LL (H)) and r K χ a SLR 1 (D + K − ) (normalized by −λ LR (H)), as an example, is shown in Fig. 18, with (3.17) A similar behavior is also observed in the chirality-flipped sectors (SRR and SRL). As in the case for the charged gauge boson, we shall also split the discussions into three different scenarios. Firstly, in scenario-I where only one nonzero effective coefficient is present in eq. (3.15), it is found that all the deviations observed inB 0 (s) → D All the other cases in the presence of only a single effective coefficient are, however, ruled out already by the combined constraints from the ratios R ( * ) (s)L collected in Table 6 at the  Table 6. The other captions are the same as in Fig. 5.
2σ level. As an explicit example, we show in Fig. 20 the individual constraint on the two effective coefficients λ LL (H) and λ LR (H) from the ratios R π , R * π , R ρ , R K , R * K , R K * , R sπ , R sK , R * sπ , as well as R * sK , respectively. Secondly, we show in the left and the right panel of Fig. 21 the individual constraint on the effective coefficient λ LL (H) from the ten ratios R  Table 6. The other captions are the same as in Fig. 5.

Conclusions
In this paper, motivated by the deviations observed between the updated SM predictions and the current experimental measurements of the branching ratios ofB 0 (s) → D ( * )+ (s) L − decays with L ∈ {π, ρ, K ( * ) }, we have investigated possible NP effects in these class-I non-leptonic B-meson decays. In order to facilitate a full NLO analysis, we have also calculated the one-loop vertex corrections to the hadronic matrix elements of the NP fourquark operators involved in these decays, within the QCDF framework.
Firstly, we have performed a model-independent analysis of the effects from twenty linearly independent four-quark operators that can contribute, either directly or through operator mixing, to the quark-level b → cūd(s) transitions. Under the combined constraints from the ten ratios R ( * ) (s)L collected in Table 6, we found that the deviations observed could be well explained at the 1σ level by the NP four-quark operators with γ µ (1−γ 5 )⊗γ µ (1−γ 5 ) structure, and also at the 2σ level by the operators with (1+γ 5 )⊗(1−γ 5 ) and (1+γ 5 )⊗(1+γ 5 ) structures. However, the NP operators with other Dirac structures fail to provide a consistent interpretation, even at the 2σ level. In the case where only a single nonzero NP Wilson coefficient is present in the effective weak Hamiltonian given by eq. (2.1), the resulting allowed ranges for the corresponding NP Wilson coefficients are obtained both at the low-energy scale µ b = m b and at the electroweak scale µ W = m W . In the case where two NP four-quark operators with the same Dirac but different color structures are present in eq. (2.1), with the corresponding two NP Wilson coefficients varied simultane-   Table 7: Allowed ranges for the NP Wilson coefficients Ci(mb) under the individual and combined constraints (last column) from the ten ratios R ( * ) (s)L varied within 1σ and 2σ error bars, respectively. Here "∅" represents an empty set and "R" the set of all real numbers within the plot ranges for Ci(mb).