Completing $1/m_b^3$ corrections to non-leptonic bottom-to-up-quark decays

We compute the contributions of dimension six two-quark operators to the non-leptonic decay width of heavy hadrons due to the flavor changing bottom-to-up-quark transition in the heavy quark expansion. Analytical expressions for the Darwin term $\rho_D$ and the spin-orbit term $\rho_{\rm LS}$ are obtained with leading order accuracy.


Introduction
The HQE [1][2][3] provides a well-established theoretical tool to study hadrons containing a single heavy quark Q. In particular, it has been used intensively over the last decades to study inclusive decays of heavy hadrons [4][5][6]. The study of these systems is specially interesting because it requires understanding the interplay between electroweak and strong interactions. Within this framework, the decay width of the heavy hadron is written as a combined expansion in the strong coupling constant α s (m Q ) [7] and the inverse heavy quark mass m Q [8]. The leading term of that expansion predicts that the decay of the heavy hadron is given by the weak decay of the "free" heavy quark it is made of [9][10][11], whereas it is blind to the spectator quark it contains. As a consequence, lifetimes are predicted to be the same for all hadrons which contain the same heavy quark. The first correction to this picture appears at O(1/m 2 Q ), which introduces lifetime differences between hadrons with different spin configurations. Lifetime differences caused by a different spectator quark flavour content appear for the first time at O(1/m 3 Q ) through four-quark operators, which explicitly contain spectator quarks of a particular flavour. These contributions are suppressed by the heavy quark mass and, as a consequence, they are expected to be small corrections to the "free" heavy quark decay.
Thus, the HQE provides a systematic way of improving theoretical predictions by either including matrix elements (ME) of higher order operators (including more terms in the 1/m Q expansion) [12] or improving the precision in α s (m Q ) of the perturbative coefficients in front of them.
When it comes to test it, the HQE has proven to be extremely successful to describe lifetimes and lifetime differences of bottomed hadrons, whose experimental determination has reached an impressive precision. The current (2019) experimental averages obtained by the Heavy Flavor Averaging Group (HFLAV) are [13] τ (B s ) τ (B d ) = 0.969 ± 0.006 , and even higher precision seems to be achievable from the most recent results from LHCb [14] and ATLAS [15], so B-physics is in its precision era.
The contribution from the second power correction (dimension six two-quark operators, also called Darwin and spin-orbit terms) is known at NLO-QCD [43][44][45] for semi-leptonic decays and at LO-QCD for non-leptonic decays [46,47].
Finally, the contribution from dimension six four-quark operators, which induces lifetime differences and which is enhanced by a phase space factor 16π 2 , is known at NLO-QCD [11,48,49].
In this work, we compute the coefficients of the dimension six two-quark operators (Darwin and spin-orbit terms) for the Cabibbo suppressed channels of O(λ 3 ) stemming from b → u transitions. These coefficients were already computed in Ref. [46] but not in Ref. [47]. Our aim is to perform a check of the results presented there while using a completely different approach following the lines of Ref. [47]. The coefficients presented here may be relevant to improve the theoretical precision of B-hadron lifetimes and lifetime differences.
The results obtained here, especially for the b → uūd channel, can be easily applied to charm decays, where their contribution is expected to be more important because the HQE for charm has a slower convergence than for bottom. These coefficients might be helpful to clarify the status of the HQE for charm [49,50], whose validity has been often questioned due to the smallness of the charm quark mass.
The paper is organized as follows. In Sec. 2 we present the outline of the calculation. We give definitions in Sec. 2.1 and compute the matching coefficients of the dimension six twoand four-quark operators at tree level in Secs. 2.2 and 2.3, respectively. Renormalization of the coefficients is discussed in Sec. 2.4. In Sec. 3 we present the results and discuss the role of evanescent operators in Sec. 3.1. Finally, we perform a numerical analysis in Sec. 3.2. We also give some technical results in the Appendix.
2 Outline of the Calculation

