NLO QCD corrections to inclusive semitauonic weak decays of heavy hadrons up to $1/m_b^3$

This paper presents $\alpha_s$ corrections to the inclusive semitauonic $B \to X_c \tau \bar{\nu}_\tau$ decay width and leptonic invariant mass spectrum up to order $1/m_b^3$ in the heavy quark expansion. Analytical results with full dependence on the final-state quark and tau lepton masses are obtained. The results up to $1/m_b^2$ can be extrapolated to the $B \to X_u \tau \bar{\nu}_\tau$ decay by taking the limit of vanishing final-state quark mass.


Introduction
to discuss the impact of the new results.

HQE for inclusive decays of heavy flavoured hadrons 2.1 The total rate
This section briefly describes the theoretical framework used for the calculation of inclusive semileptonic decays of heavy hadrons and provides the main definitions. At low momentum transfer compared to M W the heavy quark decay b → qτν τ can be described by an effective Fermi Lagrangian where the subscript L denotes left-handed fermionic fields, G F is the Fermi constant and V qb is the corresponding CKM matrix element. By using the optical theorem, the inclusive decay rate of B → X q τν τ is obtained from the imaginary part of the forward hadronic matrix element of the transition operator T where M B is the heavy hadron mass. Since the heavy quark mass m b (hard scale) is much larger than the hadronization scale of QCD Λ QCD (soft scale) the equation above contains perturbatively calculable contributions which can be factorized from the non-perturbative ones by using an operator product expansion, the so-called HQE (see e. g. [32][33][34][35] for more details). Then the Im T , computed as an expansion near the heavy quark mass shell in QCD, is matched to an expansion in Λ QCD /m b by using local operators in HQET [36,37].
where Γ 0 = G 2 F m 5 b /(192π 3 ), C i are the matching coefficients which have a perturbative expansion in the strong coupling constant α s (m b ), and O i are the HQET operators listed below }h v (Spin-orbit operator) .
Here π µ = iD µ = i∂ µ + g s A a µ T a is the QCD covariant derivative, a µ ⊥ = a µ − v µ (v · a) with v = p B /M B being the heavy hadron velocity, and h v is the HQET field whose dynamics is determined by the HQET Lagrangian [37]. Operators which are of higher dimension on shell have been neglected.
It is convenient to trade the leading term operator O 0 in Eq. (3) by the local QCD operatorb/ vb, since its forward hadronic matrix element is completely normalized. To that purpose, the HQE of theb/ vb operator is needed up to the desired order whereC i are the matching coefficients which are pure numbers. In addition, the equation of motion (EOM) of the HQET Lagrangian is used for the operator O v in Eq. (3).
where c F (µ) = 1 + α s 2π are the coefficients of the chromomagnetic and Darwin operators in the HQET Lagrangian at NLO [37]. Note that, due to reparametrization invariance [38], the coefficient of the spin-orbit operator is linked to the one of the chromomagnetic operator c S = 2c F − 1. The parameter µ is the renormalization scale and C F = 4/3, C A = 3 are colour factors. Finally, the HQE for inclusive semitauonic weak decays is written as whereC i ≡ C i − C 0Ci are defined as the difference between the coefficients C i of the HQE of the transition operator in Eq. (3) and the current in Eq. (10) multiplied by C 0 . In the b → cτν τ transition, the matching coefficients C i (i = 0, µ π , µ G , ρ D , ρ LS ) depend on two ratios ρ = m 2 c /m 2 b and η = m 2 τ /m 2 b , where m c is the charm quark mass and m τ is the tau lepton mass. In the b → uτν τ transition the matching coefficients only depend on η since the up quark mass is neglected. Note that reparametrization invariance also links the coefficients of the HQE of the rate, C 0 = C µπ and C µ G = C ρ LS [39][40][41]. The HQE hadronic parameters µ 2 π , µ 2 G , ρ 3 D and ρ 3 LS are defined as the following forward matrix elements of HQET local operators taken between full QCD states [39].

