Towards the NNLO theory prediction for the width differ-ence ∆Γ s

. The width di ↵ erence ∆Γ s that can be extracted from lifetime measurements of the two mass eigenstates of the B 0 s − ¯ B 0 s system is one of the key ﬂavor precision observables and has been experimentally measured at per cent level accuracy. The current theory prediction is much less accurate and a sizable reduction of scale uncertainties can only be achieved by means of eval-uating the uncalculated 2- and 3-loop QCD corrections. This is precisely the issue addressed in this work where we report on the results that have been obtained so far and explain some of the technical and conceptual challenges that we encountered in the course of our calculations.


Introduction
Precision flavor observables constitute an important tool in the task of probing our understanding of the Standard Model (SM) and thus constraining discrepancies between theory and experiments that may be interpreted as manifestations of beyond Standard Model (BSM) physics. The B 0 s −B 0 s system happens to be an ideal laboratory for such investigations, where we can extract three relevant physical quantities that are equally well accessible to theoretical calculations and experimental measurements. These are the oscillation frequency ∆m s , the width di↵erence ∆Γ s and CP asymmetry in flavor-specific decays a s fs . Theoretical investigations of ∆Γ s allow us to quantify how well our perturbative calculations describe the physics that arises from the SM flavor sector. In principle, light BSM particles that feebly interact with SM model fields could induce a shift in the SM value of the width di↵erence. However, one normally tends to regard ∆Γ s as a SM precision probe, since its sensitivity to BSM physics (in particular the e↵ects of heavy particles) is rather limited.
is much more precise than the most up-to-date theory prediction [5][6][7][8][9][10] ∆Γ MS s = (0.088 ± 0.011 pert. ± 0.002 B,B S ± 0.014 ⇤ QCD /m b ) ps −1 , ⇤ e-mail: v.shtabovenko@kit.edu where we would like to attract the attention of the reader to the first uncertainty (denoted as "pert.") that stems from the uncalculated perturbative QCD corrections to the B s -meson mixing at 2 and 3 loops. Our project aims at calculating these corrections analytically and providing new theory predictions for ∆Γ s that should hopefully reduce the existing perturbative uncertainties and allow for more meaningful comparisons between experimental and theoretical values.

