Three-loop b → sγ vertex with current-current operators

We compute three-loop vertex corrections to b → sγ induced by current-current operators. The results are presented as expansions in mc/mb with numerical coefficients which allow to cover all relevant values for the heavy quark masses in different renormalization schemes. Moreover we provide for the first time analytic results for the next-to-leading order contribution. Our results present an important building block to the next-to-next-to-leading order interference contributions of the current-current operators Q1 and Q2 with the electric dipole operator Q7.


Introduction
The inclusive rare decay B → X s γ, where X s is any charmless hadronic state of strangeness S = −1, is an important probe to look for phenomena beyond the Standard Model.It is thus of prime importance to both measure it with highest accuracy, and to provide precise predictions.
The current combination of various experimental measurements [1][2][3][4][5][6] leads to the CPand isospin-averaged branching ratio [7] B(B → X s γ) Eγ >1.6 GeV where E 0 = 1.6 GeV is a lower cut on the energy of the photon.The current uncertainty of about 5% is expected to be reduced by a factor two once the full Belle-II data set is analysed [8].
The most precise theory prediction from Refs.[9,10] also has an uncertainty of 5% and is given by where the updates from Ref. [11] are taken into account.
The decay B → X s γ is well approximated by the decay of a free quark b → X s γ.Perturbative QCD corrections are calculated within the framework of an effective theory obtained after decoupling the top quark, the Z, W and Higgs bosons.The weak effective theory then contains four-quark and dipole-type operators (see e.g.Refs.[12,13] and Section 2).The prediction in Eq. ( 2) includes next-to-next-to-leading order (NNLO) QCD corrections.However, for the contribution arising from the interference between the four-quark operators (Q 1 and Q 2 ) and the electromagnetic dipole operator (Q 7 ) only an interpolation from a large charm quark mass [14,15] to a massless charm quark [10] is available.This is responsible for 3% of the uncertainty cited in Eq. (2).The remaining theoretical uncertainties come from unknown higher-order corrections (3%) and the input and nonperturbative parameters (2.5%).In this paper we provide important contributions that are necessary to eliminate the uncertainty due to the charm mass interpolation.
In order to calculate the interference of Q 1,2 and Q 7 , one often applies the method of reverse unitarity [16] to map the calculation of the interference into the evaluation of cuts of two-point functions.Such contributions occur for the first time at next-to-leading order (NLO) where three-loop diagrams, as the one on the left of Fig. 1, have to be considered.These diagrams have both two-and three-particle cuts as indicated by the dashed lines.Correspondingly, at next-to-next-to-leading order (NNLO) one has to consider four-loop diagrams which in general have two-, three-and four-particle cuts.
In this paper we concentrate on the two-particle cut contribution.It can be obtained by calculating QCD corrections to the b → sγ vertex and subsequently performing the integration over the two-particle phase space.Of course, in the physical decay rate all cuts going through the photon line have to be considered.In fact, individual contributions contain divergences which cancel only when all cuts are taken together and after including the counterterm from renormalization of the ultraviolet divergences.The latter is conveniently done for the complete contributions and available from Ref. [17] which is why we only provide results for the bare amplitude.
There are a number of NNLO results available for the interference of Q 1,2 and Q 7 .Analytic results for the light-fermion two-particle cut contributions are available from Ref. [18] in an expansion for m c /m b → 0. The corresponding four-particle cut terms are available from Ref. [19].The contributions with closed charm and bottom quark loops have been computed in Refs.[11,20].A numerical approach has been used to solve the differential equations for the master integrals with boundary conditions computed for large values of m c .The same approach has been used in Ref. [17] to extend the NLO results by one order in ϵ and to obtain all counterterm contributions.The non-fermionic contribution has been computed in the large-m c limit in Refs.[14,15] and for massless charm quarks in Ref. [10].These results are used in Refs.[9,10] to construct an interpolation which induces the 3% uncertainty mentioned above.
Recently, Ref. [21] has considered the subclass of diagrams for the b → sγ vertex where no bottom quark propagator is present in the loop.Here we provide an independent check of these results and also compute the more involved contributions with internal bottom quarks.
The paper is organized as follows.In Section 2 we introduce the notation and present the setup of our calculation.We then discuss in Section 3 the calculation of the master integrals at two and three loops.Results are presented in Section 4 and we conclude in Section 5.In Appendix A we provide further details on the analytic calculation of the two-loop master integrals and in Appendix B analytic expressions for the individual NLO contributions are given.

