Two-body non-leptonic heavy-to-heavy decays at NNLO in QCD factorization

We evaluate in the framework of QCD factorization the two-loop vertex corrections to the decays $\bar{B}_{(s)}\to D_{(s)}^{(\ast)+} \, L^-$ and $\Lambda_b \to \Lambda_c^+ \, L^-$, where $L$ is a light meson from the set $\{\pi,\rho,K^{(\ast)},a_1\}$. These decays are paradigms of the QCD factorization approach since only the colour-allowed tree amplitude contributes at leading power. Hence they are sensitive to the size of power corrections once their leading-power perturbative expansion is under control. Here we compute the two-loop ${\cal O}(\alpha_s^2)$ correction to the leading-power hard scattering kernels, and give the results for the convoluted kernels almost completely analytically. Our newly computed contribution amounts to a positive shift of the magnitude of the tree amplitude by $\sim 2$\%. We then perform an extensive phenomenological analysis to NNLO in QCD factorization, using the most recent values for non-perturbative input parameters. Given the fact that the NNLO perturbative correction and updated values for form factors increase the theory prediction for branching ratios, while experimental central values have at the same time decreased, we reanalyze the role and potential size of power corrections by means of appropriately chosen ratios of decay channels.


Introduction
Non-leptonic two-body decays of bottom mesons and baryons are interesting for phenomenological studies of the quark flavour sector of the Standard Model (SM) of particle physics. They yield observables like branching ratios and CP asymmetries that are relevant for studying the CKM mechanism of quark flavour mixing and allow access to the quantities of the unitarity triangle (cf. refs. [1][2][3]).
Oscillations and decays of B-mesons received considerable attention for the first time in the 1980s and 90s when the experiments ARGUS at DESY and CLEO at Cornell started to collect a lot of statistics. In the last decade, non-leptonic two-body B (s)decays have been extensively measured at the asymmetric e + e − colliders (B-factories) at SLAC and KEK, but also in hadronic environments such as the Tevatron, and the results obtained by the Babar, Belle, D0 and CDF collaborations have reached a high level of precision (see, e.g. [4]). In recent years the LHCb experiment at the LHC at CERN has become the main player as far as experimental physics of the bottom quark is concerned. A large data set on bottom mesons and baryons has been accumulated, and results related to non-leptonic decays have been published (cf. [5,6]) and further analyses are ongoing. In the near future also Belle II will contribute significantly to further improve the measurements [7].
With the plethora of precise experimental data on non-leptonic decays at hand, theoretical predictions at the same level of accuracy are very much desired. However, the theoretical description of non-leptonic two-body B (s) decays is notoriously complicated. A straightforward computation of the hadronic matrix elements which describe the weak transition is not feasible due to the presence of the strong interaction in the purely hadronic initial and final states. This circumstance entails QCD effects from many different scales which are, moreover, largely separated. In a first approach, known as naïve factorization, the hadronic transition matrix elements were factorized into a product of a form factor and a decay constant [8]. Subsequent studies built on flavour symmetries of the light quarks [9] and on factorization frameworks such as perturbative QCD (pQCD) [10,11] and QCD factorization (QCDF) [12][13][14], to mention the most prominent ones. Certain combinations of these approaches can also be found (see e.g. [15]).
In the present work we adopt the QCDF framework and consider non-leptonic heavyto-heavy transitions, which at the quark-level are mediated by the weak decay b → cūd(s), where we treat the bottom and the charm quark as massive and the light quarks as massless. Performing an expansion of the amplitude in powers of Λ QCD /m b , where Λ QCD is the typical hadronic scale, a systematic separation of QCD effects from different scales can be achieved and corrections to naïve factorization be systematically included. Taking the decayB → D + L − as an example, the transition amplitude in the heavy-mass limit is then given by [13] where the local four-fermion operators Q i describe the underlying weak decay. The F B→D j form factors and the light-cone distribution amplitude (LCDA) Φ L of the light meson contain long-distance effects and can be obtained from non-perturbative methods like QCD sum rules and lattice QCD. The hard-scattering kernels T ij , on the other hand, only receive contributions from scales of O(m b ) and are accessible in a perturbative expansion in the strong coupling α s . After the convolution over the momentum fraction u of the valence quark inside the light meson, they yield a perturbative contribution to the topological tree amplitude a 1 (D + L − ). Taking the decayB → D + π − as a specific example, the latter is defined via [13] The leading-power hard-scattering kernels have been known to next-to-leading order (NLO) accuracy for more than a decade for both heavy-to-light [12,14,16] and heavy-toheavy [13] decays. In the latter case, expanding the LCDA in Gegenbauer moments up to the first moment α L 1 , the topological tree amplitude a 1 to NLO reads [13] For the light meson being π or ρ we have α π(ρ) 1 = 0 and for the kaon |α K 1 | < 1 is assumed [13]. With this mild dependence on the light meson LCDA we encounter a quasi-universal value |a 1 | 1.05 for heavy-to-heavy decays in QCDF to NLO accuracy. A quasi-universality was also found upon extracting a 1 from experimental data [17]. However, the favoured central value |a 1 | 0.95 for the decaysB → D ( * )+ L − (L = π, K), with errors in the individual channels at the 10 -20% level, is considerably lower.
In recent years next-to-next-to-leading order (NNLO) corrections to heavy-to-light decays have become available [18][19][20][21][22][23][24][25][26][27][28][29][30][31], and besides the prospects of increasing precision on the experimental side, there is multiple motivation to go beyond NLO in heavy-toheavy transitions as well: First, the NLO correction is small since it is proportional to a small Wilson coefficient and, in addition, is colour-suppressed. At NNLO the colour suppression gets lifted and the large Wilson coefficient re-enters, and therefore the NNLO correction could be comparable in size to the NLO term. Moreover, it is interesting to see whether the quasi-universality of a 1 persists at NNLO. At leading power the de-caysB (s) → D ( * )+ (s) L − receive only vertex corrections to the colour-allowed tree topology. Interactions with the spectator quark as well as the weak annihilation topology are powersuppressed [13] and there are neither contributions from penguin operators nor is there a colour-suppressed tree topology. Therefore, a precise knowledge of the colour-allowed tree amplitude a 1 allows to reliably estimate the size of power corrections to eq. (1) by comparison to experimental data, and at the same time provides a test of the QCDF framework. This requires that the perturbative expansion of the hard scattering kernel is under control, and that also the uncertainties of the non-perturbative input parameters (form factors, decay constants, LCDAs) can be minimized. In the present work we therefore calculate the two-loop vertex correction to the leading-power hard scattering kernels in the framework of QCDF. Parts of the computational procedure were already presented in [32,33]. Here, we give the full result of the technically challenging two-loop calculation. Besides, we present an updated phenomenological analysis ofB (s) → D ( * )+ (s) L − decays, with a light meson from the set L = {π, ρ, K ( * ) , a 1 } 1 , using the most recent values for non-perturbative input parameters (for another recent analysis, see [34]).
Recently, non-leptonic Λ b decays have received considerable attention as well. Data on Λ b → Λ + c L − with L being π or K [35] and on baryonic form factors have become available [36]. Therefore, we extend our study to these decays. Factorization has not yet been systematically established for baryonic decays, but was discussed in ref. [37]. As a systematic derivation of the baryonic factorization formula is beyond the scope of this work we adopt the factorization formula eq. (4) of ref. [37], with appropriate modifications to take perturbative corrections into account.
This article is organized as follows: In section 2 we present our theoretical framework by specifying our operator basis in the effective weak Hamiltonian. Subsequently, we derive the master formulas for the hard scattering kernels by performing a matching onto Soft-Collinear Effective Theory. In section 3 we discuss the calculation of the two-loop Feynman diagrams and specify the input to the master formulas. The analytical results of the hard scattering kernels after the convolution with the LCDAs are presented in section 4. In section 5 we give the formulas for converting from the pole to the MS scheme for the band c-quark masses. We present the results of our extensive phenomenological analysis in section 6, and conclude in section 7.
2 Theoretical framework

