Four-body contributions to B ->Xs gamma at NLO

Ongoing efforts to reduce the perturbative uncertainty in the B ->Xs gamma decay rate have resulted in a theory estimate to NNLO in QCD. However, a few contributions from multi-parton final states which are formally NLO are still unknown. These are parametrically small and included in the estimated error from higher order corrections, but must be computed if one is to claim complete knowledge of the B ->Xs gamma rate to NLO. A major part of these unknown pieces are four-body contributions corresponding to the partonic process b ->s qbar q gamma. We compute these NLO four-body contributions to B ->Xs gamma, and confirm the corresponding tree-level leading-order results. While the NLO contributions arise from tree-level and one-loop Feynman diagrams, the four-body phase-space integrations make the computation non-trivial. The decay rate contains collinear logarithms arising from the mass regularization of collinear divergences. We perform an exhaustive numerical analysis, and find that these contributions are positive and amount to no more than 1% of the total rate in the Standard Model, thus confirming previous estimates of the perturbative uncertainty.


Introduction
The inclusive radiative B meson decayB → X s γ is the paradigm for precision tests of the Standard Model (SM) in quark flavor physics. Its branching ratio is measured with a precision of ∼ 7% [1][2][3][4]  To match this experimental precision, a theory calculation to next-to-next-to-leading order (NNLO) accuracy is necessary. This program has been carried out during the last two decades and it is almost -but not quite -finished. The current theory estimate results in the following value [6]: SM Eγ >1.6 GeV = (3.15 ± 0.23) · 10 −4 , (1.2) 1 The semi-inclusive measurement in the first reference in [2] has recently been superseded by a new, more precise one -see Ref. [5]. However, this new measurement is not yet taken into account in the HFAG average of Eq. (1.1).
TheB → X s γ rate is given by the inclusive partonic rate of the b quark, up to nonperturbative effects that, for a photon energy cut E 0 = 1.6 GeV, are estimated at the level of ∼ 5% [48], where b → X parton s γ denotes the partonic decay of the b quark into any number of particles with an overall strangeness S = −1, plus a hard photon with E γ > E 0 , and excluding charm: with q = u, d, s. Each individual rate is an interference of different amplitudes, computed as matrix elements of local operators in the effective Lagrangian: (1.5) where the operators P i are defined as [19]: P q 1 = (s L γ µ T a q L )(q L γ µ T a b L ) , P q 2 = (s L γ µ q L )(q L γ µ b L ) , P 3 = (s L γ µ b L ) q (qγ µ q) , P 4 = (s L γ µ T a b L ) q (qγ µ T a q) , P 5 = (s L γ µ γ ν γ ρ b L ) q (qγ µ γ ν γ ρ q) , P 6 = (s L γ µ γ ν γ ρ T a b L ) q (qγ µ γ ν γ ρ T a q) , (1.6) With this notation, C q 1,2 contain CKM phases: C q 1,2 = −λ q C 1,2 , with λ q ≡ V * qs V qb /V * ts V tb and C 1,2 defined in the usual way, e.g. Ref. [19]. We will also use the notation C 1u ≡ C u 1 etc. The other Wilson coefficients are simply C 3,..,8 = C 3,..,8 as in Ref. [19].
The objects G ij arise from the interference of operators P i and P j in the squared matrix elements, integrated over phase space. They depend on the photon energy cut through the variable δ ≡ 1 − 2E 0 /m b . In the notation of Ref. [38], where normalization to the inclusive semileptonic b → u rate is used, K ij = G ij /G u , with G u = 1+C F ( 25 2 −12ζ 2 ) αs 4π +O(α 2 s ) [49]. In this paper we focus on the four-body contributions to Γ(b → X parton s γ) Eγ >E 0 , corresponding to Γ(b → sqqγ): where we define G ij as the b → sqqγ contribution to G ij . In addition, we expand G ij as: There is a hierarchy in the size of the Wilson coefficients of the operators in Eq. (1.6), which can be divided into two classes: 1 , P c 2 , P 7 , P 8 } ; B = {P u 1 , P u 2 , P 3 , P 4 , P 5 , P 6 } . (1.11) Operators in class A have large Wilson coefficients, whereas class-B operators have either very small Wilson coefficients or are CKM-suppressed. Among the four-body leading and next-to leading contributions we distinguish four groups: • Tree-level (B, B) interference ( Fig. 1.a). These are the leading-order (LO) contributions and provide the complete matrix G (0) (δ). This matrix has been computed in Ref. [45].
ij , with i, j = 1u, 2u, 3, .., 6. The ones in panel (f) include five-particle cuts since the one-loop four-body diagrams must be complemented with the five-body gluonbremsstrahlung correction b → sqqγ + g. We therefore leave the contributions from panel (f) aside from the present four-body calculation.
We note that NLO (A, B) interference terms are presumably as large as the (B, B) interference at the LO since C 1u,2u,3,..,6 ∼ α s /π C 1c,2c,7,8 . For the same reason, the partly neglected (B, B) interference terms at the NLO are expected to be much smaller than the (A, B) interferences that we calculate in a complete manner. As a final comment, we note that four-body NNLO contributions of the type b → sggγ are part of the calculation in Ref. [46], and do not interfere with our results.
The structure of the paper is the following. In Section 2 we discuss the details of our calculation and the structure of the different contributions, including the set of diagrams needed, the UV renormalization, and the treatment of collinear divergences. In Section 3 we collect the final results. In Section 4 we perform a numerical study of all the evaluated interferences. Section 5 contains our conclusions. In Appendix A we collect some intermediate results of the calculation, where analytic cancellation of UV and collinear divergences can be explicitly checked.
3. Renormalization: This requires the evaluation of the tree-level diagrams with counterterms shown in Fig. 5, and the corresponding phase-space integration. As a bonus, this step allows one to check the LO result for G (0) ij of Ref. [45]. This step is described in Section 2.5.
4. Collinear logarithms: Having regularized collinear divergences in dimensional regularization, we use the method of splitting functions [45,[50][51][52][53] to transform 1/ coll poles into collinear logarithms of the form log(m q /m b ). This requires the computation of the corresponding b → sqq corrections with subsequent photon emission q → q γ (with q = q, s), by evaluation of the diagrams shown in Fig. 6, the convolution with the splitting function, and the three-particle phase-space integration. This step is described in Section 2.6.
We note that every diagram has to be computed inserting all the operators P 1u,2u,3,..,6 to the right of the cut, and P 1u,2u,1c,2c,3,.., 8 to the left (see e.g. Fig. 1), leading to all the different interference terms in the matrix G (1) ij . In Section 2.1 we argue that most of the elements of this matrix can be obtained from a reduced set by the use of different operator relations and Fierz identities. In addition, this spells out transparently the dependence of the full result on the Wilson coefficients before any calculation is performed. We will see that -with one exception discussed at the end of Section 2.1.3 -only diagrams with P 7,8,1c to the left of the cut and P 4 to the right must be considered. This simplifies the calculation considerably.
Finally, for each diagram in Figs. 2−6, there is the corresponding mirror image. These "mirror" contributions are just obtained by complex conjugation, and ensure that the rate is real, i.e. G (1) ij is hermitian. We disregard these mirror contributions in the calculation, but at the end we substitute G

Color
Diagrams with insertion of the color octet operators P 4,6 are related to the ones with insertion of color singlet operators P 3,5 by a simple color factor, which can be obtained just from the color structure of the gluon penguin: This leads to the rule that P 3,5 can always be replaced by: For the same reason, one can always put P q 2 → −6P q 1 , meaning that C q 1,2 always come in the combination (C q 1 − 6 C q 2 ).

Insertions to the right of the cut
We restrict ourselves here to the insertion of operators to the right of the cut. Using the 4D identity γ µ γ ν γ λ = g µν γ λ + g νλ γ µ − g µλ γ ν + i µνλα γ α γ 5 we find: . We now consider the following "crossed" insertion of P 4 into a massless fermion loop with an arbitrary number of vector currents: There is always an even number of gamma matrices to the left of γ 5 , which can be moved besides the projector P L . After that, the relationship P L γ 5 = −P L provides the given negative sign.
The "straight" insertion of P 4 does not vanish in general but does not contribute in our case: In the situation where one vector current is attached to the loop (Fig. 2), the result is proportional to Tr[γ µ γ ν γ ρ γ λ γ 5 ] ∼ µνρλ , but there are not enough independent momenta to be contracted with the antisymmetric tensor, so this contribution vanishes. This is true also in the case where the photon couples twice to the quark loop (Fig. 4). In the case of two current insertions (Fig. 3) the result is non-zero, but summing over u, d, s quarks in the loop the result is proportional to Q u + Q d + Q s = 0.
Summing up, in the diagrams with a P 6 insertion one can always substitute: The replacement rules (2.7) combined with Eqs. (2.3) and (2.4) imply that the full contribution from P 3,4,5,6 can be obtained from the terms proportional to C * 4 : Result(P 3,4,5,6 ,straight) = (C * 4 + 10 C * 6 ) × Result(P 4 ,straight), (2.8) Result(P 3,4,5,6 ,crossed) = (−6 C * 3 + C * 4 − 96 C * 5 + 16 C * 6 ) × Result(P 4 ,crossed). (2.9) Since this is based on a 4D identity, it is in principle only true up to evanescent terms. Below we show that up to the needed order in these terms do not contribute. We have also checked this by direct computation. The (crossed) insertion of the operators P u 1,2 can also be obtained from P 4 by an argument almost identical to that of Eq. (2.6). In this case one must pay attention to the case where the photon couples to the loop, where the P 4 and P u 1,2 contributions are proportional to different charge factors.

Insertions to the left of the cut
We have shown that we only need to compute diagrams with an insertion of P 4 to the right of the cut. To the left of the cut we must insert P 7,8 as well as P 3,4,5,6 and P u,c 1,2 . As before, P 2,3,5 contributions are related to P 1,4,6 by a simple color factor. In addition, the contribution from P u 1 is obtained from the P c 1 insertion with the replacement m c → 0. We now show that insertions of P 4,6 are also known from the insertions of P 7 , P 8 and P c 1 . First we consider the case of the photon penguin, where the gluon does not couple to the fermion loop to the left of the cut. By direct inspection be find that: where C eff 7 = C 7 − C 3 /3 − 4 C 4 /3 − 20 C 5 /3 − 80 C 6 /9 is the usual "effective" Wilson coefficient [20], which includes the contributions from b-quark loops. The O( ) corrections are irrelevant to our calculation as the contributions from P 7 are finite. The term Xq /q µ denotes the contribution from four-quark operators proportional to the structure [q /q µ − q 2 γ µ ]. This last term does not contribute in our case. To see this, consider the insertion of P 1,2 into the full diagram: Here the gluon is attached to either '×'. Since we cut the photon propagator, the photon is on-shell but there is no q 2 denominator, and therefore the q 2 term cancels. The term q /q µ also cancels due to the Ward identity q µ Γ µν = 0. Note that non-zero contributions to the Ward identity vanish since they either involve a massless fermion tadpole, or if the gluon couples to the loop then it does not depend on the incoming/outgoing quark momenta. We have also checked this result by explicit computation, and indeed the different sets of diagrams satisfying the Ward identity vanish identically.
To summarize: all contributions from photon penguins to the left of the cut are obtained from the diagrams with insertion of P 7 by the replacement C 7 → C eff 7 . Next, we consider the case of the gluon penguin, where the photon does not couple to the fermion loop to the left of the cut. We find: where C eff 8 = C 8 + C 3 − C 4 /6 + 20 C 5 − 10 C 6 /3, as usual (e.g. Ref. [32]). As before, O( ) corrections to the chromomagnetic contribution are irrelevant for our calculation because contributions from P 8 are UV-finite (collinear divergences are inconsequential here). In the last term, the quantity X is given by: where n = 3 is the number of light flavors and is the loop integral to all orders in . Therefore, all contributions from gluon penguins to the left of the cut are known from the contribution of P 8 and P c 1 . Now we consider the case in which both the photon and the gluon couple to the loop: When inserted into the full diagrams, these contributions are both UV-finite and collinear safe 2 , and we can use 4D identities. The first term in the right-hand side of Eq. (2.15) can be written as Q u (C c 1 − 6C c 2 )h µν (m c ), which defines the function h µν (m). The second term can be obtained from the first term by the replacement m c → 0: . The third term (q = u, d, s, c, b) contains only the insertion of P 6 , because P 3,5 insertions are zero due to color, while the insertion of P 4 vanishes due to Furry's theorem. For P 6 we can make use of the Fierz identity 3 : which implies that we can trade the straight insertion of P q 6 with the crossed insertion of −36P q 1 . Note also that the contributions from q = u, d, s add up to zero in the massless limit due to electric charge: Q u +Q d +Q s = 0. This means that the third term in Eq. (2.15) is given by −36 The fourth term can be obtained from the first one using the identities in Eqs. (2.5) and (2.6), leading to: Q d (−6 C 3 + C 4 − 96 C 5 + 16 C 6 )h µν (0). The fifth term cannot be completely determined from the insertion of P 1 due to the chirality structure. Using the Fierz identity (2.5) we can trade P b The second piece, together with the corresponding contribution from P 5 , results in Q d (−72 Altogether, the right-hand side of Eq. (2.15) can be written as: Therefore only diagrams with insertions of the operators P 7,8 and P c 1 to the left of the cut need to be calculated, plus the extra contribution from h µν (m b ). All these relationships shape the structure of the full results displayed below in the following sections.

Set of Diagrams
There are three types of diagrams: 3 Here we use the notation P u 3 = (sLγµbL)(ūγ µ u), etc. Type (i): The photon does not couple to the cut fermion loop ( Fig. 2): In this case crossed and straight diagrams contribute. In addition, straight diagrams contain a factor n . All in all, the contribution from these diagrams is: (J) and P 4 ×,k (J) denote the contributions to (P k , P 4 ) interference terms from straight and crossed insertions of the operator P 4 to the right of the cut, to diagrams of type (J), respectively.

Type (ii):
The photon couples to the cut fermion loop once ( Fig. 3): In this case the straight diagrams are proportional to Q u + Q d + Q s = 0 and need not be considered. The P u 1,2 contributions are proportional to Q u . Therefore the total contribution from these diagrams is: The photon couples to the cut fermion loop twice ( Fig. 4): In this case crossed diagrams are proportional to Q 2 s = Q 2 d (or Q 2 u in the case of P u 1,2 ) and straight diagrams to We assume that the objects P 4 I,k (J) are phase-space-integrated matrix elements containing no prefactors or couplings, such that: i=3,..,8,1q,2q j=3,..,6,1u,2u in the notation of Eq. (1.9). In Eq. (2.21), D C (J) denotes the UV counterterms, and both D C (J) , D k (J) include the relevant collinear counterterms. Both are discussed below in Sections 2.5 and 2.6. The structure of the objects P 4 can be deduced from the discussion in Section 2.1.3. In the case of P 7,8 we have: The contributions with penguin operators to the left of the cut are given by: The functionsF I,k (J) = F I,k (J) + F I,k coll(J) include the collinear regulators discussed in Section 2.6. The functions F I,1 (J) and F I,4 (J) are related to diagrams where the photon couples to the left-hand quark loop, corresponding respectively to the terms with h µν and h µν in Eq. (2.17). Explicit results for all these functions are collected in Appendix A.

Other details 2.3.1 Irrelevance of evanescent terms to the right of the cut
In the case of (P 7 , P i ) interference, there are no UV or collinear divergences, and therefore evanescent structures are irrelevant for the O( 0 ) result.
In the case of (P 8 , P i ) interference, collinear divergences appear which combined with evanescent terms give finite contributions in the dimensionally regularized result. However these finite terms cancel when we express the dimensional regulators in terms of logarithms of masses, via the splitting-function approach (see Section 2.6): since the 1/ terms cancel in both 4D and evanescent terms separately.
In the case of (P 1,2 , P i ) interference, UV and collinear divergences are nested inside dimensionally regularized expressions. However all UV divergences cancel against coun-terterm diagrams, including finite terms from evanescent operators: All "finite" terms from collinear divergences now disappear when going to mass regularization, as in the case with P 8 .

Cancellation of
terms Traces with γ 5 will introduce terms proportional to i µνρσ k µ 1 k ν 2 k ρ 3 k σ 4 in the differential decay rate. Here we show that these terms always cancel if we perform a full angular integration over phase-space.
Consider fixing all double invariants k i · k j . Then all k i are fixed only up to an Euler rotation and an orientation. To see this go to the rest frame of the b-quark. Momentum conservation fixes all the energies (since k i · p b are fixed). This implies that k i · k j are also fixed. We can rotate the frame to put k 1 along the positive z axis, and k 2 in the (y, z) plane. Then k 3 is fixed only up to a two-fold ambiguity (an orientation), given by the sign of its x component. Once this sign is chosen k 4 is also fixed. This proves that i µνρσ k µ 1 k ν 2 k ρ 3 k σ 4 is fixed by k i · k j up to a sign, which is given by the orientation of ( k 1 , k 2 , k 3 ). Now consider phase-space integration. Terms in the integrand of the form F (k i · k j ) do not depend on the Euler rotation nor the orientation, and the angular integral over dΩ 3 dΩ 2 dΩ 1 can always be performed trivially, giving a factor 16π 2 . Terms of the form , however, change sign under change of orientation, and vanish upon integration over dΩ 1 . Obviously parity-odd terms cancel out in parity-even observables. Therefore we drop these terms from the beginning in the calculation of the integrated decay rate.

Phase-space integration
The phase-space measure for a (1 → 4) decay of a particle of mass M into four massless particles with momenta k 1,2,3,4 is given in terms of kinematic invariants by [54]: The unpolarized decay rate is given by the phase-space integral: where the sum runs over the spins and color of all particles (we assume the parent is a color triplet). |M| 2 depends only on s ij : |M| 2 ≡ K(s ij ), so the angular integrations can be performed trivially: , (2.31) and the general formula for the decay rate becomes This integral might contain soft and/or collinear divergences associated to regions of phase space where some particles are soft or collinear. These divergences can be regularized in dimensional regularization by setting d = 4 − 2 . If we insist on integrating over these regions, one must include virtual corrections to cancel the divergences. Otherwise, the regulator must be traded by a physical cutoff at a later stage.
We now specify to the b → q(k 1 )q(k 2 )s(k 3 )γ(k 4 ) case. We consider a cut on the photon , which defines the parameter δ. This translates into the constraint s 14 +s 24 +s 34 > 1−δ, which can be included in the phasespace integral in the following way. We include a delta function δ(1 − z − s 14 − s 24 − s 34 ) in the integrand, and we integrate over the new variable z from 0 to δ: The delta functions can be used to integrate over two invariants, e.g. s 13 and s 24 : , where a ± are complicated functions of the rest of the invariants: Now it is convenient to perform the following changes of variables: where u, v, w, x are integrated independently from 0 to 1, and This gives In the following we must consider the kernel K(u, v, w, x, z). As mentioned above, K is polynomial in s 23 : K = f n (v, w, x, z)s n 23 . Expanding s 23 according to (2.38)-(2.40) will provide a sum of terms of the form 42) and the integral over the variable u gives a factor β d−3 for each term. The next steps depend on the diagram at hand. Consider the diagrams with P 7,8 . In this case, for some a, b, c, ... ∈ R. The integral over x gives again a β-function: β d−2 2 + f, d−2 2 + g . Because of the (1 −zv) q and (1 − zw) r factors, the next steps will introduce hypergeometric functions. The integral over v gives 44) and the integral over w, gives: The next step is to expand around → 0 (with d = 4 − 2 ). The expansion of hypergeometric functions is performed automatically by the package HypExp [56]. This will give finite results in the case of P 7 , but 1/ poles in the case of P 8 , corresponding to collinear divergences. The integration over the photon energy z ∈ (0, δ) can then be performed, also analytically, for all terms. The case of loop diagrams is in principle more complicated, as the function K contains already a hypergeometric function. For instance, in the case of diagrams (i) with the photon not attaching to the quark loop, the variable s 12 appears in the function g(m q ) ∼ 2 F 1 ( , 2; 5 2 ; s 12 4zq ) (cf. Eq. (2.14)). However, by a suitable choice in the order of integration, analytic results can be obtained as before. In the case of diagrams such as (iii), the hypergeometric function depends on the triple invariant s 124 = s 12 + s 14 + s 24 , and the sequential-integration procedure described above does not seem to work up to finite order in . In this case we extract all the 1/ 2 and 1/ poles analytically and leave the finite terms differential in one of the variables, which we integrate numerically afterwards. This is also the case for the diagrams where the photon couples to the charm loop, which are both UV and collinear finite. In general, for the loop contributions, some finite terms turn out to be complicated functions of δ and z c ≡ m 2 c /m 2 b . We give these results as polynomial expansions in δ around the physical value δ = 0.316. The coefficients of this expansion are presented as numerical interpolations in the variable z c , reproducing the exact results to enough precision for all practical purposes. We have checked that the interpolated expressions in the appendix reproduce the exact results with high precision in the full range z c ∈ (0, 1) for values of δ near 0.316.

