Two-loop QCD penguin contribution to the width difference in $B_s-\bar{B}_s$ mixing

We consider two-loop QCD corrections to the element $\Gamma_{12}^q$ of the decay matrix in $B_q-\bar{B}_q$ mixing, $q=d,s$, in the leading power of the Heavy Quark Expansion. The calculated contributions involve one current-current and one penguin operator and constitute the next step towards a theory prediction for the width difference $\Delta\Gamma_s$ matching the precise experimental data. We present compact analytic results for all matching coefficients in an expansion in $m_c/m_b$ up to second order. Our new corrections are comparable in size to the current experimental error and slightly increase $\Delta\Gamma_s$.


Introduction
In particle collisions B q mesons, where q = d, s labels the flavour of the light valence quark, are produced as flavour eigenstates. This means that they are either meson or antimeson, with beauty quantum number B = 1 or B = −1, respectively. Subsequently, this pure B q orB q state evolves into a quantum-mechanical superposition of B q andB q following the time evolution of damped oscillations. Two accidental features of the Standard Model (SM) permit the precise study of B q −B q oscillation in modern experiments: First, the smallness of the element V cb of the Cabibbo-Kobayashi-Maskawa (CKM) matrix implies a large B q lifetime of around 1.5 ps, which makes decay-time dependences experimentally observable. Second, the heaviness of the top quark enhances the B q −B q mixing box diagram, which governs the B q −B q mixing amplitude, to a level that the oscillation frequency is in the same ballpark as the B q lifetime.
B q −B q mixing is described by the 2 × 2 matrix M q − iΓ q /2 with the hermitian mass and decay matrices M q and Γ q , respectively. Diagonalizing M q − iΓ q /2 leads to a "heavy" (H) and a "light" (L) mass eigenstate which are commonly denoted by B q H and B q L , respectively, and have masses M H,L and widths Γ H,L . The oscillation phenomena involve the three quantities |M q 12 |, |Γ q 12 | and arg(−M q 12 /Γ q 12 ) which are related to the experimentally accessible quantities where the CP asymmetry in flavour-specific decays, a q fs , is typically measured in semileptonic decays. ∆M q and ∆Γ q are related to the elements of the mass and decay matrices as In the Standard Model (SM) the phase between −Γ q 12 and M q 12 is small, so that a q fs is much smaller than ∆Γ q /∆M q and ∆Γ q ≃ 2|Γ q 12 |. M q 12 is a ∆B = 2 amplitude probing virtual effects of physics beyond the SM (BSM physics) up to mass scales of several 100 TeV. By contrast, Γ q 12 is sensitive to new physics in ∆B = 1 transitions. While Γ q 12 probes much lower scales than M q 12 , it is instead sensitive to effects of feebly coupled BSM particles which are light enough to be produced in B q decays. Such particles are predicted in theories addressing the strong CP problem [1,2] or as members of the dark sector, see e.g. Ref. [3] for a baryogenesis mechanism utilising B q −B q mixing and B q decays into dark matter.
In this paper we calculate QCD corrections to Γ q 12 in the SM needed to better predict both ∆Γ q /∆M q and a q fs . Currently, better theory predictions are needed in the case of [5] Furthermore, there is steady progress with measurements of ∆Γ s at LHCb [6], CMS [7], and ATLAS [8].
For the calculation of Γ q 12 one employs a special operator product expansion, the Heavy Quark Expansion (HQE), which treats the b quark mass m b as a hard scale. In this way one expresses Γ q 12 as a simultaneous expansion in Λ QCD /m b and α s (m b ). Each term of the 1/m b expansion involves perturbative coefficients multiplying hadronic matrix elements of local ∆B = 2 operators. Next-to-leading logarithmic order (NLO) QCD corrections at leading power in 1/m b have been computed in Refs. [9][10][11][12]. The 1/m b contribution is known to leading order in α s [13]. The uncertainty resulting from the truncation of the perturbative series of the currently known SM prediction for ∆Γ s is larger than the experimental error in Eq. (3), which calls for the calculation of higher-order QCD contributions.
First steps towards next-to-next-to-leading order (NNLO) have been undertaken in Ref. [14] where the fermionic corrections of order α 2 s N f , where N f = 5 is the number of active quark flavours, have been computed including linear terms in the expansion in m c /m b . Note that this calculation cannot be used to obtain a q fs , which is proportional to m 2 c /m 2 b . In this paper we denote any O(α s ) contribution to Γ q 12 as "NLO", irrespective of the Wilson coefficients involved. This complies with the commonly used notation in connection with higher-order QCD calculations, but differs from the language used in previous papers on Γ q 12 , in which the small ∆B = 1 penguin Wilson coefficients C 3−6 are counted as O(α s ). In order to match the precision of the experimental value in Eq. (3) one needs the yet unknown complete NNLO corrections proportional to two factors of the current-current Wilson coefficients C 1,2 , while the contributions proportional to C 1,2 C 3−6 are only needed at NLO. In Ref. [15] for the first time penguin contributions have been considered beyond LO, presenting the terms proportional to C 1,2 C 3−6 α s N f .
In this paper we present the QCD corrections to all penguin contributions proportional to the product of C 1,2 with one of C 3−6 in an expansion in where the superscript "OS" refers to the on-shell (or pole) scheme, i.e. two-loop contributions of order O(α s ). Thus this is a step towards the completion of the NLO prediction of Γ q 12 , which is a necessary preparation for NNLO. This calculation is more convenient in the "CMM" operator basis of Ref. [16], which avoids problems in connection to γ 5 . We also adopt this basis in the calculation presented in this paper. As a byproduct we reproduce the NLO result for the contribution with two copies of C 1,2 of Refs. [9][10][11][12] (expanded in z) after transforming the ∆B = 1 Wilson coefficients to the CMM basis, which is a powerful check of our calculational set-up.
The paper is organised as follows: Sec. 2 introduces the ∆B = 1 and ∆B = 2 operator bases employed by us, Sec. 3 and Appendix A present the methodology of our calculation, Sec. 4 contains the results, and we conclude in Sec. 5.

