Master Integrals for double real radiation emission in heavy-to-light quark decay

We evaluate analytically the master integrals for double real radiation emission in the b -->u W* decay, where b and u are a massive and massless quark, respectively, while W* is an off-shell charged weak boson. Since the W boson can subsequently decay in a lepton anti-neutrino pair, the results of the present paper constitute a further step toward a fully analytic computation of differential distributions for the semileptonic decay of a b quark at NNLO in QCD. The latter partonic process plays a crucial role in the study of inclusive semileptonic charmless decays of B mesons. Our results are expressed in terms of multiple polylogarithms of maximum weight four.

1 Introduction The study of inclusive semileptonic B meson decays is important for the determination of the Cabibbo-Kobayashi-Maskawa (CKM) matrix elements |V ub | and |V cb | and, therefore, it constitutes a stringent test on the unitarity of the CKM matrix. Inclusive determinations of |V ub | rely on an Operator Product Expansion (OPE) [1][2][3][4], according to which the total B meson semileptonic decay rate and various kinematic distributions can be described, at leading order in a power expansion with respect to the inverse b-quark mass, in terms of the partonic decay rate of an on-shell b quark into a lepton-neutrino pair and a u quark. Within this framework, theoretical predictions for the partonic decay rates are obtained by means of perturbation theory. Phenomenological predictions for observables related to the semileptonic B-meson decays are then obtained by combining perturbative calculations for the semileptonic b-quark decays with a finite number of non-perturbative parameters.
However, the measurements of the B → X u eν decay are affected by large backgrounds due to the B → X c eν decay. In order to suppress this background, experiments impose sharp cuts (for example, cuts on the final state hadronic invariant mass). This in turn leads to problems with the convergence of the OPE in theoretical predictions.
These issues can be addressed by parameterizing the residual motion of a b-quark in the B meson by means of the shape function [5][6][7]. Further studies [8,9] showed that with a combination of cuts on the hadronic and leptonic invariant masses the impact of the shape function can be suppressed and the OPE can be used to describe the B → X u eν. In both approaches the QCD corrections to the partonic process b → uW * (where W * indicates an off shell W boson and the mass of the up quark is set to zero) play a crucial role.
An analytic result for the next-to-leading order (NLO) QCD corrections to the b → lν l u differential decay rate was obtained in [10]. A fully differential calculation of this decay rate at next-to-next-to-leading order (NNLO) was carried out in [11,12]. A related result for the differential top-quark semileptonic decay at NNLO was presented in [13]. These NNLO calculations are based on numerical techniques. Analytic results at NNLO are also known; however so far these studies were carried out in the shapefunction region by using Soft-Collinear Effective Theory [14][15][16]. The contribution of these corrections to |V ub | were then considered in [17].
In this paper, we focus on the analytic calculation of the Master Integrals (MIs) necessary for the determination of the contribution of the QCD double-real radiation to the b → uW * decay. This process is one of the three elements that are necessary for the evaluation of the triple-differential distribution in the charged lepton energy, leptonic invariant mass and final-state hadronic invariant mass at NNLO in QCD. For what concerns the other two elements, the two-loop corrections to the b → uW * decay were evaluated analytically in [14][15][16]18], while the one-loop real-virtual contribution will be the subject of future work.
In order to carry out the calculation, we employed the method of reverse unitarity. This approach was introduced in the context of the evaluation of the NNLO QCD corrections to the production of a Higgs boson in gluon fusion [19] and then applied to several other processes (see for example [20][21][22][23][24][25]). The method consists in applying Cutkosky rules [26] in order to map the calculation of the interference between two leading order (LO) 2 → 3 diagrams integrated over the final state phase-space into the evaluation of "cuts" of two-loop 2 → 2 diagrams. In this way, the on-shell condition for the real particles in the final state of the 2 → 3 process is converted into the difference of two propagators with opposite i0 + prescription. Subsequently, one can calculate the cut diagrams by means of techniques for the analytic evaluation of multi-loop diagrams which were developed starting from the late '90s. In particular, dimensionally regularized scalar integrals are reduced to MIs by using Integration by-Parts Identities (IBPs) [27][28][29]. Reduction algorithms based on IBPs are implemented in publicly available computer programs [30][31][32][33][34][35][36][37][38]. The MIs are analytically evaluated by means of the Differential Equations (DE) method [39][40][41][42][43] and expressed in terms of generalized polylogarithms (GPLs) [44][45][46] of two variables: t, which is connected to the invariant mass of the hadronic final state, and z, which is related to the leptonic invariant mass. The evaluation of the expansion of the MIs was carried out up to terms that include GPLs of maximum weight four.
Our results are also relevant for the determination of the total width of a top quark that decays into a massless bottom quark and a lepton-neutrino pair [13].
The paper is structured as follows. In Section 2, we discuss the calculation, by introducing the notation and the kinematics of the process, and by identifying the MIs and the range of validity of their analytic expressions. In Section 3, we present the results for the MIs. Section 4 contains our conclusions. Appendix A collects the explicit expressions of the -poles of the MIs.
The analytic expressions of the twelve MIs evaluated in this work are collected in an ancillary file included in the arXiv submission.

