QCD corrections to inclusive heavy hadron weak decays at $\Lambda_{\rm QCD}^3 /m_Q^3$

We present an analytical calculation of the $\alpha_s$ corrections for the coefficient of $\rho_D/m_Q^3$ term in the heavy quark expansion for the inclusive semileptonic decays of heavy hadrons, such as $B \to X_c \ell \bar{\nu}$. The full dependence of the coefficient on the final-state quark mass is taken into account. Our result leads to further improvement of the theoretical predictions for the precision determination of CKM matrix element $|V_{cb}|$.


Introduction
The discovery of a Higgs boson a few years ago completed the the Standard Model of particle physics (SM), which thereby became a highly predictive framework, i.e. it allows us to perform very precise calculations. On the experimental side there is currently no hint at any particle or interaction which is not described by the SM, even at the highest possible energies. This implies that particle physics is about to enter an era of precision measurements of the SM parameters.
In particular, accurate measurements accompanied by precise theoretical calculations in the flavour sector of the SM have already proven to have an enormous reach at scales that are much larger than the center-of-mass of any existing or projected colliders [1]. Aside from large-scale experimental efforts, this strategy also requires accurate theoretical computations.
The need in obtaining a high precision of theoretical predictions in particular in the flavour sector is urgent, since the structure of the quark mixing is expected to be rather sensitive to possible the effects from physics beyond the SM (BSM). While the SM has successfully passed a variety of tests within current precision (as a review, see e.g. [2,3,4]), any further insights will require the use of even more accurate theoretical predictions.
The weak decays of quarks mediated by charged currents occur at a tree level and are believed to not have sizable contributions from BSM Physics. However, the study of such decays is importance for the precise determination of the numerical values of the SM parameters, in particular CKM matrix elements. For heavy quarks (i.e. for heavy hadrons) a reliable theoretical treatment of weak decays is possible, because the mass m Q of the decaying heavy quark constitutes a perturbative scale that is much larger than the QCD infrared scale Λ QCD , m Q Λ QCD .
Heavy quark expansion (HQE) techniques provide a systematic expansion of physical observables in powers of the small parameter Λ QCD /m Q . Quantitatively, the techniques work well for bottom quarks, since a typical hadronic scale associated with binding effects in QCD is Λ QCD ∼ 400 − 800 MeV. With some reservations, the HQE has been used for charmed quarks as well though the analysis is expected to be more of qualitative nature, since the charm-quark mass is not sufficiently large. Thus, the HQE and the corresponding effective theory of heavy quarks and soft gluons (HQET) have become the major tools of modern precision analyses in heavy quark flavor physics [5,6,7,8].
In particular for the determination of V cb from inclusive b → c semileptonic transitions the HQE has brought an enormous progress. The HQE expansion for the total rate and for spectral moments have been driven to such a high accuracy that the theoretical uncertainty in the determination of V cb is now believed to be at the order of about one percent. However, this assumes that higher order terms in Λ QCD /m Q and α s (m Q ) are of the expected size. In fact, the leading order terms (i.e. the partonic rate) has been fully computed to order α 2 s (m Q ), the first subleading terms of order (Λ QCD /m Q ) 2 are known to O(α s (m Q )) while all higher-order term in the (Λ QCD /m Q ) n (n = 3,4,5) are known only at tree level.
In the present paper we analytically compute parts of the QCD corrections to the contributions of order (Λ QCD /m Q ) 3 , which is one of the not yet known pieces. We point out that the size of these terms is expected to be of the same order as the partonic α 3 s (m Q ) contributions, likewise the terms of order α 2 s (m Q )(Λ QCD /m Q ) 2 . However, at the level of the current precision these terms turn out to be small and hence their calculation is to validate the assumption that they are of the expected size.
Specifically we compute the coefficient of the power suppressed dimension six Darwin term ρ D at next-to-leading order (NLO) of the strong coupling perturbation theory with the full dependence on the final state quark mass. Compared to the calculations at lower orders in the HQE, it has some new features since the mixing of operators of different dimensionality in HQET has to be taken into account for the proper renormalization of the Darwin term coefficient.

