The width difference in B−B¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ B-\overline{B} $$\end{document} beyond mixing at order αs and

We complete the calculation of the element Γ12q\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\Gamma}_{12}^q $$\end{document} of the decay matrix in Bq−B¯q\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {B}_q-{\overline{B}}_q $$\end{document} mixing, q = d, s, to order αs in the leading power of the Heavy Quark Expansion. To this end we compute one- and two-loop contributions involving two four-quark penguin operators. Furthermore, we present two-loop QCD corrections involving a chromomagnetic operator and either a current-current or four-quark penguin operator. Such contributions are of order αs2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\alpha}_s^2 $$\end{document}, i.e. next-to-next-to-leading-order. We also present one-loop and two-loop results involving two chromomagnetic operators which are formally of next-to-next-to-leading and next-to-next-to-next-to-leading-order, respectively. With our new corrections we obtain the Standard-Model prediction ∆Γs/∆Ms = (5.20 ± 0.69) · 10−3 if Γ12s\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\Gamma}_{12}^s $$\end{document} is expressed in terms of the MS¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{\mathrm{MS}} $$\end{document} b-quark mass, while we find ∆Γs/∆Ms = (4.70 ± 0.96) · 10−3 instead for the use of the pole mass.


Introduction
The description of B q −B q mixing, where q = d or s, involves two hermitian 2 × 2 matrices, the mass matrix M and the decay matrix Γ. Their off-diagonal elements M q 12 and Γ q 12 enter the observables related to B q −B q mixing, namely and a q fs = Im (1.2) Here M L,H and Γ L,H denote the masses and widths of the two eigenstates found by diagonalizing M − iΓ/2. The mass difference ∆M q and the width difference ∆Γ q are related to M q 12 and Γ q 12 as ∆M q 2|M q 12 |, (1. 3) In the Standard Model (SM) the phase between −Γ q 12 and M q 12 is small, so that the CP asymmetry in flavor-specific decays, a q fs , is much smaller than ∆Γ q /∆M q and further ∆Γ q 2|Γ q 12 |. All results calculated in this paper equally apply to the B s and B d systems. For definiteness, we quote all formulae for the case of B s −B s mixing, the generalization to B d −B d mixing is found by replacing the elements V qs , q = u, c, t, of the Cabibbo-Kobayashi-Maskawa (CKM) matrix by V qd .

JHEP04(2022)006
Currently, better theory predictions are needed for the case of B s −B s mixing to be competitive with the precise experimental values ∆M exp s = (17.7656 ± 0.0057) ps −1 [1] , where the quoted number for ∆Γ exp s is derived from data of LHCb [3], CMS [4], ATLAS [5], CDF [6], and DØ [7]. M s 12 probes virtual contributions of very heavy particles, while Γ s 12 is mainly sensitive to new physics mediated by particles with masses below the electroweak scale. Nevertheless, a better theory prediction of |Γ s 12 | also helps to quantify new physics in M s 12 : both ∆M s and ∆Γ s are proportional to |V ts | 2 , where V ts is an element of the CKM matrix. |V ts | is calculated from (and is essentially identical to) |V cb | extracted from measured b → c ν, = e, µ, branching ratios. The unfortunate discrepancy between the values for |V cb | found from inclusive and exclusive decays inflicts an uncertainty of order 15% on the predicted ∆M s . Now V ts cancels from ∆Γ s /∆M s in eq. (1. 3), so that the SM prediction of this ratio is not affected by the V cb controversy.
In this paper we address Γ s 12 at leading order of the heavy-quark expansion ("leading power"), which expresses Γ s 12 as a series in powers of Λ QCD /m b . At this order one encounters only two physical ∆B = 2 operators, whose hadronic matrix elements have been calculated with high precision with lattice QCD [8]. These matrix elements are multiplied with Wilson coefficients which are calculated in perturbative QCD. The insufficient accuracy of the Wilson coefficients dominates the uncertainty of the SM prediction of ∆Γ q [9][10][11][12][13][14][15], which exceeds the experimental error in eq. (1.4). The perturbative calculation of power corrections to Γ s 12 has been carried out to order α 0 s [16] and first lattice results for the associated hadronic matrix elements are also available [17] (for results from sum rules, see ref. [18] and references cited therein).
The |∆B| = 1 Hamiltonian H |∆B|=1 eff comprises current-current operators with large coefficients C 1,2 and the four-quark penguin operators whose coefficients C 3−6 are small, with magnitudes well below 0.1, at the scale µ 1 = O(m b ) at which they enter Γ s 12 . At order α 0 s , Γ s 12 is composed of one-loop contributions proportional to C j C k with j, k ≤ 6. H |∆B|=1 eff further involves the chromomagnetic penguin operator with coefficient C 8 ∼ −0.16, whose leading contribution is of order α s and enters Γ s 12 as products C 8 C k with k ≤ 6. In refs. [9][10][11][12][13][14] the small coefficients C 3−6 have been formally treated as O(α s ). With this counting the one-loop terms with C 1,2 C 3−6 contribute to next-to-leading order (NLO) and those involving two factors of C 3−6 are already part of the next-to-next-to-leading order (NNLO). First steps towards NNLO accuracy have been done in refs. [13,14] by calculating contribution proportional to the number of active quark flavors, i.e. loop diagrams with a closed fermion line.
As in ref. [15] we use the conventional notion of "NLO" and "NNLO" in this paper and treat C 3−6 on the same footing as C 1,2 . With this counting the NLO prediction of Γ s 12 requires the calculation of the yet unknown two-loop contributions with one or two four-quark penguin operators. In this paper present several two-loop calculations, namely:

JHEP04(2022)006
• penguin contributions proportional to the product of two C 3−6 coefficients. This contribution completes the prediction of Γ s 12 to order α s , which is NLO in the abovementioned conventional power counting. The corresponding one-loop corrections have been computed in ref. [16]. Two-loop contributions with one current-current and one four-quark penguin operator have been calculated in ref. [15].
• the contribution proportional to the product of C 8 and one of C 1−6 . The calculated one-loop and two-loop terms contribute to NLO and NNLO, respectively. The piece of the one-loop correction proportional to the number N f of active quark flavors (stemming from diagrams with closed quark loops) has been computed in ref. [14].
• the contribution proportional to C 2 8 . Here the one-loop contribution is already of NNLO and not yet available in the literature, except for the α 2 s N f part [14]. We further provide results for the two-loop term which is N 3 LO.
As in ref. [15] we use the CMM basis [19] for the |∆B| = 1 operators and calculate the two-loop QCD corrections as an expansion in up to linear order. The paper is organized as follows: in the next section we briefly discuss the operator bases of the |∆B| = 1 and |∆B| = 2 theories. Afterwards, in section 3 we provide some details of our calculation and in particular describe the matching procedure for the case of dimensionally regularized infra-red singularities. Analytic result for all new matching coefficients are listed in section 4 and we present our numerical result for ∆Γ s in section 5. Section 6 contains our conclusions. In the appendix we provide results for the renormalization constants relevant for the operator mixing in the |∆B| = 2 theory.

Operator bases
The framework of our calculation is identical to the one used in ref. [15] and thus in the following we repeat only the essential formulae needed to compute the width difference. The new contributions considered in this paper require an extension of the |∆B| = 2 operator basis which is discussed in more detail.
For the effective |∆B| = 1 theory we use the weak Hamiltonian where explicit expressions for the (physical and evanescent) operators can be found in ref. [19]. Q 1 , Q penguin operators. Q 8 is the chromomagnetic penguin operator. In eq. (2.1) we have introduced the quantities λ s a = V * as V ab , a = u, c, t, which contain the CKM matrix elements. Furthermore, we have used λ s t = −λ s c − λ s u and G F is the Fermi constant. Our two-loop calculations involve one-loop diagrams with counterterms to the physical operators in eq. (2.1) and these counterterms comprise both physical and evanescent operators.
As mentioned in the Introduction, we specify our discussion to b → s decays relevant for B s −B s mixing. The corresponding expressions for B d −B d mixing are obtained by replacing V as with V ad . Using the optical theorem we can relate Γ s 12 to theB s → B s forward scattering amplitude: where "Abs" stands for the absorptive part and T is the time ordering operator. Note that Γ s 12 encodes the information of the inclusive decay rate into final states common to B s and B s . Following ref. [9] we decompose Γ s 12 as Let us now discuss the effective |∆B| = 2 theory. To leading power in 1/m b it is convenient to introduce the following four operators where i, j are color indices. In four space-time dimensions there are only two independent operators which we choose as Q and Q S since we have (for D = 4) Q = Q and where R 0 describes 1/m b -suppressed contributions to Γ s 12 [16]. α 1,2 are QCD correction factors which ensure that the MS renormalized matrix element R 0 has the desired power suppression [9,12].
Using the Heavy Quark Expansion (HQE) it is thus possible to write Γ ab 12 in eq. (2.3) as ) 2 and the ellipses denoting higher-order terms in Λ QCD /m b . Here z is defined in terms of pole quark masses. Later we will trade z for the ratio of MS masses which leads to a better behavior of the perturbative series. H ab and H ab S are ultra-violet and infra-red finite matching coefficients which we decompose as follows where for our calculation the values of e j,l , which parametrize the O( ) terms in the definition of the evanescent operators are irrelevant since the operators in eq. (2.9) do not appear in one-loop counterterm contributions. These numbers, however, become important at NNLO to fully specify the renormalization scheme at this order.

