The width difference in the $B-\bar{B}$ system at next-to-next-to-leading order of QCD

We extend the theoretical prediction for the width difference $\Delta \Gamma_q$ in the mixing of neutral $B$ mesons in the Standard Model to next-to-next-to-leading order in $\alpha_s$. To this aim we calculate three-loop diagrams with two $|\Delta B|=1$ current-current operators analytically. In the matching between $|\Delta B|=1$ and $|\Delta B|=2$ effective theories we regularize the infrared divergences dimensionally and take into account all relevant evanescent operators. Further elements of the calculation are the two-loop renormalization matrix $Z_{ij}$ for the $|\Delta B|=2$ operators and the $\mathcal{O}(\alpha_s^2)$ corrections to the finite renormalization that ensures the $1/m_b$ suppression of the operator $R_0$ at two-loop order. Our theoretical prediction reads $\Delta\Gamma_s/\Delta M_s = {(4.33\pm 0.93)\cdot 10^{-3}}$ if expressed in terms of the bottom mass in the $\overline{\rm MS}$ scheme and $\Delta\Gamma_s/\Delta M_s = {(4.20\pm 0.95)\cdot 10^{-3}}$ for the use of the potential-subtracted mass. While the controversy on $|V_{cb}|$ affects both $\Delta\Gamma_s$ and $\Delta M_s$, the ratio $\Delta\Gamma_s/\Delta M_s$ is not affected by the uncertainty in $|V_{cb}|$ .