Technicalities
The relevant weak interaction Lagrangian for the calculation of radiative B meson decays is given by a linear combination of four-quark and dipole-type operators and can be written as where the three operators considered in this paper are given by Here, s, c and b are the fields of the strange, charm and bottom quarks, respectively.The subscripts L and R denote the projection to left-and right-handed states and T a = λ a /2, where λ a are the Gell-Mann matrices.F µν is the field strength tensor of the photon field, e is the electric charge, m b the bottom quark mass, and We can write the amplitude for b → sγ as where ε µ is the polarization vector of the photon.We parameterize M µ as where all momenta are incoming and on-shell.We construct projectors to extract the three scalar coefficients t 1 , t 2 and t 3 , which depend on the ratio m c /m b .The function t 1 does not contribute to b → sγ since the corresponding tensor structure vanishes once we contract M µ with the photon polarization vector in Eq. ( 5).Furthermore, from the Ward identity q µ γ M µ = 0 we have which we use as a cross-check for our results.Note that this relation holds for renormalized quantities, i.e., in particular for the two-loop amplitudes.Therefore all counterterm contributions obtained from multiplicative renormalization fulfil the Ward identity.However, the relation (7) does not hold for the NNLO contributions from the bottom quark mass counterterm since they are not obtained from a global renormalization factor.As a consequence, we do not expect that the bare three-loop amplitude fulfils Eq. ( 7) per se, but only in combination with the bottom mass counterterm contributions.
Since the scalar coefficients t 1 , t 2 and t 3 depend only on the ratio m c /m b and not on the external momenta, we can straightforwardly obtain the two-particle cut contribution to the decay rate: where the sum includes the operators in Eq. ( 4), the penguin operators and the chromomagnetic dipole operators.We write the perturbative expansion of the interference terms as In this paper we compute the two-particle cut contributions to Ĝ17 and Ĝ27 taking into account only tree-level contributions on the Q 7 side.The only other non-zero contribution with two-loop corrections on the Q 1,2 and one-loop corrections on the Q 7 side can be obtained from lower-order results.Sample diagrams for the three-loop corrections on the Q 1,2 side are given in Fig. 2. Ĝ17 and Ĝ27 are obtained by taking the interference between the b → sγ amplitude in Eq. ( 6) and the tree-level one mediated by the operator Q 7 .The subsequent integration over the d-dimensional two-particle phase space leads to (i=1,2) Ĝ2P,Q tree where γ E is the Euler's constant.The superscripts on the left-hand side indicate that this formula provides the two-particle cut contributions where no loop corrections are considered for Q 7 .The last factor comes from the d-dimensional two-particle phase space.Note that in case we use the relation t 3 = −t 2 /2, the first factor becomes proportional to (1 − ϵ) and we recover Eq. (3.5) of Ref. [17].
For the computation of the b → sγ vertex at two and three loops we use qgraf [22] to generate all 30 and 591 diagrams, respectively; see Fig. 2 for sample Feynman diagrams.We process them with tapir [23] and use exp [24,25] to prepare FORM [26][27][28] code for their evaluation.We then apply the projectors and perform the Dirac and colour algebra [29] to express the form factors as linear combination of Feynman integrals belonging to 10 and 181 integral families and two and three loops, respectively.The scalar integrals are reduced to master integrals by employing integration-by-parts [30,31] as well as Lorentzinvariance relations [32] and the Laporta algorithm [33] as implemented in Kira [34,35].We use Fermat [36] as back-end to process the coefficients.Before the actual reduction we reduce sets of simpler sample integrals for each integral family and feed them to an improved version of ImproveMasters.m[37] to find a good basis of master integrals in which the denominators of the coefficients factorize in the kinematic and the space-time variable [37,38].The integrals in the amplitude are then reduced to the good basis.As last step we exploit symmetry relations between the integral families with Kira to reduce No QCD corrections to external lines are considered.
the number of master integrals at three loops from 3975 to 479.At two loops we have 14 master integrals.We perform our calculation for general QCD gauge parameter ξ and check that ξ drops out in the final result.