Calculation
Let us begin by explaining the procedure needed to obtain ∆Γ s . The main ingredient is a matching calculation between two e↵ective theories, where on the |∆B| = 1 side we have [11] H |∆B|=1 with λ s a = V ⇤ as V ab , where V i j are CKM matrix elements. The Wilson coefficients C i are calculated from matching SM to the |∆B| = 1 e↵ective Hamiltonian and are known at NNLO accuracy [12]. Eq. (4) does not include evanescent operators E[Q i ] [13,14], although they are of course properly taken into account in our calculation. Evanescent operators formally vanish in 4 dimensions, as they are of O("), where " is the dimensional regularization parameter from d = 4 − 2". Their presence encodes the ambiguous nature of 4-dimensional identities such as the Chisholm identity γ µ γ ⌫ γ λ = g µ⌫ γ λ + g ⌫λ γ µ − g µλ γ ⌫ + i" µ⌫λσ γ σ γ 5 (5) or Fierz relations in d dimensions. The fact that a matrix element of an evanescent operator multiplying an " pole yields something of O(" 0 ) explains why the treatment of such operators in higher order perturbative calculations requires great care. This is also one of the reasons why the present calculation turns out to be highly nontrivial from the conceptual point of view. We refer to our upcoming publication for the complete list of the |∆B| = 1 operators relevant for our project at 3 loops. The operator basis at 2 loops is given in [15]. As far as the |∆B| = 2 side of the matching is concerned, the quantity we need to consider is Γ 12 from the relation ∆Γ s ⇡ 2|Γ 12 | which depends on the absorptive part of a bilocal matrix element constructed from time-ordered insertions of two H |∆B|=1 e↵ . Involving Heavy Quark Expansion (HQE) [16][17][18][19][20][21][22][23][24] one can rewrite Γ 12 as [8] with where z ⌘ m 2 c /m 2 b . Our main interest is devoted to the matching coefficients H ab (z) and e H ab S (z) which we want to determine at NNLO accuracy by calculating contributions from perturbative QCD corrections to the relevant operators on both sides of the matching. The 4-fermion operators Q andQ S are defined as  where i and j denote the fundamental color indices of the b and s quarks. At intermediate stages of our calculation tree-level matrix elements of further |∆B| = 2 arise so that it is useful to mention them here Apart from that we also must properly account for the 1/m b suppressed operator R 0 since, as has already been observed in [5], the unrenormalized matrix element of this operator does not exhibit the 1/m b suppression. We would like to stress that the proper treatment of these contributions from R 0 is crucial to obtain correct matching coefficients at NLO and beyond. Finally, it is worth mentioning that the |∆B| = 2 NLO operator basis also comprises a suitable set of evanescent operators which can be found in [15]. The NNLO operator basis will be given in [25].
We choose to carry out our matching calculation on-shell, i. e. with p 2 b = m 2 b and p 2 s = m 2 s . Since the s quark is taken massless we can choose a kinematics that allows us to set p s = 0, which reduces our calculation to the evaluation of 2-point functions. Even though we do not have external charm legs, m c still enters our calculation due to internal charm lines induced by |∆B| = 1 operators and gluon self-energy subloops. From the calculational point of view, the simplest situation arises in the massless charm limit which makes our integrals depend only on a single scale m b . This is indeed the path we choose at 3 loops, where our matching coefficients are accurate up to O(z 0 ). At 2 loops we, however, perform an asymptotic expansion in z up to first order.
At 2 loops we consider all possible insertions of |∆B| = 1 operators, i. e. permutations of Q 1,2 , Q 3−6 and Q 8 in either of the two vertices. This also comprises the 2-loop double insertion of Q 8 that, strictly speaking, contributes to NNNLO. As far as the 3 loop calculation is concerned, only the double insertion Q 1,2 (as the dominant contribution) it taken into account. We then match the so-obtained 2-and 3-loop |∆B| = 1 amplitudes to the 1-and 2-loop |∆B| = 1 amplitudes respectively, as can be inferred from Fig. 1.

Technical details
Let us briefly comment on the calculational setup employed in this computation. All |∆B| = 1 and |∆B| = 2 Feynman diagrams are generated with Qgraf [26], while the insertion of Feynman rules and the topology identification are done using our in-house codes q2e/exp 3 EPJ Web of Conferences 258, 04002 (2022) https://doi.org/10.1051/epjconf/202225804002 A Virtual Tribute to Quark Confinement and the Hadron Spectrum (vConf21) [27,28] or tapir [29]. Then, our well-tested calc framework written in FORM [30] takes care of the amplitude evaluation, asymptotic expansions and further simplifications.
We follow two di↵erent strategies to deal with tensor multi-loop integrals appearing in this calculation. In the one case we construct a set of suitable Dirac and color projectors (cf. appendix of [15]) that leave us with scalar loop integrals only. In the other case we calculate all required tensor reduction formulas [31] using FeynCalc [32][33][34] and Fermat [35] and export them to FORM. We explicitly verify that both approaches lead to the same results, including cases of diagrams with Dirac structures featuring a high number of Dirac matrices (complicated projectors) and diagrams with multiple gluon propagators with a generic gauge parameter (high rank tensor integrals). The implementation of the Feynman rules for various 4-fermion operators was also cross checked (on the level of selected diagrams) using Feyn-Rules [36], FeynArts [37] and FeynCalc. The IBP reduction [38,39] of the occurring loop integrals is carried out using FIRE [40] and LiteRed [41].
As far as the master integrals are concerned, we mainly need to deal with single scale propagator-type on-shell integrals. Up to 2 loops the corresponding analytic results can be readily extracted from the literature [42,43]. At 3 loops we encountered many integrals that do not seem to be known in the literature. Luckily, all of them can be evaluated by directly integrating the Feynman parametric representation (derived using Feyncalc) with the aid of HyperInt [44]. To further simplify the HyperInt-results written in terms of Goncharov Polylogarithms (GPLs) [45] containing 6th root of unity, we made use of HyperLogProcedures [46] and PolyLogTools [47]. Finally, every analytic result has been extensively checked numerically with the aid of pySecDec [48][49][50] and FIESTA [51,52].
Some new FeynCalc functions that were purposely developed during author's work on this project will become an official part of the upcoming FeynCalc 10 [53] but are already publicly available and documented in the development version of the package. These include FCFeynmanParametrize (Feynman parametrization of loop integrals), FCLoopIntegralToGraph (graphical representation of loop integrals), FCLoopFindIntegralMappings (mappings between master integrals similar to FindRules in FIRE) or FCMatchSolve (extraction of matching coefficients or renormalization constants).

