Towards next-to-next-to-leading-log accuracy for the width difference in the $B_s-\bar{B}_s$ system: fermionic contributions to order $(m_c/m_b)^0$ and $(m_c/m_b)^1$

We calculate a class of three-loop Feynman diagrams which contribute to the next-to-next-to-leading logarithmic approximation for the width difference $\Delta\Gamma_{s}$ in the $B_s-\bar{B}_s$ system. The considered diagrams contain a closed fermion loop in a gluon propagator and constitute the order $\alpha_s^2 N_f$, where $N_f$ is the number of light quarks. Our results entail a considerable correction in that order, if $\Delta\Gamma_{s}$ is expressed in terms of the pole mass of the bottom quark. If the $\overline{MS}$ scheme is used instead, the correction is much smaller. As a result, we find a decrease of the scheme dependence. Our result also indicates that the usually quoted value of the NLO renormalization scale dependence underestimates the perturbative error.

For the calculation of Γ 12 one employs an operator product expansion, the heavy quark expansion (HQE) [13]- [16], which results in a systematic expansion of Γ 12 in powers of Λ QCD /m b ∼ 0.1 and α s (m b ) ∼ 0.2. Γ 12 has been calculated to next-to-leading order (NLO) in both Λ QCD /m b [17] and α s (m b ) [9,10,18,19]. The leading-power (i.e. (Λ QCD /m b ) 0 ) term involves two |∆B| = 2 operators (B denotes the beauty quantum number) (1. 4) Here the i, j are colour indices and V ±A means γ µ (1±γ 5 ) while S ±P stands for (1±γ 5 ). The hadronic matrix elements, which must be calculated with non-perturbative methods, are usually parameterized as (1.5) Here f Bs is the B s decay constant and µ 2 = O(m b ) is the renormalization scale at which the matrix elements are calculated. In a lattice-gauge theory calculation µ 2 is the scale at which the latticecontinuum matching is performed. In the expression for Γ 12 the matrix elements of Eq. (1.5) are multiplied by perturbative Wilson coefficients which also depend on µ 2 such that the dependence on the unphysical scale µ 2 cancels from Γ 12 . In the same way the dependence on the renormalization scheme cancels between the Wilson coefficients and B(µ 2 ), B ′ S (µ 2 ). In this paper we use the scheme of Ref. [18].
∆Γ is proportional to m 2 b and the theoretical prediction depends on the renormalization scheme chosen for m b (for a detailed discussion see Ref. [10]) and further on the scale µ 1 = O(m b ) at which the |∆B| = 1 Wilson coefficients are evaluated. Both dependences are unphysical and diminish order-byorder in perturbation theory. At NLO the scheme and scale dependence is still sizable and indicates that higher orders of α s should be calculated. With up-to-date values for quark masses and the elements of the Cabibbo-Kobayashi-Maskawa (CKM) matrix (stated below in Sec. 5) one finds ∆Γ = (1.74 ± 0.24) f 2 Bs B + (0.40 ± 0.05)f 2 Bs B ′ S + (−0.65 ± 0.35) f 2