The spectrum on the dilepton invariant mass q 2
This section briefly describes the construction of the HQE for the spectrum on the dilepton invariant mass q 2 starting from the expression for total rate Eq. (2). The approach described in [42] is followed, which takes advantage of the fact that the leptonic tensor is unaffected by QCD corrections, so it always appears factorized from the hadronic tensor. Therefore, one can use a dispersion representation defined in dimensional regularization [43] to write the leptonic tensor as an integral differential in q 2 . For the case of a massive lepton and a massless antineutrino the leptonic tensor reads , and is the fourmomentum flowing through the leptons. Note that, as long as the dispersion representation is defined in dimensional regularization, there is no need for subtractions since singularities are properly regularized. Indeed only the imaginary part of Eq. (20) is required in the present calculation. Since taking the imaginary part makes the integral non-singular, commutation of the integral over q 2 with expansion in and imaginary part is justified and it can be conveniently taken Let us stress that, whereas this is technically not necessary, it is useful as it simplifies the calculation. Note that it is enough to expand the right hand side of Eq. (21) to O( 0 ) since renormalization can be performed for the differential rate, as it holds at the level of the hadronic tensor. Also note that, after using the dispersion representation, the leptonic tensor becomes an "effective massive propagator" of mass q and the dependence on m τ factorizes from the hadronic part. Therefore, the master integrals necessary for the computation of the differential rate are the same that appear in the b → ceν e transition, where e is considered to be massless. These master integrals were computed analytically in [42] and they are used for the calculation of the q 2 spectrum. The HQE for the dilepton invariant mass spectrum can be written in analogy to the HQE for the total rate are defined as the difference between the coefficients C i of the HQE of the transition operator (in differential form) and the current Eq. (10) multiplied by C 0 . Note that the q 2 spectrum also satisfies reparametrization invariance, so C 0 = C µπ and C µ G = C ρ LS [44]. The coefficients C i (i = 0, µ π , µ G , ρ D , ρ LS ) of the differential rate depend on three ratios r, ρ and η and they are related to the corresponding coefficients C i of the total rate by However, as it happens with the lepton energy spectrum, the q 2 spectrum cannot be interpreted point-by-point. In order to compare to experiment, one computes moments of the spectrum for which a clean HQE exists. In analogy to the total rate, the HQE for q 2 moments is written as , whose coefficients are related to the coefficients of the differential rate by Note that the coefficients M 0,i of the zeroth moment correspond the coefficients C i of the total rate. Finally, note that the minimum dilepton invariant mass in inclusive semileptonic decays is q 2 min = m 2 . In the case of decays to electron or muon the masses can be safely disregarded and q 2 min = 0. Such a low q 2 is difficult to detect and experimentalists use cuts while integrating up to the available q 2 [45]. However, for decays to tau q 2 min = m 2 τ ∼ 3.17 GeV 2 and it might not be necessary to introduce cuts. For example, the lowest q 2 cut introduced in the measurement of B → X c eν e moments by [45] was q 2 min = 3.0 GeV 2 , which suggests that for semitauonic decays the whole q 2 range can be measured. Despite of this, it might be useful to study moments as a function of low q 2 cuts if one pursues an analysis similar to the one in [46]. Cut moments can be easily obtained from Eqs. (24) and (25) by integrating in the desired range.