Calculation
The calculation of double emission corrections to the b → uW * process is first mapped into the problem of calculating three-particle cuts in two-loop bW * → bW * forward box-diagrams, using the method proposed in [19]. Three auxiliary topologies which encompass all of the combinations of denominators which can appear in the cut diagrams were subsequently identified. The MIs belonging to each topology were identified by using IBPs as implemented in LiteRed [31,32]. The MIs, which depend on two dimensionless parameters (defined below), are then calculated by employing the DE method. The technique employed in this work is by now a standard method in the analytic calculation of Feynman diagrams. In this section we describe the way in which we parameterized the kinematics of the process, we define the MIs which were identified and we discuss the way in which the integration constants arising in the DE method were fixed.

Kinematics
At tree-level, the kinematics of the decay we are interested in is where p 1 , p 2 and q are the four-momenta of the bottom quark, up quark and W boson, respectively. Consequently, where m b indicates the mass of the bottom quark and the mass of the up quark is neglected. The W boson is off-shell. At tree level the energy of the up quark in the bottom-quark rest frame is smaller that m b /2. Therefore, we introduce the dimensionless parameter z defined as Beyond leading order, it is necessary to consider the process where X indicates an inclusive state involving light quarks and gluons, so that p 2 X = 0 in general. As it was done in [10], the invariant mass of the state X is parameterized by introducing the variable t defined through the relation where z remains defined as in Eq. (2.3), provided that p 2 is replaced by p X . Tree-level kinematics, i.e. p 2 X = 0 with 0 ≤ z ≤ 1, is recovered in the t → 1 limit. Beyond tree level, the available phase space in the {p 2 , z} and {t, z} plane is shown in Figure 1. The physical region in the {p 2 , z} plane is shown in the left panel of Figure 1 and it is delimited by the conditions In the {t, z} plane the physical region (shown in the right panel of Figure 1) is given by In the calculation of the MIs which we carry out in this paper, we keep t = 1. In fact, the differential distribution which we ultimately want to obtain by employing the integrals which we evaluate here is divergent for t → 1 and includes "star" distributions of the form 1 [10] ln np2 However, it is sufficient to calculate the MIs for real radiation corrections by keeping t = 1 since the expression of the partonic double differential distribution in the t → 1 limit was evaluated by using Soft Collinear Effective Theory (SCET) in [17]. The partonic real radiation corrections in the t → 1 (p 2 → 0) limit were recalculated also by us. These corrections, which receive contributions from double real and real virtual diagrams, involve poles in and finite terms proportional to δ(p 2 ) as well as the terms proportional to the star distributions of Eq. (2.1). If one neglects the star distribution bracket, the star distributions arising in the SCET calculation must match exactly the singular terms arising from the calculation carried out at t = 1. The latter depends on the MIs evaluated here as well as on the two cut MIs which will be the subject of a future work. Therefore, at the stage in which the differential distribution will be assembled, the SCET calculation of the real radiation will serve as a further cross check of the calculation carried out at t = 1. At the same time, it will provide the complete contribution of the real corrections proportional to δ(p 2 ). We already tested with success this procedure by recalculating at NLO the differential distribution originally derived in [10].

