Master integrals for the two-loop penguin contribution in non-leptonic B-decays

We compute the master integrals that arise in the calculation of the leading penguin amplitudes in non-leptonic B-decays at two-loop order. The application of differential equations in a canonical basis enables us to give analytic results for all master integrals in terms of iterated integrals with rational weight functions. It is the first application of this method to the case of two different internal masses.


Introduction
The study of flavour-changing quark transitions provides an important indirect probe to search for new heavy particles as well as to test the CKM mechanism of flavour mixing and CP violation. One prominent class of such transitions are non-leptonic B-meson decays, which offer a rich and interesting phenomenology including many CP-violating asymmetries. Non-leptonic two-body decays therefore play a central role at current and future B-physics experiments. The extraction of the underlying decay amplitudes is, however, complicated by the strong-interaction dynamics of the purely hadronic environment. A systematic formalism to compute the hadronic matrix elements arises in the heavy-quark limit [1][2][3]. Schematically, where M 1,2 are light (charmless) pseudo-scalar or vector mesons and Q i is a generic operator of the effective weak Hamiltonian. The hadronic dynamics in the above factorisation formula is encoded in a form factor F and in light-cone distribution amplitudes φ. The hard-scattering kernels T , on the other hand, can be computed to all orders in perturbation theory in a partonic calculation. In the last few years, the perturbative corrections have been worked out to next-to-next-to-leading order (NNLO) accuracy. While the full set of O(α 2 s ) corrections to the spectator-scattering kernels T II i is known [4][5][6][7][8], NNLO corrections to the kernels T I i have to date only been determined for the topological tree amplitudes [9][10][11].
The missing NNLO ingredient consists of a two-loop calculation of the hard-scattering kernels T I i in the penguin sector. The calculation involves various types of operator insertions, for details we refer to a future publication [12]. The one-loop contribution of the magnetic dipole operator has been computed in [13]. The most difficult part of the calculation consists in the computation of massive two-loop penguin diagrams like the ones shown in Fig. 1. Whereas the integrals that entered the two-loop tree calculation [14,15] can be expressed in terms of Harmonic Polylogarithms (HPLs) [16], the massive propagator in the penguin loop introduces an additional scale and complicates the calculation. In the present paper we give analytic results for the master integrals that arise in this calculation.
A convenient technique for the calculation of multi-scale integrals is the method of differential equations [17][18][19]. In combination with integration-by-parts identities [20,21] and Laporta's reduction algorithm [22], the master integrals are computed by solving a set of differential equations where the derivatives are taken with respect to the external scales of the process. It has recently been pointed out that the solution simplifies considerably if the basis of master integrals is chosen appropriately [23]. We will discuss the properties of such a canonical basis in detail below. The method has been successfully applied to compute various massless as well as massive two-loop and three-loop integrals [24][25][26][27][28][29][30][31][32][33]. The present calculation is the first application of the method in which the integrals have two different internal masses. Our paper is organised as follows. We first discuss the kinematics of the process and introduce a generalisation of the HPLs in Section 2. The canonical basis of master integrals is defined in Section 3, and analytic results for all master integrals are given in Section 4. We comment on several cross-checks of our calculation in Section 5, before we conclude in Section 6. The paper is complemented by three appendices with various technical details, as well as an electronic file that contains the analytic results of all master integrals and is attached to the arXiv submission of the present work.

Kinematics
The kinematics of the process is depicted in Fig. 1. We write p b = q + p with p 2 b = m 2 b and p 2 = q 2 = 0. The momentum q of the emitted final state meson is split up into two parallel momenta q 1 = uq and q 2 = (1 − u)q ≡ūq of the quark and anti-quark, respectively, where u ∈ [0, 1] is the convolution variable that enters the first term of Eq. (1.1). The quark in the penguin loop can either be massless in the case of up, down and strange quarks, or massive of mass m c or m b in the case of charm or bottom. For massless quarks, the master integrals are already known from the calculation of the two-loop tree amplitudes in [14,15]. We therefore only consider the situation with a massive quark in the penguin loop in the following. The problem then depends on two dimensionless variables, which we choose as the momentum fractionū of the anti-quark and the mass ratio In order to express the solution to the master integrals in terms of iterated integrals with rational weights, it will be convenient to trade the variablesū and z f for other sets of variables. Our default choice is the set (r, s) with which, when solved for the original variables, implies Let us have a look at the possible values of s. Whenū runs from 0 → 1, the variable s for 4z f > 1 runs from +i∞ → r along the imaginary axis. For 4z f < 1, s runs from +i∞ → 0 along the imaginary axis, followed by 0 → r along the real axis. In this case the threshold atū = 4z f is mapped onto s = 0. Another convenient choice of variables will be the set (r, s 1 ), with and z b = 1 − iη. The variable s 1 runs from +i∞ → +i √ 3 along the imaginary axis once we letū run from 0 → 1.
A third choice of variables consists of the set (r, p) with When solved for the original variableū one obtains

