Two-loop master integrals for non-leptonic heavy-to-heavy decays

We compute the two-loop master integrals for non-leptonic heavy-to-heavy decays analytically in a recently-proposed canonical basis. For this genuine two-loop, two-scale problem we first derive a basis for the master integrals that disentangles the kinematics from the space-time dimension in the differential equations, and subsequently solve the latter in terms of iterated integrals up to weight four. The solution constitutes another valuable example of the finding of a canonical basis for two-loop master integrals that have two different internal masses, and assumes a form that is ideally suited for a subsequent convolution with the light-cone distribution amplitude in the framework of QCD factorisation.


Introduction
Non-leptonic B-decays are interesting for a number of phenomenological applications like the extraction of CKM elements and the study of CP asymmetries. Their study has already entered the area of precision physics, both on the experimental [1] and on the theoretical side. However, their theoretical description is complicated by the purely hadronic environment, entailing QCD effects from many widely separated scales. The two main approaches to non-leptonic B-decays are flavour symmetries of the light quarks [2] and factorisation frameworks such as pQCD [3] and QCD factorisation (QCDF) [4][5][6]. In the latter framework, next-to-leading order (NLO) corrections to both, heavy-to-heavy [5] and heavy-to-light [4,7] transitions have been known since more than a decade. More recently, also next-to-next-to-leading order (NNLO) results for heavy-to-light decays have become available [8][9][10][11][12]. In the present article, we consider NNLO corrections also to the heavy-toheavy decays such as B → Dπ in the framework of QCDF [13]. In the heavy-quark limit, the decay amplitude forB 0 → D + π − is given by [5] where O i are the operators from the effective Hamiltonian that describe the underlying weak decay. The F B→D j form factors and the pion light-cone distribution amplitude (LCDA) Φ π (u), with momentum fractions u and 1 − u shared among the pion constituents, are the non-perturbative inputs. The hard-scattering kernels T ij (u), on the other hand, can be evaluated in a perturbative expansion in the strong coupling, and are known in QCD to NLO accuracy [5]. Yet it is interesting to go beyond NLO in B → Dπ transitions: Since the contribution at NLO is colour suppressed and appears alongside small Wilson coefficients, the NNLO corrections may be significant in size. Moreover, since there is neither a colour-suppressed tree amplitude nor penguin contributions, and spectator scattering and weak annihilation are power-suppressed [5], we have only the vertex kernels to the colour-allowed tree amplitude. A precise theory prediction of this single contribution, together with comparison to experimental data, might give a reliable estimate of the size of power corrections in the QCDF framework.
The evaluation of Feynman diagrams that contribute to the NNLO hard-scattering kernel amounts to the computation of ∼ 70 two-loop diagrams. By using contemporary techniques to evaluate multi-loop integrals, the two-loop Feynman diagrams are reduced to a small set of a few dozens of master integrals. A powerful method to evaluate the latter analytically are differential equations [14][15][16]. This method was recently refined by Henn [17]. Considering that the basis of master integrals is not unique, Henn discovered that in a suitably chosen basis -denoted as canonical basis -the differential equations can be cast into a form that factorises the dependence on the kinematic variables from that on the number of space-time dimensions D. In this case, the solution is expressed in terms of iterated integrals. This method was recently applied to a number of problems for loop [11,[18][19][20][21][22][23][24][25][26] and phase-space [27,28] integrals.
To the present day, the construction of the canonical basis is mostly based on experience or experimentation, rather than on a systematic procedure, although developments in this direction have recently become available [21,29,30]. In the future it would be most desirable to have a general algorithm for finding a canonical basis for arbitrary external kinematics and numbers of loops, legs, scales, and space-time dimensions. Therefore, every non-trivial example of a canonical basis is most valuable, and our results contribute towards finding a general algorithm for constructing the canonical basis.
Last but not least, if the master integrals that enter the hard-scattering kernels T ij (u) are written in terms of iterated integrals, the convolution with the pion LCDA in (1.1) simplifies to a large extent. Our results therefore catalyse the steps necessary to obtain the decay amplitudes considerably, and constitute an important step towards the phenomenology of B → Dπ decays at NNLO in QCDF. This paper is organized as follows. In section 2 we introduce the kinematics of the two-body decay and present the generic form of the differential equations with respect to the kinematic variables. We proceed by defining Goncharov polylogarithm in section 3, which are a class of iterated integrals suited to describe the solutions to the differential equations. In section 4 the canonical basis is defined and the expressions for the master integrals in this basis are presented. We also elaborate on strategies to find a canonical basis. The boundary conditions for the integrals are discussed in section 5 and the results are presented in section 6. In section 7 we comment on the performed cross-checks before concluding in section 8. In appendix A we collect the matrices that contain all relevant information on the differential equations. The analytic results of all master integrals are also available electronically [31].

Kinematics
We consider the kinematics of the decayB 0 → D + π − , which emerges from the underlying weak transition b → cūd. A sample of Feynman diagrams contributing to the two-loop hard-scattering kernels is given in figure 1. The complete set of diagrams consists of those shown in figures 15 and 16 of [5], supplemented by gluon self-energy insertions in one-loop diagrams. All external momenta are taken to be incoming throughout this work. q 4 and q 3 denote the external momenta of the b and the c quark, respectively, which fulfill the on-shell constraints q 2 4,3 = m 2 b,c . The constituents of the pion share the momentum q with q 1 = uq and q 2 = (1 − u)q ≡ūq, where u ∈ [0, 1] is the momentum fraction of the quarks inside the pion entering eq. (1.1) in a convolution of the hard-scattering kernel with the pion LCDA. We consider the pion to be massless, i.e. q 2 = q 2 1,2 = 0. Due to the linear dependence of the momenta, q 1 + q 2 = q = −q 3 − q 4 , the kinematics is completely determined by two of the on-shell conditions and one additional kinematic invariant, for instance We apply commonly used multi-loop techniques which include integration-by-parts identities [32,33] and the Laporta algorithm [34], and reduce the two-loop Feynman diagrams to master integrals [35,36]. Furthermore, we construct the differential equation of the latter with respect to kinematic variables. In the derivation of eq. (1.1) the charm quark was assumed to be heavy. Hence, the ratio m c /m b remains fixed in the heavy-quark limit and our master integrals depend on two scales: the momentum fraction u and the ratio of the heavy quark masses z ≡ m 2 c /m 2 b . They are further functions of the kinematic invariants (2.1) Thus, the total derivative of a generic master integral C with respect to u is given by whereas the one in z reads Figure 1: Sample of Feynman diagrams: q 4 and q 3 are the momenta of the quark lines with masses m b and m c , respectively. q 1 = uq and q 2 =ūq are the momenta of the light quark and anti-quark, respectively. q = q 1 + q 2 is the momentum of the pion. All momenta are incoming. The black square denotes an operator insertion from the weak effective Hamiltonian.
The computation of ∂C/∂z is straightforward. The partial derivatives of C with respect to the kinematics on the r.h.s. of eq. (2.4) can be expressed in terms of partial derivatives with respect to the momenta q 3,µ and q 4,µ [37], which can be easily carried out. Note that the last term on the r.h.s. vanishes since dq 2 4 /dz = 0. We finally obtain This is the differential equation with respect to z valid for a generic master integral C(u, z).

Iterated integrals and Goncharov polylogarithms
The classical example of iterated integrals is given by the harmonic polylogarithms (HPLs) [38]. They generalise the ordinary polylogarithms and are defined by H a 1 ,a 2 ,...,an (x) = x 0 dt f a 1 (t) H a 2 ,...,an (t) , where the parameters a i can be 0 or ±1, and n is the weight of the HPL. The integral (3.1) diverges for HPLs with trailing zeroes. In order to handle HPLs in such cases, one defines H 0n (x) = 1 n! ln n (x). The weight functions f a i (x) are simply The HPLs fulfil a Hopf algebra according to where a b are all possibilities of arranging the elements of a and b such that the internal order of the elements of a and b is preserved individually (cf. also [39]). Hence the product of two HPLs of weights w 1 and w 2 has weight w 1 + w 2 . The Hopf algebra can also be used to extract singular behaviour near x = 0 or x = 1. Due to the relation with k − 1 zeroes and k > 1, one also assigns the weight k to numbers like ζ k and π k . A generalisation of the HPLs are the Goncharov polylogarithms [40], whose definition reads G a 2 ,...,an (t) (3.5) and G 0n (x) = H 0n (x). They fulfil a Hopf algebra that has the same structure as (3.3), and allow for more general weights a i than just 0 or ±1. In particular, in multi-scale problems the argument x can be represented by one scale, and the remaining scales are comprised in the weights a i . In our problem at hand, it is most convenient to choose u as the argument of the Goncharov polylogarithm whenever there is a dependence on this scale, bearing in mind that this choice simplifies a subsequent convolution with the light-cone distribution amplitude, which in a Gegenbauer expansion is a u-dependent polynomial. In this case the weights are either integer (0, ±1) or one of the following six z-dependent weights 1 , Goncharov polylogarithms that do not depend on u are written in terms of integer weights and argument z or √ z. Products of Goncharov polylogarithms of the same argument are expanded by means of the Hopf algebra.

The canonical basis
We work in dimensional regularisation with D = 4−2 and evaluate the two-loop, two-scale master integrals by applying the method proposed by Henn [17]. Considering a specific power in the -expansion of a master integral, the associated function is called uniform if each summand has the same weight. Moreover, a uniform function is called pure, if its derivative with respect to any one of its arguments yields a uniform function whose weight is lowered by one unit.
The proposal in [17] now states that a basis C of master integrals can be found such that the system of differential equations in the kinematic variables x j is given by where d i ≡ d/dx i . The C(x j , ) denote the N master integrals and A i (x j ) are N × N matrices which are independent of . It turns out that eq. (4.1) can be expressed in a compact form with the functionÃ determined by the differential d iÃ = A i . We note thatÃ, together with the boundary conditions, completely determines the solution to a master integral. The master integrals in such a basis have in turn several pleasant features: First, the solution decouples order-by-order in the -expansion. Second, it is given by pure functions to all orders in . Consequently, assigning a weight −1 to each power of the expansion parameter and multiplying each master integral by an appropriate power of renders the total weight of the master integral to be zero to all orders. Third, the solution can be expressed in terms of iterated integrals. If the coefficients A i (x j ) are rational functions of the x j , the Goncharov polylogarithms discussed above represent a suitable class of iterated integrals to describe the master integrals. We will refer to such a basis as a canonical basis.
In the absence of a completely general algorithm for the systematic construction of the canonical basis, the procedure of finding such a basis requires a certain amount of experience and experimentation. In our case, we start from a "traditional" basis that consists of undotted and singly-dotted integrals, and compute them up to terms that involve functions of weight two. For this task, alternative approaches like Feynman parameters or Mellin-Barnes representations [41,42] have to be used. Afterwards one plugs these expressions into seemingly more complicated integrals like the ones in figures 2 and 3 and investigates if the resulting expressions are uniform or even pure. This method is mostly based on trial and error, but has proven to be successful as we show below.
In the case at hand, many master integrals can be adopted from several B → ππ calculations [8][9][10]43]. In order to describe the yet unknown ones in the canonical basis, a set of 39 integrals is needed. We obtain the following expressions for the canonical master integrals C 1−39 in terms of the integrals I 1−42 , which are defined in figures 2 and 3 (x = 1 − x).