Five-flavour theory
We work in the effective five-flavour theory where the top quark, the heavy gauge bosons W ± , Z 0 and the Higgs boson are integrated out and their effects are absorbed into shortdistance Wilson coefficients. The decaysB (s) → D ( * )+ (s) L − and Λ b → Λ + c L − are mediated at parton level by a b → cūd(s) transition for L = π, ρ, a 1 (K, K * ). The corresponding QCD amplitude is computed in the framework of the effective weak Hamiltonian [14,38], which for the problem at hand simply reads We restrict our notation to the case of a b → cūd transition. The expressions for a strange quark in the final state are obtained by obvious replacements. The local current-current operators in the Chetyrkin-Misiak-Münz (CMM) basis [39,40] read Figure 1: The tree-level Feynman diagram for the b → cūd transition in full (five flavour) QCD: the black square represents the vertex of the effective weak interaction. The momenta q 4 and q 3 belong to the quark lines with masses m b and m c , respectively, and q 1 + q 2 = q is the momentum of the light meson. All momenta are taken to be incoming.
where Q 1 and Q 2 are referred to as colour-octet and colour-singlet operator, respectively. The use of the CMM basis allows for a consistent treatment of γ 5 in the naïve dimensional regularization scheme with fully anti-commuting γ 5 . Moreover, as the computation will be performed in dimensional regularization, we have to augment our physical operators Q 1,2 by a set of evanescent operators, for which we adopt the convention [41,42] These unphysical operators vanish in D = 4 dimensions but contribute if D = 4 since they mix under renormalization with the physical operators. At two-loop accuracy the set of operators (5) -(10) closes under renormalization.

