Towards B ->X_s gamma at the NNLO in QCD without interpolation in m_c

Strengthening constraints on new physics from the B ->X_s gamma branching ratio requires improving accuracy in the measurements and the Standard Model predictions. To match the expected Belle-II accuracy, Next-to-Next-to-Leading Order (NNLO) QCD corrections must be calculated without the so-far employed interpolation in the charm-quark mass m_c. In the process of evaluating such corrections at the physical value of m_c, we have finalized the part coming from diagrams with closed fermion loops on the gluon lines that contribute to the interference of the current-current and photonic dipole operators. We confirm several published results for corrections of this type, and supplement them with a previously uncalculated piece. Taking into account the recently improved estimates of non-perturbative contributions, we find B_{s gamma} = (3.40 \pm 0.17) * 10^{-4} and R_gamma = B_{(s+d) gamma}/B_{c l nu} = (3.35 \pm 0.16) * 10^{-3} for E_gamma>1.6 GeV in the decaying meson rest frame.


Introduction
Flavour Changing Neutral Current (FCNC) processes receive the leading Standard Model (SM) contributions from one-loop diagrams only, often with additional suppression factors originating from the Glashow-Iliopoulos-Maiani (GIM) mechanism [1]. It makes them sensitive to possible existence of new weakly-interacting particles with masses ranging up to O(100 TeV). Significant deviations from the SM predictions are observed in the GIM-unsuppressed FCNC processes mediated by the b → sµ + µ − transition (see, e.g., the recent summary in Ref. [2]). On the other hand, no deviations are seen in the closely related b → sγ transition, despite higher accuracy of both the measurements and the SM predictions in its case.
The physical observable giving the strongest constraints on the b → sγ amplitude is the inclusive B sγ branching ratio, i.e. the CP-and isospin-averaged branching ratio ofB → X s γ and B → Xsγ decays, withB and B denoting (B 0 or B − ) and (B 0 or B + ), respectively. The states X s and Xs are assumed to contain no charmed hadrons. B sγ is being measured [3][4][5][6][7][8] with E γ > E 0 for E 0 ∈ [1.7, 2.0] GeV, and then extrapolated to the conventionally chosen value of E 0 = 1.6 GeV to compare with the theoretical predictions (that would be less accurate at higher E 0 ). The current experimental world average for B sγ at E 0 = 1.6 GeV reads (3.32 ± 0.15) × 10 −4 [9], which corresponds to an uncertainty of around ±4.5%. With the full Belle-II dataset, the world average uncertainty at the level of ±2.6% is expected [10,11]. Achieving a similar accuracy in the SM predictions is essential for improving the power of B sγ as a constraint on beyond-SM theories. It is the goal of the calculations we describe in what follows.
The SM prediction for B sγ (see Refs. [12,13]), is based on the formula where α em = α on shell em , while the so-called semileptonic phase-space factor C is given by Its numerical value is determined [14] using the Heavy Quark Effective Theory (HQET) methods from measurements of theB → X c ℓν decay spectra. The quantity P (E 0 ) is defined through the following ratio of perturbative inclusive decay rates of the b quark: with X p s and X p u denoting all the possible charmless partonic final states in the respective decays (X p s = s, sg, sqq, . . .). The non-perturbative contribution from N(E 0 ) in Eq. (1.1) is estimated 1 at the level of around 4% of B sγ . To achieve O(3%) precision in P (E 0 ), evaluation of the Next-to-Next-to Leading (NNLO) QCD corrections to this quantity is necessary.  Perturbative calculations of P (E 0 ) are most conveniently performed in the framework of an effective theory obtained from the SM via decoupling of the W boson and all the heavier particles. The relevant weak interactions are then given by the following Lagrangian density 2 Evaluation of the Wilson coefficients C i to the NNLO accuracy (O(α 2 s )) at the renormalization scale µ b ∼ m b required computing electroweak-scale matching up to three loops [15], and QCD anomalous dimensions up to four loops [16]. Since C i in the SM have no imaginary parts, one can write the perturbative decay rate as whereĜ ij come from interferences of amplitudes with insertions of the operators Q i and Q j . The dominant NNLO effects come fromĜ 17 ,Ĝ 27 andĜ 77 that originate from the operators WhereasĜ 77 has been known up to O(α 2 s ) since a long time [17][18][19][20][21], no complete NNLO calculation ofĜ 17 andĜ 27 at the physical value of the charm quark mass m c has been finalized to date. Instead, calculations of these quantities at m c ≫ m b [22,23] and m c = 0 [13] gave a basis for estimating their physical values using interpolation [13]. The related uncertainty in B sγ (due to the m c -interpolation only) has been estimated at the level of ±3%, which places it among the dominant contributions to the overall theoretical uncertainty (see Sec. 3). By analogy to what has been done in theĜ 77 case [17][18][19][20][21], evaluation of O(α 2 s ) contributions toĜ 27 is performed in two steps. First, no restriction on the photon energy E γ is assumed. Next, one performs the calculation for E γ < E 0 , which requires considering diagrams with three-and four-body cuts only. The desired resultĜ is then obtained without necessity of determining the differential photon spectrum close to the endpoint E max γ = 1 2 m b . In the present paper, we describe our calculation ofĜ (2) 27 in at the physical value of m c , and with no restriction on E γ . Final results are presented for contributions originating from diagrams with closed fermion loops on the gluon lines, like those in the first row of Fig. 1. They undergo separate renormalization and are gauge invariant on their own, so they serve as a useful test case for our calculation of the completeĜ 27 . Most of such contributions have already been determined in the past [24][25][26][27] and implemented in the phenomenological analysis [12,13]. We confirm the published results, and supplement them with a previously uncalculated piece. Some of the previous results have been obtained by a single group only, which makes our verification relevant.
The article is organized as follows. In the next section, our algorithm for evaluation of the completeĜ (2) 27 is sketched, and the current status of the calculation is summarized. Next, we focus on the closed fermionic loop contributions, displaying our numerical results and comparing them with the literature wherever possible. In Sec. 3, the SM prediction for the branching ratio is updated, taking into account the recently improved estimates of non-perturbative effects [28]. We conclude in Sec. 4. In the Appendix, large-z expansions of our final results are presented, and one of the counterterm contributions is discussed.
27 from diagrams with closed loops of massless fermions -see the text. They have already been multiplied by n l = 3, i.e. the number of flavours we treat as massless.