Calculation and matching
The setup which we use for our calculation has already been described in ref. [15]. For convenience of the reader we repeat the essential steps and stress the differences in the following. Figure 1 shows typical one-and two-loop Feynman diagrams for the new contributions considered in this paper. The displayed diagrams correspond to the |∆B| = 1 side of the matching equation. In addition, one needs the one-loop diagrams with a gluon dressing the ∆B = 2 operators Q andQ S to determine the desired Wilson coefficients H ab and H ab S in eqs. (2.6) and (2.7). We perform the calculation for a generic QCD gauge parameter which drops out in the final result for each matching coefficient and thereby provides a non-trivial check of our calculation. The counterterms to the ∆B = 1 operators and the gauge coupling g s in the Feynman diagrams exemplified in figure 1 are all evaluated at the renormalization scale µ 1 . Conversely, operators and couplings on the |∆B| = 2 side are evaluated at the scale µ 2 . The unphysical µ 1 dependence of H ab (z) and H ab S (z) diminishes order-by-order in perturbation theory and can be used to assess the accuracy of the calculated result. The µ 2 dependence of H ab (z) and H ab S (z), however, cancels with the µ 2 dependence of the JHEP04(2022)006 hadronic matrix element, which enters the lattice-continuum matching. For calculational convenience we first choose µ 1 = µ 2 and implement the separation µ 1 = µ 2 with the help of renormalization group techniques.
We pursue two different approaches to treat the four-quark amplitudes. The first one is based on tensor integrals combined with various manipulations of the Dirac structures and relies on FeynCalc [22][23][24] and Fermat [25]. The so-obtained formulae are then exported to FORM [26].
For the contribution Q 3−6 × Q 3−6 routines are needed which can handle tensor integrals up to rank 6. The second approach is based on projectors (see appendix of ref. [15]) which allows taking traces. Thus, one only has to deal with scalar expressions. However, one needs to calculate products of two traces with up to 18 γ matrices 1 in each trace. We find that both approaches lead to the same expressions for the amplitude with two ∆B = 1 operators once the latter is expressed in terms of tree-level ∆B = 2 matrix elements.
For the reduction of the ∆B = 1 amplitude we use FIRE [27] with symmetries from LiteRed [28,29] and obtain four two-loop master integrals. Their evaluation as an expansion in is straightforward. The amplitudes in the |∆B| = 1 and |∆B| = 2 theories contain both ultra-violet and infra-red singularities. The former are cured with parameter, quark field, and operator renormalization. We use the one-loop counterterms for α s in the MS scheme and renormalize the charm quark in the one-loop expression in the on-shell (or pole) scheme. The renormalization of the bottom quark, which we also renormalize on-shell, is only needed for the contributions involving Q 8 . We also perform the renormalization of the external quark fields in the MS scheme. The counterterms needed for the renormalization of the ∆B = 1 operator mixing can be taken from the literature [30]. The renormalization constants of the ∆B = 2 part are given in appendix A.

