Contribution of the Darwin operator to non-leptonic decays of heavy quarks

We compute the Darwin operator contribution ($1/m_b^3$ correction) to the width of the inclusive non-leptonic decay of a $B$ meson ($B^+$, $B_d$ or $B_s$), stemming from the quark flavour-changing transition $b \to q_1 \bar q_2 q_3$, where $q_1,q_2 = u, c$ and $q_3 = d, s$. The key ideas of the computation are the local expansion of the quark propagator in the external gluon field including terms with covariant derivative of the gluon field strength tensor and the standard technique of the Heavy Quark Expansion (HQE). We confirm the previously known expressions of the $1/m_b^3$ contributions to the semi-leptonic decay $b \to q_1 \ell \bar \nu_\ell$, with $\ell = e, \mu, \tau$ and of the $1/m_b^2$ contributions to the non-leptonic modes. We find that this new term can give a sizeable correction of about $- 4 \, \%$ to the non-leptonic decay width of a $B$ meson. For $B_d$ and $B_s$ mesons this turns out to be the dominant correction to the free b-quark decay, while for the $B^+$ meson the Darwin term gives the second most important correction - roughly 1/2 to 1/3 of the phase space enhanced Pauli interference contribution. Due to the tiny experimental uncertainties in lifetime measurements the incorporation of the Darwin term contribution is crucial for precision tests of the Standard Model.