Bs
(1. 6) in the scheme using the pole mass definition of m b in the prefactor of ∆Γ. Here and in the following the hadronic parameters are understood at µ 2 = m b . The last term in Eq. (1.6) is the Λ QCD /m b correction. If instead the MS scheme is used for m b one finds The errors quoted in the brackets in Eqs. (1.6) and (1.7) are found by varying µ 1 between m b /2 and 2m b . Ref. [11] has quoted all results for the scheme of Eq. (1.7), while in Ref. [12] the average of results in the two schemes has been given. A recent lattice calculation [20] has found Here we have added two errors from different sources in quadrature. Ref. [20] has also calculated some of the matrix elements appearing at order Λ QCD /m b and these results went into the last terms of Eqs. (1.6) and (1.7). With f Bs = 0.224 GeV and neglecting the correlation of the uncertainties in B and B ′ S we find From Eq. (1.9) we observe that the both scale and scheme dependences exceed the uncertainties from the hadronic parameters B and B S . Furthermore, the theoretical uncertainty inferred from these dependences is larger than the present experimental error. This calls for a NNLO calculation of the perturbative coefficients multiplying Q and Q S . In this paper we present the first step in this direction, the calculation of the terms of order α 2 s N f , where N f is the number of quark flavours, neglecting quadratic and higher powers of m c /m b . Eq. (1.9) will further improve from a future calculation of the NLO corrections to the Λ QCD /m b part and progress in the lattice calculations of the hadronic matrix elements appearing in this order. The contributions of order (Λ QCD /m b ) 2 , however, have been estimated to be small [10,21]. The theoretical prediction can be further refined, if ∆Γ is predicted from the ratio ∆Γ/∆M and the experimental value of ∆M , which is proportional to f 2 Bs B. This procedure eliminates the uncertainty associated with B altogether, at the price of making the prediction sensitive to possible new physics in ∆M . From Eqs. (1.6) and (1.7) one realises that the numerically dominant term in ∆Γ/∆M will not contain any hadronic parameter [10]. This feature also alleviates the problem that the lattice-continuum matching is currently only known to NLO. This paper is organized as follows: In the following section we summarize the theoretical framework of the calculation. In Sec. 3 we describe details of the renormalization procedure and the regularization of infrared singularities. We present our analytical results in Sec. 4 and perform a phenomenological analysis in Sec. 5. Finally we conclude. Results for matrix elements and master integrals needed for the calculation are relegated to the appendix.