Renormalization
Tree-level four-body contributions from four-quark operators arise at LO in α s and have been computed in Ref. [45]. At NLO the corresponding counterterm contributions must be included, which cancel the UV divergences from the loop diagrams. One must consider the insertion of the bare operators P ,..,6,1u,2u i=3,..,6,1q,2q j= 3,..,6,1u,2u The first term leads to the LO contributions in Ref. [45], while the second term contributes to the NLO result and takes care of the UV divergences. For this we need, a priori, the tree level results with P 3,..,6,1u,2u including O( ) terms, and the renormalization factors δZ ij . The relevant renormalization factors are simple to compute. Using the relationships developed in Section 2.1, and expressing the result in terms of tree-level matrix elements of four-quark operators, we find that: (2.48) We also see that we need only tree level diagrams with insertion of P 4 to the left of the cut. All the diagrams needed are shown in Fig. 5.
For the operator insertions to the right of the cut we can (and must) use the 4D identities derived in Section 2.1, noting that evanescent terms cancel in the renormalization process by virtue of Eq. (2.27). This leads to exactly the same structure as Eqs.

Collinear divergences and splitting functions
The region of phase space in which the photon is collinear to one of the light quarks gives rise to collinear divergences. These divergences are regulated dimensionally in our computation. However, these are just artifacts of the massless limit used for light quarks, and there is a more natural regulator: a physical cut-off given by the light meson masses. A suitable parametrization of such (near-) collinear effects consists in keeping the light quarks massive and perform a massive phase-space integration. This is quite complicated from the practical point of view, taking into account that the massless phase-space integrals computed here are already rather challenging. Fortunately, one may resort to the factorization properties of the amplitudes in the quasi-collinear limit (see e.g. [45]). The idea is that close to the collinear region, the b → q 1q2 q 3 γ amplitude may be expressed as a b → q 1q2 q 3 amplitude times a splitting function f i , describing the quasi-collinear emission of a photon from q i , summed over i = 1, 2, 3. The splitting functions encode the collinear divergences, and can themselves be regulated by quark masses or in dimensional regularization. Both approaches are rather simple, since in this limit the four-body phase space factorizes into a convolution of the three-body phase space of the non-radiative process and the phase space of the radiative process alone: dΦ 4 = dΦ 3 ⊗ dΦ. By comparing the splitting functions regulated in these two different schemes, one can write a formula to switch from one to the other at the level of the decay rate [45]: (2.51) Here K 3 = |M 3 | 2 is the spin-summed squared matrix element of the b → q 1q2 q 3 decay obtained by evaluating the diagrams in Fig. 6, and dP S 3 is the three-particle phase-space measure in d = 4 − 2 dimensions [54]: (2.52) Integrating Eq. (2.51) over z ∈ [0, δ] provides the terms F coll contained in the functionŝ F. The contributions from the chromomagnetic operator (Fig. 6, left) enter into Eq. (2.22). The contributions from four-quark operators (Fig. 6, center) go into Eqs. (2.23), (2.24) and (2.25). The counterterm contributions (Fig. 6, right), enter into Eq. (2.49). The functions F coll (δ) are collected in Appendix A.
One can check that adding the collinear contribution removes the 1/ terms that survive the renormalization process, trading them for collinear logarithms of quark-mass ratios. These collinear logarithms are of the form log(m q /m b ), with q = u, d, s. The quark masses are collinear regulators and it is difficult to relate them to physical masses. In our numerical analysis we will take a common constituent-quark mass m q ∼ 100 − 250 MeV for all three light flavors, and use the notation L q = log(m q /m b ) ∼ log(m u /m b ) ∼ log(m d /m b ) ∼ log(m s /m b ). This should provide a reasonable estimate of the effect of collinear logarithms.