Computation of the master integrals
In this Section we describe the computation of the master integrals which we encounter in the calculation of the b → sγ vertex at two-and three-loop order.We use LiteRed [39,40] and a subsequent reduction with Kira to establish differential equations in for the master integrals [41][42][43][44].At two-loop order we observe poles in the differential equations at the physical cuts x = 0 and 1/2.At three-loop order we have additional spurious poles at various values of x.
At two loops we manage to obtain analytic expressions for all master integrals appearing at NLO, and in the bottom mass counterterm at NNLO.This extends the results of Ref. [17] where the master integrals for the two-and three-particle cut contributions have been obtained only by solving the differential equations numerically with boundary conditions fixed in the large charm mass limit.Furthermore, analytic expansions around m c → 0 have been presented including terms only up to order (m 2 c /m 2 b ) 3 .Here, instead, we provide for the master integrals contributing to the two-particle cuts their analytic expressions with exact dependence on m c /m b , which allows to obtain precise numerical results for an arbitrary ratio of the heavy quark masses.
For the analytic calculation of the two-loop master integrals we use the algorithm outlined in Ref. [45].We do not transform the system of differential equations into a so-called canonical or ϵ form (see e.g.Refs.[46,47]), but decouple coupled blocks of the differential system with the help of OreSys [48] and Sigma [49].This leads to a higher order differential equation for a single integral of the block which can be solved via the factorization of the differential operator as implemented in HarmonicSums [50].The boundary conditions of the integrals are provided in the limit m c → ∞ via a large mass expansion.For the matching of the boundary conditions, it is convenient to first solve the differential equations in the variable y = m b /m c since expansions for y → 0 are easier to compute.Afterwards we perform an analytic continuation to the variable x = m c /m b (further details are given in the Appendix).
To solve the integrals in the variable y we introduce generalized iterated integrals over the alphabet 1 y , Iterated integrals containing the first three letters yield Harmonic Polylogarithms (HPLs) [51].A representation of the iterated integrals in terms of rational letters is beneficial in order to evaluate them to high precision with public programs like ginac [52].
Thus we have to rationalize the occurring square root valued letter with the variable change on the integrals which contain the new letters.
For the HPLs with argument x we do not apply the variable change.In case the new letters come together with the HPL letters, we have to introduce the additional letters √ y/(1 + y) and √ y/(1 − y) for individual iterated integrals after changing from argument x to w.However, we see that contributions containing them cancel in the final amplitude and we can express the results purely in terms of HPLs with arguments x and w, which enables a fast and precise numerical evaluation for an arbitrary mass ratio.Explicit analytic results for two-loop quantities originating from the interference of the operators Q 1,2 and Q 7 can be found in Section 4.1 and are provided in electronic format in the supplementary material [53].
For the master integrals at three loops we apply the "expand and match" approach developed in Refs.[54][55][56] to obtain a semi-analytic expansion which covers all physically relevant values for x that may arise when using different mass schemes for the charm and the bottom quark.The numerical values m c ∼ 1 GeV, . . ., 1.7 GeV and The basic idea of this approach is to use the differential equations in order to construct deep series expansions of all master integrals around several values of x 0 .We make an ansatz for the master integrals as Taylor expansion in x − x 0 (or a power-log expansion if x 0 is a singular point) and insert the ansatz into the differential equations.With Kira we solve the resulting system of linear equations for the expansion coefficients and express them in terms of a minimal set of initial values.
In our case we choose x 0 = 1/5 as a first expansion point.The boundary conditions are fixed by evaluating the master integrals numerically with a precision of 60 digits at this point with AMFlow [57], which implements the auxiliary-mass-flow method [58][59][60][61][62].We use Kira with Fermat as back-end for the reduction.The precise AMFlow results allow to fix all coefficients of the expansions around x 0 = 1/5.We also perform an AMFlow run at x = 1/10 to cross check the expansion constructed around x 0 = 1/5.
This expansion alone already covers the physically relevant x region mentioned above.
In practical applications it is convenient to work with an expansion around x 0 = 0. To obtain a good precision of the expansion coefficients we introduce new expansions.First we construct an expansion around x 0 = 1/10 that we match to the x 0 = 1/5 expansion at x = 0.15.In a second step we produce a power-log expansion around x 0 = 0 that is matched at x = 1/15 to the previous one.
4 Two-and three-loop results In the following we present results for the bare NLO and NNLO contributions.Sample two-and three-loop Feynman diagrams are shown in Fig. 2.

