O(\alpha_s^2) corrections to semileptonic decay b \to c l \bar \nu_l

We compute the next-to-next-to-leading order (NNLO) QCD corrections to b \to c l \bar \nu_l decay rate at a fully differential level. Arbitrary cuts on kinematic variables of the decay products can be imposed. Our computation can be used to study the NNLO QCD corrections to the total decay rate as well as to the lepton energy, hadronic invariant mass and hadronic energy moments and to incorporate those corrections into global fits of inclusive semileptonic B-decays.


I. INTRODUCTION
Inclusive semileptonic decays of B-mesons into charmed final states are benchmark processes at Bfactories. Because of relatively large rates and clean experimental signatures, these decays can be studied with great precision. On the other hand, theoretical description of semileptonic B decays is robust thanks to the Operator Product Expansion (OPE) in inverse powers of the b-quark mass m b . The application of the OPE to semileptonic decays of B-mesons leads to the conclusion that both the total decay rate and various kinematic distributions can be described by power series in Λ QCD /m b [1]. For infinitely heavy b-quark, the decay rate B → X c lν l coincides with the rate computed at the quark level. For realistic values of bottom and charm masses, a few non-perturbative matrix elements that enter at order (Λ QCD /m b ) n , n = 2, 3 are accounted for in existing theoretical predictions.
In recent years, many measurements of moments of charged lepton energy and hadronic invariant mass in B → X c lν l decays have been performed by BABAR, BELLE, CLEO, CDF and DELPHI [2,3,4,5,6]. Comparison of these experimental results with theoretical predictions for corresponding observables leads to the determination of the Cabibbo-Kobayashi-Maskawa (CKM) matrix element |V cb |, bottom and charm quark masses and a number of non-perturbative parameters such as µ 2 π and µ 2 G [7,8]. Typical precision claimed in these analyses is about one percent for |V cb | and m b and a few percent for m c and non-perturbative matrix elements [7,8].
To achieve such precision, advances in theoretical understanding of semileptonic B-decays were necessary, including subtle interplay between perturbative and nonperturbative physics and significant developments in technology of multi-loop computations. While one-loop corrections to both b → clν l total decay rate [9] and a number of important differential distributions are known since long [10], it is interesting to remark that phenomenologically relevant triple differential distribution in charged lepton energy, leptonic invariant mass and hadronic invariant mass was computed through O(α s ) only a few years ago [11,12]. This fact illustrates the complexity of perturbative calculations, when massive particles are involved, at a fully differential level.
Given the precision of available experimental measurements, good understanding of non-perturbative effects and a fairly large value of the strong coupling constant α s (m b ) ≈ 0.24, it is expected that O(α 2 s ) corrections to b → X c lν l decays are required for a consistent theoretical description. However, as was realized long ago, technical complexity of such an endeavor is daunting. To simplify the problem, the O(α 2 s ) corrections to b → X c lν l were computed in three specific kinematic points [13,15,16]. These results were used in Ref. [13] to estimate the NNLO QCD corrections to Γ(b → X c lν l ). Unfortunately, such a description is necessarily limited in its scope even for the total rate and a generalization of such an approach to more differential quantities, such as lepton energy and hadronic invariant mass moments, is clearly out of question.
On the other hand, a subset of the NNLO QCD corrections, the BLM corrections [17], received significant attention recently. The BLM corrections are associated with the running of the strong coupling constant; they are potentially important since the QCD β-function is large. For B decays, however, the BLM corrections are known to be modest if proper definition of quark masses is adopted and judicious choice of the renormalization scale in the strong coupling constant is made. The BLM effects are the easiest NNLO effects to calculate since they can be obtained from a one-loop computation if the latter is performed with a non-vanishing gluon mass [18]. For this reason, in the past, the BLM corrections to b → clν l were calculated for the total rate and various kinematic moments [19,20]. However, the NNLO QCD corrections beyond the BLM approximation, for which genuine two-loop computations are required, remained missing.
Calculation of these two-loop corrections became possible recently thanks to developments in numerical approaches to multi-loop computations [21]. These numerical methods benefit from the absence of mass hierarchy in the problem which is the case for b → c decays, since masses of bottom and charm quarks are close. The possibility to use the approach of Ref. [21] to describe decays of charged particles was recently pointed out in [22] where electron energy spectrum in muon decay was computed through second order in the perturbative expansion in QED.
The goal of this Letter is to present the computation of O(α 2 s ) corrections to b → X c lν l decay rate at a fully differential level. Our results can be used to calculate arbitrary observables related to inclusive b → c transition through NNLO in QCD. For example, second order QCD corrections to such popular observables as lepton energy, hadronic invariant mass and hadronic energy moments can be studied in dependence of the cut on charged lepton energy. Inclusion of the results of our computation into global fits, should lead to a reduction of the theoretical uncertainty in the determination of |V cb |, the bottom and charm quark masses and the non-perturbative parameters that contribute to the decay rate.