Heavy Quark Expansion for Heavy Flavour Decays
In this section we set the stage by giving the very basics for the theoretical descriptions of semileptonic decays, in particular for the decay B → X c lν. A more detailed description can be found in e.g. [9].
The low-energy effective Lagrangian L eff for the semileptonic b → clν l transitions reads where the subscript L denotes the left-handed projection of the fermion fields and V cb is the relevant CKM matrix element.
Using optical theorem one obtains the inclusive decay rate B → X c ν from taking an absorptive part of the forward matrix element of the leading order transition operator T (see e.g. [9]) The transition operator T is a non-local functional of the quantum fields participating in the decay process. Since the quark mass is a large scale compared to the scale Λ QCD of QCD m Q Λ QCD , the relevant forward matrix element still contains perturbatively calculable contributions. These can be separated form the non-perturbative pieces by employing effective field theory tools, which allows us an efficient separation of the kinematical m Q and dynamical Λ QCD scales involved in the decay process.
For a heavy hadron with the momentum p H and the mass M H , a large part of the heavy-quark momentum p Q is due to a pure kinematical contribution due to its large mass p Q = m Q v + ∆ with v = p H /M H being the velocity of the heavy hadron. The momentum ∆ describes the soft-scale fluctuations of the heavy quark field near its mass shell originating from the interaction with light quarks and gluons in the hadron. This is implemented by re-defining the heavy quark field Q(x) by separating a "hard" oscillating phase and a "soft" field b v (x) with a typical momentum of order Inserting this into (2) we get where L is the same expression as L with the replacement Q(x) → b v (x). This makes the dependence of the decay rate on the heavy quark mass m Q explicit and allows us to build up an expansion in Λ QCD /m Q by matching the transition operator T in QCD onto an expansion in terms of Heavy Quark Effective Theory (HQET) operators [11,12].
Generally, the HQE for semileptonic weak decays is written as (e.g. [13]) The coefficients a i , i = 0, 2, D, LS depend on the ratio m 2 c /m 2 b , while µ 2 π , µ 2 G , ρ D and ρ LS are forward matrix elements of local operators, usually called HQE parameters.
The precise definition of the appropriate mass parameter for the heavy quark field is of utmost importance for the precision of the predictions of the HQE and is thus extensively discussed in the literature, see e.g. [10]. The HQE parameter µ 2 π is the kinetic energy parameter for the B-meson in HQE, µ 2 G is the chromo-magnetic parameter. The term ρ LS contains the spin-orbital interaction and ρ D is the Darwin term which is of our main interest in the present paper.
The power suppressed terms are becoming important phenomenologically as the precision of experimental data continues to improve. The coefficients a i have a perturbative expansion in the strong coupling constant α s (m Q ). The leading coefficient a 0 is known analytically to O (α 2 s ) precision in the massless limit for the final state quark [14]. At this order the mass corrections have been analytically accounted for the total width as an expansion in the mass of the final fermion [15] and for the differential distribution in [16].
The coefficient of the kinetic energy parameter is linked to the coefficient a 0 by reparametrization invariance (e.g. [17]). The NLO correction to the coefficient of the chromo-magnetic parameter a 2 has been investigated in [18] where the hadronic tensor has been computed analytically and the total decay rate has been then obtained by direct numerical integration over the phase space.
This calculation allows for the application of different energy/momentum cuts in the phase space necessary for the accurate comparison with experimental data.
The NLO strong interaction α s correction to the chromo-magnetic coefficient a 2 in the total decay rate has been analytically computed in [19,20]. The techniques of ref. [19,20] allow also for an analytical computation of various moments in the hadronic invariant mass or/and that of the lepton pair. In the present paper we give the NLO result for the coefficient a D of ρ D in the analytical form retaining the full dependence on the charm quark mass.