Introduction. The weak interaction of the Standard Model (SM) permits transitions between a neutral B q meson and its antiparticleB q , where q = d or s. The corresponding transition amplitude is mediated by box diagrams with W bosons and up-type quarks u, c, or t on the internal lines. The time evolution of the two-state system (|B q , |B q ) is governed by two hermitian 2×2 matrices, the mass matrix M q and the decay matrix Γ q . By diagonalizing M q − iΓ q /2 one finds the mass eigenstates |B q L and |B q H expressed in terms of the flavour eigenstates |B q , |B q . The mass eigenstates differ in their masses M q H,L and decay widths Γ q H,L with "L" and "H" denoting "light" and "heavy". There are three observables, the mass and width differences ∆M q = M q H − M q L and ∆Γ s = Γ q L − Γ q H as well as the CP asymmetry in flavor-specific decays, a q fs . Experimentally, ∆M q is read off from the B q −B q oscillation frequency, ∆Γ q is found by measuring lifetimes in different decay modes, and a q fs is usually measured through the time-dependent CP asymmetry in semileptonic B q decays. These observables are related to the off-diagonal elements of M q and Γ q as follows: , a q fs = Im with |∆Γ q | 2|Γ q 12 |. M q 12 is sensitive to new physics mediated by particles with masses well beyond 100 TeV. On the contrary, Γ q 12 , probes effects of light new particles with feeble couplings to quarks (see e.g. Refs. [1,2]). While this is one motivation for a more precise SM prediction of Γ q 12 , a better knowledge of Γ q 12 will also help to reveal new physics in M q 12 : Inclusive and exclusive semileptonic B decays give different values for the element |V cb | of the Cabibbo-Kobayashi-Maskawa (CKM) matrix and * gerlach. marvin@protonmail while ∆Γ exp d is still consistent with zero. The precise value in Eq. (2) calls for a better SM prediction of ∆Γ s , which is the topic of this Letter. We specify to q = s from now on.
At one-loop level the SM predictions for Γ s 12 is calculated from the dispersive part of the B s ↔B s amplitude. One must therefore only consider diagrams with light internal u, c quarks; i.e. diagrams with one or two t quarks only contribute to M s 12 . To properly accomodate strong interaction effects associated with different energy scales one employs two operator product expansions (OPE). In the first step one matches the SM to an effective theory with |∆B| = 1 operators [9], where B is the beauty quantum number. The most important operators, i.e. those with the largest coefficients, are the current-current operators Q 1,2 describing tree-level b decays. The effective |∆B| = 1 hamiltonian is known to next-to-leading (NLO) [10-12] and next-to-next-to-leading order (NNLO) [13][14][15] of Quantum Chromodynamics (QCD). The OPE employed in the second step of the calculation is the Heavy Quark Expansion (HQE) [16][17][18][19][20][21][22][23][24] (cf. also [25] for a review), which expresses the B s ↔B s transition amplitude as a series in Λ QCD /m b , where Λ QCD ∼ 400 MeV is the fundamental scale of QCD and m b is the b quark mass. The HQE involves local |∆B| = 2 operators; to find the corresponding Wilson coefficients one must calculate the ∆B = 2 amplitude in both the |∆B| = 1 and |∆B| = 2 theories to the desired order in α s .
The state of the art is as follows: QCD corrections to Γ s 12 are only known for the leading term of the Λ QCD /m b expansion ("leading power"). These include NLO QCD corrections to the contributions with current-current and chromomagnetic penguin operators [26][27][28][29], the corresponding NNLO corrections (and NLO corrections involving four-quark penguin operators) enhanced by the number N f of active quark flavours [30][31][32] as well as NLO results with one current-current and one penguin operator [33] or two penguin operators [34]. The latter paper also presents two-loop results with one or two chromomagnetic penguin operators which are part of the NNLO and N 3 LO contributions. (The four-quark penguin operators Q 3−6 have Wilson coefficients which are much smaller than those of Q 1,2 and the chromomagnetic penguin operator contributes with a suppression factor of α s .) The corrections of Refs. [30] and Refs. [33,34] have been calculated in an expansion in m c /m b to first and second order, respectively. ∆Γ s /∆M s further involves a well-computed ratio of two hadronic matrix elements [35][36][37]. The contribution to Γ s 12 being sub-leading in Λ QCD /m b is only known to LO of QCD [38] and the hadronic matrix elements still have large errors [39].
Both the described perturbative contribution and the power-suppressed term have theoretical uncertainties exceeding the experimental error in Eq. (2). In this Letter we present NNLO QCD corrections to the numerically dominant contribution with two current-current operators and reduce the perturbative uncertainty of the leading-power term to the level of the experimental error.
Calculation. To obtain ∆Γ s /∆M s we use the known two-loop QCD corrections to M s 12 from Ref. [11]. It is convenient to decompose Γ s 12 according to the CKM structures where λ s a = V * as V ab with a = u, c. Γ s 12 is obtained with the help of a tower of effective theories. In a first step we construct a theory where all degrees of freedom heavier than the bottom quark mass m b are integrated out and the dynamical degrees of freedom are given by the five lightest quarks and the gluons. We adopt the operator basis of the |∆B| = 1 theory from Ref. [40]. The matching to the Standard Model happens at the scale µ 0 ≈ 2m W ≈ m t (m t ). Afterwards, renormalization group running is used to obtain the couplings of the effective operators at the scale µ 1 which is of the order m b .
In a next step we perform a HQE which allows us to write Γ s 12 as an expansion in 1/m b . At each order Γ s 12 is expressed as a sum of Wilson coefficients multiplying respective operator matrix elements. The latter has to be computed using lattice gauge theory [35] or QCD sum rules [36,37]. To leading order in the 1/m b expansion we have where ab ∈ {cc, uc, uu}. G F is the Fermi constant and M Bs is the mass of the B s meson. The main purpose of this Letter is the computation of the matching coefficients H ab and H ab S to next-to-next-to-leading order (NNLO) in the strong coupling constant α s . They depend on z = m 2 c /m 2 b . For the ∆B = 1 theory one distinguishes current-current and penguin operators. At leading and next-to-leading orders the current-current operators provide about 90% of the total contribution to Γ ab 12 [34]. Thus, in this work we restrict ourselves to the current-current contributions.
For the calculation of the NNLO corrections one has to overcome several challenges. First, it is necessary to perform a three-loop calculation of the amplitude bs →bs in the ∆B = 1 theory. Sample Feynman diagrams are shown in Fig. 1. In total about 20,000 three-loop diagrams have to be considered which requires an automated setup for the computation. In our case the combination of qgraf [41], tapir [42] and q2e/exp [43,44] turned out to be useful. For the leading term in the HQE we are allowed to set the momentum of the strange quark to zero. Furthermore, we expand in the charm quark mass up to second order, 1 which reduces the integrals to on-shell two-point functions with external momentum q 2 = m 2 b . The propagators inside the loop diagrams are either massless or carry the mass m b . We use FIRE [45] combined with LiteRed [46,47] to reduce all occurring integrals to 23 genuine three-loop master integrals. For the latter analytic results have been obtained with the help of FeynCalc [48][49][50][51], HyperInt [52], PolyLogTools [53] and HyperlogProcedures [54].
On the ∆B = 2 side a two-loop calculation is necessary; sample Feynman diagrams are shown in Fig. 2. From the technical point of view the calculation is significantly simpler. However, in the practical calculation one has to consider three physical and 17 evanescent operators, cf. Ref. [34]). It is necessary to compute the corresponding renormalization constants for the operator mixing up to two-loop order.
The calculation of the ∆B = 2 matrix elements entails a field theoretical subtlety. In fact, in four dimensions there are only two physical operators whereas for the calculation in d dimensions three have to be taken into account. For our calculation it is convenient to choose Q, Q S and R 0 where (i and j are colour indices) and with Note that at lowest order in α s we have α 1 = α 2 = 1 and the matrix element of R 0 is 1/m b suppressed in four dimensions. At higher orders the quantities α 1 and α 2 are chosen such, that the 1/m b -suppression is maintained. The one-loop corrections are known since more than twenty years [26] and the fermionic two-loop terms are available from Ref. [30]. For the NNLO calculation performed in this Letter the α 2 s corrections to α 1 and α 2 are needed.
The 1/m b -suppression of R 0 beyond tree-level is manifest only if one is able to distinguish between ultraviolet (UV) and infrared (IR) divergences, e.g., by regularizing the latter using a gluon mass m g . Otherwise, R 0 develops an unphysical evanescent piece E R0 that scales as m 0 b [34] and hence must be included into the definition of R 0 to obtain correct matching coefficients. One cannot isolate E R0 from R 0 at the operator level, but one can distinguish evanescent and physical pieces in the matrix elements: We use R 0 from Eq. (6) including the finite UV renormalization encoded in α 1 and α 2 in our matching calculation. To this end we have first calculated the linear combination of the renormalized two-loop matrix elements Q (2) , Q S (2) and Q S (2) as given in Eq. (6). After introducing a gluon mass along the lines of Ref. [55] and using Feynman gauge we observe that each of the individual matrix elements becomes manifestly finite upon UV renormalization. α 1 and α 2 to order α 2 s are extracted from the requirement that the linear combination must vanish in the limit m b → ∞.   The matching between the |∆B| = 1 and |∆B| = 2 effective theories is conceptually simple in case IR divergences are not regularized dimensionally. In this case the UV renormalization renders amplitudes of both theories manifestly finite, allowing us to take the limit d → 4, where all matrix elements of evanescent operators vanish. However, for technical reasons we prefer to use = UV = IR , which simplifies the evaluation of the amplitudes but complicates the matching. Following [56] we need to extend the leading order (LO) matching to O( 2 ) and the NLO matching to O( ) in order the determine the NNLO matching coefficients. Furthermore, we need to determine the matching coefficients of both physical and evanescent operators: Since the UV-renormalized amplitudes still contain IR poles, we must keep all matrix elements of evanescent operators until the very end. A powerful cross check of this procedure is the explicit cancellation of the remaining IR poles and of the QCD gauge parameter ξ in the matching.
Results. For our numerical analysis we use the input values listed in Tab. I and the |∆B| = 1 Wilson coefficients from Refs. [13][14][15] and calculate the running and decoupling of quark masses and α s with RunDec [57].
In the following we present the NNLO predictions in three different renormalization schemes for the overall factor m 2 b (cf. Eq. (4)) whereas the quantity z and the strong coupling constant are defined in the MS scheme. The overall factor m 2 b is defined in the MS scheme, as a pole mass, or as a potential-subtracted (PS) mass [61]. The latter is an example of a so-called threshold mass, with similar properties as the pole mass, but is nevertheless of short-distance nature. H ab (z) and H ab S (z) are adapted accordingly, so that the scheme dependence of Γ s 12 is O(α 3 s ). Several renormalization and matching scales enter the prediction for the width difference. We choose µ 0 = 165 GeV for the matching scale between the SM and the |∆B| = 1 theory. Varying µ 0 barely affects ∆Γ s /∆M s and we can keep it fixed. In our nu-merical analysis we identify the matching scale µ 1 and µ b and µ c , the renormalization scales at which m b and m c are defined. We simultaneously vary µ 1 = µ b = µ c between 2.1 GeV and 8.4 GeV with a central scale µ 1 = 4.2 GeV. That is, z enters the coefficients as z = (m c (µ 1 )/m b (µ 1 )) 2 . The |∆B| = 2 operators are defined at the scale µ 2 which has to be kept fixed, because the µ 2 dependence only cancels in the products of H ab (z) and H ab S (z) with their respective matrix elements. In our analysis we set µ 2 = 4.75 GeV which is the bottom quark pole mass m pole b obtained from m b (m b ) with two-loop accuracy. The terms of order Λ QCD /m b in Γ s 12 are only known to LO, so that the µ 1 -dependence of these terms is non-negligible.
We now discuss the results for ∆Γ s /∆M s . In our three schemes we have where the subscripts indicate the source of the various uncertainties. The dominant uncertainty comes from the matrix elements of the power-suppressed corrections ("1/m b ") [35,39]) followed by the renormalization scale uncertainty from the variation of µ 1 in the leading-power term ("scale"). The uncertainties from the leading-power bag parameters ("BB S ") and from the scale variation in the 1/m b piece ("scale, 1/m b ") are much smaller and the variation of the remaining input parameters ("input") is of minor relevance. In Fig. 3 we show the dependence of ∆Γ s /∆M s on the simultaneously varied renormalization scales µ 1 = µ b = µ c for the MS and PS schemes. The small contributions involving four-quark penguin operators are only included at NLO in both the NLO and NNLO curves. Dotted, dashed, and solid curves correspond to the LO, NLO, and NNLO results, respectively. In both schemes one observes a clear stabilization of the µ 1 dependence after including higher orders. Furthermore, we observe that the NNLO predictions (solid lines) in both schemes are close together which demonstrates the expected reduction of the scheme dependence. In the MS scheme we observe that the LO and NLO curves intersect close to the central scale. As a consequence the NLO corrections are relatively small and the NNLO contributions are of comparable size. Close to 9 GeV the NNLO contribution is zero and the NLO corrections amount to about +21%. At the same time the NNLO predictions for µ 1 = 4.2 GeV and µ 1 = 9 GeV differ only by +5% and +9% in the MS and PS schemes, respectively. Note that in the MS scheme the scale dependence of the leading-power term drops from +0 −29 % at NLO to +5 −10 % at NNLO and is now of the same order of magnitude as the ±6% experimental error in Eq. (2). In the PS scheme the scale uncertainty is of the same order of magnitude as in the MS scheme. Note that the scheme dependence inferred from the MS and PS central values in Eq. (8) is only 3%. Eq. (8) clearly shows that one needs better results for the 1/m b matrix elements. A meaningful lattice-continuum matching calls for NLO corrections to the power-suppressed terms, which will further reduce the uncertainty labeled with "scale, 1/m b ".
For the pole scheme we only show the NNLO prediction in Fig. 3. While we also see a relatively mild dependence on µ 1 , the corresponding solid curve lies significantly below the predictions in the MS and PS schemes. This feature can be traced back to the large two-loop corrections in the relation between the MS and the pole bottom quark mass affecting NNLO contributions as much as the genuine NNLO corrections, underpinning the well-known issues with quark pole masses [62][63][64]. For this reason we recommend to not use the pole scheme for the prediction of ∆Γ s .
The most precise prediction for ∆Γ s is obtained from the results in Eq.
The comparison to Eq.
(2) shows that the uncertainty is only about a factor three bigger than from experiment and dominated by the 1/m b corrections. With our NNLO result for Γ q 12 we can also improve the predictions for width difference in the B d −B d system and the CP asymmetries a s fs and a d fs in Eq. (1), whose experimental results are still consistent with zero. We postpone this to a future publication [66].
In Fig. 4 we confront our predictions for the ratio ∆Γ s /∆M s in the MS scheme (green band) with the individual predictions of ∆Γ s and ∆M s . The latter are dominated by the uncertainty in the CKM matrix element V ts which is obtained from V cb through CKM unitarity and cancels in the ratio. Fig. 4 illustrates this feature with |V incl cb | = 42.16(51)10 −3 from [67] and |V excl cb | = 39.36 (68)10 −3 from [68]. The current experimental results for ∆Γ s and ∆M s are indicated by the black bar. Once the prediction of ∆Γ s /∆M s is improved further, it will be possible to test the SM without CKM uncertainty, and with progress on |V cb | one will be able to constrain new physics in ∆M s and ∆Γ s individually.
Conclusions. The SM prediction of ∆Γ s /∆M s based on the long-standing NLO calculation has two sources of uncertainty which exceed the experimental error: the hadronic matrix elements of the power-suppressed operators and the perturbative coefficients, as inferred from the scale and scheme dependences of the calculated result. With the NNLO calculation presented here we have brought the latter uncertainty to the level of the accuracy of the experimental result. For this we had to calculate 20,000 three-loop diagrams and to solve subtle problems related to the interplay of infrared divergences and evanescent oprators. We have pointed out that ∆Γ s adds information to the usual study of ∆M s , because both quantities probe different new-physics scenarios and |V cb | drops out in the ratio ∆Γ s /∆M s .
Acknowledgements. We thank Artyom Hovhannisyan and Matthew Wingate for useful discussions and we are grateful to Erik Panzer and Oliver Schnetz for helpful advice regarding the calculation of the master integrals and the simplification of the obtained results with HyperInt and HyperLogProcedures. VS thanks David Broadhurst for enlightening discussions on iterated integrals. This research was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under grant 396021762 -TRR 257 "Particle Physics Phenomenology after the Higgs Discovery". The Feynman diagrams were drawn with the help of Axodraw [69] and JaxoDraw [70].