Iterated integrals
One of the classical examples of iterated integrals are HPLs [16]. They are generalisations of ordinary polylogarithms and appear in many calculations of higher-order corrections in perturbative Quantum Field Theory. The HPLs are defined by H a 1 ,a 2 ,...,an (x) = x 0 dt f a 1 (t) H a 2 ,...,an (t) , (2.6) where the parameters a i can take the values 0 or ±1, and n is called the weight of the HPL.
In the special case that all indices are zero, one defines H 0n (x) = 1 n! ln n (x). The weight functions f a i (x) are given by In addition one assigns the weight k to numbers like π k , ln k (2) and the Riemann zeta function ζ k , and one uses that the product of two expressions of weights k 1 and k 2 has weight k 1 + k 2 . These definitions were generalised in [34] by introducing linear combinations of f 1 (x) and f −1 (x), the so-called "+" and "−"-weights, according to (2.9) In the present work we further generalise the weights by allowing more generic expressions to appear in the weight functions. For any expression w = 0 we define (2.10) and accordingly (2.12) Also with these newly introduced weight functions we define a general HPL by means of Eq. (2.6), but we also allow the weights (2.10) -(2.12) to enter the integrand. In the current calculation, we encounter the following expressions for w, We will refer to w 1 -w 5 as rational weights, since any of the w i is rational either in r or m f , given that As a matter of fact, the generalised HPLs are closely related to Goncharov polylogarithms [35], which are defined by G a 1 ,a 2 ,...,an (x) = x 0 dt t − a 1 G a 2 ,...,an (t) (2.14) and G 0n (x) = H 0n (x). We can therefore always write a generalised HPL as a linear combination of Goncharov polylogarithms, for example and similarly for higher weights. The structure of the differential equations in the subsequent sections reveals that the results of the master integrals are most compactly written in terms of HPLs with generalised weights. For their numerical evaluation described in Section 5, however, we prefer the notation in terms of Goncharov polylogarithms.

Canonical Basis
Within dimensional regularisation where space-time is analytically continued to D = 4 − 2 dimensions, integration-by-parts identities [20,21] provide non-trivial relations between different loop integrals. It has now become a standard tool to use automated reduction algorithms to express complicated multi-loop calculations in terms of a much smaller set of irreducible master integrals. The choice of the master integrals is, however, not unique.
Henn recently conjectured that the set M of master integrals can always be chosen in a way such that the set of differential equations assumes the form [23] where x n are dimensionless kinematic variables and A m (x n ) is a matrix which does not depend on . In this form the system of differential equations decouples order-by-order in the -expansion. The system (3.1) can be written as a total differential, The matrixÃ contains the relevant information about the structure of the occurring weight functions. Together with suitably chosen boundary conditions, this entirely fixes the solution. As an additional feature, the solutions to the master integrals contain functions that are of uniform weight at each order in , and the weight increases by unit steps as one goes from one power to the next one in the -expansion. As a consequence, by assigning the weight −1 to and multiplying the master integrals by an appropriate power of , one can achieve that the total weight of each master integral is zero to all orders in . Integrals with the latter property and a system of differential equations of the form (3.2) will be referred to as a canonical basis. At present there does not exist a systematic algorithm to find a canonical basis of master integrals. The construction therefore requires some level of experimentation, for some guidelines cf. the discussions in [24,27,29,31,32]. In the current calculation we mainly used explicit integral representations to find the canonical basis. The basis consists of 29 master integrals which we denote by M 1−29 . In terms of the integrals I 1−34 defined in Fig. 2, they are given by

11)
M 10 (r, s) = 3 u I 11 (ū, z f ) , (3.12) (3.31) The variables r, s, s 1 and p have been introduced in Section 2.1, and the definition of the integrals I 4,5 (z f ) can be found in Appendix B. In addition there are seven auxiliary integrals, labeled M 1−7 , which are already known from previous calculations but which are needed in order to close the system of differential equations. In the given integral basis the system of differential equations takes the form (3.2). Instead of one large matrixÃ, we solve each topology separately and in turn get several smaller matricesÃ k . We give the solution to the basis integrals M 1−29 in the next section, together with the relevant boundary conditions. The solution to the auxiliary integrals M 1−7 can be found in Appendix B.

Results
We write the results for the master integrals in the form with the number of loops L and an integer N which denotes the sum of all propagator powers. The integralM is therefore dimensionless. Our integration measure per loop is d D k/(2π) D and the pre-factor S Γ reads Once the differential equations are set up, the only missing ingredient are the boundary conditions. It turns out that the following conditions -almost all of which describe the vanishing of an integral in a particular kinematic point -are sufficient to write down the entire solution to an integral. We find that M 1,3,4,6,7,9,11,14−17,20−22 (r, s) and M 27 (s 1 ) vanish inū = 0, corresponding to s = +i∞ or s 1 = +i∞. Furthermore, M 8,10,12,13,18,19 (r, s), M 2 (ū), M 23 (r, s 1 ), M 26 (s 1 ), and M 28,29 (r, p) vanish inū = 1, corresponding to s = r, which can be derived using the Laporta reduction algorithm [22]. All these considerations lead to the full set of solutions which we list below.