Boundary conditions
Before we present the differential equations, we specify the boundary conditions that are used to completely fix the solution. In the simplest cases, the master integrals vanish in a specific kinematic point. This is the case for C 13,38,39 , which vanish in u = 0, whereas C 3,4,5,6,16,17,23,24,28,29,30,35 vanish in u = 1. Moreover, C 19,31,32 vanish in z = 0, whereas C 7,14,33,34 vanish in z = 1. In other cases we find special relations between integrals, that hold either in general, or in certain kinematic points, and can be used as boundary conditions. Examples are the relation C 26 = z − C 21 , or the following relations that hold in u = 1, where the symbol "↔" is used for the corresponding "mass-flipped" integral, in which m c ↔ m b and q 3 ↔ q 4 , see section 6 for more details. Hence, the integrals C ↔ 19,20 can be easily obtained from C 19,20 or from [31]. Relations that have a similar structure than (5.1) hold in z = 1 for For the remaining integrals we either use that they assume simple, closed forms that are valid to all orders in the -expansion, or asymptotic forms as u → 0 or z → 0. Examples of the former type are (see below in section 6 for the precise definition ofC i ) where for C 36 we give the result for each loop separately, such that also the boundary conditions for C 21,37 can be read off. Asymptotic expansions as u → 0 or z → 0 were derived by means of MBasymptotics.m [44] for