II. COMPUTATION
In this Section, we set up our notation and briefly describe technical aspects of the computation. A detailed description of the method can be found in [21,22].
Consider the decay b → X c lν l where the final state lepton is massless. The differential decay rate can be written as where G F is the Fermi constant, m b is the b-quark pole mass, a = α s /π and α s is the MS strong coupling constant defined in the theory with five active flavors and renormalized at the scale m b . For numerical computations, we use m b = 4.6 GeV and m c = 1.15 GeV. While these numerical values for the quark masses can not be justified in the pole scheme, our choice is motivated by an eventual necessity to transform the pole scheme computation to a more suitable scheme. The values of the quark masses that we employ in this Letter correspond to the central values of m b,c in the "kinetic scheme" [23], derived in recent fits to inclusive semileptonic B-decays [7,8].
To calculate the functions dF 0−2 , we have to account for different processes. At leading order, dF 0 is computed by squaring the matrix element of the process b → clν l and summing or averaging over spins and colors, as appropriate. At next-to-leading order, dF 1 receives contributions from virtual O(α s ) corrections to b → clν l and from the real-emission process b → clν l + g. To compute dF 2 , we require two-loop O(α 2 s ) corrections to b → clν l , one-loop O(α s ) corrections to b → clν l +g and the double real-emission corrections b → clν l + X, where X refers to two gluons or a quark-antiquark pair or a ghost-antighost pair. We will refer to these corrections as double-virtual, real-virtual and double-real, respectively. In addition, we have to account for a variety of renormalization constants, when computing higher order corrections. We do not include the process b → clν l + cc in our calculation since the energy release in this process is so small that it can not be treated perturbatively.
To calculate the NNLO QCD corrections, the method for multiloop computations developed in [21,22] is employed; in those references a detailed discussion of many technical issues relevant for the current computation can be found. One technical aspect that we improve upon relative to Refs. [21,22] is how virtual corrections to single gluon emission process b → clν l + g are treated. In Refs. [21,22] these corrections were dealt with by an analytic reduction to master integrals followed by a numerical evaluation of those. This method, however, becomes impractical quite rapidly, once the number of external particles or the number of massive particles in the problem increases. In principle, the real-virtual corrections can be computed numerically, but for heavy-to-light decays this is complicated because some Feynman diagrams develop imaginary parts. To handle these imaginary parts, we proceed as follows. For all Feynman diagrams that contribute to real-virtual corrections, it turns out possible to identify a Feynman parameter that enters the denominator of the integrand linearly. Let us call this Feynman parameter x 1 . Then, a typical integral that has to be computed reads Here n ≥ −1, b > a > 0 and both a and b depend on other Feynman parameters and the kinematic variables. The two arguments of the function I refer to lower and upper limits for x 1 integration. To calculate I(0, 1), we note that by extending upper integration boundary to infinity, a solvable integral for arbitrary a, b and n is obtained. On the other hand, since and because the denominator of the integrand in I(1, ∞) is sign-definite, I(1, ∞) can be computed numerically in a straightforward way. It turns out that, up to minor modifications, this trick can be used to avoid dealing with the imaginary parts for all Feynman diagrams that contribute to one-loop corrections to single gluon emission process in b → c decays.
Because couplings of quarks and leptons to the charged current are chiral, proper treatment of the Dirac matrix γ 5 in d = 4 − 2ǫ dimensions is important. While this problem can be avoided in the computation of the total decay rate, for more differential quantities it becomes an issue. We use the approach of Ref. [24] where a consistent framework for extending the axial vector current to ddimensions is given.
Our computation can be checked in a number of ways. First, the double-virtual, real-virtual and double-real corrections are divergent when taken separately but these divergences must cancel in physical observables. We have checked these cancellations for a variety of observables, from the inclusive rate to various moments with cuts on both charged lepton energy and the hadronic invariant mass. Second, in the limit m c → m b , the NNLO QCD corrections to the decay rate b → clν l are described by the so-called zero-recoil form factors computed through O(α 2 s ) long time ago [16]. We have checked that in the limit m c → m b our computation reproduces the zerorecoil form factors. Third, we can use published results for the BLM corrections to the total rate and charged lepton energy, hadronic invariant mass and hadronic energy moments [20] to check parts of our computation related to massless quark contributions to gluon vacuum polarization. Finally, considering the limit m c ≪ m b , we reproduce the NNLO QCD corrections to b → ulν l decay rate reported in Ref. [25].