Theoretical framework
The effective ∆B = 1 weak Hamiltonian, relevant for our calculation, is the following [22] H ∆B=1 Here the i, j are colour indices and summation over q = u, d, s, c, b is implied. V ± A refers to γ µ (1 ± γ 5 ) and S ± P (which we need below) to (1 ± γ 5 ). C 1 , . . . , C 6 and C 8 are the corresponding Wilson coefficient functions. G F is the Fermi constant and V jk denotes an element of the CKM matrix. Cabbibo-suppressed contributions proportional to V * ub V us are neglected in (2.1). To find ∆Γ ≃ 2|Γ 12 | we must calculate where 'Abs' denotes the absorptive part of the matrix element and T denotes time ordering. The HQE expresses Eq. (2.3) in terms of matrix elements of local operators. The leading term (in powers of Λ QCD /m b ) reads Using the notation of Refs. [9,10,18], the coefficients G and G S are further decomposed as Here F and F S are the contributions from the current-current operators Q 1,2 while the small coefficients P and P S stem from the penguin operators Q 3−6 and Q 8 . The coefficients G, G S are calculated by expressing the bilocal matrix elements ("full theory") in terms of the local matrix elements Q , Q S ("effective theory"), the coefficients of the latter are the desired coefficients. Since G, G S are short-distance quantities, this matching calculation can be done order-by-order in perturbation theory, with quarks instead of mesons as external states in Eq. (2.6). The NLO result of Refs. [9,10,18,19] involves Eq. (2.6) at the two-loop level for i, j = 1, 2. The chromomagnetic operator O 8 is proportional to the strong coupling g s , so that for i = 8 or j = 8 a one-loop calculation is sufficient for NLO accuracy. It is further customary to count the small penguin Wilson coefficients C 3−6 as O(α s ) and only one-loop diagrams are considered for i ≥ 3 or j ≥ 3.
The first ingredient of an NNLO result are the Wilson coefficients of the ∆B = 1 weak Hamiltonian in Eq. (2.1). The NNLO Wilson coefficients involve the three-loop anomalous dimension matrix governing the renormalization-group evolution of C 1−6,8 from the electroweak scale down to the scale µ 1 ∼ m b , at which the matrix elements in Eq. (2.6) are evaluated. The NNLO effective hamiltonian has been calculated in Refs. [23,24], albeit in a different operator basis than the one in Eq. (2.2), which is used in the NLO calculations of Refs. [9,10,18,19] and in this paper.
The NNLO contributions presented in this paper all involve a closed quark loop and would be dominant in the case of a large number N f of light quarks. However, the limit N f → ∞ is in conflict with asymptotic freedom of QCD, as the first term β 0 of the QCD β function would change sign. It has been suggested to trade N f for β 0 , so that the α 2 s N f term is replaced by a term of order α 2 s β 0 (naive non-abelianization [25,26]). In some applications this procedure gives a good approximation to the full α 2 s term. However, in quantities involving effective four-quark operators, it is pure speculation whether the original α 2 s N f term or its naively non-abelianized version ∝ α 2 s β 0 approximates the full result in a better way, because neither term cancels the scheme dependence of the operator renormalization. That is, in one scheme the α 2 s N f term may be a good approximant, while in another one the α 2 s β 0 term does better, or neither of them is sensible. For the standard NDR renormalization scheme used by us, e.g. the calculation in Ref. [27] revealed that the α 2 s β 0 term is not a good approximation to the full result. In light of this finding we do not advocate the use of naive non-abelianization in our case. Nonetheless, the α 2 s N f portion of the full NNLO result is gauge invariant and therefore a meaningful quantity. One can also overcome the scheme-dependence issue by only keeping the α 2 s N f terms of the NNLO correction to the RG-improved Wilson coefficients. However, we find that applying this procedure to the known NLO result gives a poor approximation, so that we refrain from using it.
The desired α 2 s N f contribution requires the calculation of the diagrams in Fig. 1. We formally distinguish the charm mass in the lines attached to an effective operator (i.e. to a weak vertex) from that in the charm loop correcting the gluon propagator: The latter give rise to corrections which are linear in m c /m b and we keep a non-zero charm mass in these loops. On the contrary, the dependence on the charm mass arising from the lines in which the charm originates from a weak vertex is only quadratic and we use m c = 0 for these lines. Denoting the MS-renormalized mass of the quark q with m q (µ), where µ is the renormalization scale, we define If the LO and NLO terms are expressed in terms of z, the error associated with the above approximation is of order α 2 s N f z log 2 z. If, however, one usesz instead, the approximation only inflicts an error of order α 2 s N fz and the logarithmic terms α n s z log n z, z = 1, . . . are summed to all orders. This feature has been studied in Ref. [10,28]. The NLO result for ∆Γ expressed throughz is numerically very well reproduced ifz is set to zero in the NLO correction. Since we discard terms of order α 2 sz , one may also expand the z-dependence from the charm quark loop to order z log z and neglect terms of order z and higher. We calculate the tree-loop diagrams with charm loop indeed as an expansion in z, but keep all terms to order z 3 , to check whether the expansion is numerically under good control. Furthermore, a future NNLO calculation keeping higher powers of z terms will benefit from these results.