Outline of the calculation
This section sketches the outline of the calculation. First one computes the matching coefficients for the q 2 spectrum by using the spectral representation given in Eq. (21), and latter one integrates over q 2 to obtain the coefficients of the total rate. Since Eqs. (3), (10) and (11) hold a the operator level, the matching calculation can be performed by computing on shell matrix elements with partonic states i.e. with quarks and gluons. For power corrections it is possible to compute the coefficients by using a small momentum expansion near the heavy quark mass shell [35,47].
The LO-QCD and NLO-QCD contributions to the differential rate are given by one-loop and twoloop heavy quark to heavy quark (b → b) or heavy quark to gluon heavy quark (b → gb) scattering Feynman diagrams which are computed in dimensional regularization. The Feynman gauge is chosen. The scattering is computed in the external gluonic field by using the background field method. The amplitudes are reduced to a combination of master integrals by using LiteRed [48,49]. The corresponding master integrals were computed analytically in [42]. Dirac algebra is handled by using Tracer [50].
As for renormalization, the MS scheme is used for the strong coupling α s (µ) and the HQET Lagrangian. The bottom and charm quarks are renormalized on-shell. That is where quantities with the subscript B stand for bare, quantities which do not have subscript stand for renormalized, and Z OS 2 = Z OS m b to this order. Therefore, results are conveniently presented in the the on-shell scheme for both quark masses m c and m b .
Note that in order to achieve more precise theoretical predictions one typically uses a short distance mass for the bottom quark, like the 1S mass [51][52][53][54][55] or the kinetic mass [56,57]. For the charm quark mass one often uses the MS mass. A change in the mass scheme can be implemented by using the known one-loop relation between the different mass schemes.
For the presentation of results, the coefficients of the q 2 spectrum and the total rate are split as it follows (i = 0, v, µ G , ρ D ) Analytical results for coefficients of the right-hand side are provided in the Mathematica notebook "cobqtv.nb", which is supplied as an ancillary file. Following [42], the coefficients of the q 2 spectrum are given in terms of the convenient variables together with η. The coefficients of the differential rate have a simple dependence on the tau lepton mass. They are polynomials of degree three in η, which follows from Eq. (21). In other words, it is a consequence of the factorization of the leptonic and hadronic tensors. Contrarily, the dependence on the the charm quark mass and the invariant mass of the leptons is more complicated and it is efficiently expressed in terms of x ± . Analytical results require the use of dilogarithms Li 2 (x). The coefficients of the total rate are obtained by analytical integration of Eq. (23) which, as pointed out in [42], it can be achieved after changing variables to x − by using x + = ρ/x − . After the change of variables, the integral in Eq. (23) becomes The coefficients of the total rate have a complicated dependence on ρ, η. It is convenient to express them in terms of the variables z ± defined by since they simplify the results drastically. Analytical results require the use of dilogarithms Li 2 (x) and trilogarithms Li 3 (x). It is instructive to see how results in the B → X c τν τ decay channel can be compactly written in terms of x ± and η for the differential rate, and in terms of z ± for the total rate. For this reason, the LO expressions will be explicitly written in the text. The NLO expressions are too lengthy to be explicitly written and they are provided only in the ancillary file. For the B → X u τν τ decay channel expressions are reasonable in length and they will be explicitly shown up to NLO.

Leading power coefficient
To leading power, the heavy hadron decay corresponds to the free heavy quark decay. The NLO QCD corrections to the inclusive b → qτν τ decay rate and distributions were computed analytically to this order almost thirty years ago [21,22]. This section reproduces the known results for the decay rate and the q 2 spectrum as a check of the calculation.
For the computation, one takes the heavy quark to heavy quark scattering amplitude with onshell external momentum p 2 = m 2 b and projects it to the O 0 operator. The diagrams that contribute to the coefficient of the differential rate at leading power are shown in Fig. [1]. The coefficient of For the B → X c τν τ case, the partonic coefficient of the differential rate at LO reads At NLO the expression is too lengthy to be explicitly written in the text and it is provided in the ancillary file. After integration over r in the whole range one obtains the partonic coefficient of the total rate at NLO. The LO expression is obtained after integrating Eq. (35). It reads Note that, in the expression above, the overall functions multiplying the (z − − z + ) and logarithmic structures are written in terms of ρ and η instead of z ± . To LO this compacts the results significantly. However, this is not the case to NLO where writing the expression in terms of z ± simplifies the result drastically. The coefficients given in the ancillary file are written only in terms of z ± . Again the NLO the expression is too lengthy to be explicitly written in the text and it is provided in the ancillary file. For the leading power coefficient the limit ρ = 0 can be taken, which allows to find expressions for the B → X u τν τ case. The coefficient of the differential rate reads After integration of Eqs. (37) and (38) over r in the whole range one obtain the partonic coefficient of the total rate. It can be also obtained by taking the limit ρ = 0 in the coefficient of the partonic rate for the B → X c τν τ case. It reads where ζ(x) is the Riemann zeta function. These expressions are in agreement with the known results [21,58]. By taking the limit η = 0 one recovers the known expressions for the B → X q eν e case, where leptons are massless [32,42,59]. These two independent limits are a strong check of the calculation.
If desired, the leading power coefficient of moments can be easily obtained by numerical integration of Eq.(25).