M 1
As a warm-up exercise and to demonstrate how the method of differential equations in the canonical basis works, we consider the one-loop integral (4.4) The auxiliary integral appears as a subtopology and has to be taken into account in order to make the system of differential equations complete. The solution to the auxiliary integral M 1 (z f ) is elementary and can be found in Appendix B.
In terms of the variables r and s, the system of differential equations becomes and The system of differential equations can be brought into the canonical form (3.2), with Solving Eqs. (4.6) and (4.8) together with the aforementioned boundary condition gives The solution can also be obtained from the following closed form, by expanding the hypergeometric function e.g. with HypExp [36,37].

M 2
From now on, we will not give the explicit form of the differential equations anymore, but only the corresponding matricesÃ i and the final solution to the integrals. The integral M 2 only depends on one kinematic variable, which we choose to beū. The set of integrals is now given by M = M 2 (ū),M 1 (z f = 1),M 2 (ū) , and we havẽ Also here the topology consists of six integrals, namely and the matrixÃ 8,9 (r, s). Owing to simple boundary conditions, the result is quite short, (4.28) and the matrixÃ 10,11 (r, s). The result is rather long since we need functions up to weight four in M 10 (r, s), Again we need seven integrals to complete the system of differential equations. They are (4.34) The integrals in this topology only depend on one non-trivial scale ratio, and their solution can be written in terms of ordinary HPLs. The topology involves five integrals, 35) and the matrixÃ 15−17 (r, s). The result reads The results up to functions of weight three arẽ  where r = i √ 3 corresponds to z f = 1. Here we choose the set of variables (r, s 1 ). The fact that the number of integrals is large is not the only complication of this topology. As can be seen from the matrixÃ 23−25 (r, s 1 ) in Eq. (A.11), many factors appear in the differential equations which are irrational in both r and s 1 . For example, (4.47) Fortunately, we can still find a form of the differential equations which allows us to apply the formulas for iterated integrals from Section 2. There are two reasons why this is possible. First, there exist variable transformations which rationalise either of the square roots, namely For later convenience we also define 50) which correspond to the limit s 1 → +i √ 3 of t and v, respectively. Second, it turns out that we only need the lowest order in the -expansion for each of the integrals M 23−25 . This ensures that M 24 appears only in combination with t, whereas M 25 appears only with v, without any admixture of the respective other variable. This does not hold at higher orders in , which can be concluded for instance from the appearance of the logarithm L 15 inÃ 23−25 (r, s 1 ) in Eq. (A.11) which contains both t and v. Having said this, we find   The integrals in this topology already appeared in the two-loop calculation of the tree amplitudes [9][10][11], where explicit Mellin-Barnes (MB) representations have been used for their numerical evaluation (for a convenient parameterisation cf. also the appendix of [38]). With the current techniques, we are now in the position to compute these integrals analytically. For this topology, it will be convenient to use the variables (r, p) defined in Section 2. We need seven integrals, (4.63) and the matrixÃ 28,29 (r, p). The integral M 28 is required up to functions of weight three, but M 29 is only needed up to weight two. The solution is again lengthy, and we introduce a short-hand notation for p 0 = 1 − 2 √ z f . We find in (4.60). Another purely numerical method is sector decomposition, where we used both the SecDec-package [43,44] as well a Mathematica-based in-house routine. For the most complicated coefficients the numerical evaluations confirm the analytic results at the level of 10 −4 , and for the simpler coefficients the agreement is of several orders of magnitude better.

Conclusion and Outlook
We computed the master integrals that arise in the computation of the two-loop correction to the vertex kernel of the leading penguin amplitudes in non-leptonic B-decays. The calculation is complicated by the presence of two non-trivial scales (ū and z f = m 2 f /m 2 b ), as well as the kinematic threshold atū = 4z f . We computed the master integrals in a recently advocated canonical basis, which enabled us to derive analytic results for all master integrals in terms of generalised HPLs. The results are given up to the relevant order in the -expansion that is needed to obtain the finite terms of the penguin amplitudes. Our calculation is the first application of a canonical basis to integrals with two different internal masses. Apart from the integral basis, we find that the choice of the kinematic variables is of utmost importance since it renders the logarithms in the matricesÃ k rational and therefore makes the formulas for iterated integrals applicable.
The results of this paper form the basis to derive fully analytic expressions for the hardscattering kernels T I i in the factorisation formula (1.1). In phenomenological applications, one has to integrate over the product of the kernels and the Gegenbauer expansion of the light-cone distribution amplitudes. The presence of the charm threshold makes the numerical evaluation of the convolutions delicate. The threshold is much easier to handle in an analytic approach, and the convolutions can now be computed to very high precision.
The integrals presented here are also relevant for other applications such as rare or radiative B-meson decays. For example, the two-loop QCD correction to the matrix elements of current-current operators in inclusiveB → X s + − decays have to date only been computed numerically [45] or as expansions in the lepton-invariant mass q 2 [46,47]. With the present results, one can now obtain completely analytical expressions for any value of q 2 . In exclusiveB → K ( * ) + − decays, one can study non-factorisable corrections to charm-loop effects.

A. MatricesÃ k
In this appendix we list the matricesÃ k for the different subtopologies. To this end, we define the following logarithms, The matricesÃ k now assume a compact form,