Renormalization and infrared regularization
In this section we specify our renormalization scheme, present the various counterterms, and clarify the regularization procedure used to isolate infrared (IR) divergences. The latter factorize between the full-theory and effective-theory diagrams (see Fig. 1) and render the desired Wilson coefficients IR-finite.
For the three-loop diagrams involving two insertions of O 1,2 we need C 1,2 at NNLO (i.e. calculated with three-loop anomalous dimensions). The result of Ref. [23] has been transformed to the traditional operator basis in Eq. (2.2) in Ref. [29] and we use the result of this paper. We renormalize the operators in the usual naive dimensional regularization (NDR) scheme. To fully specify the scheme one must further define the evanescent operators [30]. In Ref. [29] the usual NLO definition of these operators has been extended to NNLO in such a way that the diagonal RG evolution of O 2 ± O 1 is maintained at NNLO. For our calculation we must specify the evanescent operators related to Q andQ S : Their operational definition involves the following replacements in the D-dimensional Diracstructures (D = 4 − 2ǫ): These relations, as well as their colour-flipped counterparts, extend the result of Ref. [18] to order ǫ 2 . Formally, the evanescent operator E 1 [Q] (see Ref. [30]) is defined as the difference between the expression on the left and on the right of the arrow in Eq. (3.2), supplemented with the quark field operators on the left and right of the Dirac structures, and analogously Eq.
At NNLO the ǫ 2 terms matter, and these are chosen to preserve the Fierz symmetry, i.e. the two-loop matrix elements of Q andQ S are equal to the matrix elements of the operators obtained from Q and Q S by 4-dimensional Fierz transformations.
In a first step of the calculation the diagrams contributing to Eq. (2.6) generate three effective operators, Q,Q S and Q S = (s i b i ) S−P (s j b j ) S−P . However, one linear combination of Q, Q S , andQ S is 1/m b suppressed [17], so that one can choose any two of them in the leading-power result addressed in this paper. The 1/m b -suppressed operator reads with α 1,2 = 1 at LO. In Ref. [18] it was found that α 1,2 receive corrections of order α s . To our order α 2 s N f and in the scheme defined by Eqs. (3.1) and (3.2) these coefficients read: Here C f = 4/3 is a colour factor and µ 2 is the scale at which the operators in Eq.(3.3) are defined. N H = 1, N V = 1, and N L = 3 count the numbers of b, c, and light (u, d, s) quarks, respectively. The redundant parameters N H,V are introduced for an easier recognition of the various contributions in the formulae for the coefficients. The results for α 1,2 are further expanded in z to the third order. Later we will have to express α s (µ 2 ) in terms of α s (µ 1 ), which occurs in the Wilson coefficients. To this end one can use the following formula: One may freely choose two of the three operators Q, Q S , andQ S . The choice of the basis Q,Q S leads to numerically more stable results [10] than the choice Q, Q S and renders the unknown NLO corrections proportional to R 0 color-suppressed. Nevertheless the NNLO calculation is more convenient in the latter basis and one may easily transform the result between the bases by using Eqs. (3.3) to (3.5).
We next discuss the infrared regularization. For the gluon propagator we use the following expression (similar to the W boson propagator in an R ξ gauge with ξ = 0) where m g is a gluon mass. Our choice of a gluon mass as IR regulator instead of using dimensional regularization has two advantages: In the matching procedure we do not need the ǫ and ǫ 2 parts of NLO and LO Wilson coefficients and the disapperance of m g from the Wilson coefficients provides a non-trivial check of the calculation.
The NLO renormalization constants of the gluon mass and g s in MS scheme read [31,32] δZ For the NNLO calculation we need NLO diagrams with counterterms, so that the full-theory NLO diagrams are needed up to order O(ǫ). For this reason we have extended the calculation of Ref. [18] to order ǫ 1 for m c = 0. Since the two-loop counterterms have 1/ǫ 2 poles, we further need the full-theory LO diagrams to order ǫ 2 . The results of these diagrams can be found in Appendix A.
The NNLO-large-N f piece of the field renormalization constant for the external quark lines is in the 't Hooft-Feynman gauge.
We now turn to the counterterms for the ∆B = 1 operators. The hamiltonian in Eq. (2.1) reads (3.10) The last lines illustrates that one can view Z jk as either renormalising the operator O k or the Wilson coefficient C j . Traditionally the renormalization is attributed to the operator, but we adopt the latter viewpoint, with C j ≡ C ren Writing Z jk = δ jk + δZ jk and expanding δZ jk = αs 4π δZ jk + O(α 3 s ) we find the following counterterms (first calculated in Ref. [33]) at order α 2 s N f : For the penguin-diagram contributions we need the counterterms δZ 2k related to the mixing of O 2 into the four-fermion operators O 3−6 , necessary to renormalize the penguin diagram D 11 . There are two types of contributions. The first type induces the mixing between O 2 and O 3−6 . The non-zero contributions are: In the result for ∆Γ the counterterms in the first line multiply the matrix elements M (1) 62 noted above. In the effective theory the counterterms for gluon mass, strong coupling constant g s , and external fields (b and s) are treated as in the full theory. For the counterterms of the ∆B = 2 operators note that here only the NNLO renormalization constants can contain parts proportional to N f , while the NLO renormalization constants have no pieces proportional to N f . Thus the MS renormalization of the ∆B = 2 operators at order α 2 s N f is trivial, one just has to drop the divergence from the considered two-loop diagrams with quark loop.
4 Results for the coefficients G, G S at order α 2 s N f We first discuss the contributions F ,F S to G,G S with two insertions of O 1,2 (see Eq. (2.5)). We decompose F defined as with an analogous definition of F S,ij . We further write and similarly for F S (z). N H,V,L are defined after Eq. (3.5). The argument of F (2),N H,V,L ij is the ratio z q = m 2 q /m 2 b , where m q is the mass of the quark running in the loop in the gluon propagator, i.e. z q equals 1,z, or 0.
The NNLO functions F (2),N f ij and F (2),N f S,ij for the b quark loop read: (1) = 704 The result for the charm loop quark is expanded in  where P NLO (z) and P NLO S (z) are the NLO results of Ref. [18], while ∆P (z) and ∆P S (z) are the NNLO corrections with z . Since we treat C 3−6 as O(α s ), the latter contain terms of order C 3−6 C 3−6 , α s C 2 C 3−6 , and terms of order α 2 s C 2 2 . The large-N f part of ∆P NNLO (z) is decomposed as with an analogous formula for ∆P NNLO S (z). In the penguin contributions the charm mass on all lines touching O 2 are set to zero, while all other charm loops are kept massive. These include not only the loop in D 11−13 , but also the loops connecting two penguin operators O 3−6 or one penguin operator and a charm-gluon vertex. The latter two contributions appear in counterterm diagrams (to e.g. D 11−13 ) and must be treated in the same way as the diagrams which they renormalize. Consequently, the argument z q (with z q = 1, z, or 0) in ∆P NNLO,N H,V,L (z q ) refers to the mass in the loop of any of these three situations. (At NNLO there are no diagrams with more than one loop.) The results are: In the matching procedure one has to take into account that the operators, couplings and masses on the full-theory side are defined at the scale µ 1 , while the effective operators are defined at the scale m b (m b ) = (4.18 ± 0.03) GeV [37]m c (m c ) = (1.286 ± 0.013 stat ± 0.040 syst ) GeV [38][39][40] m s (m b ) = (0.079 ± 0.002) GeV [20,41]  is done with f Bs = (0.224 ± 0.05) GeV [41]. Subsequently this result is used to rescale B R0 /B from the value in Ref. [20] to the one in the table. B S is larger than B ′ S by a factor of M 2 Bs /(m b +m s ) 2 = 1.588, so that B ′ S /B = 1.83 ± 0.21.
To compare both sides one must choose the same expansion parameter on both sides, e.g. α s (µ 1 ), and use Eq. (3.6) for this. Therefore the α 2 s N f results quoted in this section also contain contributions from the α 1 s parts through Eq. (3.6).