Results
We write the four-body contribution to theB → X s γ rate as: where Γ 0 is the absolute normalization of the decay rate: The sum runs over i, j = 1u, 2u, 3, .., 6, 1c, 2c, 7, 8. The Wilson coefficients C 3,..,8 are real, but C 1u,2u,1c,2c contain CKM phases: with C i the Wilson coefficients in the notation of Ref. [19]. They are needed here to NLO: their numerical values are given below. The matrix elements G ij (µ, δ) depend on the renormalization scale and the photonenergy cut and can be split into LO and NLO components: The LO matrix G (0) is real and symmetric and was computed in Ref. [45]: we reproduce and confirm these results (after the 2014 update of that paper). We write (here i, j = 1u, 2u, 3, .., 6): where: The NLO matrix G (1) contains perturbative strong phases from on-shell contributions from light quarks, as well as from charm quarks when the photon-energy cut is low enough. The matrix G (1) is the main result of this paper. It has the following structure: ij in the file "Gij.m" attached to the arXiv submission of the present manuscript. The first is given by the 6 × 6 matrix "GijLO" (i, j = 1u, 2u, 3, .., 6) and the second by the 10 × 10 matrix "GijNLO" (with i, j = 1u, 2u, 3, .., 6, 1c, 2c, 7, 8).