Auxiliary topologies and Master Integrals
All of the Feynman diagrams contributing to the double emission corrections to the b → uW * process can be calculated once the scalar integrals belonging to the (three cuts of the) three auxiliary topologies are known.
1 Star distributions can be defined through the relation where f is a smooth test function.
The integrals belonging to the first auxiliary topology considered, which is referred to as topology A, are defined as follows The seven propagators P i in Eq. (2.8) are where the last three propagators in the list are cut propagators. For example (2.10) Equivalent relations hold for P 6 and P 7 . As a consequence of the presence of cut propagators, the integrals I A in Eq. (2.8) are zero when at least one among the powers α 5 , α 6 , α 7 is zero or negative. The eight MIs belonging to topology A are shown in Figure 2. We introduce the following notation in order to label the MIs of topology A: The integrals belonging to topology B are defined in analogy to Eqs. (2.8,2.9): with where again Q 5 , Q 6 , Q 7 are cut propagators. Topology B involves eleven MIs, which include the eight MIs already needed for topology A, plus the three non planar MIs shown in Figure 3. The latter are labeled as follows I 9 ≡ I B (0, 1, 1, 0, 1, 1, 1) , I 10 ≡ I B (0, 2, 1, 0, 1, 1, 1) , (1, 0, 1, 1, 1, 1, 1) . (2.14) Topology C is defined by the integrals Topology C involves five MIs: I 1 and I 4 , which are present also in topologies A and B, I 9 and I 10 , which are already needed for topology B and one additional non-planar integral which is shown in Figure 4.

Differential equations
Each MI satisfies differential equations with respect to the two dimensionless parameters z and t which we use in order to parameterize the phase space. The differential equations can be derived starting from the IBPs and are directly obtained from LiteRed. Only the integrals of topology A involving the propagators P 2 , P 3 , P 5 , P 6 , P 7 are found to involve more than two MIs. Consequently, since the solutions of a system of three first-order differential equations cannot be found with standard methods, the evaluation of the three MIs with the five propagators listed above could be problematic. However, the integrals I 5 , I 6 and I 7 , which we chose as the MIs for this set of propagators, satisfy a system of three differential equations in t and y which decouple order by order in . For this reason, in the case at hand it was possible to evaluate the twelve MIs in the problem without employing a canonical basis [47]. In order to eliminate square-root weights in the GPLs which appear in the solution of differential equations satisfied by the MIs, we traded the variable z with the variable y defined through the relation For 0 ≤ z ≤ 2 the variable y is a pure phase which we choose to parameterize as The differential equations with respect to t and y satisfied by the MIs are solved order by order in using iterated integration. The solutions depend on several integration constants. Most of these constants can be fixed by imposing the regularity of the Numerical Constants MIs in the t → 0 limit. However, a subset of seven constants is left undetermined once the regularity in t → 0 has been required. MIs are in general not regular in the t → 1 limit; indeed, one expects the differential distribution to be singular in the tree-level kinematic limit. The singular part of the distribution in the tree-level limit was evaluated by using SCET. However, all of the MIs which are finite in (i.e. I 1 , I 2 , I 4 , I 5 , I 9 ) vanish in the t → 1 limit. In particular, the behavior of I 2 and I 5 in the tree level limit is sufficient to overconstrain the seven remaining constants. The analytic expression of the MIs in terms of GPLs of argument y and t are very long but they can be evaluated to arbitrary precision by means of the GiNaC routines of [48]. Therefore, six of the seven constants were fixed by requiring that integrals I 2 and I 5 vanish in the t → 1 limit. They are given in numeric form in  Table 2. Values of the constants determined by evaluating the analytic result for the MIs with GiNaC. The condition imposed in order to fix the constants is that integrals I 2 and I 5 vanish in the t → 1 limit. The second column indicates the MIs used to fix a given constant and the power n at which the constant first appears as a cofactor of n . In the last column one can find the complex value of the constant, which can in principle be determined at arbitrary precision.
digits. In general, we evaluated MIs up to the order in the expansion where GPLs of weight four first appear in the result, since one does not expect GPLs of weight five to be present in the NNLO differential distributions we are ultimately interested in. (Table 1 summarizes the order in at which the various MIs were evaluated.) A notable exception is represented by the MI I 5 . Indeed, one of the integration constants which appears alongside GPLs of weight four in I 7 , appears only at order 2 in I 5 . The term of 2 in I 5 also involves GPLs of weight five. This is not surprising since the set of MIs that was chosen does not have uniform transcendentality. Therefore, we evaluated I 5 up to order 2 and required that it vanishes in the t → 1 limit in order to fix this residual constant. In the following section we present analytic results for the MIs.

