Heavy Quark Expansion for Heavy Hadron Lifetimes: Completing the $1/m_b^3$ Corrections

We complete the calculation of the contributions from the dimension six operators in the heavy quark expansion for the total lifetime of heavy hadrons. We give the expressions for the Wilson coefficients of the Darwin term $\rho_D$ and the spin-orbit term $\rho_{\rm LS}$ at tree level.


Introduction
Inclusive weak decays of hadrons with a single heavy quark Q have been intensively studied over the last decades [1][2][3].The most inclusive quantity is the lifetime of the ground state heavy hadrons which is determined by their weak decay [4][5][6].The theoretical method is the heavy quark expansion (HQE) [7][8][9], which is a combined expansion in the strong coupling α s (m Q ) [10] and inverse powers of the heavy quark mass [11].The leading term of this expansion is free of any hadronic parameter and is given by the decay rate of the "free" heavy quark.The corrections to this statement appear only at order 1/m 2 Q and are given in terms of the residual kinetic energy µ 2 π and the chromomagnetic moment µ 2 G which are both of order Λ 2  QCD .Consequently these corrections should be at the level of a few percent, since the leading order result implies that all lifetimes of hadrons with a single heavy quark Q should be identical up to corrections of order Λ 2 QCD /m 2 Q .In the early days of the HQE this was taken as an embarrassment, since the lifetimes of the bottom hadrons had not been measured yet, and the lifetimes between charmed hadrons differ by factors of two to five, which is related to the fact that the c-quark is too light for the HQE to be a good approximation for these observables [12].
Since then the methods have been refined and the HQE makes quite precise predictions for the lifetime pattern of bottom hadrons and qualitatively describes the pattern of charmed hadrons.In fact, assuming SU (2) flavour symmetry for the light quarks the lifetime differences between the three ground state mesons are driven by the terms of 1/m 3 Q and higher, in particular by the four-quark operators appearing at tree level, which involve light quarks of a particular flavour.
The progress in the HQE for lifetimes rests on two pillars.On the one hand, there are refinements in the HQE by including higher order terms in the 1/m Q expansion [13], on the other hand there are perturbative calculations improving the Wilson coefficients appearing in the HQE.The higher orders in the HQE contain hadronic matrix elements, for which precise lattice predictions became available recently.Based on this, we are entering the precision era for these observables, in particular for the lifetime differences.
However, a few ingredients have not yet been worked out in detail, since they were believed to be irrelevant.For this reason, the full calculation of all terms appearing at 1/m 3  Q has not yet been done, not even at tree level, since it was assumed that such terms will be small and mainly independent of the light-quark flavour.In the present paper we complete the tree/level calculation of the 1/m 3 Q terms for the lifetime of a heavy hadron with a single heavy quark Q.While these contributions are known since some time for the inclusive semi-leptonic case, the full calculation of the terms at order 1/m 3 Q for the non-leptonic width was still missing.In section 2 we describe the current status of bottom-hadron lifetimes.In section 3 we give a short description of the method of the calculation of the non-leptonic width.In section 4 we present our results, and discuss their implications in section 5.

Synopsis on the status of bottom-hadron lifetimes
The measurements of lifetimes and lifetime rations for bottom hadrons have become very precise of the last decades.The current (2019) experimental averages obtained by the Heavy which shows that the HQE technique can be successfully applied to bottom hadron decays, allowing us to make precision predictions.Therefore, B-physics is entering in its precision era.To arrive at such precise theoretical values, several advancements have been made: The leading term in the total decay rate, i.e. with the absence of power corrections, which describes the free b-quark decay and does not contain non-perturbative corrections, is currently known at NLO-QCD [20][21][22][23][24][25][26][27] and at NNLO-QCD in the massless case [28] for non-leptonic decays.For semi-leptonic decays the current precision is NNLO-QCD [29][30][31][32][33][34][35][36][37][38].
The contribution from the second power correction due to the dimension six Darwin and spin-orbit operators is known at LO-QCD [46] and NLO-QCD [47] only for the semi-leptonic case.However, the ρ D contribution for inclusive non-leptonic decays, is still missing.That is precisely the task we address in this this work.In fact, previous studies focus on the four-quark operators appearing at this order which induce lifetime differences at tree level and which are parametrically enhanced by a phase space factor 16π 2 .However, our explicit calculation shows that the coefficient in front of ρ D turns out to be enhanced and thus needs to be taken into account.The contribution from dimension six four-quark operators is known at NLO-QCD [48][49][50].