Chromomagnetic operator coefficient
The NLO corrections to the 1/m 2 b terms in inclusive semileptonic decays of massless leptons are known [32-35, 40, 60, 61], but to the best of my knowledge they have never been computed for the case where one of the leptons is massive. However, to LO the 1/m 2 b corrections to inclusive semitauonic decays have been extensively studied [26][27][28]. This section addresses the computation of the missing α s /m 2 b corrections for the coefficient of the dilepton invariant mass spectrum and the total width.
Before addressing the calculation of the 1/m 2 b corrections, one needs to compute the auxiliary coefficientC v , which appears in the HQE of the transition operator at order 1/m b . This coefficient contributes to higher orders in the 1/m b expansion after using the EOM and, in particular, shifts the coefficients of the 1/m 2 b and 1/m 3 b terms. For its calculation one takes the amplitude of the quark to quark-gluon scattering diagrams given in Fig. [2] with vanishing soft momenta k 1 ⊥ = k 2 ⊥ = 0 and projects it to the O v operator. This is achieved by taking a longitudinally polarized gluon exchange (v · ) without momentum transfer. In practice, one directly computes the difference between the HQE of the transition operator and the current The explicit expression for theC v coefficient is not written here because it appears only at intermediate steps. However, the analytical expression is provided in the ancillary file.
To order 1/m 2 b one only needs to compute the coefficient of µ 2 G , since the coefficient of µ 2 π is related to the leading power coefficient due to reparametrization invariance [38,41]. Note that at dimension five the operator O I =h v (v · π) 2 h v also appears, but it contributes to higher orders in the 1/m b expansion after using the EOM. The computation follows the lines of [32], where one takes the amplitude of quark to quark-gluon scattering, expands to linear order in the small momentum and projects it to the chromomagnetic operator O G .
The diagrams that contribute are shown in Fig. [2]. For the determination of the chromomagnetic operator coefficient one takes k 1 ⊥ = 0 and k 2 ⊥ = k ⊥ . In other words, one considers a single small Like in the previous case, one directly computes the difference between the HQE of the transition operator and the currentC again with with C i,B = C i,B (m c,B = Z OS mc m c ) and where is the renormalization factor of the chromomagnetic operator. The coefficientC G is finite and the cancellation of poles provides a solid check of the computation. The coefficient of the differential rate C µ G in front of µ 2 G is finally obtained from For the B → X c τν τ case, the chromomagnetic operator coefficient of the differential rate at LO reads At NLO the expression is too lengthy to be explicitly written in the text and it is provided in the ancillary file.
For η = 0 one recovers the known results for the chromomagnetic operator coefficient of the leptonic invariant mass spectrum [32] and the total rate [33][34][35] in B → X c eν e , where leptons are massless. Note that a typo in the latter references was identified and corrected in [32].
If desired, the chromomagnetic operator coefficient of moments can be easily obtained by numerical integration of Eq.(25).