Results
In order to facilitate the presentation of the results we write the master integrals as with an integer n that denotes the sum of all propagator powers, such that the integralC is dimensionless. Our integration measure is d D k/(2π) D per loop and we use the pre-factor Besides the integrals defined in section 4, the QCD amplitude also contains the same set of integrals but with m c ↔ m b and q 3 ↔ q 4 . We will refer to these as "mass-flipped" integrals and denote them as C ↔ , see section 5. However, we note here that in order to defineC ↔ we factor out an appropriate power of m b , rather than m c . As stated earlier the QCD amplitude requires terms of order O( 4 ) for most of the integrals. However, in order to keep the paper at a reasonable length, we only give terms up to order O( 3 ) explicitly below. If desired, terms of weight four can be derived from theÃ and the boundary condition, which we actually give to weight four. Moreover, we refrain from presenting the "mass-flipped" integrals explicitly. They can be obtained by letting z → 1/z, keeping in mind that analytic continuation is done via z → z − iη, with infinitesimal η > 0. We provide the results to all integrals, including the "mass-flipped" ones, to order O( 4 ) in electonic form in [31].
Last but not least, instead of dealing with one large 39 × 39 system of equations, we solve each topology separately and therefore deal with several, smaller matricesÃ i which we collect in appendix A. This finally puts us in the position to present the analytic results to the C 1−39 .