NLO for Darwin term ρ D : Calculation and Results
In this section we describe the actual computation of the coefficient a D of the Darwin term. The present calculation follows the techniques used earlier for the determination of NLO corrections to the chromo-magnetic operator coefficient in the total width [19,21,20]. Here we give a brief outline of the calculational setup, for details of the techniques, see [20].
We consider a normalized transition operatorT defined by The heavy quark expansion for the rate is constructed by using a direct matching from QCD to where we retain only the Darwin term in the 1/m 3 Q order. The local operators O i in the expansion (6) are ordered by their dimensionality O Here the field h v is the heavy quark field the dynamics of which is given by the QCD Lagrangian expanded to order 1/m 3 Q . Furthermore, π µ = iD µ is the covariant derivative of QCD and π µ = v µ (vπ) + π µ ⊥ . The coefficients C 0 , C v , C G , C D of the operators are obtained by matching the appropriate matrix elements between QCD and HQET.
After taking the forward matrix element with the B-meson state one can use the HQET equations of motion for the field h v in order to eliminate the operator O v . We note that, in general, there is an additional operator O 5 =h v (vπ) 2 h v in the complete basis at dimension five, however it will be of higher order in the HQE after using equations of motion of HQET.
Note that one can use the full QCD fields for the heavy quark expansion afterwards [22]. It is convenient to choose the local operatorb/ vb defined in full QCD as a leading term of the heavy quark expansion [23] as it is absolutely normalised and provides a direct correspondence to the quark parton model as the leading order of the HQE.
We note that (6)  The leading order result is given by two-loop Feynman integrals of a simple topology -the so called sunset-type diagrams [24,25], while at the NLO level one has to evaluate three-loop integrals with massive lines due to the massive c-and b-quarks. In Fig. 1 we show a typical three-loop diagram for the power corrections in the heavy quark expansion.
We use dimensional regularization for both ultraviolet and infrared singularities. We used the systems of symbolic manipulations REDUCE [26] and Mathematica [27] with special codes written for the calculation. For reduction of integrals to master integrals the program LiteRed [28] was used. The master integrals have been then computed directly. Some of them have checked with the program HypExp [29].
For the Darwin term one takes an amplitude of quark to quark-gluon scattering and projects it to an HQET operator. We choose a momentum k gluon and take the structure ( v)k 2 ⊥ . There are several operators in HQET that can have such a structure, for instance,h v (πv)π 2 ⊥ h v . This operator is irrelevant because it is of higher power on shell. One disentangles the mixing of such operators with the Darwin term by using two quark momenta k 1 and k 2 and pick up the structure (k 1 k 2 ) that emerges in the coefficient of the Darwin term. The other operators can have k 2 1 or k 2 2 structure. The coefficient a D is defined in front of the meson matrix element. After taking the matrix element one can use equation of motion of HQET to reduce the number of the operators in the basis. One more conventional step is to trade the leading order operatorh v h v for the QCD operatorb/ vb that provides correspondence to the parton model.

The final expression for the coefficient a D is then
where C HQET D is the NLO coefficient of the operator O D in HQET Lagrangian, and C bvb D is the NLO coefficient of the operator O D in the expansion ofb/ vb. At the LO we find a LO D = −5r 4 − 8r 3 + 24r 2 + 36r 2 log(r) − 88r + 48 log(r) + 77 (13) where r = m 2 c /m 2 b that agrees with [30]. The coefficient contains a logarithmic singularity log(r) at small r. This singularity reflects the mixing to hidden/intrinsic charm contribution [31]. At higher powers even more singular terms (like 1/r) can appear [32]. The matching is performed by integrating out the charm quark simultaneously with the hard modes of the b-quark. This means that we treat m 2 c /m 2 b as a number fixed in the limit m b → ∞, and therefore our results cannot be used to extrapolate to the limit m c → 0.
An important check of a loop computation consists in verifying the cancellation of poles after performing the appropriate renormalization of the physical quantity in question. Since in the case at hand this is quite delicate, we briefly discuss the renormalization of the ρ D coefficient at NLO within our computation.
We single out the pole contribution to the NLO coefficient in the form The contribution to the coefficient C N LO−pol D from one-particle irreducible diagrams reads The proper cancellation of these poles requires to consider the mixing between HQE operators of different dimensionality which is known to be possible in HQET [33,34,35,36,37,38]. The anomalous dimensions of the operators are numbers independent of r while the functions of r appearing in C D should cancel in the renormalization of the coefficient. This is a rather strong restriction because only the coefficient functions of the lower power operators (which are basically C 0 (r) and C v (r)) can be used for the pole cancellation in C D . The leading order C G coefficient is proportional to C LO 0 (r). Indeed, the explicit expression reads These properties of HQE are important for the implementation of the renormalization procedure of the Darwin-term coefficient.
The operator O π from HQE after the insertion of one more O π from the Lagrangian can mix with O D that produces the pole structure proportional to C 0 (r) The relation (16) means that the ghh vertex computed in perturbation theory within HQET with one insertion of O π gets a contribution from higher powers of the HQET Lagrangian (see, e.g. [33]). By the same token the operator O G from HQE after the insertion of one more O G from the Lagrangian that produces the pole structure proportional to C 0 (r) again because of eq. (15). The cross-insertions (O G from HQE to O π from the Lagrangian and vice versa) renormalize the spin-orbit operator at the order 1/m 2 b . This type of mixing is known for a long time from computation of 1/m 2 b running of coefficients of HQET Lagrangian.
In the literature the renormalization is considered often for the static heavy fields when the contributions of reiterated terms in the HQET Lagrangian are accounted for through the bi-local operators [36,37]. We consider the standard approach and treat higher order terms as perturbations (see, [35,39]). In our case the inclusion of these mixings does not suffice to cancel all the poles in C D as there are other structures than C 0 (r) necessary. The operator O v can mix with double insertions of higher dimensional terms, i.e.
and the counterterm proportional to O D emerges from two insertions of the operators O π . The mixing matrix γ vD is unknown. But the effect of such a mixing leads to the appearance of the coefficient C v (r) in the expression for the poles. One can now fit the pole function with two entries C 0 (r) and C v (r).
Thus we infer the corresponding mixing anomalous dimensions and find that the combination the NLO contribution shifts the ρ D coefficient by 10%.