Definitions
Weak decays of hadrons are mediated by the weak interactions of their quark constituents. These flavor changing transitions are described very efficiently through an operator product expansion (OPE) approach which gives rise to an effective Lagrangian [51]. Despite of the plethora of effective operators necessary to describe the weak interactions we only consider here the tree-level four-quark operators of current-current type, since their Wilson coefficients are the largest. The effective Lagrangian describing b → q 3q1 q 2 transitions reads where G F is the Fermi constant, O 1 , O 2 are four-quark operators of the current-current type, and V q 1 q 2 , V * q 3 b are the corresponding Cabibbo-Kobayashi-Maskawa (CKM) matrix elements which describe the mixing of quark generations under weak interactions. However, not all the operators of this kind have the same importance. Their weight is determined by the size of the corresponding CKM matrix elements L eff ∼ λ 2 (bcūd +bccs) (2) +λ 3 (buūd +bucs +bcūs +bccd) (3) +λ 4 (buūs +bucd) , where λ = V us = 0.2257 +0.0009 −0.0010 is a Wolfenstein expansion parameter. The contribution to the total non-leptonic width coming from the dominant O(λ 2 ) operators was already considered in Refs. [46,47]. In this work, we consider the Cabibbo suppressed contribution of O(λ 3 ) stemming from the b → u transitions (q 3 = u), which was already considered in Ref. [46]. Our aim is to perform a check of the results presented there while using a completely different approach following the lines of Ref. [47]. It is worth mentioning that the remaining Cabibbo suppressed contributions of O(λ 3 ) and O(λ 4 ) can be directly obtained from the results presented here and in Refs. [46,47] by exchanging the d and s quarks 1 .
In summary, we are interested in weak decays of B-hadrons mediated by transitions of the type b → uq 1 q 2 with (q 1 , q 2 ) = (u, d) and (q 1 , q 2 ) = (c, s). In that case, the operator basis reads [51] O with Γ µ = γ µ (1 − γ 5 )/2. However, for the purpose of the computation it is convenient to use an operator basis diagonal in the color space [11] which is obtained after applying a Fierz tranformation in the operator O 2 . The computation is carried out using the standard techniques of dimensional regularization (D = 4 − 2 ) [52] with anticommuting γ 5 [53,54] and renormalization. Therefore the Dirac algebra of γ-matrices usually defined in D = 4 needs to be extended to D-dimensional spacetime [55][56][57][58]. A consequence of it is that using Fierz transformations can lead to non-trivial dependences [59][60][61][62] and in general a change in the operator basis valid in four-dimensional spacetime is not allowed in D-dimensional spacetime 2 . In our particular case, the transformation of Eq. (5) into Eq. (6) is allowed because at the order we are working in, only tree level expressions for the Wilson coefficients C 1 and C 2 are needed.
According to the optical theorem the B-hadron decay rate for the inclusive non-leptonic decays can be related to the discontinuity of the forward scattering matrix element which we compute in the HQE up to 1/m 3 , V ub and V q 1 q 2 are the corresponding CKM matrix elements, q stands for a massless quark and are HQET local two-quark 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. 2.3. Note that the QCD spinor b of the bottom quark is replaced by the HQET fermion field h v . They are related as follows

Matching of two-quark operators: computation of C i
In this section we discuss the computation of the matching coefficients of Eq. (8) associated to the two-quark operators Eqs. (9)(10)(11)(12)(13)(14). To that purpose, we compute the Im T in the full theory and then match it with the HQE given in Eq. (8). The Feynman diagrams representing the different contributions are shown in Fig. [1]. We take the c-quark to have mass m c and the u, d, s-quarks as massless. As a consequence, the matching coefficients will be just numbers for b → uūd and will depend on the mass ratio r = m 2 c /m 2 b for b → ucs. For the explicit calculation we find expressions for the quark propagator S F in an external gluon field [63], which automatically generates the proper ordering of the covariant derivatives, which is really convenient and optimal. We proceed as in Ref. [47] where a more detailed explanation is available, as well as explicit expressions for the quark propagators.
Therefore, the problem requires the computation of the imaginary part of two-loop diagrams of the sunset type with zero and one massive lines. We use LiteRed [64,65] to reduce integrals to a combination of a small set (indeed only one) of master integrals which we compute analytically (see Appendix). The Mathematica packages Tracer [66] and HypExp [67,68] are also used to deal with gamma matrices and to perform the expansion of Hypergeometric functions, respectively.
Let's us discuss the peculiarities of each contribution: The color structure only allows the radiation of a single gluon from the u-quark in thebS u b line. Therefore we only need to expand this particular u-quark propagator. The computation is analogous to the semi-leptonic b → u¯ ν case. Due to the gluon emission from (or expansion of) a massless quark propagator the coefficient of ρ D becomes IR divergent, which signals the mixing between the local operators of the HQE, in particular with the four-quark operators.
After renormalization, such IR divergences are canceled by UV divergences developed by the local operators in the effective theory.
• O 2 ⊗ O 2 contribution: This case is identical to the previous one after the exchange of the u and q 2 quarks. Since both quarks are massless, the coefficients are invariant under the exchange of C 1 and C 2 .
• O 1 ⊗ O 2 contribution: It is precisely in this case where we observe the peculiarities of the non-leptonic decay. Unlike in the previous cases, the gluon emission from (or expansion of) the quark propagators have to be taken into account from all internal quark lines. Again, the gluon emission from (or expansion of) massless quark propagators will lead to IR divergences in the coefficient of ρ D which properly cancel after renormalization. It is remarkable that in the b → uūd case, where all internal quark lines are massless, IR divergences cancel giving a finite coefficient.
Overall, the coefficient of ρ LS is IR safe even for massless quarks, so it does not mix with the four quark operators. Note that it can be computed in four dimensions.