Preliminaries
The effective |∆B| = 1 weak Hamiltonian in the CMM operator basis [16] reads: where λ s a = V * as V ab , a = u, c, t, contains the CKM matrix elements and λ t = −λ c − λ u . For definiteness we specify to b → s decays relevant for B s −B s mixing. The corresponding expressions for B d −B d mixing are trivially found by replacing V as with V ad . G F is the Fermi constant and the dimension-six ∆B = 1 operators are given by where q L = P L q with P L = (1 − γ 5 )/2. Q (u) 1 and Q (u) 2 are the current-current operators describing the W -mediated tree-level decay of the b quark including QCD effects. Q 3 , . . . , Q 8 are four-quark penguin operators. We list the operator Q 8 (with σ µν = i[γ µ , γ ν ]/2) for completeness; it does not enter the calculations in this paper. g s is the strong coupling constant and G a µν denotes the gluon field strength tensor. In Eq. (6) the sum over q runs over all five quark fields u, d, s, c or b. For our calculation we also need the following evanescent operators [16] and the counterparts of E 1 [Q 1,2 ] with one or both c replaced by u. The ∆B = 1 operators in Eqs. (6) and (7) destroy a b ands quark while creating ab and s quark and thereby describe the transition of aB s ∼ bs into a B s ∼bs meson. The corresponding Feynman diagrams have incoming b quark and outgoing s quark lines.
Using the Hamiltonian in Eq. (5) the width difference ∆Γ ≈ 2|Γ 12 | is obtained from where "Abs" stands for the absorptive part and T is the time ordering operator. Γ s 12 encodes the information of the inclusive decay rate into final states common to B s andB s and Eq. (8) employs the optical theorem to relate Γ s 12 to theB s → B s forward scattering amplitude.
It is convenient to decompose Γ s 12 as [9] where in the practical calculation the quantities Γ ab 12 are considered. The Heavy Quark Expansion (HQE) allows us to express the quantities Γ ab 12 in Eq. (9) in terms of infrared-safe perturbative coefficients and hadronic matrix elements of ∆B = 2 operators. To leading power in 1/m b one only needs two ∆B = 2 operators, which are conveniently chosen as with colour indices i, j. At intermediate steps of the calculation one also encounters and operators with more than two Dirac matrices on both quark lines. Q S can be traded for Q, Q S , and an operator R 0 describing 1/m b -suppressed contributions to Γ s 12 [13], with the usual ǫ = (4 − D)/2 of dimensional regularisation. The operators on the RHS are understood to be expressed in terms of the minimal physical basis Q , Q S , e.g. Q is to be read as Q + E 2 . The choice of the O(ǫ) terms in the coefficients affect the expressions of the renormalised coefficients H ab , H ab S of Q , Q S [18]. That is, their specification is part of the definition of the renormalisation scheme of the operators (along with the MS prescription and the use of anticommuting γ 5 ). Our definitions in Eq. (13) ensure that the coefficients do not depend on the Fierz arrangement [17,18], i.e. a fourdimensional Fierz transformation of Q, Q S does not change C and C S .
It is thus possible to write Γ ab 12 in Eq. (9) as with z defined in Eq. (4). The ellipses denote higher-order terms in Λ QCD /m b . The matching coefficients H ab and H ab S are related to the functions G ab and G ab S defined in Refs. [9] via (see, e.g., Eq. (21) of Ref. [12]) with where where the superscript "(c)" denotes the contributions with two current-current operators Q 1,2 , while "(cp)" refers to those with one operator Q 1,2 and one penguin operator Q 3−6 and "(p)" labels the terms involving two penguin operators. The functions H (cp) ab (z) and H (cp) ab S (z) are the main focus of this paper.