Alphabet
The analytic expressions for the MIs which we evaluated are written in terms of Harmonic Polylogarithms [46] and (two-dimensional) GPLs [48][49][50][51] of arguments t and y. GPLs can be defined recursively, for n ≥ 0, via the iterated integral G(a 1 , . . . , a n ; z) = z 0 dt t − a 1 G(a 2 , . . . , a n ; t) , for a generic argument z and weights {a 1 , · · · , a n }, assuming G(z) = G(; z) = 1. The case in which a i = 0 for all i needs to be considered separately: The weights of the GPLs or argument t can depend on y. For convenience we define the combinations Consequently, w 2 is real while w 1 is imaginary. The list of weights, i.e. the "alphabet", appearing in the GPLs of argument t is

Integrals
The simplest MI in the list is the two-loop phase space diagram shown in the top left panel of Figure 2, for which one can find an expression which is exact in d = 4 − 2 dimensions [52] . (3.7) The first two terms of the expansion ofÎ 1 read (3.8) In analogy with Eq (3.6), we introduce the notation The analytic expression of allÎ i up to terms involving GPLs of weight four can be found in the ancillary file MasterIntegrals.txt. We cross-checked nine out of the twelve MIs which we calculated analytically in this work by comparing the numerical evaluation of their analytic expressions (carried out by means of GiNaC) to the direct numerical integration of the MIs, carried out with the package SecDec [53][54][55]. We found agreement within the SecDec numerical integration error in all points tested. The remaining three MIs (I 8 , I 11 , I 12 ) belong to subtopologies which also admit at least one two-particle cut. Consequently, these MIs cannot be evaluated directly by comparing them to the imaginary part of a 2 → 2 forward box calculated with SecDec, since that imaginary part receives a contribution from twoparticle and three-particle cuts. Once the MIs corresponding to the two-particle cuts are known, they can be combined with the three-particle cuts calculated here and the sum of two-and three-particle cuts can finally be compared with the imaginary part of the corresponding 2 → 2 box integrals.
Finally, we employed GiNaC to evaluate the analytic results we found in all of the phase space. In Figures 5 and 6 we plot, for all of the integralsÎ i , the order in at which GPLs of order four appear first.

Conclusions
In this paper we evaluated analytically the MIs needed for the calculation of the double emission corrections to the b → uW * decay at tree level. The problem was mapped into the calculation of three-particle cuts in two-loop bW * → bW * forward box diagrams. We identified a set of twelve MIs belonging to three auxiliary topologies by using IBPs. The MIs depend on two dimensionless parameters. Their analytic expression in terms of GPLs was found by means of the DE method. The integrals can be evaluated with arbitrary numerical precision by means of the GPLs functions implemented in GiNaC. The complete analytic expression of all of the MIs can be found in an ancillary file included in the arXiv version of this paper. The result were cross checked against a direct numerical integration of the MIs carried out by means of the program SecDec. The results obtained here are needed for the analytic evaluation of the b → uW * decay double differential distribution to NNLO in QCD. Since the two-loop virtual corrections to the b → uW * decay are already known analytically [14][15][16]18], the only missing element at this stage are the one-loop, single-emission diagrams contributing to b → uW * decay to NNLO.

A MIs Poles
In this Appendix we collect the explicit expression of the poles of the MIs which are divergent in the → 0 limit. int 5 order ϵ int 6 order ϵ 2 Figure 5. Plots of the first order in which contains GPLs of weight four for the MIŝ I 1 , · · · ,Î 6 . int 7 order ϵ 3 int 8 order ϵ int 9 order ϵ 2 int 10 order ϵ 2 int 11 order ϵ int 12 order ϵ Figure 6. Plots of the first order in which contains GPLs of weight four for the MIŝ I 7 , · · · ,Î 12 .