Numerical analysis
We briefly discuss here the numerical impact of the four-body contributions to the total B → X s γ rate. We consider for convenience the following quantity: given by Eq. (3.1) and normalized to the leading contribution to the decay rate. The Wilson coefficients are given by: which are computed following Ref. [12]. For the reference matching and renormalization scales µ 0 = 160 GeV, µ = µ b = 2.5 GeV, we have 4 : for i = 1q, 2q, 3, .., 8. However, the µ-dependence of the Wilson coefficients is important and we will analyze it here. In addition, λ q ≡ V * qs V qb /V * ts V tb denote the appropriate CKM factors, given by [57]: The quantity ∆Γ can be expanded in α s : We begin with a discussion of the µ-dependence of our results. To leading order, the µdependence is given purely by the LL (leading-log) running in the effective theory. Note that to this order, only C u 1,2 and C 3,4,5,6 contribute. At NLO, three new contributions arise: (i) the contribution from NLO Wilson coefficients, (ii) NLO matrix elements and (iii) contributions from C c 1,2 , C 7,8 , absent at LO. The µ-dependence should cancel up to a residual scale-dependence from higher orders, and up to the neglected contributions shown in Fig. 1.f (note that the Z factors in Eq. (2.48) are not the full renormalization constants).
In Fig. 7 we show the µ-dependence of the LO result, and LO+NLO excluding C c 1,2 , C 7,8 , LO+NLO excluding C c 1,2 and LO+ full NLO. We also gauge the impact of collinear logarithms, showing the result for two different choices of L q , corresponding to m q = m b /50 (m q ∼ 100 MeV) and m q = m b /20 (m q ∼ 250 MeV). Collinear logarithms are, as expected, numerically important.
Contributions from P c 1,2 and P 7,8 arise only at NLO and therefore introduce at this order a novel µ-dependence. Although, as we will see, certain cancellations make the NLO contribution small, there is a considerable reduction in the renormalization-scale dependence of the full LO+NLO result as compared to the LO contribution alone. This is due to the fact that the main µ-dependence of the leading order contribution arises from the mixing of P c 1,2 onto penguin operators, which is compensated at NLO by the matrix elements of P c 1,2 . This can be seen in Fig. 7: the reduction in the µ-dependence is achieved only after including C c 1,2 contributions. In the left plot of Fig. 7 one can see that for the value µ 4 GeV strong cancellations make the NLO contribution very small. More concretely, for the inputs µ 0 = 160 GeV, where the term labeled 'NLO WCs' corresponds to the second term in Eq. (4.6). This cancellation is very efficient for µ 3.8 GeV, but depends strongly on m q and z c . However, it is a general feature of our results that the contribution from C c 1,2 is of the same order as the rest of the NLO contribution, but with opposite sign, leading always to some level of cancellation. Note also that the (NLO) C c 1,2 contribution is also as large as the LO result. In the following we fix the renormalization scale to µ = 2.5 GeV and study the dependence on the charm mass and the photon-energy cut. This is shown in Fig. 8. In general the full LO+NLO result increases with δ and decreases with m c , always remaining below the 1% level for δ 0.45. We note that these results are only valid for δ not far from 0.316 as some of the functions are expanded up to second order in (δ − 0.316).
Finally, we provide some results for two different values of E 0 of interest: E 0 = 1.6 GeV, corresponding to δ = 0.316, and E 0 = 1.9 GeV, corresponding to δ = 0.188. For the input parameters and their uncertainties we take: µ 0 = 160 +90 −80 , µ = 2.5 +2.5 −0.5 and z c = 0.07 ± 0.02, which captures the different values for m c within different schemes. For the value δ = 0.188 (E 0 = 1.9 GeV) we have used the exact results (not the expanded ones), as for this value of δ the quadratic expansion is not expected to be accurate enough.