Matching onto SCET and master formulas
We construct the master formulas for the hard scattering kernels by performing a matching from the effective weak Hamiltonian 2 onto Soft-Collinear Effective Theory (SCET) with three light flavours. The procedure follows similar lines than the derivation of the master formulas for the hard kernels in heavy-to-light transitions [28]. The kinematics of the b → cūd transition is shown in the tree-level Feynman diagram depicted in figure 1. The b and the c quark are considered to be massive and carry momenta q 4 and q 3 , respectively. The massless d andū quarks share the momentum q with q 1 = uq and q 2 = (1 − u)q ≡ūq, where u ∈ [0, 1] is the momentum fraction of the valence quark inside the light meson. All external momenta are taken to be incoming and are subject to the on-shell constraints q 2 4 = m 2 b , q 2 3 = m 2 c , and q 2 = 0. We consider a reference frame in which the b quark within the B meson moves with momentum q b = m b v +k, where k is a residual momentum of order of the typical hadronic scale Λ QCD , and v is the velocity of the B meson. The b quark can then be described by a heavy-quark field h v which satisfies the equation of motion / vh v = h v . We further choose a reference frame such that the energetic light meson moves in the light-cone direction n + . The light-like vectors n + and n − = 2v − n + then fulfill the constraints n 2 ± = 0 and n + n − = 2. As the quark and the anti-quark in the light meson nearly move in the same direction we can describe them by the same type of collinear SCET field χ, which satisfies the equations of motion / n + χ = 0 andχ/ n + = 0. In the derivation of the factorization formula (1) the power counting m c /m b ∼ O(1) was adopted. Hence, we treat the charm quark as a heavy quark and consequently describe it -in analogy to the b quark -by another heavy-quark field h v with velocity v and equation of motion / v h v = h v . The amplitudes in full QCD and in SCET are made equal by adjusting the corresponding hard coefficients at the matching scale. We express the renormalized matrix elements of the QCD operators (5) and (6) as a linear combination of a basis of SCET operators, where H ia and H ia are the matching coefficients. The basis of SCET operators is given by Here, the perpendicular component of a Dirac matrix is defined by Moreover, we have omitted the Wilson lines which render the non-local light currents χ(tn − )[. . . ]χ(0) gauge invariant. One therefore has to keep in mind that the coefficients H ia are also functions of the variable t, and the products H in eq. (11) are in fact convolutions. We also remark that the SCET operator basis is chosen such that all operators with index a > 1 are evanescent, and we have the two physical SCET operators O 1 and O 1 . The operators (12) - (14) have the same structure as those in [28] for heavy-to-light transitions, but with a heavy-quark field h v instead of an anti-collinear SCET field ξ in direction n − . For heavy-to-heavy transitions this set of operators has to be extended by those in (15) - (17) which have a different chirality structure, to take into account the non-vanishing mass of the charm quark. For technical details on the operators see [19,43].
We first consider the expansion of the left-hand side of eq. (11) in terms of on-shell QCD amplitudes. The expression for the renormalized matrix elements reads Here, a sum over a = 1, 2, 3 is understood, and α s is the MS strong coupling constant with five active flavours. The index i = 1, 2 denotes the physical operators from (5) and (6) only, whereas j includes physical as well as evanescent operators from (5) -(10), hence j = 1, . . . , 6. The A (l) are the bare l-loop on-shell amplitudes and A * (1) (A * * (1) ) is the oneloop bare on-shell amplitude with a b (c) quark mass insertion on the heavy b (c) line. The primed amplitudes are defined analogously. The renormalization factors Z ij , Z ext and Z α stem from operator renormalization, wave-function renormalization of all external legs and coupling renormalization, respectively. They are defined in a perturbative expansion The operator renormalization is performed in the MS scheme, whereas for the mass and the wave-function renormalization the on-shell scheme is applied. Renormalized matrix elements of evanescent operators vanish also beyond tree level. Nevertheless, these operators cannot be neglected right from the beginning as they yield physical contributions to the products Z ja . Similarly, we can write down the expression for the renormalized matrix elements of the SCET operators that enter the right-hand side of eq. (11) and obtain Here, a = 1, 2, 3 and a sum over b = 1, 2, 3 is understood. The MS strong coupling constantα s has three light flavours and M (l) are the bare l-loop SCET amplitudes. The Y α are the l-loop wave-function, operator and coupling renormalization constants, respectively. They are defined in a perturbative expansion analogous to eq. (20) except that the strong coupling has only three light flavours. The corresponding expression for the primed operators from eqs. Eq. (21) can be simplified to a large extent. In dimensional regularization the on-shell renormalization constants Y ext are equal to unity. Moreover, the bare on-shell amplitudes only contain scaleless integrals, which vanish in dimensional regularization. We thus arrive at the following simplified expression of eq. (21) which for the primed operators takes a similar form.
For relating the matching coefficients H ia and H ia in eq. (11) to the hard scattering kernels we introduce two factorized QCD operators which are by definition the products of the two currents in brackets. The upper sign corresponds to the un-primed operator. The renormalized operators Q ( )QCD are then matched onto the renormalized SCET operators O 1 and O 1 by adjusting the corresponding hard coefficients. This can be done separately for the light-to-light and heavy-to-heavy currents. For the renormalized light-to-light current we make the ansatz The heavy-to-heavy currents with different chiralities mix in the matching. Thus, we make the ansatz Since these equations are symmetric under interchanging P L ↔ P R we have C LL Since by construction Q QCD and Q QCD factorize into a light-to-light and a heavy-toheavy current, the matrix element of these operators is the product of an LCDA and the full QCD form factor with the corresponding helicity structure. We now consider the two hard scattering kernelsT i andT i that are defined by the following expression Comparing eqs. (11) and (29),T i andT i can be related to the matching coefficients as follows Plugging in the matching coefficients as expansions in the five-flavour coupling α s , the matrix can be inverted order-by-order in α s . We remark that Cq q = 1 + O(α 2 s ), i.e. it receives a correction at two loops only since at one loop only scaleless integrals contribute. The explicit one-loop expressions for the heavy-to-heavy coefficients will be derived in section 3. For the diagonal coefficients we have C D F F = 1 + O(α s ). In contrast, the nondiagonal matching coefficients C ND F F that induce the chirality mixing only arise beyond tree level, C ND F F = O(α s ). Putting everything together, the master formulas for the hard scattering kernels read The expression for the primed kernelsT i is given by eq. (31) with the replacement A ↔ A , H ↔ H ,T ↔T . Note that the quantities H (l) , A (l) and the hard kernelsT (l) depend on the quark mass ratio z c = m 2 c /m 2 b and the momentum fraction u of the quark inside the light meson (as do the corresponding primed quantities). Whenever they appear alongside a renormalization factor Y (l) such as H b1 we must keep in mind that these expressions must be interpreted as a convolution product in eq. (31) are termed "non-factorizable". At one-loop the corresponding amplitudes are given by all Feynman diagrams with one gluon connecting the heavy and the light current. The one-loop Feynman diagrams where the gluon is attached solely to either the light or the heavy current are part of the LCDA and the form factor, respectively. The Feynman diagrams contributing to A can be found in figures 15 and 16 of ref. [13] and in addition include the one-loop self-energy insertions to the "non-factorizable" one-loop amplitudes. A sample of two-loop diagrams is shown in is technically the most challenging contribution to the two-loop kernels. Therefore, we briefly describe their evaluation in the next section and, moreover, specify the remaining input to eq. (31). The final expression of the hard scattering kernels must be free of ultraviolet and infrared divergences. We comment on this at the end of the next section.
Finally, we remark that eq. (31) has a structure similar to the corresponding expressions for the two-loop hard scattering kernel in the right-insertion contribution to the decay B → ππ, which is given in eq. (24) in [28]. The main difference is three-fold: First, we find two contributionsT andT to the hard scattering kernel as a result of the extended operator basis. Second, we encounter the off-diagonal element C due to the mixing of the heavy-to-heavy currents with different chirality structures. Finally, we have a mass counterterm for the massive charm quark in eq. (31).