Outline of the calculation
In this section we give a brief outline of the calculation which in fact contains a few subtleties.A detailed description will be deferred to a more technical publication.We start from the effective Lagrangian for flavor changing transitions due to charged hadronic currents [51] where G F is the Fermi constant, O 1,2 are four-quark operators, and V CKM , V ′ CKM are the corresponding Cabibbo-Kobayashi-Maskawa matrix elements describing weak mixing of quark generations.We consider only the tree-level four-quark operators of current-current type since the Wilson coefficients of these operators are the largest.The numerical values of the Wilson coefficients C 1,2 (µ) at the scale µ = m b , where m b is the value of the b-quark mass, are known in the SM with high precision mainly thanks to using the high order renormalization group improved QCD perturbation theory at the scales between m b and M W .
We are interested in weak decays of beauty hadrons mediated by the CKM leading transitions with the flavour structure b → cūd and b → ccs.The latter decay is additionally slightly suppressed by the phase space available for the decay products due to the mass of the c-quarks.The canonical choice of the operator basis for the decays b → cq 1 q 2 reads [51] where q L denotes the left-handed quark.It is the basis (4) that is used for the computation of the Wilson coefficients.However, for the purposes of the present computation we use a different operator basis (cf.[6]), which is obtained after applying a four-dimensional Fierz transformation to O 2 .The operators of the new basis are diagonal in the color space and have the form with Γ µ = γ µ (1 − γ 5 )/2.We consider two Cabibbo favoured decay channels, (q 1 , q 2 ) = (u, d) and (q 1 , q 2 ) = (c, s) .The main technical tool for our computation is dimensional regularization (D = 4−2ϵ) [52].The Dirac algebra of γ-matrices is usually defined in D = 4 and needs to be properly extended to D-dimensional space time [53][54][55][56].In particular, using Fierz transformations can lead to a non-trivial ϵ dependence [57][58][59][60].With this in mind, an arbitrary change of the operator basis valid in four-dimensional space is not allowed if perturbative corrections of higher order are to be included: the change will require the corresponding change of the set of evanescent operators associated with a given basis.For our computation however the required accuracy is such that one can use the new basis without changing the coefficients C 1,2 .A review of the relevant techniques can be found in, e.g.[61,62].
The B meson decay rate for the inclusive non-leptonic decays can be computed from the discontinuity of the forward scattering matrix element which is computed in the HQE.The property of unitarity of the S-matrix and the optical theorem lead to an expression for the decay width in the form The HQE of the transition operator up to 1/m 3 b is given by Here Γ 0 q1 ), V cb and V q 1 q 2 are the corresponding CKM matrix elements, q stands for a massless quark and are HQET local operators with π µ = iD µ = i∂ µ + g s A a µ T a and π µ = v µ (vπ) + π µ ⊥ .The four-quark operators O (q) 4F i will be defined in Sec.3.1.Finally, the QCD spinor b of the bottom quark is replaced by the HQET fermion field h v .They are related as follows 3.1 Matching of four-fermion operators: computation of C (q) 4F i In this subsection we give some sketch of the matching calculation of four-quark operators relevant for the renormalization of the coefficient C D .We chose the version of the HQE where the bottom and the charm quarks are integrated out at the same scale m c ≤ µ ≤ m b such that only the u, d, s quarks remain as soft (massless) dynamical quarks.As a consequence, the matching coefficients will depend on the mass ratio r = m 2 c /m 2 b ∼ O(1).We only compute those pieces which are relevant for the renormalization of C D .Such contributions are diagramatically represented in Fig. [1].

The channel b → cūd
The relevant operators in the HQE are with the matching coefficients in D = 4 − 2ϵ dimensions These expressions coincide with the results of ref. [6].
The one-loop matrix elements of these four-fermion operators also contribute to the coefficient C D (see Fig. 2).We find that Im T = . . .
and we can determine the counterterm of C D in the MS renormalization scheme.We obtain where
and we can determine the counterterm of C D in the MS renormalization scheme, for which we obtain where