Calculation
The Wilson coefficients H ab (z) and H ab S (z) encode the short-distance physics and are independent of the external states in the matrix elements in Eqs. (8) and (14). Thus one may replace the mesons by free quarks, i.e. calculate the forward-scattering amplitude b +s →b + s in perturbation theory and apply the optical theorem in order to extract the desired absorptive part. By equating Eq. (8) with Eq. (14) one determines H ab (z) and H ab S (z). The infrared singularities present in both sides of this matching equation factorise, which makes the desired coefficients meaningful infrared-safe perturbative quantities. The external quarks are on-shell, i.e. we have p 2 b = m 2 b and may choose p s = 0 since we use m s = 0 and terms proportional to p b · p s match onto power-suppressed matrix elements. Thus we must evaluate two-point loop integrals with external momentum q 2 = m 2 b . In our calculation we regulate the infrared divergences with a gluon mass which introduces a further mass scale, m g . We introduce the gluon propagator as It is possible to expand the Feynman integrals for m g ≪ m b . We perform this expansion at the level of the master integrals as described below. We further employ an arbitrary QCD gauge parameter ξ and use its cancellation as a check of our calculation.
In the following we describe our methodology for the dominant contribution encoded in H cc (z) and H cc S (z). The calculational steps for the CKM-suppressed contributions involving H uc (z) and H uc S (z) are the same. Our practical calculation proceeds as follows: We consider the bilocal matrix elements where O i and O j are operators from Eqs. (6) and (7). At one-loop order we have to consider 4 }. The matrix elements with evanescent operators enter via the renormalisation procedure. One may formulate this procedure in terms of either bare and renormalised Wilson coefficients or bare and renormalised operators. With the former choice we have where Z QQ and Z QE are 6×6 and 6×4 matrices, respectively. They can be extracted from Ref. [19]. The entries in Eq. (20) represented by a * are irrelevant for our calculation. The UV poles contained in Z QQ and Z QE force us to include O(ǫ) terms in the one-loop transform the result to the traditional operator basis [20,21], and compare to the literature [9] in order to have a non-trivial cross check for the implementation of the CMM operator basis. New results are obtained for i d 4 x T Q 1−2 (x)Q 3−6 (0) .
For our calculation we use a well-tested program chain including qgraf [22] for the generation of the amplitudes, q2e and exp [23,24] for the identification of the integral families and FORM [25] for the algebraic manipulations and the traces of the γ matrices. As an alternative to q2e we also use the program tapir [26] which automatically generates FORM code, in which scalar products in the numerator are re-written in denominator factors and relations implementing a partial fraction decomposition are applied, if necessary. Furthermore, the input files for FIRE [27] are automatically generated. The Feynman rules involving the ∆B = 1 and ∆B = 2 operators have been obtained with the help of FeynRules [28] and FeynCalc [29,30].
At one-loop order only the type of diagrams shown in Fig. 1(a) contribute. At two loops we can distinguish four different classes of Feynman diagrams, see also Fig. 1. Figures 1(b) and (c) show the type of diagrams which contribute to the matrix element i d 4 x T Q 1,2 (x)Q 1,2 (0) . These topologies are also present if one of the operators is replaced by a penguin operator. Note that in Fig. 1(c) one of the closed quark loops contains charm or up quarks whereas the other may contain all five active flavours. In Fig. 1(d) and (e) we show sample diagrams which require the presence of a penguin operator. In Fig. 1(d) it is the left operator whereas in (e) it is the one on the external quark line.
We have implemented two approaches for the manipulation of the fermion spinor lines. In the first approach we concentrate on tensor integrals and various manipulations of Dirac structures. We use FeynCalc [29][30][31] together with Fermat [32] to obtain formulae for tensor reduction which we then implement in FORM. To this end the tensor reduction algorithm of FeynCalc was improved using ideas from [33]. In the second approach we construct projectors to all Dirac structures. This has the advantage that we can take traces and afterwards only scalar expressions have to be manipulated. More details can be found in Appendix A.
At this point a comment concerning the expansion in m c is in order. Since we restrict ourselves to quadratic terms in m c , i.e. linear terms in z, all loop integrals with both bottom and charm quark lines present in the same loop can be naively Taylor-expanded in m c before performing the loop integrations. 1 All such diagrams are contained in the class which is represented by Fig. 1(a) and (b). In all other cases we can apply the so-called large-momentum expansion [35] as implemented in exp [23,24]. However, our explicit calculation shows that up to order z indeed a naive expansion in m c is sufficient.
For the reduction to master integrals we use FIRE [27] and LiteRed [36,37]. For all infrared contributions the reduction is performed for general gluon mass m g . Afterwards we consider the limit of small m g and perform an asymptotic expansion [35] for m g ≪ m b at the level of the master integrals. We have performed numerical cross-checks of the expansions with the help of FIESTA [38]. After the asymptotic expansion we have to compute single-scale one-and two-loop integrals, most of which are available in the literature (see, e.g., Ref. [35]). The remaining ones are straightforward to compute.
We multiply the matrix element on both sides of the matching equation with Z 2 ψ , where Z ψ is the quark field renormalisation constant in the MS scheme. This renders both expressions UV-finite. Note, that they still depend on the gauge parameter which is due to the gluon mass used as infrared regulator. For the renormalisation of the charm quark mass we use both the MS and on-shell scheme, see also Section 4. No renormalisation of the bottom or strange quark mass is needed since in the considered order there are no corresponding self-energy diagrams.
For the ∆B = 2 theory we calculate one-loop QCD corrections for the matrix elements of the minimal operator basis in Eq. (10). Sample Feynman diagrams, which have to be considered at NLO, are shown in Fig. 2. The results of the matrix elements in both the ∆B = 1 and ∆B = 2 theories can be expressed as a linear combination of the tree-level matrix elements of Q, Q S , R 0 and the unphysical operators in Eq. (13). Since both results are UV-finite we can take the limit ǫ → 0 and then read off the desired NLO corrections to the ∆B = 2 Wilson coefficients H ab and H ab S . We observe that the infrared regulator m g and the gauge parameter cancel from these coefficients, providing a non-trivial check of the calculation. H ab and H ab S depend on the renormalisation scales µ 1 and µ 2 , at which the renormalised operators are defined in the ∆B = 1 and ∆B = 2 theories, respectively. The µ 1 -dependence of H ab and H ab S diminishes order-by-order in perturbation theory and is commonly used as a means to estimate the accuracy of the truncated perturbative series. The µ 2 -dependence cancels in the matching procedure of the perturbative ∆B = 2 matrix elements with their non-perturbative counterparts.