Computational details 3.1 Technical aspects of the two-loop computation
We work in dimensional regularization with D = 4 − 2 and expand the amplitudes in the parameter . The Feynman diagrams contributing to the bare two-loop amplitude A (2)nf then contain up to 1/ 4 poles stemming from ultraviolet (UV) and infrared (IR) regions. We calculated them by applying commonly-used multi-loop techniques, including a new method for evaluating the master integrals. The procedure goes as follows: First, we decompose all tensor integrals into scalar ones by applying the Passarino-Veltman decomposition [44]. We then perform the reduction of the Dirac structures to the SCET operator basis given in eqs. (12) - (17) in Mathematica by using simple algebraic transformations. The number of remaining scalar two-loop integrals exceeds several thousands and can be simplified by using the Laporta algorithm [45,46], which is based on integration-by-parts identities [47]. Here, we apply the implementations AIR [48] (in Maple) and FIRE [49] (in Mathematica) of this algorithm to reduce the large number of integrals to a small set of master integrals. Many of the latter are already known from several B → ππ calculations [25,26,28]. In addition, we find 23 yet unknown two-loop master integrals. Since most of them depend on two scales (the momentum fraction u and the quark-mass ratio z c = m 2 c /m 2 b ), an analytic solution by common techniques is hardly feasible. We therefore evaluate them by applying the approach of differential equations in a canonical basis recently advocated in [50]. The solution is given by iterated integrals and falls into the class of Goncharov polylogarithms [51]. We obtain analytic results for all 23 master integrals. Details on their calculation and the result of all master integrals can be found in [33].

Input to the master formulas
Here we give the explicit expressions for the renormalization factors and matching coefficients that enter the master formula, and in the end comment on the cancellation of the poles in once all pieces of the master formula are plugged in.
The operator renormalization factors Z ij of the effective weak Hamiltonian were calculated to two-loop accuracy in the MS scheme in [41,42]. The explicit one-and two-loop expressions read Here, n f = 5 is the total number of active quark flavours and T f = 1/2. The row index of these matrices corresponds to (Q 1 , Q 2 , E 1 , E 2 , E 1 , E 2 ) and the column index to (Q 1 , Q 2 ). The strong coupling constant is renormalized in the MS scheme as well, whereas the renormalization of the masses and the wave-functions is performed in the on-shell scheme. The corresponding renormalization factors are well known and shall not be repeated here.
In eq. (31) we further encounter the SCET operator renormalization factor Y 11 that can be split into the following two parts Here, Z Jh and Z BL are the renormalization factors for the HQET heavy-to-heavy and the SCET light-to-light current, respectively. Since one collinear sector in SCET is equivalent to full QCD, the renormalization constant Z BL coincides with the ERBL kernel in QCD [52,53]. We take Z BL from [54], which for pseudoscalar and longitudinally polarized vector mesons reads The plus-distribution for symmetric kernels f is defined as follows, The renormalization factor Z Jh can be obtained in a matching of the heavy-to-heavy In this process also the matching coefficients C F F can be determined. Beyond tree-level the QCD current also mixes into the chirality-flipped HQET currenth v / n + (1 + γ 5 )h v . Hence, we make the following ansatz for the renormalized currents where we have already made use of the fact that both equations are symmetric under interchanging P L ↔ P R . The renormalization factor Z Jh is defined via the on-shell oneloop matrix element of the HQET currents The one-loop renormalized matrix elements of the QCD currents can be calculated straightforwardly. Inserting their explicit expressions in eq. (37) we can identify Z (1) Jh as the pole term in , that is which correctly reproduces the IR behaviour of QCD currents in the effective theory. The C F F , on the other hand, are given by the coefficients that are finite in . Their explicit As a last step the contribution of the b1 in eq. (31) needs to be further specified (the primed quantities are obtained by obvious substitutions). We find that only H (1) ib with i = 1 and b = 2 yields a non-vanishing contribution. A straightforward calculation yields H has already been used in the NNLO calculation of the vertex corrections to the decay B → ππ and is given in eq. (45) of [25]. Its explicit expression reads With this we have specified all input to the master formulas and are now ready to produce an expression for the hard scattering kernels. The final expressions for the hard scattering kernels are free of poles in , even though most of the individual terms in eq. (31) contain divergences. At the one-loop level we checked the cancellation of all poles analytically. We find that our expressions for the finite pieces of the one-loop kernels agree with the results given in eqs. (89) and (90) in [13] 3 . Some of the one-loop quantities that enter the two-loop master formula (last equation in (31)) have to be evaluated to higher orders in the -expansion since they multiply poles in contained in the renormalization factors. We checked that in the limit m c → 0 the O( ) piece of the one-loop hard scattering kernel coincides with the one used in [28].
At two loops we could check the pole cancellation numerically to an accuracy of 1 × 10 −10 or better for 12 different points in the u-z c plane. To this end, we evaluate the Goncharov polylogarithms and the harmonic polylogarithms [55] that are contained in A (2)nf numerically with the C++ routine GiNaC [56] and the Mathematica program HPL [57,58], respectively. The explicit results for the two-loop hard scattering kernels are lengthy, not very illuminating, and enter the physical quantities only after convolution with the LCDAs. For these reasons we refrain from presenting them explicitly here, but they can be obtained from the authors upon request. However, after the convolution of the hard scattering kernels with an expansion of the LCDAs in terms of Gegenbauer polynomials up to the second moment, the expressions simplify considerably and we can express the result almost entirely in terms of harmonic polylogarithms. At this level we convolute also the pole terms in and checked that for the convoluted kernels all poles cancel analytically. We give the corresponding finite parts in the next section.