The NNLO contribution toĜ 27
The quantityĜ (2) 27 is given by a few hundreds of four-loop propagator diagrams with unitarity cuts, as those presented in Fig. 1. We generate them using QGRAF [29] and/or FeynArts [30,31]. After performing the Dirac algebra with the help of FORM [32], we express the fullĜ (2) 27 in terms of several hundred thousands scalar integrals grouped in O(500) families. 3 Next, the Integration By Parts (IBP) identities [33][34][35] for each family are generated and applied using KIRA [36,37], as well as FIRE [38,39] and LiteRed [40,41]. In effect,Ĝ where the rational functions R kl (z, ǫ) on the r.h.s. are determined [42][43][44] from the IBP, too. 4 Similar equations are explicitly displayed in Eq. (3.6) of Ref. [45] where ultraviolet counterterm contributions toĜ (2) 27 have been determined. We solve the DEs using the same method as in Refs. [26,45,46]. The MIs are expanded in ǫ to appropriate powers, with the expansion coefficients being functions of z only. Boundary conditions for these functions at large z are found using asymptotic expansions [47]. Next, the variable z is treated as complex, and the DEs are numerically solved along half-ellipses in the z-plane, to bypass singularities on the real axis.
In practice, the codes q2e and exp [48,49] are used to determine the asymptotic expansions at large z. Coefficients at subsequent powers of 1/z are given in terms of one-, two-and threeloop single-scale integrals, either massive tadpoles or propagator-type ones with unitarity cuts (see Fig. 2). Only at the level of the latter integrals, we perform cross-family identification, which gives us O(50) essentially different and non-vanishing integrals. They are evaluated [50] using various techniques, in particular the Mellin-Barnes one. Once the large-z expansions are found, numerical solutions of the DEs starting from the boundary at z = 20 are worked out using the code ZVODE [51] upgraded to quadrupole-double precision with the help of the QD [52] computation package. Half-ellipses of various sizes are considered to test the numerical stability.
At present, our IBP reduction for the fullĜ (2) 27 is (almost) completed, and the evaluation of the boundary conditions is being finalized [50]. However, for the diagrams with closed fermionic loops (as the ones in the first row of Fig. 1   The displayed results correspond to various contributions toĜ . The renormalization has been performed with the help of the counterterm contributions evaluated 5 in Refs. [45,46]. In all the plots, the black dots correspond to numerical solutions that we have obtained using the DEs. Dots corresponding to the physical value of z are bigger and highlighted in red. Blue dots of similar size on the left boundaries of each plot indicate the z → 0 limits for each contribution, known from the calculation in Ref. [13]. Thin dashed curves continuing to large values of z describe our large-z expansions evaluated up to O(1/z 2 ) (see the Appendix). The dash-dotted vertical lines indicate the cc production threshold at z = 1/4, in the vicinity of which neither the large-z nor the small-z expansions are expected to converge well.
In Fig. 3, three distinct contributions from diagrams with closed massless fermion loops are presented. The first (upper left) plot corresponds to diagrams with two-body cuts. The thin dashed line in the small-z region shows the analytic expansion in powers of z evaluated in Ref. [25]. It is the only case for which such an expansion is known. The solid blue curve shows the numerical fit corresponding to Eq. (3.2) of Ref. [26] where a numerical method (identical to ours) has been used.
The second (upper right) plot of Fig. 3 shows all the four-body-cut contributions except the diagrams displayed in Fig. 5. The latter diagrams have been skipped in evaluating the photon spectrum in the Brodsky-Lepage-Mackenzie (BLM) [53] approximation by the authors of Refs. [24,27]. The solid blue curve is based on the numerical fit from Eq. (3.6) of Ref. [13] that corresponds to no restriction on E γ , and has been obtained as a by-product of the calculation in Ref. [27].
The third (bottom) plot in Fig. 3 corresponds to the very diagrams from Fig. 5. In this case, no numerical result valid for arbitrary m c has existed prior to our present calculation. For z < 1 4 , we can describe our findings by the following fit:  5 In the charm loop case (the right plot in Fig. 4), we had to rely on our so-far unpublished results for the UV counterterms -see the Appendix.
It is shown as a solid blue curve in the considered plot. A quick look at Fig. 5 is sufficient to realize that ∆ 4-b ✘ ✘ BLM m=0Ĝ (2) (2) 27 , due to the identity T a T b T a = − 1 6 T b for the SU(3) c generators. The same relative colour factor is valid for all the plots in Figs. 3 and 4. Fig. 4 shows contributions toĜ (2) 27 from diagrams with closed loops of quarks with masses m b (left) and m c (right). Only the two-body cuts are included. The solid blue lines correspond to the numerical fits from Eqs. (3.3) and (3.4) of Ref. [26]. In these cases, no four-body cuts are allowed, as the state X p s in Eq. (1.5) is assumed to contain no charm quarks. We do not consider three-body cuts here, as their effect can be included by multiplying the well-known three-body contribution toĜ 27 . As evident from the plots, our results are in perfect agreement with all the previously available expansions and fits. It is particularly important in the massive case (Fig. 4) where our verification comes as the first one from an independent group. The massless results from the upper two plots of Fig. 3 have already been cross-checked before.
As far as the new contribution (the third plot in Fig. 3) is concerned, it has so far been included in the interpolated part of the NNLO correction, and resulted in a tiny effect, around one per-mille of the decay rate only. Now we remove it from the interpolated part and replace by the fit in Eq. (2.9). It turns out that the interpolation estimate was correct within ∼10% of the considered contribution, so the effect remains tiny.

Updated SM predictions for B sγ and R γ
In the present section, we work out updated SM predictions for B sγ , as well as for the ratio R γ ≡ B (s+d)γ /B cℓν , where B cℓν is the CP-and isospin-averaged branching ratio of the inclusive semileptonic decay. Our main motivation for performing an update right now is not due to the NNLO corrections evaluated in the previous section. The new contribution is tiny, while the sizeable ones (that we have confirmed) were already included in the phenomenological analysis of Ref. [13]. However, there has been an important progress in estimating non-perturbative effects (see below). An update of the SM prediction should thus be performed right now, even though the m c -interpolation uncertainty remains essentially unchanged.
The first improvement in estimating the non-perturbative effects becomes possible thanks to the new Belle measurement of the isospin asymmetry In the SM, the dominant contribution to this asymmetry arises from a process where no hard photon but rather a hard 7 gluon is emitted in the b-quark decay [55]. Next, the gluon scatters on the valence quark, which results in emission of a hard photon. Instead of the valence quark, 6 Z OS G stands for the on-shell renormalization constant of the gluon wave-function, while Z g renormalizes the QCD gauge coupling in the MS scheme. 7 with momentum of order m b but possibly much smaller virtuality also a sea quark (u, d or s) can participate in such a Compton-like scattering. Taking this fact into account, one can write the decay rates as where Q u,d,s denote electric charges of the quarks participating in the Compton-like scattering, while the quantities A, . . . , D are given by interferences of various quantum amplitudes whose explicit form is inessential here. Since the considered effect gives only a small correction to the decay rate (B, C, D ≪ A), quadratic terms in Q u,d,s have been neglected above. We have also neglected isospin violation in the quark masses (m u = m d ) and in the electromagnetic corrections to theB-meson wave functions (suppressed by extra powers of α em ). The leading term A contains the dominant contribution originating from the operator Q 7 . The corrections B, C, D are suppressed w.r.t. A both by g 2 s (as the gluon is hard) and by Λ/m b , with Λ ∼ Λ QCD . The latter suppression can be intuitively understood by realizing that the gluon scatters on remnants of theB meson, i.e. on a diluted target whose size scales like 1/Λ. Such a suppression is confirmed in Refs. [55,56] where the Soft-Collinear Effective Theory (SCET) has been applied to analyze non-perturbative corrections to B sγ . From Eq. (3.11), one easily obtains the isospin-averaged decay rate 12) and the isospin asymmetry It follows that the relative correction to the isospin-averaged decay rate that arises due to the considered effect reads where, in the last step, Q s = −Q u − Q d has been used. The second term in the square bracket vanishes in the SU(3) F limit, i.e. when the three lightest quarks are treated as mass-degenerate. In this limit, as observed in Ref. [57], δΓ c /Γ and ∆ 0− are related to each other in a simple manner that is free from non-perturbative uncertainties. The authors of Ref. [56] suggested ±30% as an uncertainty estimate stemming from the SU(3) F -violating effect in Eq. (3.14). Following this suggestion, we find (3.15) where the experimental errors from Eq. (3.10) were combined in quadrature, giving ∆ 0− = (−0.48 ± 2.12)%; next, the multiplicative factor was taken into account as follows [58]: (1 ± 0.3)(−0.48 ± 2.12)% = −0.48 ± 2.12 2 + (0.3 · 0.48) 2 + (0.3 · 2.12) 2 %. (3.16) In the above estimates, we have treated the measured ∆ 0− in Eq. (3.10) as already extrapolated from the experimental cutoff of E 0 = 1.9 GeV down to our default E 0 = 1.6 GeV, even though no such extrapolation has actually been done in Ref. [54]. However, we expect the extrapolation uncertainty in this case to be negligible w.r.t. the remaining ones in Eq. (3.10).
If the uncertainty on the r.h.s. of Eq. (3.15) is treated as 1σ of a Gaussian distribution, then the 95% C.L. range is [−1.3, +1.6]%. The corresponding 8 range [−1.4, +2.0]% in Sec. 3.5 of Ref. [28] is somewhat wider due to a different method of combining uncertainties and using the PDG [59] central value of −0.6% for ∆ 0− . When determining our SM predictions below, we calculate B sγ without including the photon emission from the valence/sea quarks and, in the final step, we multiply by 1 + δΓc Γ , employing the number from the r.h.s. of Eq. (3.15). Another important non-perturbative correction to be considered arises in the interference of Q 1,2 and Q 7 . Its presence in the inclusiveB → X s γ rate was first pointed out in Ref. [60]. It amounts to around +3% of B sγ , as established in Ref. [61] at the leading order of an expansion in powers of m b Λ/m 2 c . The corresponding leading contribution to N(E 0 ) in Eq. (1.1) reads where µ 2 G ≃ 0.3 GeV 2 is one of the HQET parameters that matter in the determination of C in Eq. (1.2). Since m b Λ/m 2 c is not a small parameter, the authors of Ref. [56] argued that no expansion in its powers can be used at all. Instead, they estimated the considered correction in the framework of SCET, where essential constraints on models of the relevant soft function came from moments of the semileptonicB → X c ℓν decay spectra. A recent update of these estimates in Ref. [28] implies that δN V (3.17) needs to be multiplied by (3.18) The final numerical value above has been derived by us from ranges for Λ 17 given in Ref. [28], assuming that these ranges can be interpreted as 1σ ones. The remaining parameters on which κ V depends were set to the values corresponding to the widest range for Λ 17 in Ref. [28]. Since the expression for δN V (3.17) is calculated at the leading order in QCD only, the renormalization scheme for m 2 c in the denominator is unspecified. We assume that the corresponding uncertainty is included in the overall ±3% higher-order one that is being retained the same as in Ref. [13]. As the total effect from δN V amounts to around 3% in B sγ , uncertainties due to scheme-dependence of m c in δN V can safely be treated this way. In our numerical calculations, the quark masses and HQET parameters are included with a full correlation matrix (see Appendix D of Ref. [13]), except for the very m c in δN V that is now fixed to 1.17 GeV. The parameter κ V (3.18) will be treated as uncorrelated.
Apart from the two effects we have discussed above, the authors of Ref. [56] identified a third source of uncertain contributions to N(E 0 ) that arise at the order O(Λ/m b ). They come proportional to |C 8 (µ b )| 2 , where C 8 is the Wilson coefficient of the gluonic dipole operator Previous estimates of these corrections in Refs. [62,63] focused on large collinear logarithms ln m b ms that are present in the corresponding contributions to P (E 0 ). In Ref. [13], such logarithms were varied in the range [ln 10, ln 50] ≃ ln m B m K , ln m B mπ , which served as a crude estimate of the very uncertain but otherwise small contributions to B sγ where light hadron masses are the physical collinear regulators. However, according to Ref. [56], possible non-perturbative effects that come multiplied by |C 8 (µ b )| 2 can be unrelated to collinear logarithms, and affect B sγ by relative corrections in the range [−0.3, 1.9]% with respect to the m b ms = 50 case, for µ b = 1.5 GeV and E 0 = 1.6 GeV. Numerically, we can reproduce this range by performing a replacement ln m b m s → κ 88 ln 50 with κ 88 = 1.7 ± 1.1 (3.20) in all the perturbative contributions proportional to |C 8 (µ b )| 2 .
In the following, we shall treat the quantities δΓc Γ (3.15), κ V (3.18) and κ 88 (3.20) on equal footing with all the other parameters that B sγ depends on. Since they account for all the non-perturbative effects estimated in Refs. [28,56], we shall no longer include the overall ±5% non-perturbative uncertainty that entered the analysis of Ref. [13] as an input from Ref. [56]. This way we determine our updated SM predictions for B sγ and R γ in the SM, namely for E 0 = 1.6 GeV. The overall uncertainties have been obtained by combining in quadrature the ones stemming from higher-order effects (±3%), interpolation in m c (±3%), as well as the parametric uncertainty where all the non-perturbative ones are now contained. Not only δΓc Γ , κ V and κ 88 but several other inputs parameterize non-perturbative effects, too, namely the collinear regulators (see above), as well as the HQET parameters that enter either directly or via the semileptonic phase-space factor C (1.2). In the B sγ case, our parametric uncertainty amounts to ±2.5% at present. All the input parameters listed in Appendix D of Ref. [13] have been retained unchanged.
The overall uncertainty in R γ (3.21) amounts to ±4.8%, noticeably improved w.r.t. to ±6.7% in Ref. [12]. The main reason for the improvement comes from Ref. [28] where the dominant non-perturbative uncertainty, stemming from Λ 17 in Eq. (3.18), was reduced by a factor of three. Further improvement requires removing the m c -interpolation, and re-considering the higher-order and parametric uncertainties. If they remain unchanged, the expected future accuracy in the SM prediction for B sγ amounts to √ 3 2 + 2.5 2 % ≃ 3.9%, still somewhat behind the experimental expectation of ±2.6% that was mentioned above Eq. (1.1).