JHEP04(2022)006
In order to regulate the infra-red singularities two possibilities come to mind: one can either introduce a (small) gluon mass, m g , or instead use dimensional regularization.
The choice m g = 0 is conceptually simpler and has the advantage that after renormalization the ∆B = 1 and ∆B = 2 amplitudes are separately finite and one can take the limit D → 4, which eliminates all evanescent operators before matching the two theories. Furthermore, it is possible to use four-dimensional relations in order to arrive at a minimal operator basis. However, a finite gluon mass breaks gauge invariance and thus, in general, additional counterterms have to be introduced for its restoration. In our application the two-loop ∆B = 1 amplitudes with four-quark operators do not involve three-gluon vertices, and thus it is safe to regulate the infra-red divergences with m g = 0. However, at three-loop level this is not the case. Furthermore, the two-loop corrections with two Q 8 operators also contain infra-red divergences in the non-abelian part.
Regulating the infra-red divergences dimensionally using the same regulator as for the ultra-violet divergences has the advantage that the loop integrals are simpler. However, the matching has to be performed with divergent quantities in D = 4 dimensions. As a consequence lower-order corrections have to be computed to higher order in , meaning that also the evanescent operators have to be taken into account.
In our calculation we proceeded as follows: we have computed the contributions Q 1−6 × Q 1−6 and Q 1−2 × Q 8 both for m g = 0 and m g = 0 and have obtained identical results for the matching coefficients, which provides sufficient confidence that the conceptually more involved approach where the infra-red divergences are regularized dimensionally is understood. Thus, the calculation of the Q 3−6 × Q 8 and Q 8 × Q 8 have only been performed for m g = 0.
In the following we provide some details to the matching procedure. In this context we also refer to ref. [31] where the contribution Q 1,2 × Q 1,2 is discussed. We introduce the |∆B| = 1 and |∆B| = 2 amplitude in a schematic way as where the H X , A X and B XY have an expansion both in α s and with B QQ = B EE = 1 and B EQ = B QE = 0 at LO. · 0 denote tree-level matrix elements. Starting from two-loop order 2 A X and B XY contain infrared 1/ poles. The presence of these poles force us to calculate the LO coefficients H X to order in order to obtain the correct finite 0 piece on the right-hand side of eq. (3.1). Thus, the desired finite matching coefficients H X have the following expansion in α s and : the difference A ∆B=1 − A ∆B=2 is finite, which constitutes an important consistency check. In a next step we concentrate on the part of A ∆B=1 − A ∆B=2 proportional to Q 0 , which contains H For definiteness, we now consider the LO expression of the Q 3−6 × Q 3−6 contribution, where for simplicity we set the matching coefficients C 4 , C 5 and C 6 to zero and display only the terms proportional to C 2 3 . Then the LO ∆B = 1 amplitude including terms of O( ) is given by where we set the number of colors to N c = 3. At the same order the ∆B = 2 amplitude reads and from the matching procedure we obtain we observe an explicit cancellation of all 1/ poles multiplying C 2 3 which allows us to take the limit D → 4. We also find the coefficients independent of the gauge parameter.
The presence of R 0 in eqs. (3.3) to (3.5) requires some explanation: is ensured by a finite renormalization). To derive this result one employs four-dimensional Dirac algebra (such as using the Fierz identity from [16]) and for D = 4 the definition of R 0 in eq. (2.5) thus includes an evanescent piece. One may write Clearly, if one uses a gluon mass as infra-red regulator, this subtlety does not occur, because the matching is done in D = 4 dimensions. In our case of dimensional infra-red regularization, however, E R 0 must be included in the LO matching just as any other evanescent operator. If we were interested in the C 2 3 contributions to the 1/m b -suppressed part (which is beyond the scope of this paper), we would have to provide different coefficients for the physical operator R phys 0 and the unphysical E R 0 . For our choice of external states, namely zero momenta p s for the light strange quarks, we cannot determine the coefficient of R phys 0 , because R phys The O( ) terms of the coefficients of evanescent operators, i.e. H (0,1) are not needed for the NLO calculation presented in this paper. However, they will be relevant at NNLO and beyond.