Phenomenology of ∆Γ
In this section we show the impact of the new α 2 s N f terms on ∆Γ s . Our input parameters are collected in Tab. 1. We use the complete NNLO ∆B = 1 Wilson coefficients C 1 , C 2 [29] and the complete NLO expressions for C 3 , ...C 6 , with the numerical values listed in Tab. 2. The α 2 s N 0 f terms of the coefficients inflict a scheme dependence on ∆Γ, which will only be cancelled once the full NNLO calculation is performed. Nevertheless we can study whether the new large-N f terms help to reduce scale and scheme dependences.
The coefficients G = F + P and G S = −F S − P S correspond to the pole scheme for ∆Γ. For the MS scheme we must multiply these coefficients withm 2 b /m pole b and expand this ratio to the order in α s to which G, G s are calculated [10], in our case this is O(α 2 s N f ). In both schemes we usez defined in Eq. (2.7); the transformation from z toz in the NLO formula can be found in Eq. (18) of Ref. [28]. Since we have set z = 0 in the charm lines attached to weak vertices, no NNLO corrections to the transformation occur.
We further must calculate m pole b fromm b and we use the full 2-loop result for this [34][35][36]. This is a reasonable approach, if the missing α 18 GeV using α s (m b ) = 0.226 (implementing the formula of Ref. [42] with QED effects set to zero) and the matching scale µ 0 = M W . We have used Ref. [29] to compute C  is taken from the calculation in a different basis [24] and the quoted value therefore neglects a numerically small contribution from an evanescent operator. in the MS scheme, because the sizable µ 1 dependence of the prefactorm b (µ 1 ) 2 cancels nicely with the µ 1 dependence of G,G S . In our partial NNLO result this efficient cancellation is less pronounced than in the NLO result of Eq. (1.9). To be conservative, we therefore use a different approach in this section: We keepm b (m b ) 2 fixed and, for consistency, also eliminate the log(µ 1 /m b ) terms related to the running ofm b from G,G S . This leads to a larger µ 1 dependence at NLO. where the scale dependence is calculated by varying µ 1 between m b /2 and 2m b and for the quoted central values of ∆Γ we took µ 1 = m pole b and µ 1 =m b for the pole and MS schemes, respectively. Unlike in Eq. (1.9) other sources of error are neglected here. The µ 1 dependence is plotted in Fig. 2.
We observe that the partial NNLO corrections calculated in this section decrease the scheme dependence and give preference to the NLO result in the MS scheme. The result also suggests that in Eq. (1.9) the µ 1 dependence is underestimated and that the partial NNLO calculation does not reduce the scale dependence to a satisfactory level.
We have discussed the naive non-abelianization approach (NNA) in Sec.