Conclusions
The inclusive radiative decayB → X s γ has beyond any doubt reached the era of precision physics, with the total uncertainties on both the experimental and theoretical side being at the ±7% level. The foreseen improvement in precision on the experimental side -the envisaged uncertainty with 50/ab at Belle II is of O(6%) [58], although this might even be a conservative estimate -justifies every effort to reduce the theoretical error to at least the same level.
The present article aims at addressing a particular higher-order perturbative contribution, namely the four-body contributions b → sqqγ to Γ(B → X s γ) at NLO. The smallness of the Wilson coefficients of penguin operators and CKM-suppression of current-current operators suggests that this contribution should be small. However, only an explicit calculation can turn this estimate into a firm statement. The calculation arises from tree and one-loop amplitudes, but it involves the four-body phase-space integration in dimensional regularization, which makes the calculation non-trivial owing to the appearance of higher transcendental functions. Moreover, the cancellation of poles in the dimensional regularization parameter is only achieved after proper UV and IR renormalization. The latter gives rise to logarithms ln(m q /m b ) when turning the dimensional into a mass regulator. These logarithms stem from the phase space region of energetic collinear photon radiation off light quarks in the final state. They are computed with the splitting function technique and treated in the same way as in [45,55].
We find indeed that the contribution of our four-body NLO correction to the total rate is below the 1% level, as expected. This statement even holds true once we vary the input parameters such as the charm mass m c , the photon energy cut (parameterized by δ), the masses m q of the light quarks, or the renormalization and matching scales, as can be seen by the numbers and the plots in Section 4. We also confirm the LO results presented previously in Ref. [45].
Yet the NLO calculation ofB → X s γ is still not complete. There are certain yet unpublished three-particle cuts contributing to Γ(b → sgγ), mainly interferences of P u 1,2 with P c 1,2 , which are also of the (A, B)-interference type. These contributions can be extracted from the results of Ref. [27]. The only missing pieces are given by the diagrams in Fig. 1.f. These are NLO interferences of the type (B, B) and are expected to be negligible with respect to the (A, B) ones that we have calculated in a complete manner. While these contributions can be calculated with the techniques described in this paper, they are left for future work.
Up to subleading terms in , we have always The functions F ×,8 (δ) are given by: where L µ = log(µ/m b ) and L q = log(m q /m b ), and: Li 2 (δ) ; (A.14) For F I,1 (J) (δ) we give analytical results for m c -independent functions, but the m c -dependence is given as interpolated functions. Up to subleading terms in , we have always The functions F × (z q , δ) are given by: Counterterms are given by: where again L µ = log(µ/m b ) and L q = log(m q /m b ). From these expressions one can check the pattern of cancellation of UV and collinear divergences. The z q -independent functions are given by (with the notationδ ≡ 1 − δ, Lδ ≡ log(1 − δ), L δ ≡ log δ), Lδ ; (A.28)  While our calculation provides exactly all the functions B ( ) (z q , δ), C ( , ) (z q , δ) and E ( ) (z q , δ), the corresponding expressions depend on z q and the photon energy E γ through complex functions of harmonic polylogarithms of various weights, which must be integrated in the region 2E γ /m B ∈ [1 − δ, 1]. Solving these integrals analytically is highly non-trivial, and even the numerical integration is computationally demanding. We have performed a numerical evaluation of such integrals and find it more convenient to present the results as an expansion in δ around the value δ = 0.316, and as an interpolation in z q . These interpolations coincide with the exact results in the region z ∈ [0, 1] to a very good precision. The relevant functions are written as:  Table 3. Padé coefficients for C(z q , δ).