Summary
We reported on our calculation of the NNLO QCD corrections to B sγ without interpolation in m c , and presented final results for contributions originating from propagator diagrams with closed fermion loops on the gluon lines. They correspond either to the two-body (sγ) or fourbody (sqqγ) final states. In all the previously investigated cases, we confirmed the results from the literature, some of which had been obtained by a single group only. The new part comes from four diagrams with four-particle cuts that had not been determined before, as they are not included in the BLM approximation. Their contribution turns out to be tiny (∼ 0.1% of the decay rate) and quite well reproduced by our former interpolation algorithm.
In view of the recent progress in estimating the non-perturbative contributions, we performed an update of the phenomenological analysis within the SM. The obtained results yield B sγ = (3.40±0.17)×10 −4 and R γ ≡ B (s+d)γ /B cℓν = (3.35±0.16)×10 −3 for E 0 = 1.6 GeV. The main improvement in the uncertainty came from the analysis in Ref. [28] where nonperturbative effects in the Q 1,2 -Q 7 interference were re-analyzed including up-to-date constraints from moments of the semileptonic decay spectra. In effect, the corresponding uncertainty was reduced by a factor of three.
The next contribution to suppressing the overall theoretical uncertainty is expected from the calculation ofĜ (2) 17 andĜ (2) 27 for E 0 = 0 and at the physical value of m c , thereby removing the need for m c -interpolation in these quantities. where S 2 = 4 9 √ 3 Im Li 2 e iπ/3 . Finally, for the closed charm loops (the right plot in Fig. 4), the large-z expansion reads Ref. [26]. Analytical expressions for the leading terms agree with the findings of Ref. [23].
Determining the renormalized results plotted in Figs. 3 and 4 required taking into account three-loop counterterm diagrams with vertices proportional to Q 4 = (s L γ µ T a b L ) q (qγ µ T a q), namelyĜ (1)bare 47 . An expression for this quantity in Eq. (2.4) of Ref. [13] contains no contributions from closed loops of charm quarks, as all the other results in Sec. 2 of that paper. Such contributions arise in the two-body channel only. They take the form Small-z expansion of the function b(z) has been given in Eq. (3.9) of Ref. [64], while the large-z expansion of Re b(z) can be found Eq. (5.2) of Ref. [22]. As far asb(z) is concerned, we obtain the following expansions: No explicit expressions for the expansions ofb(z) have so far been published, even though this function must have been used for UV renormalization in Ref. [26].