Analytic results
In this section we present analytic results for the new contributions to H (p) ab and H (p) ab S introduced in eq. (2.7). For this purpose it is convenient to decompose these quantities according to the |∆B| = 1 matching coefficients as follows and to write the perturbative expansion of the coefficients to two-loop contributions. We define the strong coupling constant with five active quark flavors at the renormalization scale µ 1 , i.e. we have α s (µ 1 ) ≡ α (5) s (µ 1 ). Both the charm and bottom quark masses are defined in the on-shell scheme. Furthermore, we fix the number of colors to N c = 3. Computer-readable expressions for all results for generic N c are available as supplementary material to this paper and can be downloaded from [32].

Four-quark penguin operators
We start with the Q 3−6 × Q 3−6 contribution. Both at one-and two-loop order, which contribute to LO and NLO, respectively, the "cc", "uc" and "uu" contributions agree, because penguin operators come with the CKM factor −λ s At one-loop order exact results are available [16], which we repeat for convenience The symbols N L and N V label closed fermion loops with mass 0 and m c , respectively. In the numerical evaluation we set N L = 3 and N V = 1.
The two-loop results are new. Their expansions up to linear order in z are given by (4.24) Here N H = 1 labels closed fermion loops with mass m b and As a novel feature compared to the NLO calculation with two current-current operators [9], the penguin operator contributions involve Feynman diagrams with an FCNC b → s self-energy in an external leg, cf. figure 1. Owing to p 2 b = m 2 b = p 2 s = 0 these diagrams contribute to the result in the same way as all other diagrams [33]. Indeed, we find that their omission would lead to a divergent result.

Chromomagnetic and four-quark operators
In this subsection we present results for all contributions involving one chromomagnetic and one of the four-quark operators Q 1 , . . . , Q 6 . Here the one-and two-loop corrections correspond to NLO and NNLO contributions.
We start with Q 1,2 × Q 8 where the (exact) one-loop result is given by [9] p cc,(0) 18 At two-loop order the results are new. The "cc" contribution is given by (4.28) Note that the (uu) contribution is not simply obtained by taking the limit z → 0 in the expressions of eq. (4.28) since there are charm quark loops not connected to the external operators. We thus have For the contribution Q 3−6 × Q 8 we observe that both at one-and two-loop order we obtain the same results for the "cc", "uu" and "uc" contributions and thus we have At two-loop order our results read

Two chromomagnetic operators
Finally, we come to the Q 8 × Q 8 contribution, where the one-loop corrections are already of NNLO. The one-loop result, for which only the N f -piece has been known in the literature, is given by (4.44)

Numerical results
In this section we present the numerical effect of the new corrections to ∆Γ s and a s fs . We start with discussing the relative size of the contributions from the various operators and consider afterwards the ratio ∆Γ s /∆M s , from which |V ts | and the ballpark of the hadronic uncertainties cancel. Finally, we use the measured result for ∆M s and present updated results for ∆Γ s in two different renormalization schemes. We also present updated results for a s fs . The calculations described in the previous sections and the analytic results presented in section 4 use the MS scheme for the strong coupling constant and the operator mixing and the on-shell scheme for the charm and bottom quark masses. It is well known that the latter choice leads to large perturbative corrections. Thus, we choose as our default renormalization scheme the one where all parameters are defined in the MS scheme. It is obtained with the help of the one-loop relations between the on-shell and MS charm and bottom quark masses. We define a second renormalization scheme where the overall factor m 2 b (see, e.g., eq. (2.6)) is defined in the on-shell scheme, but H ab andH ab S depend on the quark masses in the MS scheme. In the following we refer to this scheme as the "pole" scheme [13,14]. Note that after each scheme change, which adds z-exact expressions to the two-loop term, we re-expand the latter in z up to linear order to be consistent with our genuine two-loop calculation.
For convenience, we summarize in table 2 the input parameters needed for our numerical analysis. In addition we have (see ref. [14]) B Bs = 0.813 ± 0.034 [8] B S,Bs = 1.31 ± 0.09 [8] f Bs = 0.2307 ± 0.0013 GeV [36]   For the matrix elements of the 1/m b suppressed corrections we have The results for B s |R 2 |B s , B s |R 2 |B s , B s |R 3 |B s , and B s |R 3 |B s can be found in ref. [17] and we extract the remaining three matrix elements from [8]. For B s |R 1 |B s and B s |R 1 |B s the ratio of the bottom and strange quark masses is needed m b (µ)/m s (µ) = 52.55 ± 0.55 [37]. Let us next discuss our choices for the various renormalization schemes. We fix the high scale in the ∆B = 1 theory to µ 0 = 165 GeV ≈ 2m W ≈ m t (m t ). Since µ 2 is closely connected to the lattice results for B Bq ,B S,Bs and the 1/m b matrix elements of eq.

