NNLO QCD counterterm contributions to $B$ \to $X_s \gamma$ for the physical value of $m_c$

One of the most important O($\alpha_s^2$) corrections to the $B \to X_s \gamma$ branching ratio originates from interference of contributions from the current-current and photonic dipole operators. Its value has been estimated using an interpolation in the charm quark mass between the known results at $m_c=0$ and for $m_c \ll m_b/2$. An explicit calculation for the physical value of $m_c$ is necessary to remove the associated uncertainty. In the present work, we evaluate all the ultraviolet counterterm contributions that are relevant for this purpose.


Introduction
The inclusive decayB → X s γ belongs to the most interesting flavour changing processes. Many popular extensions of the Standard Model (SM) receive important constraints from its branching ratio B sγ . For instance, the 95% C.L. bound it provides on the charged Higgs boson mass in the Two-Higgs-Doublet Model II is in the vicinity of 580 GeV [1,2]. Strengthening such bounds in the future is likely to be achieved with the upcoming precise measurements of B sγ at Belle II [3], as well as with more accurate theoretical calculations within the SM.
The SM prediction for the CP-and isospin-averaged branching ratio B sγ amounts to [4,5] B SM sγ = (3.36 ± 0.23) × 10 −4 , (1.1) for E γ > E 0 = 1.6 GeV. The corresponding experimental world average evaluated by HFAG [6] reads It is based on the Belle [1,7], Babar [8][9][10] and CLEO [11] measurements performed at E 0 ∈ [1.7, 2.0] GeV, and then extrapolated down to E 0 = 1.6 GeV using the method of Ref. [12]. A recent discussion of the energy extrapolation issue and its effect on the averages can be found in Ref. [2]. In that paper B exp sγ = (3.27 ± 0.14) × 10 −4 (1.3) was obtained by naively averaging the measurements at E 0 = 1.9 GeV only [1,[7][8][9][10], and then performing an extrapolation to E 0 = 1.6 GeV. As far as the SM prediction is concerned, it is based on the fact that the hadronic decay rate Γ(B → X s γ) is well approximated by the corresponding partonic one Γ(b → X p s γ), where X p s = s, sg, sgg, sqq, . . .. The difference between the two rates is due to non-perturbative effects. Such effects contribute at around +3% level [13] to the central value in Eq. (1.1), while the associated uncertainty is estimated [14] to be around ±5%. This uncertainty is combined in quadrature with the remaining three ones that are due to the input parameters (±2%), unknown higher-order O(α 3 s ) perturbative effects (±3%), as well an interpolation in the charm quark mass m c (±3%). The interpolation is used to estimate some of the most important O(α 2 s ) corrections.
The purpose of our present work is providing a contribution to removing the latter uncertainty. Perturbative corrections to Γ(b → X p s γ) are most conveniently analyzed in the framework of an effective theory that arises after decoupling the W boson and all the heavier particles. The relevant weak interaction Lagrangian is then given by a linear combination of four-quark and dipole-type operators: L weak ∼ i C i (µ)Q i (see, e.g., Eq. (1.6) of Ref. [5]). For most of our discussion here, only three of them are going to be relevant, namely The MS-renormalized Wilson coefficients C i (µ) are first evaluated at the electroweak scale µ 0 ∼ M W , m t , and then evolved down to µ b ∼ m b using the effective theory renormalization group equations. A complete O(α 2 s ) calculation of all the relevant C i (µ b ) was finalized a decade ago, with the latest contribution amounting to a determination of the necessary four-loop anomalous dimension matrix [15].
Once the Wilson coefficients C i (µ b ) are found, one writes the perturbative rate as 1 where G F is the Fermi constant, and V ij are the Cabibbo-Kobayashi-Maskawa matrix elements. Both the b-quark mass m b and the electromagnetic coupling α em are assumed to be on-shell renormalized. The interference termsĜ ij form a symmetric 8 × 8 matrix. Their perturbative expansion in The Next-to-Next-to-Leading Order (NNLO) corrections to be studied in the present work are given byĜ (2) 17 andĜ 27 . They describe interferences of decay amplitudes generated by the current-current operators Q 1,2 and the photonic dipole operator Q 7 . Such interference terms are most conveniently represented via Cutkosky rules as four-loop propagator diagrams with unitarity cuts. Sample diagrams are presented in Fig. 1. Altogether, around 850 of such four-loop diagrams need to be evaluated.
Despite the fact thatĜ (1,2)7 determine some of the most important NNLO corrections to B sγ (comparable to the present experimental and theoretical uncertainties), their explicit form for the physical value of m c remains unknown, which is due to extreme complexity of the necessary calculation. What has been found so far for arbitrary m c are contributions from diagrams with massless and massive quark loops on the gluon lines [16][17][18]. They are the basis for evaluatinĝ G (1,2)7 in the Brodsky-Lepage-Mackenzie (BLM) [19] approximation. As far as the non-BLM terms are concerned, they have been calculated for m c ≫ 1 2 m b [20] and at m c = 0 [5] only. Next, their values for the physical m c were found using an interpolation, as described in detail in Ref. [5].
The interpolated contributions affect B sγ by (+5 ± 3)%, which implies that their effects on both the central value and the overall uncertainty are sizeable. The interpolation uncertainty could be completely removed via a calculation ofĜ (1,2)7 directly at the physical value of m c . Such a calculation involves the "bare" diagrams like the ones in Fig. 1, as well as many lowerloop diagrams with insertions of ultraviolet counterterms.
In the present paper, we provide results for all the necessary counterterm contributions in the m c = 0 case. Although it is only a first step towards the complete NNLO calculation ofĜ (1,2)7 , it is by no means a trivial one given that two mass scales are present in the relevant threeloop propagator diagrams with unitarity cuts. 2 After generating the considered interference contributions and expressing them in terms of Master Integrals (MIs), we evaluate all the MIs using several methods. In the first approach, we derive a system of Differential Equations (DEs) for the MIs, and solve it numerically starting from boundary conditions at m c ≫ 1 2 m b . The latter are found using asymptotic expansions in such a limit. The relevance of this method lies in the fact that it is basically the only feasible one for the future "bare" calculation. Thus, our current task serves as a preparation for the future enterprise. Apart from the numerical approach, we also use analytical methods for determining all the necessary MIs, either in closed forms or as expansions in z = m 2 c /m 2 b . Similarly to the m c = 0 calculation of Ref. [5], we restrict ourselves to the limit E 0 = 0. Considering such a limit is a necessary step on the way to the E 0 = 0 case. Results in the latter case are going to be obtained in the future by subtracting the low-E γ tail from the total rate, in the same way as it was done forĜ 77 [21,22].
Our article is organized as follows. In Section 2, we define all the relevant counterterm contributions via an explicit formula for the renormalizedĜ (1,2)7 . Expressions for them in terms of the MIs, as well as DEs for the MIs are given in Section 3. Our numerical results for all the considered counterterm contributions are presented in Section 4. Section 5 is devoted to presenting our final results as expansions in z. We conclude in Section 6.

The NNLO renormalization formula forĜ 17 andĜ 27
The renormalizedĜ (2) i7 for i = 1, 2 are obtained from the bare onesĜ (2)bare i7 and the ultraviolet counterterms. The necessary expression in the m c = 0 case together with all its ingredients was presented in Section 2.2 of Ref. [5]. Following that article, we skip the contributions that are due to charm quark loops on the gluon lines together with the corresponding counterterms. Such NNLO contributions have been already calculated for arbitrary m c in Ref. [18], and are included in the phenomenological analysis of Ref. [5] (Eq. (3.8) there). Under such assumptions, the renormalization formula for arbitrary m c takes the following form (for i = 1, 2): (1) Figure 2: Sample diagrams forĜ (1)bare 27 with unitarity cuts indicated by the dashed lines.
It turns out to differ from the m c = 0 case only in the last line where the renormalization of m c is taken into account. Moreover, we have simplified a few details of the notation without changing the actual content of the equation. Its l.h.s. corresponds to such a renormalization scheme where the Wilson coefficients, α s and m c are MS-renormalized at the scale µ b . On the other hand, the b-quark mass and the external quark fields are on-shell renormalized, which is indicated by the superscript OS at the corresponding renormalization constants. In these constants (Eq. (2.8) of Ref. [5]), as well as in Eq. (2.1) here, one should substitute The remaining renormalization constants are the MS-scheme ones (Eq. (2.9) of Ref. [5]). In The bareĜ kl on the r.h.s. of Eq. (2.1) are assumed to have been evaluated for the renormalization scale µ 2 = e γ m 2 b /(4π), where γ is the Euler-Mascheroni constant. They involve indices corresponding to the operators listed in Eqs. (1.6) and (2.5) of Ref. [5]. In some of the cases, the superscript "bare" has been replaced by either "3P " or "m". In these cases, either only the three-particle cuts were included (3P ), or one of the b-quark propagators was squared (m). The latter operation is the simplest way to obtain the counterterms that originate from renormalization of the b-quark mass.
Explicit expressions for all the necessary bareĜ kl in the m c = 0 case have been given in Eqs. (2.3)-(2.7) of Ref. [5]. Only a few of them get modified in the m c = 0 case, namelŷ In the following, we shall calculate all of them except forĜ . Some of the considered quantities are related to each other in a very simple manner. In particular (1)bare 7 (12) . It is due to the fact that in all the non-vanishing Feynman diagrams that contribute to these quantities (e.g., those in Fig. 2), the internal gluon couples only once to the charm quark loop. In consequence, the diagrams with Q 1 and Q 2 differ only by a colour factor of − 1 6 that comes from the identity T a T b T a = − 1 6 T b for the SU(3) generators. The case is analogous for the evanescent operators Another simple identity allows us to expressĜ (1)bare 7 (12) in terms of various contributions tô stand for the two-particle cut contributions toĜ that are proportional to the down-type and up-type quark electric charges, respectively. The identity (2.4) can be recovered after expressing all the involved quantities in terms of the MIs, but without actually calculating any of the MIs.
Thus, what remains to be calculated are the six quantities on the r.h.s. of the following two equations: . (2.6) Analogously to the first of the above equations, the three terms on the r.h.s. of the second one are defined according to the quark electric charges, and to the number of final-state particles.
All the considered quantities need to be found up to O(ǫ). In the case of Eq. (2.5), it is equivalent to extending the NLO calculations of Refs. [23][24][25][26][27] to one more order in ǫ. Such an extension is actually already available from Ref. [28], and we shall compare our results to that article. As far as the counterterms in Eq. (2.6) for m c = 0 are concerned, our calculation is new.

Master integrals and differential equations
We generate the relevant Feynman diagrams using FeynArts [29]. The Dirac algebra after summing/averaging over polarizations is performed with fully anticommuting γ 5 , which is allowed because only traces with even numbers of γ 5 can give non-vanishing contributions in the considered problem. Once the results are expressed in terms of scalar integrals, we reduce them to the MIs using Integration By Parts (IBP) identities, with the help of the code FIRE [30]. This way we express the six quantities on the r.h.s. of Eqs.
where p is the external momentum, while k i are the cut propagator momenta (directed from left to right in all the MIs of Tab. 1). In the above expressions, integrals over dPS n stand for the n-body massless phase space integrals, with n = 2, 3. They can be simplified [32] to integrals overŝ ij = 2(k i k j )/m 2 b for arbitrary integrands W 2 (ŝ 12 ) and W 3 (ŝ 12 ,ŝ 13 ,ŝ 23 ): In four spacetime dimensions, the global normalization of all the considered contributions to the decay rate is fixed by Eq. (1.5). In D = 4 − 2ǫ dimensions, we follow Ref. [5] where the global normalization is fixed bŷ The factor (1−ǫ) above comes from the Dirac algebra in | sγ|Q 7 |b | 2 summed over polarizations.
To evaluate the MIs, we differentiate them with respect to z, and reduce the resulting integrals to the MIs again, using the IBP identities. This gives us a closed set of DEs which takes the following form (I ′ k = dI k /dz): The necessary boundary conditions for the above DEs are evaluated in the z ≫ 1 limit, in which the integrals become considerably simpler. Large-z expansions of the three-body MIs (I 14 -I 18 ) are obtained in a straightforward manner using the Feynman parameterization. In the two-body case (I 1 -I 13 ), we have used the code exp [33,34] that applies the asymptotic expansion method in a fully automatic manner.

Numerical results for the counterterm contributions
In our first approach, we expanded all the MIs in ǫ up to the relevant orders, which gave us a system of coupled ordinary DEs for 76 functions of z being a single variable. The boundary condition was chosen at z = 20, and evaluated using the expansions up to O(1/z 8 ). Since the coefficients of the DEs (3.6) have singularities at z ∈ {0, 1 4 , 1}, the variable z must be treated as complex, and a contour in the complex z-plane must be chosen for the numerical solution. In practice, we solved the DEs along ellipses in the lower half-plane of z, following the Feynman prescription m 2 c → m 2 c − iε in the propagators, and treating m 2 b in z = m 2 c /m 2 b as a real rescaling factor. In this way, non-vanishing imaginary parts in some of the MIs were obtained on the real axis for z < 1 4 (below the cc production threshold), despite the fact that both the boundary conditions at z = 20 and the DE coefficients were purely real. In particular, the proper imaginary parts of the functions a(z) and b(z) from Eqs. (3.3)-(3.4) of Ref. [27] were recovered. Such imaginary parts drop out from our final results in Eqs. (2.5)-(2.6) but, nevertheless, they provide convenient cross-checks of the calculation.
The numerical solution was obtained using the FORTRAN code ZVODE [35], upgraded to quadrupole-double precision with the help of the QD computation package [36]. The physical value of z is around 0.06. However, apart from the vicinity of this point, we have considered a wide range of the final values, in order to test the limits at z → 0, as well as the behaviour around the threshold at z = 1 4 . To present our numerical results for the quantities listed in Eq. (3.1), we parameterize them in terms of nine functions of z as followŝ Their dependence on z is shown in Fig. 3. Each of the (blue) dots corresponds to a particular final value of the numerical solution of the DEs along an ellipse in the complex plane. The physical point around z = 0.06 is marked by a slightly bigger (red) dot. Similar (red) dots on the vertical axes are used to indicate the limits at z → 0 in all the cases when they are finite, i.e. for all the functions except those in Eq. (4.3). They can be taken over from the calculation performed at z = 0 from the outset [5]: {f 0 , f 1 , g 0 , g 1 , r −1 , r 0 , r 1 } where ζ k ≡ ζ(k). In the case ofĜ The above expression contains a 1/ǫ divergence, contrary to Eq. (4.3). It arises due to the logarithmic divergences of the functions j i (z) at z → 0. Such non-commutation of limits is frequently encountered in the framework of dimensional regularization. Our expansions above z = 20 that have served as the boundary conditions for the DEs are displayed as solid (blue) lines in the first six plots of Fig. 3. The remaining (green) solid lines describe expansions either in z (for z < 1 4 ) or in 1 z (for z > 1 4 ) that have been obtained in our second approach -see Section 5. For f i (z), the expansion plots are terminated away from the threshold at z = 1 4 , as they become inaccurate in its vicinity. In the three-body cases (g i (z), j i (z)) the applied expansions are so accurate that their mismatch at z = 1 4 is invisible within the plot resolution. 3 We refrain from presenting here any numerical fits for the functions plotted in Fig. 3. Instead, either expansions or exact expressions for them will be given in Section 5. It is evident from the plots that our expansions around z = 0 are sufficiently precise for phenomenological purposes in the vicinity of the physical point.

Expansions in z
In this section, we collect our final results for the small-z expansions of the functions defined in Eqs. (4.1)-(4.4). They have been obtained in Ref. [37] with the help of Feynman parameterizations, Mellin-Barnes techniques [38], as well as our differential equations in Eq. (3.6) that allowed to extend the expansions to higher orders.
It is likely that fully analytical results for all the considered counterterm contributions could be found using the so-called canonical basis [40,41] of the MIs, and solving the DEs order-byorder in ǫ, in terms of iterated integrals. However, at the time the main part of our calculation was performed, no automatic tools for constructing canonical bases at three loops (see, e.g., [42,43]) were available. On the other hand, the expansions in z are certainly sufficient for all the phenomenological purposes in our case, while the numerical method is most suitable for preparing the bare calculation.

Summary
We calculated all the relevant counterterm contributions toĜ (2) 17 andĜ (2) 27 at the NNLO for the physical value of the charm quark mass m c . In combination with the future bare calculation of these quantities, the current results are necessary for removing one of the most important perturbative uncertainties in the SM prediction B sγ . In our first approach, we derived differential equations for the master integrals, and solved them numerically using the boundary conditions at large m c . In the second approach, a combination of various techniques was used to determine the master integrals in terms of power-logarithmic expansions around m c = 0. In most of the cases, logarithmic divergences at m c → 0 are absent, and our results properly converge to the m c = 0 ones previously found in Ref. [5].
acknowledge support from the National Science Centre (Poland) research project, decision no DEC-2014/13/B/ST2/03969. The research of M.S. has been supported by the BMBF grant 05H15VKCCA.