Introduction
The total decay rate of heavy hadrons can be described by the Heavy Quark Expansion (HQE), whose history goes back to the work of Shifman and Voloshin in the 1980ies [1,2] 1 , as the decay of the free heavy quark plus corrections that are suppressed by inverse powers of the heavy quark mass. Since the b-quark mass is large compared to the typical hadronic scale, the corrections are expected to be small and hence, back in 1986, the following b-hadron lifetime ratios were expected [2] τ (B s ) τ (B d ) (1. 2) There has also been considerably progress on the theory side. The total width of a B meson with mass m B and four-momentum p µ B is given by where H eff represents the effective weak Hamiltonian [5] describing all possible b-quark decays. So far the dependence on the spectator quark is only due to the different values of the matrix elements and the coefficients Γ 0,2,3 are independent of the quark content of the B meson. The four-quark operators are phase-space enhanced (as indicated by the factor 16π 2 ) and first arise at order 1/m 3 b . The possible topologies, specifically weak annihilation (WA), Pauli interference (PI) and weak exchange (WE), see Fig. 4, imply that the short distance coefficients are now dependent on the spectator quark. The dimension-seven operatorsÕ 7 contain one covariant  derivative compared toÕ 6 . Due to the larger phase space, it is typically expected thatÕ 6 gives the dominant contribution to the lifetime ratios [1-3, 6, 7]. Currently Γ 0 is known at NLO-QCD [8][9][10][11][12][13][14][15] for non-leptonic decays. NNLO-QCD corrections have been computed for semi-leptonic decays [16][17][18][19][20][21][22][23][24][25] and for non-leptonic decays the massless case was determined in full QCD (i.e. no effective Hamiltonian was used) in Ref. [26]. Γ 2 was determined at LO-QCD for both semi-leptonic and non-leptonic decays [27][28][29][30] and we confirm these results. For the semi-leptonic modes even NLO-QCD corrections are available [31][32][33]. Γ 3 was first computed at LO-QCD in Ref. [34] and recently the NLO-QCD corrections were determined in Ref. [35], both for the semi-leptonic case only.Γ 3 is known at NLO-QCD for B-meson lifetimes [36,37] and D meson lifetimes [38], whileΓ 4 is only known at LO-QCD [39]. This work presents the first determination of Γ 3 for non-leptonic decays. An interesting subtlety of the computation is mixing between four-and two-quark operators. Namely, at dimension-six the renormalised one-loop matrix elements of the operatorsÕ 6 contribute to the coefficient of the Darwin operator through the diagram in Fig. 5, ensuring the cancellation of the infrared (IR) divergences that otherwise would appear in Γ 3 . This feature has been intensively discussed for semi-leptonic decays [40][41][42][43] -also under the name of "intrinsic charm", see e.g. Refs. [44,45]. Finally, for numerical analysis, the values of the non-perturbative matrix elements O Y are needed. O 5 and O 6 can be extracted from fits to the semi-leptonic spectrum for the case of B d and B + mesons, see e.g. Ref. [46]. In the literature one can also find lattice determinations [47][48][49][50][51] and sum rule estimates [52][53][54] for these parameters. For the SU (3) F violating ratios B s |O 5,6 |B s / B d |O 5,6 |B d one can use the theory estimates from Ref. [55]. The matrix elements of the four-quark operators Õ 6 have been determined by Heavy Quark Effective Theory (HQET) sum rules [56]. Violations of SU (3) F are expected to yield visible effects and a calculation of these corrections with HQET sum rules -following Ref. [57] -is currently been performed [58]. Corresponding lattice results for the matrix elements of the four-quark operators would be highly desirable. Taking all the currently known contributions into account we arrive at significant improvements compared to the pioneering work in 1986, and the 2019 status of lifetime predictions reads [3,56] τ (B s ) τ (B d ) • It was found (see e.g. Refs. [34,41,55,[62][63][64]) that the 1/m 3 b correction in semileptonic inclusive decays B → X c ν are of a similar size as the 1/m 2 b terms due to enhanced Wilson coefficients. This could lead to a visible effect in τ (B s )/τ (B d ), in particular if all other contributions are cancelling to a high degree.
• As indicated in Eq. (1.2) the experimental precision achieved so far is very high, enabling thus precision tests and an even higher precision seems to be achievable, see e.g. the two most recent results from LHCb [65] and ATLAS [66], which interestingly differ significantly. The theory precision should of course cope with these experimental advancements.
• According to the above arguments the lifetime ratio τ (B s )/τ (B d ) might thus provide a unique opportunity to test directly higher orders in the HQE that are otherwise just noise compared to the leading, numerically dominant contributions. In that respect this might also increase our insights on the assumptions of quark hadron duality, severely questioned around 20 years ago -mostly due to a very low experimental value of the Λ b lifetime, which seemed to be in severe conflict with the HQE expectation from Eq. (1.1): the 2002 HFAG average gives e.g. τ (Λ b )/τ (B d ) = 0.798 (34), more than 5σ away from the current experimental value given by Eq. (1.2). This signal for a violation of quark hadron duality was a clear false alarm. Moreover, the first measurement of the decay rate difference of neutral B s mesons, ∆Γ s , has excluded large duality violating effects, see e.g. the discussion in Refs. [67,68]. Smaller effects of duality violation can still appear and their potential size can be constrained by comparing experiment and theory for as many HQE observable as possible with a high precision, see e.g. Ref. [69].
• Assuming the validity of quark hadron duality the τ (B s )/τ (B d ) lifetime ratio can be used to constrain invisible or hardly detectable B-meson decays, like B s → τ τ , at the per mille level, see e.g. Refs. [70,71] -see also Ref. [72] for an alternative way to constrain the potential size of the bsτ τ couplings.
We will not present updated lifetime ratio predictions in this work, but we will postpone a new numerical study until the SU (3) F violation ratio of the Bag parameters is available, see Ref. [58]. The paper is structured as follows: in Section 2 we outline the main ingredients of the calculation: in Section 2.1 we explain how the double insertion of the effective Hamiltonian and the subsequent expansion in the external soft gluon field are performed. Section 2.2 is devoted to describe mixing between two-and four-quark operators at order 1/m 3 b , which guarantees the cancellation of the IR divergences. Our results are presented in Section 3 and we conclude in Section 4. The expansion of the quark propagator in the external gluon field is discussed in more detail in Appendix A, supplement material to Section 2.2 and Section 3 can be found in Appendix B and Appendix C, respectively, and for completeness we also present the non-leptonic results for Γ 2 in Appendix D. According to the optical theorem, the total width for the inclusive non-leptonic decay of a B meson induced at the quark level by the flavour-changing transition b → q 1q2 q 3 (with q 1 , q 2 = u, c and q 3 = d, s), can be computed from the discontinuity of the forward scattering matrix element: The effective Lagrangian L eff (x) reads where G F is the Fermi constant, V qq are the elements of the Cabibbo-Kobayashi-Maskawa (CKM) matrix [73,74], C 1,2 (µ 1 ) denote the Wilson coefficients at the renormalisation scale µ 1 ∼ m b and Q 1,2 are the effective ∆B = 1 four-quark operators 2 with i, j standing for colour indices and Γ µ ≡ γ µ (1 − γ 5 )/2. In our notation, Q 1 is the colour singlet operator and Q 2 the colour rearranged one, opposite to e.g. the notation in Ref. [5]. Three combinations of operators enter the decay rate in Eq. (2.1): We use the completeness property of the colour matrices t a ij 3 t a ij t a lm = where the colour octet operator T is given by Considering only the two-quark operators contribution, the decay width takes the form Note that in the case of Γ the two colour matrices t a appear on the r.h.s of Eq. (2.9). The corresponding expression for Γ (2q) 22 is obtained from that of Γ (2q) 11 by replacing q 1 ↔ q 3 as it follows by Fierz-transforming the operator Q 2 given in Eq. (2.4). To compute dimensionsix contributions we need to expand first each of the quark propagators up to one covariant derivative of the gluon field strength tensor, respectively defined as In Appendix A we derive the Fourier transform of the quark propagator in the soft external gluon field In the previous equation,G ρµ = (1/2) ρµση G ση denotes the dual field strength tensor, with ρµση being the Levi-Civita tensor and . From Eq. (2.11) it follows that the quark propagator can be split up according to its colour structure as where S (0) (k) denotes the free quark propagator and S (1) a (k) absorbs contributions with one gluon field (for more details see Appendix A, also the original Refs. [28,29]). This representation allows a straightforward treatment of colour in Eq. (2.9). In case of Q 1 ⊗ Q 1 contribution, the colour flow factorises between the (bq 1 )-quark line and the (q 2 q 3 )-loop, therefore one can only expand the propagator of q 1 up to terms linear in t a , see l.h.s. of Fig. 2. For the Q 2 ⊗ Q 2 combination the trace over colour indices involves the q 1 -and q 2 -quark propagators and only the gluon radiation from q 3 is non-vanishing, see r.h.s. of Fig. 2. Finally, in the case of the Q 1 ⊗ T contribution, the colour flow forces the gluon to be radiated off the (q 2 q 3 )-loop. Any interference term would be O(1/m 4 b ), and it is thus sufficient to expand independently each of the two quark propagators in the loop (Fig. 3). Using translation invariance of the quark field to write b(x) = e −ipx b(0), with p µ being the b-quark four-momentum, and after performing x-and momentum integration, we obtain twoloop tensor integrals of possible rank r = 1, . . . , 4. These are decomposed in terms of rank r tensors built with the metric tensor g µν and the external momentum p µ . The corresponding coefficients represent scalar integrals of the type (2.14) In the previous equation, is the Källen function. One can easily compute analytically the integral in Eq. (2.14) for two non-vanishing masses, while in the case of three non-vanishing masses the complexity highly increases and the solution involves elliptic functions, see e.g. Refs. [78][79][80][81]. We emphasize that we set d = 4 from the beginning since the discontinuity of the diagrams at LO-QCD does not develop any ultraviolet divergences. On the other side, the gluon emission off a light quark propagator (q = u, d or s) gives rise, at dimension-six, to infrared logarithmic divergences of the type log(m 2 q /m 2 b ) which are cancelled by the renormalised one-loop matrix element of the corresponding four-quark operators, as shown in detail in Section 2.2. Here we limit ourselves to state that for the computation of most of the diagrams of interest it is always possible to set one mass to zero and obtain an analytic expression for all the master integrals. This is not the case only for the gluon emission from the s quark in the b → ccs transition, where we need to keep all masses finite in order to regularise the infrared divergence.
After computing the two-loop integrals we are left with the following matrix elements where F(p), F µν (p) and F µνρ (p) are functions of the quark masses and of the external momentum p µ . Following the standard technique of the HQE we decompose the b-quark momentum as where v µ represents the four-velocity of the B meson and D µ accounts for the soft interaction with the spectator quark. We have used the phase redefinition where S n is the symmetric group of permutation of n elements. This leads to a systematic expansion of Eq. (2.9), schematically where a, b µ , c µν , d µνρ are now only functions of the quark masses and of the four-velocity v µ . Finally, the matrix elements in Eq. (2.21) admit a series expansion in powers of 1/m b given explicitly in Ref. [62] 6 , from which we can readily obtain the coefficients of µ 2 π , µ 2 G and of ρ 3 D , ρ 3 LS , which are defined as with σ µν = (i/2)[γ µ , γ ν ]. As already stressed, the numerical values of the non-perturbative parameters above depend on the spectator quark in the B meson.
We then arrive at the following form of Eq. (2.8)

24)
6 For reference we quote the expansion of the dimension-three matrix element up to order 1/m 3 b , Figure 4: One-loop diagrams corresponding, from left to right, to the WA, PI and WE topology. can be found in Appendix D. We point out that the adopted definition for the non-perturbative parameters in Eqs. (2.22), (2.23) implies that the coefficient of ρ 3 LS is found to vanish for all the ∆B = 1 operator combinations.