JHEP04(2022)006
Contribution X r X (MS) r X (pole)  Table 3. Relative contributions in percent in the MS and pole schemes. The breakdown into oneand two-loop contributions is shown inside the round brackets. In the last column we mention the corresponding perturbative order.
In table 3 we show the relative size of the individual contributions to ∆Γ s both in the MS and pole scheme. They are defined as Power-suppressed 1/m b corrections are only included in the denominator of eq. (5.4) but not in the numerator. In both renormalization schemes the dominant contribution is given by the Q 1,2 × Q 1,2 , followed about a 7% contribution from Q 1,2 × Q 3−6 . The remaining terms contribute at the 1% level or below. Note that these contributions are necessary to obtain complete NLO and NNLO corrections. It is interesting to note that the QCD corrections to Q 1,2 × Q 1,2 amount only to 9% in the MS scheme but to more than 30% in the pole scheme. Also for the contribution Q 1,2 × Q 3−6 the QCD corrections are about a factor of three larger in the pole scheme whereas for Q 3−6 × Q 3−6 the situation is vice versa. For the contributions involving Q 8 the QCD corrections in the MS scheme amount to up to about 50% of the leading order term, though their absolute contribution is small. Let us next consider ∆Γ s /∆M s . We use eq. (1.3) with Γ s 12 from eq. (2.2) and M s 12 from ref. [39] where two-loop QCD corrections have been computed. In the two renormalization schemes our results read where the subscripts indicate the source of the uncertainties: "scale" denotes the uncertainties from the variation of µ 1 , "BB S " those from the leading order bag parameters and "input" The largest uncertainty is induced by the power-suppressed 1/m b corrections. It is obtained by combining the uncertainties from the seven matrix elements of eq. (5.3) in quadrature taking into account the 100% correlation of B s |R 2 |B s and B s |R 2 |B s . Next, there is the renormalization scale uncertainty, which we use to estimate the contribution from unknown higher order corrections. We obtain the numbers in eq. (5.5) by varying µ 1 between 2.5 GeV and 10.0 GeV while keeping µ 2 , µ c and µ b at their default values. A simultaneous variation of µ 1 = µ b = µ c leads to significant larger scale uncertainties, which is expected, because the anomalous dimension of the quark mass is large and appears in the coefficient of log(µ b /m b ).
The last three uncertainties in eq. (5.5) are correlated between the two schemes. The scale dependence is plotted in figure 2(a) and leads to the asymmetric uncertainties quoted in eq. (5.5). The difference between the central values found in the pole and MS schemes is around 11%, i.e. of the expected size of an NNLO correction.
We proceed in a similar way for a s fs . We use eq. (1.2) and obtain a s fs = 2.07 +0.10 −0.11 scale ± 0.01 BB S ± 0.06 1/m b ± 0.06 input × 10 −5 (pole) , In figure 2(b) we show the dependence on µ 1 for the two renormalization schemes. Here the interval in the pole scheme is completely contained in the one from the MS scheme. The predictions in eqs. (5.5) and (5.6) are consistent with those of ref. [14], but the central values for ∆Γ s /∆M s in eq. (5.5) are larger in both schemes. In ref. [14] only partial NLO corrections to the Q 1,2 × Q 3−6 contribution and no Q 3−6 × Q 3−6 or NNLO Q 8 terms have been included. Inspecting the sources of the differences in detail, we find that almost 2/3 of these stem from the new contributions presented in ref. [15] and this paper. The remainder is due to terms, which are formally of higher order in α s . Interestingly, the µ 1 dependence of ∆Γ s /∆M s is much smaller in eq. (5.5) compared to ref. [14], while the JHEP04(2022)006 situation is vice versa for a s fs . We trace this feature back to the use of α s (µ 1 ) versus α s (µ 2 ) in certain NLO terms, both of which are allowed choices in the considered order. In view of this observation and the fact that the intervals from the scale uncertainty of ∆Γ s /∆M s in both schemes barely overlap, we conclude that the µ 1 dependence is not always a good estimate of the size of the unknown higher-order corrections.
In a next step we can use the experimental result for ∆M s [34], and obtain for ∆Γ s in the two renormalization schemes Comparing our prediction with the experimental value in eq. (1.4) we see that both the "pole" and MS results are consistent with the measured value, but the central value of the former is closer to the experimental result. One needs a better perturbative precision (which will bring the "pole" and MS results closer to each other and reduce the scale uncertainty) and more precise lattice results for the matrix elements of the 1/m b -suppressed operators to quantify new-physics contributions to ∆Γ s /∆M s .
The value for ∆Γ s /∆M s quoted in eq. (5.5) also applies to ∆Γ d /∆M d for two reasons: first, while the CKM-suppressed contribution to ∆Γ d /∆M d is a priori expected to be relevant due to |λ d u /λ d t | |λ s u /λ s t |, it merely contributes at the percent level because of a numerical cancellation in the sum of uc and uu contributions [11]. Second, the non-perturbative calculations of the B s and B d hadronic matrix elements agree well within their error bars. As a result the central values for ∆Γ d /∆M d and ∆Γ s /∆M s agree within a few percent (see e.g. [14]) and the difference is much smaller than the uncertainty in eq. (5.5

Conclusions
In this paper we have completed the calculation of the NLO contributions to the decay matrix element Γ q 12 appearing in B q −B q mixing. These new contributions involve two-loop diagrams with two four-quark penguin operators. We have further calculated two-loop contributions with one or two copies of the chromomagnetic penguin operators, which belong to NNLO or N 3 LO, respectively. All results are obtained as an expansion to first JHEP04(2022)006 order in z = m 2 c /m 2 b , except for the one-loop Q 8 × Q 8 contribution for which our result has the exact z-dependence. With our new results the theoretical uncertainties associated with the penguin sector are under full control and way below the experimental error of the width difference ∆Γ s 2|Γ s 12 | in eq. (1.4). We present updated predictions for ∆Γ s and ∆Γ d and the CP asymmetry in flavor-specific B s decays, a s fs . For the width differences we find the predictions in the pole and MS schemes to differ by 11%, which invigorates the need for a full NNLO calculation of the contributions from current-current operators.
We provide the newly obtained matching coefficients in a computer readable format with full dependence on the number of colors N c . In the same way we present the renormalization matrix Z ij of the ∆B = 2 operators including the submatrices governing the mixing of evanescent operators with physical operators and among each other.

JHEP04(2022)006
We can write the matrix of renormalization constants as a 20 × 20 matrix which is naturally decomposed into four sub-matrices where Z QQ , Z QE , Z EQ and Z EE have the dimension 3 × 3, 3 × 17, 17 × 3 and 17 × 17, respectively. We define Z ∆B=2 via the renormalization of the coefficient functions as follows where C bare and C ren are 20-dimensional vectors of the bare and renormalized |∆B| = 2 coefficient functions, respectively. The perturbative expansion of the sub-matrices is introduced as where the first superscript denotes the order in α s and the second one the order in 1/ . Note that at one-loop order the matrix Z EQ only contains finite contributions. In order to determine the matrix elements of Z ∆B=2 we compute the amplitude b +s →b + s in the kinematics described above, take into account the field renormalization of the external quarks in the MS scheme and require that the remaining poles in , which are all of ultra-violet nature, are absorbed by the operator mixing via Z ∆B=2 . This condition fixes all matrix elements but the ones in Z EQ . The latter are fixed by the requirement that the contributions of evanescent operators vanish in D = 4 dimensions [21,40]. Note that to our order we do not have to renormalize the common strange and bottom quark mass.
An important check of our calculation is the locality of the extracted renormalization constants. Furthermore, we perform the calculation for general QCD gauge parameter and observe that the matrix Z ∆B=2 is independent of ξ.
In the following we present explicit results for the one-loop corrections to Z QQ , Z QE and Z EE . For N c = 3 we have where the entries " * " are not needed for our calculation.