C 33 and C 34
This topology has seven integrals, of which only C 33 -C 34 have not yet appeared in the previous topologies. The integrals are ordered as C = C 33 ,C 34 ,C 7 ,C 22 ,C 31 ,C 32 ,C 12 . (6.40) The corresponding matrix isÃ 33,34 . The solution to O( 3 ) is very short At order O( 4 ) the solution requires also Goncharov polylogarithms of argument √ z.

C 35
This topology has three integrals, and only C 35 is new. The integrals are ordered as The corresponding matrix isÃ 35 . The solution reads The corresponding matrix isÃ 36,37 . The solution reads +4 G 1,1,0 (z) − 12 G 1,1,1 (z) + G 1,a 1 ,0 (u) + 2 G a 1 ,0,1 (u) + 2 G a 1 ,1,0 (u) + 1 6 iπ 3 ] + O( 4 ) , (6.46) 6.9 C 38 and C 39 These integrals arise from diagrams with a massive quark loop inside a gluon propagator. They appeared in a slightly different version already in the calculation of the two-loop tree amplitudes in B → ππ [9,10], and analytic results were recently derived in [11] as M 28,29 . It turns out that the results of C 38,39 can be obtained from the latter reference if one adjusts the kinematics to the present problem. To be precise, one has to replace u → u (1 − z) (6.48) in the expressions in [11]. That is, in the definition of the canonical basis (cf. (3.30) and (3.31) of [11] and (4.40), (4.41) of the present article), and also in the solution, eqs. (4.64) and (4.65) of [11]. In particular, the kinematic variable p changes to (z = 1 − z)