Matching of two-fermion operators: computation of C i
In this section we describe the matching computation of two-quark operators.The different contributions are diagramatically represented in Fig. [3].In order to optimize the computation we find expressions for the quark propagator in an external gluon field A.
In the semi-leptonic case the tree level expression for the HQE can be obtained from expanding the external field propagator for the charm or the up quark in powers of the covariant derivative which automatically generates the proper ordering of the covariant derivatives.However, in the non-leptonic case the leptonic lines are replaced by quark lines and hence we pick up additional diagrams where gluons are emitted from these quark lines and the simple trick from the semi-leptonic calculation cannot be used here.
Still such contributions can most easily be taken into account by using the expression of the quark propagator in the external gluon field in the Fock-Schwinger gauge x µ A µ (x) = 0 (see, e.g.[63]).This is especially convenient because the expansion of the propagator has then an explicitly gauge invariant form.Another important property of the Fock-Schwinger gauge is that it breaks explicitly the translation invariance of the quark propagator, namely S F (x, 0) ̸ = S F (0, x).We obtain with explicit expressions in momentum space given by where S F (p) is the free quark propagator The expressions S F (p) and SF (p) are used for the propagator of the q1 -quark and q 2 -quark respectively to compute the diagrams that do not appear in the semi-leptonic case.
Let's us discuss the peculiarities of each contribution.The computation of the O 1 ⊗ O 1 contribution goes exactly as in the semi-leptonic case.The color structure of the operator O 1 only allows for the radiation of a single gluon from the c-quark in the bS c b line.So we only need to expand the c-quark propagator.The computation is then identical to the case of semi-taonic decays and the corresponding results can be taken from ref. [64].
The computation of O 2 ⊗ O 2 proceeds as the semi-leptonic case as well after replacing c → q 2 .The color structure only allows radiation of a single gluon from the q 2 -quark in the bS q 2 b line (q 2 = d, s), so we only need to expand this q 2 -quark propagator.In this case one faces the IR divergences due to the gluon emission or the expansion of the massless quark propagator.Within the HQE (and OPE in general) the appearance of such infrared divergence signals the mixing between the local operators that constitute the basis of the expansion.The corresponding local operator develop UV divergences and should be properly renormalized.The well known advantage of using dimensional regularization is that both IR and UV divergences are dealt with simultaneously and a uniform manner.This treatment allows us to retain some vital symmetries of the theory and has technical superiority of simplicity.In fact, it is just this phenomenon of mixing that is the most essential and interesting part of the whole calculation.
The computation of O 1 ⊗ O 2 , which is found to be the same that for O 2 ⊗ O 1 , differs from the one in the semi-leptonic case.Here the gluon emission or the expansion of the quark propagators from all lines have to be taken into account.Overall, the computation of the coefficient of the ρ LS operator in HQE is infrared safe even for massless quarks, does not require considering mixing with four quark operators, and can be performed in D = 4. 4 Results for the Wilson coefficients at order 1/m 3 Q Before we give our results for the terms of order 1/m 3 b , we need to discuss the effects induced by operator mixing.The HQE as any OPE of effective theory gives an example of the general phenomenon of the separation of physics at greatly different scales.Indeed, the hadronic width in the representation depends on the heavy quark mass m b and the infrared scale of QCD Λ QCD with m b ≫ Λ QCD .
The HQE in expression (7) is organized in such a way that the coefficients are insensitive to Λ QCD while the matrix elements of the local operators are independent of the large scale m b .An explicit naive computation, however, produces at some intermediate stage both IR singularities in the coefficient functions and UV singularities in the matrix elements.The combinatorics of HQE is such that all singularities cancel.Technically the most efficient way to perform computations is to use dimensional regularization for both IR and UV divergences.In such a setup one has to take into account the mixing of local operators at UV renormalization.Thus, an infrared divergence of coefficient functions signals the UV mixing between the local operators that constitute the basis of the expansion.The corresponding local operator develop then UV divergences and should be properly renormalized.
In our computation a naive way of getting the coefficient of the ρ D operator leads to an IR singularity in the contribution of O 2 ⊗ O 2 and O 1 ⊗ O 2 correlators as one gets the radiation of a soft gluon from the light quark line.
This singularity is canceled by the UV renormalization of the four-quark operator that has the general form (quite symbolically) where the matrix Γ gives the corresponding Dirac structure of the four-quark operator and the quantity γ(Γ) is the mixing anomalous dimension depending on the Dirac structure Γ.
The UV pole in ϵ coming from the operator mixing in eq. ( 37) cancel the IR divergence in the coefficient function for the operator ρ D .
The HQE of the imaginary part of the transition operator is given by eq. ( 7).However, it is convenient to rewrite it in terms of the local operator b/ vb defined in full QCD.It can be employed to remove O 0 in the HQE by using the expansion Inserting this expression we get (omitting here and in what it follows the four-quark contributions) which has the advantage that the forward matrix element of the leading term is normalized to all orders.
For the operator O v we use the equation of motion to remove it from the expression for Im T This is the desired expression for the transition operator, from which we compute the total decay rate in terms of the HQE parameters We obtain At leading order we have c S = c D = C mag = 1.It is convenient to define new coefficients corresponding to every matrix element We find for the coefficients in the case of Γ(b → cūd) and for the case Γ(b → ccs) Here r = m 2 c /m 2 b and z = √ 1 − 4r.We note that the equalities C 0 = C µπ and C µ G = C ρ LS are a consequence of reparametrization invariance [65].We take this here as a check of our calculation.We also note that the results for the coefficient of ρ D depend on the calculational scheme.This does not only concern the use of the MS scheme, but also the treatment of the Dirac algebra in D dimensions.This is related to the fact that the Fierz-rearrangement in D dimensions generates evanescent operators which result in constants to be taken into account when comparing results [66].
After using the transformation rules (3.20-3.23) in [67], and after proper definition of evanescent operators for the b → cud channel (see Sec. 4.1), these results are in agreement with [66], where the coefficients were computed in four dimensions.We will comment in more detail on this in the next section.
Note that from these results one readily finds the coefficients of HQE in eq. ( 7) whose computation was described in Sec.3.2