Matching of four-quark operators: computation of
In this section, we define the basis of four-quark operators O

The channel b → uūd
The relevant operators in the HQE are where P L = (1 − γ 5 )/2 and P R = (1 + γ 5 )/2 are the left and right handed projectors, respectively. Their matching coefficients read These expressions coincide with the results of Ref. [11].

The channel b → ucs
The relevant four-quark operators in the HQE are O (s) with matching coefficients Again, these expressions coincide with the results of Ref. [11].

Renormalization
The expression for the decay width given in Eq. (7) formally depends on the heavy quark mass m b and the QCD hadronization scale Λ QCD with m b Λ QCD . The HQE given in Eq. (8) takes advantage of the large scale separation and it is organized in such a way that the coefficients are insensitive to the infrared scale Λ QCD , whereas the matrix elements of the local operators are independent of the ultraviolet scale m b .
However, a naive calculation produces both IR singularities in the coefficient functions and UV singularities in the matrix elements of the local operators, which cancel at the end of the computation. We use dimensional regularization to deal with both IR and UV divergences. In such a setup, IR divergences in the coefficient functions signal the UV mixing between the local operators of the HQE under renormalization.
In our computation a naive way of getting the coefficient of ρ D leads to an IR singularity due to the radiation of a soft gluon from masless quark lines. This singularity is canceled by the one-loop UV renormalization (see Fig. [3]) of the four-quark operators which has the general form (in the MS renormalization scheme) where γ (q) 4F i is the mixing anomalous dimension andμ −2 = µ −2 (e γ E /4π) − is the MS Let's define i,q C (q) whereas for b → ucs we obtain γ (q) Likewise, we can write explicitely the counterterm of C D in the MS renormalization scheme.
and for b → ucs it reads where 3 Results for the Wilson Coefficients up to order 1/m 3 b For the matching computation of the Im T we use the HQE given by Eq. (8). However, for the final presentation of the results it is convenient to write it in terms of the local operatorb/ vb defined in full QCD. Its HQĒ can be employed to remove the O 0 operator which has the advantage that the forward matrix element of the leading term is normalized to all orders. We also use the equation of motion (EOM) to remove O v from the expression for the Im T . After these changes we obtain the following expression for the decay rate in terms of the HQE parameters Note that c S = c D = C mag = 1 to the leading order. In the Γ(b → uūd) case we find the following coefficients and for the case Γ(b → ucs) It is remarkable that by explicit calculation we obtain the relations C 0 = C µπ and C µ G = C ρ LS which are a consequence of reparametrization invariance [69]. This is a strong check of the calculation. After using the transformation rules Eqs. (3.20-3.23) in Ref. [70] we can compare our results with Ref. [46], where the coefficients were computed in four dimensions. The coefficients up to O(1/m 2 b ) are in agreement with Ref. [46], where coefficients are compared to the original papers. The coefficient of ρ D for b → ucs is in agreement with Ref. [46], as well. We postpone the comparison for b → uūd to Sec. 3.1, as it requires a non-trivial change of the operator basis of four-quark operators which is related to the treatment of the Dirac algebra in D dimensions [11,47,54].

The canonical basis of four-quark operators
The results presented in Sec. 3 for the b → uūd channel are expressed in terms of the operators which constitute a basis in D-dimensions. However, one may want to use as a basisÕ (u) 4F 2 only, like in Ref. [46]. While the two sets of operators can be related straightforwardly in D = 4, the situation for arbitrary D is more involved. Doing so in D dimensions requires the addition of new operators called evanescent operators, whose choice is not unique. A particular recipe reduces to the substitution [11,54] where E QCD

1,2
are the so-called evanescent operators and a, b are arbitrary numbers, which makes clear why the choice of the evanescent operators is not unique. A conventional choice is a = 4 and b = −4, with D = 4 − 2 . This choice is motivated by the requirement of validity of Fierz transformations at one-loop order [11,56]. We will call the basis fixed by this choice to be the canonical basis of four-quark operators. The complete operator basis now reads In the new basis the imaginary part of the transition operator becomes with the corresponding change in the counterterm of the coefficient C D The operators E QCD