Checks
In order to validate the analytic results presented above, we performed several checks of analytic and numeric nature. Those integrals that possess a closed form in terms of hypergeometric functions were analytically expanded in using HypExp [45,46]. Subsequently, we re-wrote the resulting polylogarithms and HPLs in terms of Goncharov polylogarithms and compared to the results obtained by the differential equation method.
For the numerical checks we used a dozen points in the u − z plane. We first evaluated the Goncharov polylogarithms that appear in our analytic results numerically with the GiNaC-library [47,48]. We also derived Mellin-Barnes (MB) representations, partially using the AMBRE-package [49]. The analytic continuation to = 0 and subsequent numerical integration was carried out by MB.m [50]. This worked for almost all cases, even in the presence of kinematic thresholds, and yielded agreement to the GiNaC results to 5·10 −10 or better. There are, however, a few cases in which the Monte-Carlo integration implemented in MB.m failed due to highly oscillating integrands, notably for the integrals C 28−30 , and their "mass-flipped" counterparts (where m c ↔ m b and q 3 ↔ q 4 ). In these cases, we relied on the sector decomposition method implemented in SecDec [51,52], which yielded agreement with GiNaC at the level of 8·10 −7 for the highest -coefficients in C 28−30 , and at the level of 6 · 10 −4 for the highest -coefficients of their "mass-flipped" counterparts. The agreement is several orders of magnitude better for the lower coefficients in the -expansion.
Another important point to mention is the fact that the GiNaC results were obtained in the canonical basis, whereas most of the MB representations and the SecDec results were derived in an "ordinary" basis of un-dotted and singly-dotted master integrals. The change of basis was then performed using the Laporta reduction. Having calculated the numerics in two different integral bases constitutes another non-trivial check of our results.

Conclusion
We obtained analytic results to all two-loop master integrals that are necessary for the description of the non-leptonic decay B → Dπ at NNLO in QCD factorisation. They are expressed in terms of Goncharov polylogarithms of argument u and weights that are either integer numbers (0 or ±1), or contain the second kinematic varible, z. It is remarkable that six z-dependent weights are sufficient for writing down the entire set of solutions, including the "mass-flipped" integrals.
With the master integrals at hand, the bare two-loop part of the hard-scattering kernels T ij (u) in (1.1) is complete. The remaining task consists of renormalising the ultraviolet divergences and subtracting infrared divergences via matching from QCD onto soft-collinear effective theory. Steps towards this goal are outlined in [13]. Having the hard-scattering kernels T ij (u) written in terms of iterated integrals is an optimal choice for carrying out the convolution integral with the pion LCDA in (1.1), and it might be feasible to obtain the NNLO topological tree amplitude in analytic form. In any case our results constitute an important step towards the phenomenology of B → Dπ decays at NNLO in QCD factorisation.
Let us compare the integrals in the present work to those recently obtained in [11] during the evaluation of the two-loop penguin amplitude. Both are two-loop problems with scales u and z. The present integrals are a bit less involved compared to those in [11], in a sense that the linear combinations that form a canonical master integral are shorter, the occurring weights are fewer, and the choice of kinematic invariants is less complicated. The main reason for this is that in the present work the external kinematics of the final state contains also the second internal mass, notably m c . On the other hand, the only five-line integral in [11], a two-point function (M 22 ), is in fact a one-scale integral, whereas here we encountered several five-line integrals with four external legs which are genuine two-scale functions. Moreover, most of our integrals are needed to order O( 4 ), whereas in [11] all but two integrals were required only to order O( 3 ).
On more general grounds, it will be interesting to investigate how the canonical basis depends on the number of loops, legs, scales, space-time dimensions, and on the external kinematics. Every example therefore sharpens our understanding of the patterns that such bases follow, with the goal of eventually developing an algorithm for their automated construction.

Acknowledgments
We would like to thank A. Smirnov for useful correspondence on FIRE [36]. This work is supported by DFG research unit FOR 1873 "Quark Flavour Physics and Effective Field Theories".

A. The matricesÃ
Here we list the matricesÃ for the different topologies. Their entries can all be expressed in terms of the following nine logarithms, The matricesÃ now assume the following compact form, (A.9)