Comment on the basis of four-quark operators
Our results discussed above are expressed in the operator basis while one may chose as well the basis Õ(u) which has been used in ref. [66].While the two bases are equivalent in D = 4, the situation for arbitrary D is more involved.Relating the two bases in D dimensions requires the addition of new operators called evanescent operators.The choice of the evanescent operator is not unique, and a particular recipe reduces to a substitution [6,62] A conventional choice is a = 4 and b = −4, with d = 4 − 2ϵ.We will call the basis fixed by this choice to be the canonical basis of four-quark operators.The evanescent operators are thus defined as The choice of the evanescent operators E QCD 1,2 is not unique.This choice is motivated by the requirement of validity of Fierz transformation at one-loop order [6,54].
Thus the complete operator basis reads and the rule for the transformation between the two bases is In the new basis the imaginary part of the transition operator becomes The operators do not contribute to the anomalous dimension of C ρ D .However, the change of basis produces a shift in the ρ D coefficient which depends on the choice of the evanescent operators i.e. on a and b.We call the new coefficient C MS ρ D (a, b).The difference between the results obtained in the two bases is whereas the difference between our results and the ones obtained in ref. [66], where the coefficients are computed in D = 4, is Note that for the canonical choice of the evanescent operators the difference vanishes.
As we mentioned there is some freedom when choosing the evanescent operators E QCD 1,2 .The difference in the results due to the different choice of a and b corresponds to a shift in the coefficient When inserting numbers one has to keep in mind, that the coefficient is thus dependent on the scheme.This scheme dependence is compensated by the value of ρ D itself, which is a scheme-dependent quantity.

Numerical analysis
In this section we give numerical values for phenomenological applications.We will chose the canonical scheme for the evanescent operators such that the results in [66] can be directly compared to our results.We employ the MS scheme for the definition of ρ D and chose for the scale µ = m b .
For both channels we have contributions which come from the operators O 1 and O 2 , which come with the Wilson coefficients C 1 and C 2 , see ( 3).In table 1 we give the numerical values of the coefficients for the transition b → ccs.We also list the values of the coefficients for the transition b → cūd in the four-quark operator basis of Sec.3.1.1 in table 2, and in the canonical basis in table 3.
In order to get an idea about the size of the total contribution of ρ D to the non-leptonic width we insert values for the Wilson coefficients We denote 4F i |B(p B )⟩/(2M B ) and use the abbreviation Γ q1 q 2 to refer to Γ(b → cq 1 q 2 ).We obtain Assuming that ⟨E QCD 4F i ⟩/m 3 b ∼ Λ 3 QCD /m 3 b , the Darwin coefficient gives a correction to the tree level values of the coefficients of the fourquark operator of ∼ 7% for b → ccs and of ∼ 1% for b → cūd in the canonical basis (we take the largest coefficient of the four-quark operators to compare).