Conclusions
We have calculated the contributions of order α 2 s N f to the width difference in the B s −B s system in an expansion in m c /m b , neglecting terms of order (m c /m b ) 2 and higher. This calculation has involved three-loop massive master integrals with two mass scales. We find a larger correction for the decay width difference in the pole scheme and only a minuscule correction for the MS scheme. As a result, the scheme dependence reduces considerably and we advocate the use of the NLO numerical values in Eq. (5.4).

A Full-theory matrix elements
In this section we collect the needed unrenormalized LO and NLO matrix elements to order ǫ 2 and ǫ, respectively. We decompose the matrix element as where the first term denotes the contribution with two insertions of the current-current operators O 1,2 and the second term comprises the diagrams with at least one penguin operator. Recall that we count C 3−6 as order α s , so that one loop less is needed for M peng compared to M cc . We expand M cc,peng = M Here and in the following ... (0) denote tree-level matrix elements and C b k = j C j Z jk are bare Wilson coefficients (see Eq. (3.10)).
We decompose the NLO diagrams according to the diagrams in Fig. 1  The sum of the full-theory non-penguin NLO diagrams amounts to

and the Wilson coefficients as
and for the penguin diagrams

A.2 Penguin operators
For the matrix elements with two QCD penguin operators we write As usual we expand M jk as M jk = M , (A.13) . (A.14) We only need the LO contributions M (0) jk for j, k ≥ 3 : We need the coefficient functions up to ǫ 2 , findinĝ The M

B Results of master integrals
We have reduced the Feynman diagrams shown on Fig.1 to master integrals by means of the program FIRE [43]. For the full-theory diagrams we have calculated the absorptive part of master integrals, i.e. the 2-, 3-, 4-particle cuts with a massive c-, b-quarks in the closed fermion loop, with a massive gluon in infrared singular diagrams and a massless c-quark in the weak loop, using formulas for phase space integrals derived in [44]. For some integrals, with massive charm, we used a Mellin-Barnes representation [45] and expanded in terms of the small parameter z = m 2 c /m 2 b . The master integrals, which include m g , are expanded over z g = m 2 g /m 2 b . The results of master integrals have been checked numerically by means of the program SecDec-3 [46].
The results for diagrams with massless u-, d-, s-quarks in the closed fermion loop are obtained by taking the limit m c → 0 in the results with a massive c-quark in the closed fermion loop.
From the results below one can see that the first three orders in the expansion over z already exhibit a good convergence.
Our convention for the loop measure is Some of the integrals in the following subsections have more than one cut (e.g. 2, 3 and 4 particle cuts). The following subsections quote the results of the various cuts. We write Im = Im (2) + Im (3) + Im (4) to separate the contributions from these cuts.