Convoluted kernels
The light meson LCDAs are expanded in a basis of Gegenbauer polynomials C Following [13] we assume that the leading-twist LCDA is close to its asymptotic form Φ L (u, µ) = 6u(1 − u) and truncate the expansion after the second moment. The first two Gegenbauer polynomials read C . The Gegenbauer polynomials are eigenfunctions of the one-loop renormalized ERBL-kernel [59] and thus the Gegenbauer moments are multiplicatively renormalizable to leading-logarithmic (LL) accuracy [59]. The next-to-leading logarithmic (NLL) evolution was derived in [59][60][61].
The result for the hard scattering kernels after the convolution with the LCDA can be written as follows with α L 0 (µ) ≡ 1. At tree-level we obtain V In the following we use the abbreviations L ≡ log(µ 2 /m 2 b ) and H a (z c ) ≡ H a for the harmonic polylogarithms of argument z c . The one-loop results for the convoluted colouroctet kernels then read The one-loop colour-singlet kernels vanish as the corresponding colour factors are zero, 2k (µ) = 0 , for k = 0, 1, 2 .
At two loops the result is rather lengthy. Here, we only present the full result for the µ-dependent part which governs the scale dependence. For the µ-independent part we provide a fitted function in z c that agrees with the original result at the per mill level in the range of physical values 0.05 ≤ z c ≤ 0.2. The full result is attached in electronic form to the arXiv submission of the present work. For the convoluted colour-octet kernels we obtain The result for the convoluted colour-singlet kernels takes the form Finally, we checked with the full result (without interpolation in z c ) that in the limit m c → 0 eq. (44) with µ = m b coincides with the result for the vertex corrections to the colour-allowed tree topology of the decay B → ππ given in eq. (48) of [28] .

Conversion from the pole to the MS scheme
The convoluted kernels in eqs. (44) and (45) are given in the pole scheme, where m c and m b appearing in L ≡ log(µ 2 /m 2 b ) and z c = m 2 c /m 2 b denote the pole quark masses, and the renormalization scale µ ∼ m b . In order to discuss the scheme dependence of the convoluted kernels, we also give the results in the MS scheme for the quark masses. Since the LO kernels are constant and the NLO colour-singlet kernels vanish, the conversion from the pole to the MS scheme will only affect the NNLO colour-octet kernels V (2) 1k and V (2) 1k . To this end, using the one-loop relation between pole-and MS-quark mass, we find that the corresponding convoluted kernels in the MS scheme are obtained via the relation where now L ≡ log(µ 2 /m 2 b (µ)) and z c = m 2 c (µ)/m 2 b (µ), with µ ∼ m b (m b ). The tree-level and one-loop kernels will have the same functional dependence as in the pole scheme, but now depend on the above new abbreviations in the MS scheme. At two loops we explicitly give the terms that have to be added,

Phenomenological applications
In this section we perform an extensive phenomenological analysis ofB (s) → D ( * )+ (s) L − and Λ b → Λ + c L − decays in QCDF. Like before, L is a light meson from the set {π, ρ, K ( * ) , a 1 }. We take into account the expressions through to NNLO for the hard scattering kernels, and the most recent values for non-perturbative input parameters, which we specify below. We analyze the impact of the NNLO correction on the topological tree amplitude a 1 (D ( * )+ L − ), and subsequently predict the branching ratios for the mesonic decays. Afterwards, we perform tests of QCD factorization by considering suitably chosen ratios of non-leptonic to either semi-leptonic or non-leptonic channels. Finally, we give the theoretical predictions for baryonic decays.

Input parameters
Here we collect in Table 1 the theoretical input parameters entering our numerical analysis throughout this paper. They include the SM parameters such as the CKM matrix elements, quark masses, and the strong coupling constant, as well as the hadronic parameters such as meson decay constants, transition form factors, and the Gegenbauer moments of light mesons. Three-loop running is used for α s throughout this paper. Furthermore, we use a two-loop relation between pole and MS mass to convert the top-quark pole mass m pole t to the scale-invariant mass m t (m t ) [75]. For the B → D ( * ) transition form factors, we adopt the parameterization proposed by Caprini, Lellouch, and Neubert (CLN) [76], with the relevant parameters extracted from exclusive semileptonic b → c ν decays [66]. For the B s → D ( * ) s transition form factors, on the other hand, we use the results obtained by QCD sum-rule techniques, assuming a polar dependence on q 2 that is dominated by the nearest resonance [69,77]. However, to discuss the SU(3)-breaking effects in the form-factor and decay-constant ratios, we adopt the most recent lattice QCD results for the ratios [70,78] F Bs→Ds 0 (m 2 π ) F B→D 0 (m 2 π ) = 1.054 ± 0.047 stat. ± 0.017 syst. , Neither of the form-factor ratios shows significant deviation from the U-spin symmetry. For the Λ b → Λ c transition form factors, we use the most recent high-precision lattice QCD calculation with 2 + 1 dynamical flavours [36]. Here the q 2 dependence of the form factors is parameterized in a simplified z expansion [79], modified to account for pionmass and lattice-spacing dependence. All relevant formulas and input data can be found in eq. (79) and Tables VII -IX of [36]. Following the procedure recommended in [36], we calculate the central value, statistical uncertainty, and total systematic uncertainty of Table 1: Summary of theoretical input parameters. The Gegenbauer moments of light mesons are evaluated at µ = 1 GeV.
The decay constants f π and f K are averaged over the two-flavour lattice QCD results [70], while f ρ and f K * are determined from experiments [71]. The light-meson Gegenbauer moments are determined by the QCD sum rule approach [71,72] and the lattice QCD calculation [74]. For the hadronic inputs of the axial-vector meson a 1 (1260), we use the results presented in ref. [73]. It is noted that the Gegenbauer moments are evaluated at µ = 1 GeV, and are evolved to the characteristic scale µ ∼ m b [59][60][61]. We use LL running of the Gegenbauer moments for the tree-level and the one-loop amplitude, but NLL running in the two-loop amplitude. Moreover, the running of the Gegenbauer moments is performed in the four-flavour scheme.