Discussion and conclusions
We have computed the contributions of the Darwin and the spin-orbit term appearing at order 1/m 3 Q in the HQE of the non-leptonic width.Although the coefficient of the spin-orbit term is fixed by reparametrization invariance, we have explicitly computed it as a check of our methods.
The most interesting part is the computation of the coefficient of ρ D , since one has to take into account the mixing with the four-quark operators.The coefficient of the Darwin term turns out to be sizable which was also found in the semi-leptonic case.
In fact, this may become relevant for lifetime differences.To be specific, we will consider the SU (3) Flavor triplet of ground state B hadrons B = (B u , B d , B s ).It has been noticed already very early [17] that up to and including 1/m 2 b the operators appearing in Im T are SU (3) flavour singlets and hence a lifetime difference to this order can only emerge from the SU (3) breaking coming from the states, meaning that µ π and µ G differ between the three Bmeson ground states.In turn, assuming the SU (3) flavour symmetry, no lifetime differences can be induced up to this order.Since the coefficients of the HQE parameters at order 1/m 2 b are small, the effect on lifetime differences due to the SU (3) flavour breaking in µ π and µ G is very small.
The situation changes at the order 1/m 3 b where four-quark operators appear, involving light quarks.Since the weak hamiltonian is sensitive to the light-quark flavor, the resulting four-quark operators have different matching coefficients.Analyzing the SU (3) flavour structure of the four-quark operators, we may decompose them into a singlet and an octet contribution with respect to SU (3) Flavor according to hv Γ(q1q where q = (u, d, s) is the light quark triplet and 1 and T a are the generators of U (3) flavour .If we, in addition, use the equation of motion (where λ a are the Gell-Mann matrices of SU (3) color ) we can eliminate ρ D in favour of fourquark operators, contributing to the SU (3) Flavor singlet part only.Obviously the mixing of ρ D can only happen with the SU (3) Flavor singlet part of the four-quark operators.
Overall, we thus find that we can re-write ( 49) as (schematically) where the sums run over the matrix elements of the four-quark operators.
Clearly the octet part of the four-quark operators is the main source for lifetime differences, which is present also if the states are exactly SU (3) symmetric.However, the precision of the lifetime measurements has increased and thus also the SU (3) breaking through the states needs to be taken into account, which means that also the matrix elements T 4q i,singlet will contribute to lifetime differences.This effect may be important, since the complete calculation of these terms shows that the coefficients of these terms are large, even enhanced by phase space factors.However, a quantitative study of their impact on lifetime differences needs estimates of the SU (3) flavour breaking in the matrix elements T 4q i,singlet which is beyond the scope of the present paper.a combination of the following three master integrals We are only interested in the imaginary part of the corresponding integrals, related to the discontinuity across the cut.We denote J ≡ Im J. On the one hand J(0, 1, 1, 0, 0) = 0. On the other hand, we can use that d dm 2 c J(1, 1, 1, 0, 0) = J(1, 2, 1, 0, 0) + J(1, 1, 2, 0, 0) = 2J(1, 2, 1, 0, 0) , and the reduction of J(1, 2, 1, 0, 0) to a combination of the three master integrals above in order to express the master integral J(2, 1, 1, 0, 0) only in terms of J(1, 1, 1, 0, 0) and J(0, 1, 1, 0, 0) whose imaginary part is zero Therefore there is only one master integral we need to compute, which is J(1, 1, 1, 0, 0).Note that eq. ( 96) has a pole in D = 4 − 2ϵ dimensions.Therefore the O(ϵ) expansion of J(1, 1, 1, 0, 0) will be needed.We where The computation of the decay b → ccs in a scheme with a hard IR regularization The case we are considering can be used for pedagogical purposes of OPE (HQE).There is only an IR divergence in the original amplitude.One can proceed by regulating IR with a small mass and use four dim to compute the coefficients.The relevant master integral is then defined as follows J(0, 1, 1, 0, 2; m 0 )| d=4 ≡ d 4 q 1 d 4 q 2 ((p + q 1 − q 2 ) 2 − m 2 c )(q 2 2 − m 2 c )((q 2 + p) 2 − m 2 0 ) 2 . .This is a rough sketch of the procedure used in [66] (for a related discussion, see [6]).In such an approach one needs an expansion of the master integral at the limit of small m 0 /m b .
We have obtained an analytical expression for the required expansion in the form J(0, 1, 1, 0, 2;

Figure 1 :
Figure 1: One loop diagrams contributing to the matching coefficients of four-quark operators in the HQE of the forward scattering matrix element of the B meson.

Figure 2 :
Figure 2: One loop diagrams contributing the renormalization of C D .

Figure 3 :
Figure 3: Two loop diagrams contributing to the matching coefficients of two-quark operators in the HQE of the forward scattering matrix element of the B meson.

Table 3 :
Numerical values for the coefficients of b → cūd in the canonical basis (a = 4, b = −4).For illustration we take the numerical values µ = m b = 4.8 GeV and m c = 1.3 GeV.

Table 1 :
Numerical values for the coefficients of b → ccs.For illustration we take the numerical values µ = m b = 4.8 GeV and m c = 1.3 GeV.

Table 2 :
Numerical values for the coefficients of b → cūd.For illustration we take the numerical values µ = m b = 4.8 GeV and m c = 1.3 GeV.