III. RESULTS
We are now in position to discuss the results of our computation. We consider a number of observables, mostly for illustration purposes. We present the results in the pole mass scheme and use the strong coupling constant renormalized at the scale m b . While the pole mass scheme is known to be an unfortunate choice inasmuch as the convergence of the perturbative expansion is concerned, we decided to present our results in this way for clarity. However, we emphasize that the impact of the NNLO QCD corrections, computed in this paper, on the determination of |V cb |, heavy quark masses and the nonperturbative parameters, including kinetic and chromomagnetic heavy quark operators, can only be assessed once the pole mass scheme is abandoned in favor of a more suitable quark mass definition and the NNLO QCD corrections are included into the fit.
To present the results, we follow Ref. [20] and define where ... denotes average over the phase-space of all final state particles, E l,h is the energy of the charged lepton or hadronic system in the b-quark rest frame and The lepton energy moments introduced in Eq.(4) can be written as where ellipses stands for higher order terms in the perturbative expansion in QCD. Similar decomposition can be performed for the hadronic energy moments H n . In addition, we use β 0 = 11 − 2/3N f and define the non-BLM corrections L  In Tables I,II the results for lepton energy and hadronic energy moments with and without a cut on the lepton energy are displayed. The numerical accuracy of L is about 1 − 3%. It is possible to improve on the accuracy but this requires somewhat large CPU time. Nevertheless, for all practical applications the achieved numerical accuracy is sufficient.
There are a few interesting observations that follow from Tables I,II. Quite generally, the non-BLM corrections and the BLM corrections have opposite signs; given their relative magnitude and the value of β 0 , it is easy to see that the O(α 2 s ) corrections are about twenty percent smaller than what the BLM-based estimates suggest. The relative magnitude of the non-BLM and BLM corrections is largely independent of n and of whether the lepton energy cut is applied.  First row in Table I provides the NNLO QCD corrections to the total decay rate b → clν l in the pole mass scheme. Such corrections were estimated earlier in Ref. [13]. Note that in Ref. [13] the numerical results are given for the ratio of quark masses m c /m b = 0.3 and also the BLM corrections are defined with N f = 4, rather than N f = 3. Calculating the non-BLM corrections for the set of parameters employed in [13], we find L (2) 0 ≈ 1.73 which is to be compared with the estimate L (2) 0 ≈ 0.9(3), reported in [13,14]. The results of Ref. [13] were used in Ref. [26] to estimate the impact of the QCD corrections on Γ(B → X c lν l ). In Ref. [26] the perturbative corrections to b → clν l decay rate are described by a factor A pert , defined as where r = m c /m b . A pert depends on the adopted scheme for the quark masses. In the kinetic mass scheme, A pert (0.25) = 0.908 is quoted. To arrive at this result, Ref. [26] uses L (2) 0 = 1.4 which is about a factor 2.5 smaller than the corresponding entry in Table I. Correcting for this discrepancy, we derive A pert (0.25) = 0.919.
We believe that this value for the perturbative renormalization factor in the kinetic scheme for m c /m b = 0.25 should be employed in global fits of semileptonic Bdecays. Further analysis of entries in Table I suggests that the QCD corrections in general and the non-BLM corrections in particular mostly affect the overall normalization rather than shapes of kinematic distributions. This follows from the approximate independence of L (1,2) n /L (0) n of n and also of whether or not the cut on the lepton energy is imposed. It is therefore possible to speculate that the non-BLM corrections computed in this Letter will mostly affect the extraction of |V cb | whereas their influence on, e.g., the b-quark mass determination will be minor. Concerning the |V cb |, the increase of the perturbative renormalization factor A pert by 10 × 10 −3 implies the change in the value of |V cb |, extracted in Ref. [8], by about −0.25 × 10 −3 . On the other hand, since non-BLM corrections were not included in a fit of Ref. [7], the shift in the value of |V cb | derived in that reference will likely be larger ∼ −0.5 × 10 −3 . Although expected shifts in central values of |V cb | are not large, we stress that they are comparable to uncertainties in |V cb |, derived in the global fits [7,8].

IV. CONCLUSIONS
In this Letter, the computation of the NNLO QCD corrections to the fully differential b → clν l decay rate is reported. The differential nature of the computation makes it possible to apply arbitrary cuts on the kinematic variables of final state particles. This result allows to extend the existing determinations of the CKM matrix element |V cb |, the bottom and charm quark masses and the non-perturbative parameters µ 2 π and µ 2 G from global fits to semileptonic decays of B-mesons, by including the NNLO QCD corrections exactly. We note that for a consistent high-precision analysis of semileptonic B-decays, also O(α s ) corrections to Wilson coefficients of non-perturbative kinetic and chromomagnetic operators are required. Such a correction is available for the kinetic operator [27] but is still missing for the chromomagnetic.
We presented a few results for charged lepton energy moments and hadronic energy moments with and without a cut on the lepton energy in the pole mass scheme. These results suggest that the magnitude of the non-BLM corrections does not depend strongly on the kinematics; the non-BLM corrections are approximately 2% for all the moments considered. We therefore expect that the non-BLM NNLO QCD corrections will mostly affect the determination of |V cb | decreasing its central value by about one percent whereas their impact on the quark masses and the non-perturbative parameters will probably be quite mild.
As a final remark, we note that it would be interesting to extend this calculation in two ways. First, one may consider semileptonic decays of B-mesons into massive leptons. Such an extension, relevant for the description of B → X c + τ +ν τ decay is straightforward. Second, it is interesting to extend the current calculations to allow for a massless quark in the final state. This is a difficult problem but it is highly relevant for the determination of the CKM matrix element |V ub | from semileptonic b → u transitions.