Role of the four-quark operators
At order 1/m 3 b the gluon radiation from a light quark leads to IR divergences, namely where R (q 1q2 q 3 ) nm are finite functions and D (q 1q2 q 3 ) nm log(m 2 q /m 2 b ) absorb the contribution of the divergent logarithms, the latter are listed in Appendix B. As discussed already in Ref. [40], logarithmic infrared divergences signal mixing between operators of the same dimension. To see this in more detail, we start again from Eq. (2.1), 7 with the transition operatorT being the time-ordered product of the double insertion of the effective Lagrangian. The expansion in inverse powers of m b allows to expressT in terms of local new effective operators. At dimension-six one has 8 :

28)
7 For brevity we omit the indices (q 1q2 q 3 ) and (nm), although they must be always understood. 8 In principle the sum in Eq. (2.28) includes also the c-quark but for m b ∼ m c Λ QCD the four-quark operatorsÕ (c) 6,i are not relevant for the further discussion, see for more details Ref. [44].
with the corresponding coefficient C D (µ 1 ) -see Eq. (2.24). Since we are only working at LO-QCD, we do not find any explicit µ 1 -dependence at this stage. The l.h.s of Eq. (2.28) receives also contributions by the one-loop diagrams depicted in Fig. 4, corresponding to the weakannihilation, Pauli interference and weak-exchange topologies. The coefficients C (q) 6,i (µ 1 , µ 0 ) do not develop any divergences at LO-QCD [6,7] hence there is no explicit µ 1 -and µ 0 -dependence; it is however present in the NLO-QCD corrections determined in Refs. [36,37,82]. We refer to Appendix B for explicit expressions at LO-QCD. Integrating out these one-loop diagrams leads to the following ∆B = 0 four-quark operators: 10 The four-quark and the Darwin operators mix under renormalisation. Namely, the perturbative one-loop matrix elements ofÕ  instance Refs. [42,83]. Once renormalised (we adopt the MS scheme to remove the 1/ pole together with the ln(4π) − γ E factor), the above one-loop matrix elements also contribute to the coefficient C ρ D . At the matching scale µ 0 = m b one obtains where the coefficient C D has the divergent logarithmic dependence shown in Eq. (2.26). From the renormalised one-loop matrix elements of the four-quark operators obtained from Eq. (2.32), it follows that on the r.h.s of Eq. (2.33) all the logarithms log(m 2 q /m 2 b ) cancel exactly, leaving C ρ D (µ 1 , m b ) free of any IR divergences. Finally, using Eq. (2.32), one can read the anomalous dimension matrix of the dimension-six operators (at LO-QCD) and solve the corresponding renormalisation group equations (RGEs) to run the coefficients C ρ D (µ 1 , m b ) down to the scale m b ≥ µ 0 Λ QCD : The IR-finite coefficient C ρ D (µ 1 , µ 0 ) will now get an explicit µ 0 dependence from the perturbative matrix element Õ (q) 6,i (µ 0 ) ren even if we are only working to LO-QCD. Our results are presented in the next section. The remaining dimension-six contribution to Γ NL (B) (see Appendix C for the LO-QCD expressions) reads where the matrix element of the four-quark operators can be parametrised as ) 2 and f B being the decay constant of the B meson. We have separated the contribution due to the valence and non-valence quarks, where B (q) i is non vanishing only for q equal to q , the spectator quark in the B meson, and is expected to be of order one, see e.g. Ref. [56], while τ (q) i accounts for the effects of an "intrinsic" q quark [43,44] and is expected to be small. Its numerical value can be estimated via e.g. the calculation of the so-called eye contractions in the non-perturbative determination of the matrix elements, see Ref. [58].