Analytical and numerical results
In the following we discuss the results for the matching coefficients H ab and H ab S (z) introduced in Eq. (14).
We start with the analytic expressions for the penguin contributions H (cp) ab (z) and H (cp) ab S (z) (for ab = uu, uc and cc) introduced in Eq. (17). It is convenient to decompose the ∆B = 2 matching coefficients in terms of the ∆B = 1 coefficients C i of the |∆B| = 1 Hamiltonian in Eq. (5): We furthermore introduce the perturbative expansion as with the MS masses m q and the pole (on-shell) masses m OS q . While it is easier to employ on-shell masses in the calculation, their poor definition (especially of m OS c ) make them unsuited for numerical evaluations and we will always use MS values as inputs.
The one-loop coefficients p ab,(0) ij , p S,ab,(0) ij can be extracted from Ref. [13], where the full m c dependence has been taken into account, by transforming the result to the operator bases used in this paper. We can reproduce these results in an expansion in z including the linear terms. Note that p as well as p S,cc,(0) 13 The two-loop coefficients p and p S,cc,(1) 13 with L 1 = log(µ 2 1 /m 2 b ) and L 2 = log(µ 2 2 /m 2 b ). Furthermore, we introduce the symbols N L , N V and N H which label closed fermion loops with mass 0, m c and m b , respectively. In the numerical evaluation we set N L = 3, N V = 1 and N H = 1.
z and transforming to the operator basis used in [15]. We note that the NLO coefficients with i = 1, 2 and j = 8 are only one-loop quantities and can be extracted from Ref. [9].
It is interesting to note that the √ 3 in our results originate from the Feynman diagrams in Fig. 1(b) where in one of the closed loops a massive bottom quark is present. We mention that our results passes the checks mentioned above, the gauge parameter and the gluon mass vanish. As an additional check we have re-done the calculation employing dimensional regularisation of the IR divergences, which requires to do the LO matching at order ǫ, and found the same results.
The results in Eqs. (26) to (30) contain terms of order z log z which result from diagrams with charm self-energies and mass counterterms. The large coefficients of these terms, proportional to the LO term γ In the following we discuss the "cc" contribution of the quantities H ab and H ab S in the MS scheme. We refrain from showing explicit results for the "uu" and "uc" contributions which show a similar pattern. We have where "c-gb" refers to the diagrams with two current-current operators and a gluon bridge, see Fig. 1(c). The numerical values are specific to the operator renormalisation scheme chosen by us. The scheme dependence cancels in combination with the NLO Wilson coefficients C 3−6 entering the numbers label with "cp". From Eq. (34) we observe that at one-loop order the penguin contribution is about a factor 20 smaller than the terms proportional to C 1 and C 2 , which justifies to calculate penguin contributions to lower orders in α s than those with two copies of C 1,2 . However, at two loops the impact of the penguin coefficients is larger. In the case of H cc the relative factor is less than three and in the case of H cc S the penguin coefficient is even bigger than the current-current contribution. We want to remark that the numerically most important penguin contribution is the one proportional to C 4 .
We want to remark that in all cases the fermionic contributions to the the penguin coefficients, which are known from Ref. [15], are significantly smaller than the non-fermionic terms computed in this paper. Still, using N H = N V = 1 and N L = 3 we observe a screening of the non-fermionic coefficient of close to 50%.
We observe that the expansion in z is well-behaved. For example, more than 90% of the non-fermionic penguin coefficients at two-loop order in Eq. (34) are provided by the m c → 0 approximation.
We are now in the position to evaluate the shift of the new corrections to the width difference. To illustrate the numerical effect of the new corrections we omit both the fermionic NNLO contributions computed in [14] and power corrections of order Λ QCD /m b . We furthermore concentrate on ∆Γ s . In addition to the quark masses and α s given above we have the following input parameters [43][44][45] Correlator Perturbative order z-dependence O 1,2 × O 1,2 [9] 1 loop exact 1 loop exact Thus the calculated corrections increase ∆Γ s by 0.003 ps −1 , which is almost as large as today's experimental error in Eq. (3). The size of the correction is also in the ballpark of the hadronic uncertainty, if ∆Γ s is predicted from ∆Γ s /∆M s , since hadronic uncertainties largely cancel from this ratio [12,15].
Next, we discuss the relative shift of ∆Γ s due to the contribution of the penguin contribution in more detail. At one-loop order we obtain where the quantity in denominator includes all current-current and current-penguin corrections up to order α 1 s . The penguin-penguin contributions are included up to order α 0 s (one-loop order). The numerator in Eq. (38) only contains the LO current-penguin contributions (indicated by the superscript "12 × 36").
where the numerator contains the new corrections computed in this paper together with the corresponding fermion contributions from [15]. Note that the non-N f penguin contribution overcompensates the N f terms [15]. In the pole scheme this leads to tiny corrections below the percent level. In the MS scheme the non-N f contribution is about a factor three bigger than the N f terms which leads to a relative correction of −1.4%.