Discussion
The technical details of the calculation will be discussed in a more detailed paper, where we also plan to calculate moments of various distributions. However, the result presented here already have a few interesting consequences.
The first remark concerns the dependence on the mass of the charm quark which appears in the ratio r = m 2 c /m 2 b . This ratio is kept at a fixed value as m b , m c → ∞ and the behavior of the coefficients close to r = 0 is given by Note that the behaviour of the coefficient at the border of the available phase space depends on the definition of the c-quark mass. Here we use MS mass. Ii Fig. 2 (left panel) we plot the dependence on r in the full kinematically allowed region range 0 ≤ r ≤ 1. We show the ratio while the right panel of Fig. 2 focuses on the physical region around r = 0.07.
The plots show that the mass dependence in the physical region is weak, while it is sizable over the full range. As we discussed above, the massless limit cannot be taken, since in the case of a b → u transition additional operators have to be taken into account. Nevertheless at small r the NLO corrections become even smaller and have a zero at r ∼ 0.005. A similarly strong dependence has been observed also for the QCD corrections in the coefficient of µ 2 G [20]. Although the corrections are not untypically large, they will have a visible impact on the determination of V cb . This is mainly due to the fact that the coefficient in front of ρ D in the total rate is quite large, see (21). While a detailed analysis will require to repeat the combined fit as e.g. in [40] we may obtain a tendency from an approximate formula given in eq. (12) in this paper.
According to (21) the NLO correction corresponds to a reduction of the contribution of ρ D by 10%, thus eq. (12) of [40] implies a shift in the central value of which is about a third of the current theoretical uncertainty.
However, parametrically this correction is of the same size as the yet unknown corrections of order α 2 s Λ 2 QCD /m 2 b and α 3 s , which would need to be included in a full analysis up to order α s Λ 3 QCD /m 3 b . We note in passing that the corrections α s ρ 2 LS are not needed, since these are included in the known α s µ 2 G contributions [41]. Nevertheless, the contribution of ρ D is significant due to the large coefficient in front of ρ D and hence we expect that the impact of this correction is largest.

Appendix
The matching coefficient of the "operator" ρ D in the HQET Lagrangian in NLO at µ = m b gets a correction [12] 1 After using equation of motions the final ρ D -coefficient is expressed through the coefficient of the relevant operator in HQE (direct contribution) and the contributions due to HQET Lagrangian and the choice of the full QCD operator at the leading power through the relation At LO one obtains a LO D = −5r 4 − 8r 3 + 24r 2 + 36r 2 log(r) − 88r + 48 log(r) + 77 (28) that agrees with [30].
We write the coefficient of ρ D term after taking matrix elements as The C F color part reads at NLO Here Li 2 (z) is a dilogarithm.
The new master integral N m has appeared compared to our previous results. The relevant combination turns out to be (N p + N m )/2, and both N p and N m should be evaluated at LO in -expansion.