Darwin operator coefficient
The 1/m 3 b corrections in inclusive semitauonic decays have been considered to LO in [29][30][31], where expressions for the lepton energy spectrum and total rate have been obtained. To the best of my knowledge, they have never been considered for the q 2 spectrum. The NLO corrections to the 1/m 3 b terms in inclusive semileptonic decays of massless leptons have been computed quite recently [32]. However, they are unknown for inclusive semitauonic decays. This section addresses the computation of the α s /m 3 b corrections for the coefficient of the dilepton invariant mass spectrum and the total rate.
In this case, in addition to the Darwin and spin-orbit operators, there are five more operators which nevertheless contribute to higher orders in the 1/m b expansion after using the EOM. The contributions to these operators must be properly disentangled from the contributions to the Darwin and spin-orbit terms by choosing appropriate projectors. Note that reparametrization invariance relates the coefficient of ρ 3 LS and µ 2 G , so one only needs to compute the coefficient of ρ 3 D . The computation follows the lines of [32] where one takes the amplitude of quark to quark-gluon scattering, expands to quadratic order in the soft momenta and projects it to Darwin operator O D . The diagrams that contribute are shown in Fig. [2]. For the determination of the Darwin operator coefficient one takes the external gluon to have longitudinal polarization (v · ), use two soft quark momenta k 1 ⊥ and k 2 ⊥ , and pick up the structure k 1 ⊥ · k 2 ⊥ .
Like before, one directly computes the difference between the HQE of the transition operator and the currentC are the renormalization factor of the Darwin operator and the contribution to the Darwin coefficient coming from the HQET operator mixing under renormalization [62][63][64][65][66][67][68][69][70], respectively. The quantitȳ C D is finite and the cancellation of poles provides a solid check of the calculation [32]. The coefficient of the differential rate C ρ D in front of ρ 3 D is finally obtained from For the B → X c τν τ case, the Darwin operator coefficient of the differential rate at LO reads At NLO the expression is too lengthy to be explicitly written in the text and it is provided in the ancillary file. After integration over r in the whole range one obtains the Darwin operator coefficient of the total rate at NLO. The LO expression is obtained after integrating Eq. (57). It reads Again the NLO the expression is too lengthy to be explicitly written in the text and it is provided in the ancillary file. The LO expression in Eq. (58) confirms the result given in [31] and disagrees with [29,30]. In the case η = ρ, Eq. (58) can be compared to the C 2 1 structure in the nonleptonic b → ccs decay [71,72] for which we find agreement. To the best of my knowledge, the LO result for the Darwin coefficient of the q 2 spectrum has never presented explicitly. The NLO results are new.
Note that the matching is performed by integrating out the charm quark simultaneously with the hard modes of the bottom quark. That means m 2 c /m 2 b is treated as a number of order one in the limit m b → ∞, which implies that m c → ∞ as well. Therefore, the results can not be extrapolated to the ρ = 0 (B → X u τν τ ) case in general. In particular, the limit ρ → 0 happens to be non-singular for power corrections up to 1/m 2 b and it can be taken, but this is no longer true for the Darwin term. This feature shows up as a logarithmic singularity in the Darwin coefficient when taking the limit ρ → 0. In the calculation for strict ρ = 0 these singularities show up as poles which point out the mixing under renormalization of the Darwin operator with dimension six four quark operators (see e. g. [71][72][73][74]). This is an additional complication beyond the scope of this work. Therefore, α s /m 3 b corrections are presented only for the B → X c τν τ case. For η = 0 one recovers the known results for the Darwin coefficient of the leptonic invariant mass spectrum [32] and the total rate [32,75] in the case of massless leptons B → X c eν e . For the total rate, the coefficient previously computed in [47] was corrected in [32].
As discussed in [32], performing the integration over r of the Darwin coefficient to obtain the coefficients of the total rate and moments is rather subtle, since the integral is infrared (IR) singular at the upper integration limit r max = (1 − √ ρ) 2 , whereas the regularization parameter has been omitted after renormalization of the differential rate. This feature points out that expansion in and integration over r are operations which do not commute in general. Therefore, computing the coefficient of the total width and moments requires restoring the dependence in the IR singular terms. For the Darwin term there is a single IR divergent integral with singularity (r max − r) −3/2 , which should be understood as the following dimensionally regulated integral For the integration one takes to be an arbitrary complex number and computes the integral in a domain of convergence. Later on one performs an analytic continuation in through the complex plane to = 0. Note that the IR divergent integral is finite in dimensional regularization due to the Parameter Numerical value where f (r, η) is an arbitrary test function regular at r = r max and the plus distribution (r max − r) Note the close connection between IR singular integrals at the endpoint, the corresponding dimensionally regulated integrals, and delta functions sitting at the endpoint. Note that, in contrast to what it happens in the charged lepton energy spectrum, the coefficient function in front of the delta distribution is not singular. In this case, the computation of the Darwin coefficient of moments by numerical integration of Eq. (25) is more involved, since it requires computing integrals in dimensional regularization. However, the only integral which needs to be regulated is the total rate, for which analytical results are provided. In other words, all moments of the form Eq. (25) can be related to the total rate and IR safe integrals that can be evaluated numerically if desired