Conclusions
In this paper, for the first time, the ∆B = 1 operator basis from Ref. [16] has been used for the computation of NLO corrections to the decay matrix element Γ q 12 , governing the width difference between the eigenstates of the B q −B q mass matrix and the CP asymmetry in semileptonic B q decays. After reproducing known results [9-11, 14, 15] we have obtained novel two-loop contribution to Γ q 12 , namely all contributions involving one current-current operator and one four-quark penguin operator. We have computed these two-loop corrections in an expansion in m c /m b including quadratic terms. Computerreadable expressions of our results can be downloaded from [39].
The calculated NLO effects dominate over the previously known partial results which contain only fermion loop contributions. While the NLO penguin contributions are numerically less relevant than those with two large current-current coefficients C 1,2 , they are needed for the theory prediction to match the experimental precision of ∆Γ s in Eq. (3). To fully keep up with experiment one further needs the contributions involving Q 8 at the two-loop level and a full NNLO (three-loop) calculation of the contributions with two current-current operators. For the NNLO calculation it is instrumental to use the CMM operator basis [16] as we did in this paper.

A Projector methodology
In this appendix we briefly describe the approach based on the construction of projectors for the various tensor structures. In general, the scattering amplitude of the process b +s →b + s can be parametrized as where c n describe the colour and i n the spinor indices. Note that the number of colour structures Σ (m) is finite. On the other hand, the basis of the Lorentz structure is a priori not finite. For a massless s quark the Lorentz structure can be expressed as 2 Γ (n) where P R = (1 + γ 5 )/2 and the basis vectors B (n) are given by In four space-time dimensions it is possible to avoid chains with more than four Dirac matrices, which is not possible in d = 4 − 2ǫ dimensions. However, in a fixed order in the perturbative expansion only a finite number of basis vectors B (n) appear.
In general, the coefficients A (n,m) include dimensionally regularized scalar Feynman integrals. To extract A (n,m) in Eq. (40), one can apply tensor reduction to get scalar integrand expressions. Alternatively, one can make use of the composition of Eq. (40). Hence, we define projection operators for Lorentz (P (n) ) and colour space (C (m) ), acting as A (n,m) = C (m) Tr d P (n) M , where C (m) commutes with the operations applied in Lorentz space. P (n) is constructed from a linear combination of the structures introduced in Eq. (42). It is understood that the traces are evaluated in d dimensions. Note that in our case no traces including γ 5 appear since Eq.
Note that on the right-hand side one has a product of two traces.
A caveat of this approach is that the complexity of the matrix p grows considerably with the number of γ matrices in the basis elements B (n) . For our NLO calculation, we have to consider terms up to n = 9 which leads to products of two d-dimensional traces where each one contains up to 18 γ-matrices. This non-trivial computational task was done using FORM [25], where the special hints described in the manual have been used. To avoid unnecessary recomputations, we evaluate each occurring trace product separately and include the result in a lookup table.