Renormalization and matching
The renormalization of the bare amplitudes on both sides of the matching is rather involved due to the large number of sources of UV-divergences. First of all, we employ the usual QCD renormalization of quark fields, masses, ↵ s and the gauge parameter ⇠, where the latter two are relevant only at 3 loops for the |∆B| = 1 and at 2 loops for the |∆B| = 2 amplitudes. In addition to that, we must also take care of the Wilson coefficients of the e↵ective operators that, as it is customary in EFTs, happen to mix under the renormalization. More explicitly, for each of the two e↵ective theories we have where W i denote the Wilson coefficients, while the renormalization matrix Z consists of four submatrices describing the mixing of di↵erent operator types. In the case of Z QQ and Z EE we have the mixing of physical and evanescent operators among themselves, while Z EQ governs the mixing of evanescent operators into the physical ones. The submatrix Z QE is special in the sense that it contains not only poles but also finite pieces. Those are chosen such, that in the limit d ! 4 the contributions of evanescent operators to the renormalized amplitude must vanish. The O(↵ 2 s ) results for Z |∆B|=1 can be found in [54]. On the |∆B| = 2 side of the matching we work with di↵erent sets of physical operators and investigate di↵erent choices of evanescent operators. Each choice leads to a di↵erent 2-loop renormalization matrix, so that we find it more convenient to compute Z |∆B|=2 in a separate calculation for each relevant operator basis.
While the renormalization procedure renders the bare amplitudes UV-finite, in general they still contain IR divergences. Those can be regularized either dimensionally or by giving the gluon a finite mass. The former makes the evaluation of the amplitudes simpler but complicates the matching procedure. The main reason for this is the appearance of tree level matrix elements of evanescent operators multiplying IR poles, e. g. hE (1) 1 i (0) /", in the renormalized amplitudes. Such quantities are formally of O(" 0 ) so that they do not vanish in the limit d ! 4. However, if the matching procedure is carried out in a proper way, these terms will cancel in the matching as they should. To this end one has to write down the |∆B| = 2 matching coefficients as an expansion not only in ↵ s but also in ", i. e.
Notice that this applies not only to the physical but also to the evanescent operator coefficients. For an NNLO calculation, at LO the matching has to be carried out up to O(" 2 ), while at NLO the O(" 1 ) accuracy is required. Finally, when matching the O(↵ 2 s ) amplitudes, it suffices to determine the Wilson coefficients at O(" 0 ). This procedure is explained in [6] and we have explicitly verified that it works not only at 2 but also at 3 loops.
On the other hand, one might also choose to work with a finite gluon mass, which is unproblematic for diagrams of operator insertions that do not contain 3-and 4-gluon vertices i. e. for most of the 2-loop |∆B| = 1 and for all 1-loop |∆B| = 2 diagrams. In this case we need to deal with an additional scale m g appearing in our master integrals but the matching itself becomes much simpler. This is due to the fact that the UV-renormalized amplitudes are manifestly free of " poles, so that one can immediately take the limit d ! 4, in which all matrix elements of evanescent operators vanish. Furthermore, it is worth noting that it is not necessary to compute 2-loop master integrals with full the m g dependence, as we are only interested in the logarithms of m g that capture the IR divergences. Here it is fully sufficient to perform an asymptotic expansion in m g , where we make use of the program asy.m [55,56] to reveal the contributing regions of each integral under consideration.
In our calculation we have explicitly verified that for selected 2-loop |∆B| = 1 correlators the matching with m g = 0 and m g , 0 leads to exactly the same matching coefficients. In the 3-loop case we contented ourselves with the dimensional regularization of IR divergences and the fact that all IR poles explicitly cancel in the matching.

Results
An overview of the new matching coefficients that were obtained in our calculation is given in Table 1, where the first column denotes relevant operator insertions on the |∆B| = 1 side of the matching, the second column provides the existing literature result, and the third column highlights our contribution. As one can easily see, in most cases only the so-called n fpiece, i. e. the fermionic contributions, are currently known in the literature. In this sense our work constitutes an important milestone in the task of obtaining complete predictions that necessarily must incorporate both fermionic and nonfermionic contributions. We would also like to stress that up to O(z) at 2 loops and O(z 0 ) at 3 loops we explicitly confirm the correctness of all previously computed matching coefficients. This is also true for the n fpiece of the 3-loop correlator Q 1,2 ⇥ Q 1,2 computed in [9], where the check required us to use the |∆B| = 1 basis transformation formulas at O(↵ 2 s ) from [12] in order to convert our result in the operator basis of [11] to the operator basis from [58]. Due to the still ongoing analytic and numerical cross checks of the obtained results, we are not yet in the position to provide an updated theory prediction for ∆Γ s at NNLO accuracy. The 2-loop matching coefficients for the Q 1,2 ⇥ Q 3−6 insertions have been recently made public [15], while the remaining 2and 3-loop contributions are expected to appear in the near future [25,59].
In [15] we have shown that the inclusion of nonfermionic contributions to the 2-loop result for Q 1,2 ⇥ Q 3−6 induces a sizable shift of ∆Γ s as compared to the 1-loop result. To this end it is sufficient to consider the ratio between full ∆Γ s (comprising Q 1,2 ⇥ Q 1,2 and Q 1,2 ⇥ Q 3−6 contributions up to O(↵ s ) and Q 3−6 ⇥ Q 3−6 contributions up to O(↵ 0 s )) and the Q 1,2 ⇥ Q 3−6 piece only. With the 1-loop result for Q 1,2 ⇥ Q 3−6 we find More details on the numerical ingredients of these comparisons can be found in [15]. The captions "MS" and "pole" refer to the way how we treat the m 2 b prefactor in Eq. (7). Here we have the choice between regarding it as an MS or an on-shell mass. We would like to stress that in both cases all masses apart from this m 2 b prefactor are always treated in the MS scheme.

Summary
We reported on the current status of our project which aims to achieve a sizable reduction of the perturbative uncertainties in ∆Γ s . This is needed to match the current experimental precision of the width di↵erence in B 0 s −B 0 s mixing and can be accomplished by extending the available matching calculations between |∆B| = 1 and |∆B| = 2 e↵ective theories to 2 and 3 loops for the relevant operator insertions.
At the first stage of this work we have already carried out the corresponding calculations analytically by considering an expansion z ⌘ m 2 c /m 2 b up to O(z) and O(z 0 ) at 2 and 3 loops respectively. While a part of our results has already been published in [15], we still need to do more checks on the remaining contributions, especially at 3 loops. In the near future we plan to publish full analytic results [25,59] for all of the new matching coefficients obtained in the course of this work and to provide a new update of the theory prediction for ∆Γ s .