Results
The contribution of the Darwin operator to the inclusive non-leptonic decay b → q 1q2 q 3 is presented in the following form ρ D ,12 = 2 3 −9 + 12 1 − 3ρ 2 + ρ 3 log(ρ) (3.13) The dimensionless parameter η = m 2 q /m 2 b and the master integral M 112 is defined as 11 (3.14) 11 An explicit analytic expression for M 112 (ρ, η) − 1 − 4ρ log(η) η→0 has been found in Ref. [84].  Note, that because of m d = m s = 0 the following relations hold: The relative effect of the Darwin term with respect to the corresponding partonic-level contribution is given by In Fig. 6, these ratios are plotted as functions of ρ for all the colour structures and the four modes, using for reference the values µ 0 = m b , m b = 4.5 GeV and ρ 3 D = 0.2 GeV 3 . Fig. 7 shows the total relative contribution for each mode, namely (3.16) As one can see, the Darwin operator leads to sizeable corrections of the order 1 − 7 % (for ρ = 0.05) to the b → q 1q2 q 3 decay width.

Discussion and conclusion
This work presents the first determination of the Darwin term contribution to the nonleptonic decay b → q 1q2 q 3 . Using the expansion of the quark propagator in the external gluon field together with the standard technique of the HQE allows to obtain a systematic expansion in inverse powers of the heavy quark mass m b . At order 1/m 3 b operator mixing ensures that the IR divergences arising from the expansion of the light quark propagators cancel, introducing though scale-dependence in the coefficient of ρ 3 D , with the divergent log(m 2 q /m 2 b ) being in fact replaced by log(µ 2 0 /m 2 b ) 12 . Preliminary numerical analysis reveals that this contribution is sizeable. For illustration purposes, we show the total non-leptonic decay width of B + , B d up to dimension-six and just at LO-QCD, setting µ 1 = m b and µ 0 = 1 GeV for the dimension five and the Darwin term, as well as µ 0 = m b for the four quark contributions. Using inputs for the quark masses and HQE parameters from Refs. [46,56,58] (in the kinetic scheme) we get: is given by Eq. (2.25). In Eq. (4.1) the effect of the different non-perturbative parameters is shown separately. We find that the new contribution due to the Darwin operator is significant and larger than the dimension-five and the weak-exchange contributions, while in the case of B + the Pauli interference term still gives the dominant correction. The inclusion of the Darwin contribution could thus lead to a sizeable modification of the theory prediction for the lifetime ratio τ (B s )/τ (B d ) compared to the high experimental precision, which reaches the accuracy of the four per mille level. An updated theoretical analysis of the lifetime ratio τ (B s )/τ (B d ) is postponed to the near future when the SU (3) F violation ratio of the Bag parameters and the matrix elements of the Darwin operator will become available.
A Expansion of the quark propagator in the external gluon field Following Refs. [28,40], soft gluon interactions of the quark field can be accounted using the background field method. Assuming the b quark embedded in a weakly changing gluon field allows to systematically expand the quark propagator as a series in the gluon field strength tensor G µν = −i [iD µ , iD ν ]. In the case of a massive quark the expansion of the propagator up to terms linear in G µν is given in Ref. [28], while an expression including terms proportional to D ρ G µν = −[iD ρ , [iD µ , iD ν ]] can be found in Ref. [85] 13 . To compute the propagator, one starts from the Green-function equation, which admits a solution in form of the perturbative series where S (0) (x − y) is the free-quark propagator and S (1) (x, y) the first order correction: Using the the Fock-Schwinger gauge i.e. x µ A µ (x) = 0, the gluon field is expressible directly in terms of the gluon field strength tensor, see Ref. [40] for a detailed derivation: Expanding G ρµ (αz) around z = 0 and taking into account that in the Fock-Schwinger gauge z µ ∂ µ = z µ D µ , yields: Substituting the previous expression in Eq. (A.3) and setting y = 0 one obtains where the quark propagator in momentum space reads In the above,G ρµ = (1/2) ρµση G ση and ρµση is the Levi-Civita tensor, while the ellipsis stands for terms with higher derivatives as well as higher powers of G µν . In the limit m → 0, Eq. (A.7) correctly reproduces the massless expression given in Ref. [40]. Finally, we emphasise that the Fock-Schwinger gauge breaks explicitly the translation invariance of S(x, y), namely: For completeness, we also present the equivalent representation of Eq. (A.6) in coordinate space where K 0,1,2 (z) are the modified Bessel function of the second kind.

B Complementary material to Section 2.2
The divergent coefficients D  The discontinuity of the WA, PI and WE diagrams in Fig. 4 at LO-QCD, respectively is where z i = m 2 q i /m 2 b and the four-quark operatorsÕ

D Coefficients of the dimension-three and chromo-magnetic operators
Here we present the analytic expressions for the coefficients of the dimension-three and chromo-magnetic operators intoduced in Eq. (2.24). They read respectively