Numerical analysis
This section presents a numerical analysis in order to illustrate the size of the corrections which have been computed. The numerical values summarized in table 1 are taken for illustration. In order to compare the coefficients of the differential rate with LO and NLO precision, as well as the shapes of the LO and NLO contributions, their dependence on the dilpeton invariant mass is shown in Fig. 3 for the B → X c τν τ decay channel, and in Fig. 4 for the B → X u τν τ decay channel.
Note that for the leading power coefficient the shapes of the LO and NLO contributions are very similar but with opposite sign. This is no longer true for the coefficients of the power corrections. In general, the α s corrections to the coefficients of the power corrections are smaller at low q 2 than at high q 2 .  In order to see the importance of power corrections to the spectrum, Fig. [5] shows its dependence on the dilepton invariant mass by including subsequent power corrections for both decay channels, B → X c τν τ and B → X u τν τ . It is interesting to compare the spectrum of semitauonic decays with their massless lepton counterparts. For this reason, plots are also presented for the B → X c eν e and B → X u eν e decay channels. Note that the lepton mass shifts the position of the peak to the right hand side of the spectrum and it also suppresses it. The latter statement is expected from experimental observation of the rate. The spectrum at large q 2 for b → u falls softer than for b → c. In general, power corrections become larger at high q 2 and the IR behaviour at the endpoint becomes sharper for higher power corrections, pointing out the breakdown of the power expansion at large q 2 .
Some illustrative values for the different contributions to the momentsM n , normalized moments  M n , and ratios of moments R q /q n between different decay channels are also provided, wherē Numerical values are presented for the four decay channels B → X q ν with q = c, u and = e, τ . This information is complementary to the one found in plots and helps to understand better the size of the corrections that have been computed, since they are of a few percent, and their size is difficult to illustrate in plots. Moreover, normalized moments and ratios might be useful as they are ideal observables to be compared to experiment because they do not depend on the absolute normalization of the rate which is very sensitive to the heavy quark mass (proportional to m 5 b ). The results are given in the form (for normalized moments and ratios this is possible after re-expanding in 1/m b ) where and A =M n ,M n , R q /q n . Numerical values are provided in tables for the coefficients on the right hand side of Eq. (65), for the contribution of every power correction split in LO and NLO contributions (term by term sum), and for A n itself (total). Numerical values are presented in tables 2, 3, 4, and 5 for the moments, in tables 6, 7, 8, and 9 for the normalized moments, and in tables 10, 11, 12, and 13 for the ratios.
As for moments, observe that in the B → X c τν τ decay the LO 1/m 3 b corrections represent a 8-24% correction to the leading term. The α s corrections to the Chromomagnetic and Darwin  terms represent a 0.5-1.3% correction to the leading term. In general, the corrections due to the Chromomagnetic and Darwin terms are similar in size. In the B → X u τν τ decay the α s corrections to the Chromomagnetic term represent a 0.2-2.8% correction to the leading term. As for normalized moments, observe that in the B → X c τν τ decay the LO 1/m 3 b corrections represent a 4-16% correction to the leading term. The α s corrections to the Chromomagnetic and Darwin terms represent a 0.1-0.4% and 0.6-2.6% correction to the leading term, respectively. Note that, due to some numerical cancellations, the α s /m 3 b corrections are comparable or even larger than the α s corrections to the partonic term. In semitauonic decays the Darwin term corrections happen to be numerically larger than the Chromomagnetic term corrections, whereas in the case of massless leptons they are of similar size. In the B → X u τν τ decay the α s corrections to the Chromomagnetic term represent a 0.6-3.5% correction to the leading term.
As for ratios, observe that in R ce/cτ n the LO 1/m 3 b corrections represent a 5-6% correction to the leading term. The α s corrections to the Chromomagnetic and Darwin terms represent a 0.2-0.5% and 0.5-0.6% correction to the leading term, respectively. The Darwin term corrections happen to be numerically larger than the Chromomagnetic term corrections. For the other ratios, only corrections up to α s /m 2 b are available. The α s corrections to the Chromomagnetic term represent a 0.1-3.6% correction to the leading term. Note that in the ratios R uτ /cτ n and R ue/ce n numerical cancellations happen for the coefficient of the chromomagnetic term at LO which makes, in some cases, its NLO contribution to be numerically more important. Finally, note that normalized moments and ratios are insensitive to µ 2 π .