Predictions for a 1 (D ( * )+ L − )
We are now in the position to perform a numerical analysis of the coefficients a 1 (D ( * )+ L − ) according to the expressions into which eqs. (44) and (45) have to be inserted. Using the NNLO Wilson coefficients C i (µ) in the CMM basis [42], together with the input parameters collected in Table 1 where the number without bracket is the LO contribution, which has no imaginary part, and the following two numbers are the NLO and NNLO terms, respectively. The total errors comprise the uncertainties, added in quadrature, from the variation of the scales µ ∈ [m b /2, 2m b ] and µ 0 ∈ [m W /2, 2m W ], the quark masses, the Gegenbauer moments, and α s (m Z ). Unless stated otherwise, the numbers given here and below are obtained with the band c-quark masses renormalized in the pole scheme, which is set as our default scheme. It is observed that both the NLO and NNLO contributions add always constructively to the LO result. We also observe that the new two-loop correction is quite small in the real, but rather large in the imaginary part. It amounts to approximately 60% (2%) of the total imaginary (real) part of a 1 (D + K − ). We emphasize that the sizable NNLO correction to the imaginary part does not indicate a breakdown of the perturbative expansion, but is due to the fact that the imaginary part vanishes at LO, and its NLO term is colour suppressed and proportional to the small Wilson coefficient Due to the truncation of the perturbative expansion, the obtained values in eq. (58) depend on the renormalization scale µ, which is usually considered as a measure of the accuracy of the approximation at a given order in the perturbative expansion. This is shown in figure 4 for a 1 (D + K − ) up to NNLO, where results both in the pole (blue) and in the MS (red) scheme for band c-quark masses are given. We observe a pronounced stabilization of the scale dependence for the real part, but not for the imaginary part. This is again explained by the fact that the imaginary part vanishes at LO. It is also observed that the dependence on the band c-quark mass scheme is quite small, especially for the real part. We finally remark that also within a given quark-mass scheme the dependence of a 1 (D + K − ) on the value of z c is minor. The dependence of a 1 (D + K − ) on the second Gegenbauer moment is small, too.
It is also interesting to mention that, even up to NNLO, the coefficients a 1 (D ( * )+ L − ) are quasi-universal, with very small process-dependent non-factorizable corrections, a fact that was observed already at NLO in ref. [13]. This is clearly seen from the following numerical results for different final states:

Predictions for class-I decays
It is generally believed that the factorization theorem is well established in class-I decays of the formB (s) → D mesons [13,80]. We now present in Table 2 our predictions for the branching ratios of these decays through to NNLO. The explicit formulas for the branching ratios can be found in [13] and shall not be repeated here. The experimental data is taken from the Particle Data Group (PDG) [62] and/or the Heavy Flavor Averaging Group (HFAG) [66]. For the vector and axial-vector final states, the results refer to the longitudinal polarization amplitudes only, with the longitudinal polarization fractions taken from [81] forB d → D * + ρ − and [82] forB s → D * + s ρ − , respectively. From Table 2, one can see that our predictions for the branching ratios of these decays generally come out higher than the experimental data, especially forB d → D ( * )+ π − and B d → D ( * )+ ρ − decays, where the difference in central values is at the 20 -30% level. Taking into account the uncertainties, the deviation is at the level of 2 -3σ. Compared to ref. [13], which found at NLO rather good agreement between theory and experiment, essentially three things have changed: First, using the latest extraction from [66][67][68] our numerical values for the form factors are about 10% larger than the ones used in [13]. Second, the NNLO corrections add another positive shift of 2 -3% on the amplitude level. Third, the experimental central values have slightly decreased since the analysis of [13]. All three effects shift theory and experiment further apart.
Given the fact that the results show rough agreement within errors forB d → D ( * )+ K ( * )− decays, which receive only contributions from colour-allowed tree topologies, this may indicate a non-negligible impact from the W -exchange topologies appearing only in B d → D ( * )+ π − andB d → D ( * )+ ρ − decays. ForB s decays, on the other hand, since the B s → D ( * ) s transition form factors have so far received only little theoretical attention [69,[83][84][85][86][87][88], especially by the lattice QCD community [78,89], our theoretical predictions are still plagued by larger uncertainties due to these hadronic parameters. Table 2: CP-averaged branching ratios (in units of 10 −3 for b → cūd and 10 −4 for b → cūs transitions) ofB (s) → D ( * )+ (s) L − decays. The vector-and axial-vector final states refer to the longitudinal polarization amplitudes only. The theoretical errors shown correspond to the uncertainties due to renormalization scales µ and µ 0 , the CKM as well as the hadronic parameters, added in quadrature. The experimental data is taken from refs. [62,66,81,82]

Test of factorization
To further test the factorization hypothesis in class-I decays of B-mesons into heavylight final states, as well as to probe the non-factorizable corrections to the coefficients a 1 (D ( * )+ L − ), we now consider either ratios of non-leptonic to semi-leptonic decay rates [13,17,90,91], or ratios of two non-leptonic decay rates [13,91], both of which are essentially free of CKM and hadronic uncertainties. As suggested firstly by Bjorken [90], a particularly clean and direct method to test the factorization hypothesis is provided by dividing the non-leptonicB d → D ( * )+ L − decay rates by the corresponding differential semi-leptonicB d → D ( * )+ −ν decay 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 transferring to the lepton pair. In this way, the coefficients a 1 (D ( * )+ L − ) can be extracted directly from experimental data through the relation [13,91] where V ij is, depending on the constituent quark content of the meson L, the appropriate CKM matrix element. With the light lepton mass neglected, X L = X * L = 1 for a vector or axial-vector meson, whereas for a pseudoscalar X ( * ) L deviates from unity only by calculable terms of order m 2 L /m 2 B , which are numerically below the percent level; explicit expressions for X ( * ) L can be found, for example, in ref. [91]. To get the differential semi-leptonic decay rates at q 2 = m 2 L in eq. (60), we use the CLN parameterization for the B → D ( * ) transition form factors [76], with the relevant parameters summarized in Table 1. Explicitly, we get numerically (in units of 10 −3 GeV −2 ps −1 ) Together with the data on the branching ratios of non-leptonic decays given in Table 2, we arrive at the experimental values for |a 1 (D ( * )+ L − )| collected in Table 3, where, for comparison, our theoretical predictions at different orders are also shown. From Table 3, one can see clearly that our theoretical predictions based on the QCDF approach result in an essentially universal value of |a 1 (D ( * )+ L − )| 1.07 (1.05) at NNLO (NLO), being consistently higher than the central values favoured by the current experimental data. The deviation is again at the level of 2 -3σ. Similar results were obtained in [17], yet without inclusion of the NNLO correction. It would be, therefore, very encouraging to determine directly the ratios of non-leptonic and semi-leptonic decay rates at current and future experiments such as LHCb and Belle II. Compared to the NLO analysis in [13], where theory predictions for |a 1 (D ( * )+ L − )| were found to be Table 3: Theoretical predictions for |a 1 (D ( * )+ L − )| at different orders in perturbation theory. For comparison, the coefficients |a 1 (D ( * )+ L − )| determined from current data are shown in the last column. The experimental errors are estimated by adding the uncertainties of the non-leptonic branching ratios and the semi-leptonic decay rates in quadrature, while the uncertainties from the decay constants are not taken into account. in agreement with the values extracted from experiment, together with the conclusion that there was no hint for sizable power corrections, the situation has changed, owing to increased values in the theory predictions and, at the same time, decreased experimental values (see also the discussion in section 6.3). We will come back to this point below.
We now turn to discuss the ratios of non-leptonicB (s) → D ( * )+ (s) L − decay rates, following the notations used in refs. [13,91]. As a quasi-universal |a 1 (D ( * )+ L − )| is predicted in the QCDF approach, these ratios could be used to test the factorization hypothesis, as well as the SU(3) relations in B-meson decays into heavy-light final states [17]. Our results of such an analysis are presented in Table 4, where the experimental data is obtained using the corresponding branching fractions collected in Table 2.
From Table 4, one can see that, within the errors, our theoretical predictions are generally consistent with the current experimental data, indicating therefore no evidence for any significant deviation from the factorization hypothesis for these class-I B-meson decays into heavy-light final states. The last two ratios shown in Table 4 could also be used to determine the ratio of fragmentation functions f d /f s , a key quantity for precise measurements of absolute B s -meson decay rates at hadron colliders [17,92].
One possible interpretation of our findings that, on the one hand, the non-leptonic to semi-leptonic ratios come out larger in theory compared to experiment, and on the other hand the non-leptonic ratios in general agree with experiment, might be non-negligible power corrections which could be negative in sign and 10 -15% in size on the amplitude Table 4: Predictions for the ratios of non-leptonicB (s) → D ( * )+ (s) L − decay rates at different orders. The experimental data is obtained using the corresponding branching fractions collected in Table 2.
level. They would render the factorization test via non-leptonic to semi-leptonic ratios better, and at the same time could cancel out in the non-leptonic ratios, especially if they were of a certain universality. The size of power corrections stemming from spectator scattering and weak annihilation was roughly estimated in section 6.5 of [13]. Depending on the phases of the integrals over the D-meson wave function and on the value of the first inverse moment λ B of the B-meson distribution amplitude, these two contributions could in principle interfere constructively, and in this case their total effect could indeed add up to −10% in the amplitude. Another possibility which has essentially the same effect would be to reduce the values of |V cb | times the form factors by ∼ 10%. This option seems attractive in view of the fact that those non-leptonic ratios in Table 4 in which |V cb | and the form factors cancel out are in very good agreement with experiment. On the other hand, the semi-leptonic rate is measured very precisely and the current form factors times |V cb | are extracted by HFAG [66] from a global fit to all available data, whose result we quote in Table 1. Hence they are optimized to describe the shape of the semi-leptonic rate and therefore should be trustworthy. One could even conclude from this that the experimental extraction of  3.3 ± 1.2 |a 1 (D ( * )+ L − )| from eq. (60) is independent of the product of |V cb | and the form factor. We emphasize that without a rigorous treatment of power corrections in the QCDF approach nothing more can be said at the present stage. In any case, the QCDF approach per se is not invalidated.

Predictions for
While the Λ b baryons are not produced at an e + e − B-factory, they account for about 20% of the b-hadrons produced at the LHC [93]. Remarkably, the number of Λ b baryons produced is comparable to the number of B u or B d mesons, and is significantly higher than the number of B s mesons. Due to the half-integer-spin of Λ b , its decays provide complementary information compared to the corresponding mesonic ones. Therefore, this may open up a new field for flavour physics. For a review, see e.g. refs. [94][95][96][97].
Here we study the two-body non-leptonic Λ b → Λ + c L − decays, for which the factorization assumption is believed to be reliable [37,[98][99][100][101]. As demonstrated especially in ref. [37], the proof of factorization at leading order in Λ QCD /m b,c for these decays follows closely that forB d → D ( * )+ π − [80]. These decays provide, therefore, a testing ground for different QCD models and factorization assumptions used in B-meson case. It is straightforward to generalize the expressions in [37] to take radiative corrections through to NNLO into account.
Using the most recent lattice QCD results for Λ b → Λ c transition form factors [36], we present in Table 5 our predictions for the branching fractions of Λ b → Λ + c L − decays, as well as some ratios between them, where the experimental data is taken from HFAG [66]. From Table 5, one can see that, contrary to the observation made in mesonic decays, our predictions for the branching ratios of these decays now come out lower than the experimental data; especially the higher-order corrections always increase the LO predictions and shift our predictions closer to the experimental data. Our predictions for the two ratios Br(Λ b → Λ + c µ −ν )/Br(Λ b → Λ + c π − ) and Br(Λ b → Λ + c K − )/Br(Λ b → Λ + c π − ) are both consistent with the current data, indicating that the non-factorizable effects should be small in these decays. Moreover, we emphasize the fact that the non-leptonic to semileptonic ratio in the baryonic case is consistent with experiment, but shows a tension in the mesonic case (see section 6.4). The discrepancy between our prediction and the current experimental data for Br(Λ b → Λ + c π − )/Br(B d → D + π − ) makes it interesting to evaluate directly the form-factor ratios of Λ b → Λ c and B → D transitions by the lattice community.

Conclusion
We have calculated the NNLO vertex corrections to the colour-allowed tree topology in the framework of QCDF for the mesonic decaysB (s) → D ( * )+ (s) L − and the baryonic decays Λ b → Λ + c L − , with L = {π, ρ, K ( * ) , a 1 }. The calculation of the two-loop correction to the hard scattering kernels requires the evaluation of several dozens of genuine twoscale Feynman diagrams, which describe these heavy-to-heavy transitions at the quark level. We performed this calculation by means of techniques that have become standard in the business of multi-loop computations. It might be worth noting that we evaluated all master integrals analytically [33] in a so-called canonical basis [50], a result which catalyzed the convolution with the LCDA and enabled us to obtain the convoluted kernels almost completely analytically.
The NNLO contributions yield a positive shift to the colour-allowed tree amplitude a 1 , which is sizable for its imaginary part, but small for its real part and its magnitude. Moreover, the amplitude only mildly depends on the ratio of the heavy-quark masses z c = m 2 c /m 2 b . The dependence on the factorization scale gets reduced for the real part compared to the NLO result. This reduction does not occur in the imaginary parts, which is expected, as the latter only arise beyond LO. We performed our analysis using the pole scheme for the heavy-quark masses. A change to the MS scheme does not show any significant shift of the amplitude within the range of physical values for the heavyquark masses. Moreover, the results for the different final states only slightly depend on the light meson LCDA and hence, we can confirm the quasi-universality of the tree amplitude to NNLO accuracy.
In our phenomenological analysis we evaluated the branching ratios to NNLO accuracy, and with the latest values for the non-perturbative input parameters. We find that forB d decays the central values of the theoretical predictions are in general higher compared to the experimental values. Within the given uncertainties the quantities agree at the 2 -3σ level for π and ρ in the final state, and slightly better for K and K * . Compared to the analysis at NLO [13], our increased values for the form factors and the amplitude, together with decreased experimental values, have shifted theory and experiment further apart. ForB s decays, the theory predictions are still plagued by large uncertainties which are mainly due to poorly known form factors. For the baryonic decays, on the other hand, the predicted branching fractions turn out to be 20 − 30% smaller than the experimental ones. It would be interesting to understand the reason for this difference in theB d and the Λ b decays. We therefore propose a systematic analysis of factorization for Λ b decays in the future.
Moreover, we analyzed ratios of non-leptonic and semi-leptonic decay rates in order to further probe the factorization theorem to NNLO. The ratios of different non-leptonic rates turn out to be in good agreement, comparing theoretical prediction and experiment. In the case of non-leptonic to semi-leptonic ratios, on the other hand, the values for |a 1 (D ( * )+ L − )| that we extract from experiment are lower by 2 -3σ compared to the NNLO theory predictions (see also [17]).
One possibility to interpret the entity of these results could be non-negligible power corrections. Given the uncertainties of the branching ratios and a 1 they could be negative in sign and 10 -15% in size on the amplitude level. They could cure the non-leptonic to semi-leptonic ratios, without destroying the agreement in the non-leptonic ratios, especially if they were of a certain universality. It will also be very interesting to investigate what this would imply for the power corrections in charmless non-leptonic decays. Recent analyses that address weak annihilation in charmless non-leptonic decays can be found in [102][103][104].
Another, yet less favourable option would be reduced values of |V cb | times the form factors. As stated already in section 6, without a rigorous treatment of power corrections in the QCDF approach nothing more can be said at the present stage.