NLO
NLO corrections to b → sγ involving current-current operators are known since long in the literature [63][64][65][66][67].In Ref. [ The function f 0 (z) which enters the NLO prediction is usually written as For the funtions a(z) and b(z) there are no closed analytic expressions known so far but only threefold integral representations [67].
In this work we provide independent cross checks for f 0 (z), f 1 (z), r −1 (z), r 0 (z) and r 1 (z) and present for the first time their analytic expressions in terms of HPLs. 1 Analytic results for a(z) and b(z) can be found in Appendix B. For illustration we show the expressions for f 0 , r −1 and r 0 : where x and w are defined in Eqs. ( 11) and ( 13), respectively.These respresentations can be used for x ≤ 1/2 since for x > 1/2 the HPLs develop imaginary parts.Note that the HPLs in f 0 and r 0 can be expressed in terms of classical polylogarithms.This is not the case for the O(ϵ) terms f 1 and r 1 .
The comparison with Eq. ( 16) allows us to extract the function Re(a(z) + b(z)).This is the first time that analytic results for these functions are presented.Before only a three-dimensional integral representation [67] and analytic expansions for large [15] and small [67] charm quark masses were available.
Our analytic results for r 0 and r 1 are expressed in terms of HPLs up to weight 4 and 5, respectively, where as argument we have both x and w.We find agreement with the analytic expansions given in Refs.[17,64] and can determine the constant C entering r 1 in Eq. (5.5) of Ref. [17].We obtain It deviates by about 0.1% from the numerical value reported in Ref. [17]: C ≈ −292.228.

NNLO
After inserting the three-loop master integrals into the amplitude, we obtain results for the two-particle cuts at NNLO as series expansions around x 0 = 0 and x 0 = 1/5, with numerical coefficients (see Section 3).The expansions include about 40 terms and we report the x 0 → 0 expansion in electronic form in the supplementary material to this paper [53], where for convenience we keep the colour factors and the electric charges of the up-type and down-type quarks distinct.
We present in the following the real part of t 2 in an expansion up to x 2 .For the colour and charge factors we insert their numerical values.Furthermore, we provide for all coefficients six significant digits.Our results for the operator insertions where l x = log(x) and n l = 3 denotes the contribution from closed massless fermion loops while the n c = 1 and n b = 1 contributions arise from closed fermion loops with masses m c and m b , respectively.
In the following we describe the various checks which we performed on our result.We perform the calculation for generic QCD gauge parameter ξ which drops out for t Q 1 2 and t Q 2 2 at the level of the master integrals.In Ref. [18] analytic expansions for the z → 0 limit of the light-fermion contribution have been computed.We compare the expansion coefficients of the x n log k (x) terms up to n = 4 and find agreement with an accuracy of 15 digits or better.We successfully compare the heavy-fermion contributions against Refs.[11,20].
The authors of Ref. [21] have computed the subset of diagrams contributing to the b → sγ vertex at NNLO where the gluon couples only to the charm and strange quarks, i.e. there is no internal bottom quark propagator.Furthermore, all diagrams with fermion, gluon or ghost insertions in the gluon propagator have been excluded.For the comparison we isolate this subset of diagrams in our calculation and then compare t 2 to the corresponding expression in the ancillary file of Ref. [21].We specialize our result to Feynman gauge since this is the choice of Ref. [21] (the set of diagrams considered in Ref. [21] is not gauge invariant).Using our expansion around x 0 = 1/5 we observe agreement of more than 10 digits for the terms proportional to Q s , the charge of the down-type quark, and of about 5 digits for the terms proportional to Q c , the charge of the up-type quark.According to Ref. [21] the program FIESTA [70] has been used to fix parts of the boundary conditions for the second part.This may explain the reduced accuracy.Using our expansion around x 0 = 0 we can also directly compare the coefficients in front of x n log k (x).As before, we observe about 5 digits agreement for the Q c terms and at least 18 digits for the Q s terms.
In Ref. [71] the contributions induced by four-quark operators are computed via the reverse unitarity method by considering bottom quark two-point functions with an insertion of Q 1 or Q 2 and Q 7 .The three-and four-loop contributions yield the NLO and NNLO corrections, respectively, where the latter contains two-, three-and four-particle cuts.In Ref. [71] the two-particle cut contribution is computed.We compare our NNLO result for Ĝ(2),2P,Q tree and find agreement at the level of 10 −15 for x = 1/5.Let us stress that the calculations performed in Ref. [71] and the one in the current paper are completely independent and thus the comparison constitutes an important cross check for both calculations.
The Ward identity in Eq. ( 7) holds for renormalized quantities.Since it holds at NLO, all NNLO contributions obtained by multiplying the NLO amplitude by a global renormalization factor also fulfil the Ward identity.It does not hold for the mass counterterm contribution and the bare three-loop amplitude.However, in the combination of the three-loop corrections to the form factors and the contributions from the bottom mass counterterm the relation (7) has to be fulfilled.For our calculation we have checked that this is indeed the case.Using the expansion around x 0 = 0 we observe that the numerical cancellation in the combination t 3 + t 2 /2 is of order 10 −12 or smaller for x = 0.2 while it reaches the level of 10 −5 at x = 0.4.At the latter point an accuracy of at least 12 digits is obtained in case the expansion around x 0 = 1/5 is used.

Conclusions
In this paper we compute three-loop corrections to the b → sγ vertex induced by the current-current operators Q 1 and Q 2 .We apply a semi-analytic method and construct expansions around x 0 = 0, x 0 = 1/10 and x 0 = 1/5 including about 40 terms, which covers the physically relevant values of the mass ratio in different mass schemes.Furthermore, we provide for the first time analytic results for the two-loop contributions.
The current predictions for the B → X s γ branching ratio at NNLO employ only an interpolation to estimate the m c dependence in the interference term between Q 1 , Q 2 and Q 7 .Such interpolation is responsible for 3% out of the 5% (= (3%) 2 + (3%) 2 + (2.5%) 2 ) theoretical uncertainty.Our result for the b → sγ vertex at NNLO allows to calculate the two-particle cut contribution to the B → X s γ branching ratio (stemming from currentcurrent operators) without such an interpolation.Therefore it constitutes an important building block towards the reduction of the theoretical error in B → X s γ, once the missing three-and four-particle cut contributions are available as well.
the European Union's Horizon 2020 research and innovation program under the Marie Sk lodowska-Curie grant agreement No. 101065445 -PHOBIDE.We have used the program FeynGame [72] to draw the Feynman diagrams.

A Details of the two-loop master integral calculation
In this appendix we present additional details regarding the calculation of the master integrals at two loops, in particular about the analytic continuation.
As discussed in Section 3, at two loops we solve the master integrals using the method of differential equations and fix the boundary conditions in the limit m c → ∞ via a large mass expansion.To match the boundary conditions it is convenient to first solve the differential equations in the variable y = m b /m c since the expansions for y → 0 are easier to compute.It is therefore necessary at some point to perform an analytic continuation to the variable x = m c /m b .We find it beneficial to use the argument y for the master integrals and only apply the appropriate variable change once the physical amplitude is assembled.In the following we explain how the analytic continuation is performed.
We express the master integrals in terms of iterated integrals with letters from the set given in Eq. ( 12).In the end we want to evaluate the iterated integrals for values m c /m b < 1, i.e. y > 1.As we can see in Eq. ( 12) the iterated integrals therefore require an analytic continuation since the letters develop poles at y = 1 and y = 2.For iterated integrals which only contain the first three letters, i.e. harmonic polylogarithms, the analytic continuation to x = 1/y is well known; we use HarmonicSums to do this step.For the integrals which contain non-standard letters we do the analytic continuation with the help of the differential equations.We start by taking the derivative of the iterated integrals with respect to y, which leads to iterated integrals of lower weight, then do the variable change x = 1/y and integrate over x.The differential equations with iterated integrals of weight 1 can then be integrated trivially since they reduce to integrals over algebraic functions which can be performed by using the definition of iterated integrals and possibly integration-byparts identities to reduce higher powers of the letters.For the contributions with higher weight one can proceed iteratively, since the derivative of an iterated integral of weight n with respect to its argument only depends on iterated integrals of weight n − 1.After the change of variables we therefore have to perform a single integral over expressions of algebraic functions possibly multiplying iterated integrals of lower weight.These integrals can again be performed using the definition of iterated integrals and integration-by-parts identities.However, in order to find the analytic continuation of the iterated integrals in this way we have to fix the integration constants.
It is convenient to fix this integration constants at y = 2 since after the shift y → 2ỹ the non-standard letters in Eq. ( 12) transform to 1/(1 + ỹ), 1/(1 − ỹ), 1 − ỹ2 .The square root valued letter can further be rationalized which subsequently leads to cyclotomic harmonic polylogarithms [73] for which the constants at argument ỹ = 1 are known.
Note that it is crucial to consider the analytic continuation of the whole amplitude and not only a single iterated integral.The individual integrals depend on various constants for which we do not have a replacement in terms of known transcendental numbers.However, in the amplitude these cancel out and we are left with well known constants.We also observe that the analytic continuation of single master integrals is much more involved than the one of the physical amplitude.
Finally, we notice that the analytic continuation of the HPLs completely decouples from the one of the generalized iterated integrals.This can be seen by introducing y = 1/x + iδ ( ′ ) , where δ is used for the harmonic polylogarithms and δ ′ for the generalized iterated integrals.Both quantities are infinitesimally small.After performing the analytic continuation we see that the expression is independent of δ while we have to choose δ ′ as positive in order to identify the correct analytic continuation for the generalized iterated integrals.This effect is not necessarily expected.In other words: From the observation that the analytic continuation of the HPLs is independent from the one of the "non-standard" iterated integrals we conclude that they "decouple" and a separate variable transformation can be applied to the latter.Note that this discussion only affects the imaginary parts of t Q 1 2 and t Q 2 2 and is thus not relevant for the results presented in the main part of this paper.

B Analytic results for a(z) and b(z)
The functions a(z) and b(z) defined in Ref. [67] can be obtained as the contributions proportional to Q u and Q d with an additional subtraction to fulfill the normalization a(0) = b(0) = 0.In the conventions of Ref. [67] The variables x and w are defined in Eqs.(11) and (13), respectivley.Note that the separation in real and imaginary parts is correct below the two-particle threshold z < 1/4.
From our expressions we can find analytic expressions for the constants in Ref. [67].They agree with the results derived in Ref. [10] by a direct integration of the three-fold integral representation.They are given by a The results in this Appendix together with the expressions given in Eq. (3.1) of Ref. [67] provide analytic results also for all NLO penguin contributions.
In Refs.[68,69] the two-loop form factors with an off-shell photon have been calculated.Obtaining the on-shell limit from the representation in terms of Goncharov multiple poly-logarithms is, however, a non-trivial task, since many letters become singular in this limit.

Figure 1 :
Figure 1: Three-and four-loop sample diagrams for the interference of Q 1,2 and Q 7 .The dashed black lines represent possible cuts through the diagrams.

Figure 2 :
Figure 2: Two-and three-loop sample diagrams contributing to the decay vertex b → sγ.No QCD corrections to external lines are considered.