Conclusions
The current experimental value for R(D ( * ) ) in B decays mediated by the semitauonic b → cτν τ transition shows a more than 3σ deviation from the standard model prediction, which might point out the presence of new physics coupled to the τ lepton. In order to reveal the presence of new physics, the study of semitauonic decays is necessary by any means. In particular, inclusive decays may provide valuable complementary information to the one of exclusive decays, which motivates its investigation [10,11]. If new physics is present in semitauonic decays, their study at the precision level might be necessary for its future establishment. Currently, the inclusive semitauonic decays have been poorly measured as they are experimentally challenging. However, their precise measurement should be possible in the near future by Belle II. In this paper, analytical expressions for the dilepton invariant mass spectrum and total rate of inclusive semitauonic B → X c τν τ and B → X u τν τ decays have been obtained up to order α s /m 3 b and α s /m 2 b , respectively. The dependence on the final state quark and the tau lepton masses is taken into account in its full analytical form. Analytical results are provided in the ancillary Mathematica notebook "cobqtv.nb". Therefore, the largest unknown contributions are the NNNLO corrections to the leading power term and the LO corrections to 1/m 4 b terms. For the case of massless leptons, these corrections have been considered for the total rate and a number of kinematical moments [77][78][79][80].
It has been observed that for q 2 moments, normalized moments and ratios the LO 1/m 3 b corrections represent roughly a 10% correction to the leading term. The α s corrections to the Chromomagnetic and Darwin terms represent roughly a 1% correction to the leading term.
Precise theoretical predictions for the decay width, normalized moments and ratios can be derived by using our results. This requires nevertheless the use of some short distance mass like the 1S mass or the kinetic mass. The determination of the final estimates for the theory predictions are left to future publications. It might be interesting to compare the theoretical predictions for the observables studied in this paper to future experimental measurements and look for possible discrepancies.
The results of this paper can be a used as a partial check and a breakthrough for future analytic calculations of NLO corrections to power suppressed terms in inclusive non-leptonic decays, in particular for the b → ccs channel, which contains two massive particles in the final state, very much like in b → cτν τ .         Table 6: Numerical values for the coefficients of the normalized moments in the B → X c τν τ channel.    Table 8: Numerical values for the coefficients of the normalized moments in the B → X c eν e channel.