1,2
do not contribute to the anomalous dimension of C D . The difference between the results obtained in the two bases is whereas the difference with the results presented in Ref. [46], where the coefficient was computed in D = 4, is Note that for the canonical choice of the evanescent operators the difference vanishes, so our results are in agreement with Ref. [46]. As we mentioned there is some freedom when choosing the evanescent operators E QCD 1,2 . A different choice of a and b corresponds to the following shift in the coefficient

Numerical evaluation
In this section, we give numerical values for future phenomenological applications and as a first study of the expected size of the new contributions due to ρ D and ρ LS . It is important to keep in mind that the ρ D coefficient depends on the scheme. This scheme dependence is compensated by the matrix elements of the four-quark operators, which are also schemedependent. We  Wilson coefficients of the weak hamiltonian we use the numerical values C 1 (µ b ) = −1.121 and C 2 (µ b ) = 0.274 at µ b = 4.0 GeV given in Ref. [51]. We obtain b → uūd  Table 3: Numerical values for the coefficients of b → ucs.
where we used the abbreviation Γq 1 q 2 to refer to Γ(b → uq 1 q 2 ). Assuming that E QCD the Darwin coefficient gives a correction to the tree level values of the coefficients of the four-quark operators of ∼ 0.4% for b → ucs and of ∼ 0.5% for b → uūd in the canonical basis (we take the largest coefficient of the four-quark operators to compare). The coefficient of the spin-orbit term is less important since it is a factor ten smaller than the coefficient of the Darwin term. However, we have to keep in mind that the final size of the different terms is subject to the evaluation of the matrix elements, which is beyond the scope of this work. Finally, recall that not all four-quark operator contributions are included in the above expressions, but only those which were relevant for the renormalization of the Darwin term coefficient.

Conclusions
In this paper, we have computed the coefficients of the 1/m 3 b hadronic matrix elements ρ D and ρ LS present in the HQE of the non-leptonic width of B-hadrons. More precisely, we have considered the Cabibbo suppressed contributions coming from b → u transitions.
The computation of the ρ LS coefficient is not specially interesting since reparametrization invariance fixes it to be equal to the µ G coefficient, which is known. However, we compute it explicitly as a check.
The computation of the ρ D coefficient by itself is interesting from the theoretical point of view. It requires understanding the operator mixing between O D and dimension six fourquark operators under renormalization. As necessary ingredients, we have computed the leading order coefficients of the dimension six four-quark operators relevant to renormalize the ρ D coefficient. We find agreement with previously known results [11].
We find that the ρ D coefficient depends on the calculational scheme. This not only concerns the use of the MS scheme, but also the choice of the evanescent operators.
The coefficient of the Darwin term turns out to be sizable as it happened in the semileptonic and the non-leptonic (b → c) cases. It gives a correction to the tree level values of the coefficients of the four-quark operator of ∼ 0.4% for b → ucs and of ∼ 0.5% for b → uūd in the canonical basis. However, this statement requires some caveats. Since the coefficient is µ-dependent we can only talk about its size after proper cancellation with the µ-dependence of the matrix elements of the four-quark operators. Indeed, the ρ D term by itself means nothing, but only its combination with the matrix elements of dimension six four-quark operators. The coefficient of the spin-orbit term is less important since it is a factor ten smaller than the Darwin term coefficient.
The HQE has proven to be very successful to explain the lifetimes and lifetime differences of B-hadrons. Concurrently, experimental precision is becoming more impressive over the years. The results presented here may be relevant for a precise determination of B-hadron lifetimes and lifetime differences.
The results obtained here, specially for b → uūd, can also be easily applied to charm decays, where their contribution is expected to be more important due to the HQE for charm has a slower convergence than for bottom. These coefficients might be helpful to clarify the status of the HQE for charm [49], whose validity has been often questioned due to the smallness of the charm quark mass.
The ρ D term could be an important source of SU (3) violation of lifetimes (lifetime differences) [16,47]. That is because ρ D can be related to the matrix elements of fourquark operators through the gluon EOM, whose matrix elements are not SU (3) symmetric (depend on the spectator quark). However, a quantitative study of its impact on lifetime differences needs estimates of SU (3) violation in the matrix elements, which is beyond the scope of this paper.
The coefficients computed in this work were recently determined in Ref. [46] using a different approach. After the proper choice of the evanescent operators (canonical basis) we find agreement with the coefficients presented in this work.
where +iη prescriptions are assumed in the propagators. The following master integrals are needed for the computation J (0, 1, 0, 1, 1