The Pentabox Master Integrals with the Simplified Differential Equations approach

We present the calculation of massless two-loop Master Integrals relevant to five-point amplitudes with one off-shell external leg and derive the complete set of planar Master Integrals with five on-mass-shell legs, that contribute to many $2\to 3$ amplitudes of interest at the LHC, as for instance three jet production, $\gamma, V, H +2$ jets etc., based on the Simplified Differential Equations approach.


Introduction
With LHC delivering collisions at the highest energy achieved so far, 13 TeV, experiments are analysing data corresponding to an integrated luminosity of 42 pb −1 [1] and 85 pb −1 [2], as well as those already collected at an energy of 8 TeV and an integrated luminosity of 20.3 fb −1 [3] and 19.7 fb −1 [4]. In order to keep up with the increasing experimental accuracy as more data is collected, more precise theoretical predictions and higher loop calculations are required [5].
In the last ten years our understanding of the reduction of one-loop amplitudes to a set of Master Integrals (MI) based on unitarity methods [6,7] and at the integrand level via the OPP method [8,9], has drastically changed the way one-loop calculations are preformed leading to many fully automated numerical tools (some reviews on the topic are [10,11]). In the recent years, a lot of progress has been made also towards the extension of these reduction methods for two-loop amplitudes at the integral [12][13][14] as well as the integrand [15][16][17][18] level. Contrary to the one-loop case, where MI have been known for a long time already [19], a complete library of MI at two-loops is still missing. At the moment this is the main obstacle to obtain a fully automated NNLO calculation framework similar to the one-loop one, that will satisfy the anticipated precision requirements at the LHC [20].
Following the work of [21][22][23], there has been a building consensus that the so-called Goncharov Polylogarithms (GPs) form a functional basis for many MI. A very successful method for calculating MI and expressing them in terms of GPs is the differential equations (DE) approach [24][25][26][27][28][29], which has been used in the past two decades to calculate various MI at two-loops [28,[30][31][32][33][34][35][36][37][38][39][40][41][42]. In [43] a variant of the traditional DE approach to MI was presented, which was coined the Simplified Differential Equations (SDE) approach. In this paper we present a further application of this method, concerning the calculation of planar massless MI relevant to five-point amplitudes with one off-shell leg, as well as the complete set of planar MI for five-point on-shell amplitudes. This is an important step towards the calculation of the full set of MI with up to eight internal propagators needed to realise a fully automated reduction scheme, à la OPP, for NNLO QCD.
Pentabox integrals are needed in particular in order to compute NNLO QCD corrections to several processes of interest at LHC [5]. The pp → H + 2jets can be used to measure the HW W coupling to a 5% accuracy with 300 fb −1 data. The pp → 3jets to study the ratio of 3−jet to 2−jet cross sections and measure the running of the strong coupling constant. The pp → V + 2jets for PDF determination and background studies for multi-jet final states.
The paper is organized as follows. In Section 2 we set the parameterization and notation of the variables describing the two-loop MI of interest. In Section 3 we discuss the DE obtained, and the results for the pentabox MI. We conclude in Section 4 and provide an overview of the topic and some perspective for future developments. In the Appendix A we present details on the derivation of the planar pentabox MI with on-shell legs and in the Appendix B we give a few characteristic examples on how the boundary conditions are properly reproduced in our approach by the DE. Finally in the ancillary files [44], we provide our analytic results for all two-loop MI in terms of Goncharov polylogarithms together with explicit numerical results.

The pentabox integrals
The MI in this paper will be calculated with the SDE approach [43]. Assume that one is interested in calculating an l−loop Feynman integral with external momenta {p j }, considered incoming, and internal propagators that are massless. Any l−loop Feynman integral can be then written as with matrices {c ij } and {d ij } determined by the topology and the momentum flow of the graph, and the denominators are defined in such a way that all scalar product invariants can be written as a linear combination of them. The exponents a i are integers and may be negative in order to accommodate irreducible numerators. Any integral G a 1 ···an may be written as a linear combination of a finite subset of such integrals, called Master Integrals, with coefficients depending on the independent scalar products, s ij = p i · p j , and space-time dimension d, by the use of integration by parts (IBP) identities [45][46][47]. In the traditional DE method, the Master Integrals are differentiated with respect to p i · ∂ ∂p j and the resulting integrals are reduced by IBP to give a linear system of first order DE [24,27]. The invariants, s ij , are then parametrised in terms of dimensionless variables, defined on a case by case basis, so that the resulting DE can be solved in terms of GPs. Usually boundary terms corresponding to the appropriate limits of the chosen parameters have to be calculated using for instance expansion by regions techniques [48,49].
SDE approach [43] is an attempt not only to simplify, but also to systematize, as much as possible, the derivation of the appropriate system of DE satisfied by the MI.  To this end the external incoming momenta are parametrized linearly in terms of x as the parameter x captures the off-shell-ness of the external leg. The class of Feynman integrals in (2.1) are now dependent on x through the external momenta: By introducing the dimensionless parameter x, the array of MI, G({s ij }, ; x), which now depends on x, satisfies a system of differential equations in one independent variable, where H is a matrix whose elements are rational functions of the kinematics {s ij ≡ p i · p j }, of x and of . Experience up to now shows that this simple parametrization can be used universally to deal with up to six kinematical scales involved, as is the case we will present in this paper. The expected benefit is that the integration of the DE naturally captures the expressibility of MI in terms of GPs and more importantly makes the problem independent of the number of kinematical scales (independent invariants) involved. Note that as x → 1, the original configuration of the loop integrals (2.1) is reproduced, which corresponds to a simpler one with one scale less.
More specifically we are interested in calculating the MI of two-loop QCD five-point amplitudes. As it is an inherent characteristic of the SDE method to interpolate among different kinematical configurations of the external momenta the starting point is to compute  Figure 1, that are contained in the family P 1 . All external momenta are incoming.
five-point amplitudes with one off-shell leg. These amplitudes contribute to the production i.e., of one massive final state V , plus two massless final states j 1 , j 2 at the LHC: The colliding partons have massless momenta q 1 , q 2 , while the outgoing massive and the two massless particles have momenta q 3 and q 4 , q 5 respectively. Of course, by appropriately taking the limit x = 1 the pentabox MI with all external massless momenta on-shell will be obtained, that are relevant for instance to the three-jet production For the off-shell case M 2 3 = 0, there are in total three families of planar MI whose members with the maximum amount of denominators, namely eight, are graphically shown in Figure 1. Similarly, there are five non-planar families of MI as given in Figure 2. We have checked that the other five-point integrals with one massive external leg are reducible to MI in one of these eight MI families. The two-loop planar ( Fig. 1) and non-planar ( Fig.  2) diagrams contributing to (2.4) have not been calculated yet. In fact by taking the limit x = 1 all planar graphs for massless on-shell external momenta are derived as well 1 . In this paper we calculate all MI in the family P 1 as well as all the the on-shell MI as M 2 3 → 0. We use the c++ implementation of the program FIRE [51] to perform the IBP reduction to the set of MI in P 1 .
The family P 1 contains in total 74 MI. Up to five denominators integrals with doubled propagators have been used as MI whereas starting from six denominators integrals with irreducible numerators are chosen. Many of the 74 MI already appear in the families of the double box integrals discussed in [31,32,36,[38][39][40][41][42]. However, there are seventeen new fivepoint Feynman diagrams that are not contained in the double box integral families. Three of them are pentaboxes, including the scalar and two MI with irreducible numerators. There are six seven-denominator, and eight six-denominator ones, the scalar members of which are shown in Figure 3.
For the family of integrals P 1 the external momenta are parametrized in x as shown in Figure 4. The parametrization is chosen such that the double box MI with two massive external legs that is contained in the family P 1 has exactly the same parametrization as that one chosen in [41], i.e. two massless external momenta xp 1 and xp 2 and two massive external momenta p 123 − xp 12 and −p 123 . The MI in the family P 1 are therefore a function of a parameter x and the following five invariants: where the notation p i···j = p i + · · · + p j is used and p 5 := −p 1234 . As the parameter x → 1, the external momentum q 3 becomes massless, such that our parametrization (2.6) also captures the on-shell case M 2 3 → 0. The MI depend in total on 6 variables, namely the Lorentz products q i .q j with i < j < 5 and the (squared) particle mass M 2 3 = q 2 3 . The x-parameterization of the external momenta as in Figure 4 results in these variables being related to the parameter x and the five independent scalar products of our choice that are defined in (2.6). The momenta q 1 , q 2 , q 4 and q 5 of the massless external particles can correspond to either of the four massless external legs in Figure 4, while the massive particle V has an external momentum q 3 = p 123 − xp 12 with a mass: After fixing the x-parameterization as in Figure 4, the class of loop integrals describing the planar familyP 1 is now explicitly expressed in x as: where γ E is the usual Euler-Mascheroni constant.
Using the notation given in Eq. (2.8), the indices a 1 · · · a 11 for the list of MI in the planar family P 1 is as follows 2 : In the next section we discuss the DE method that we use to calculate the above 74 MI in P 1 .

Differential equations and their solution
The resulting differential equation in matrix form can be written as where G stands for the array of the 74 MI given in Eq. We found that, after absorbing the integrating factors, the resulting matrix M can be written as where the coefficients C IJ;ijk andC IJ;jk are rational functions of the kinematics {s ij } and N IJ (ε) are rational function of ε . The twenty letters l i , are given explicitly by with ∆ 1 being the usual Gram determinant. The normalization factors N IJ (ε) can be cast in the factorized form N IJ (ε) = n J (ε) /n I (ε) and can be absorbed by redefining G I → n I (ε) G I . Although the DE can be solved starting from (3.2) 3 and the result can be expressed as a sum of GPs with argument x and weights given by the letters in Eq. (3.3), it is more elegant and easier to solve, if the system of differential equations is brought into a Fuchsian form [52], where only single poles in the variable x appear. In fact the series of successive transformations where the residue matrices M a are independent of x and ε. The rationale of the above transformations is of course to bring the system in its Fuchsian form. Starting from the matrix in (3.2) and after the absorption of the ε-dependent terms N IJ (ε), the matrix at order ε 0 , with the exception of the rows 69 and 74, only contains terms proportional to (x − l i ) −2 and x 0 , where l i are letters from (3.3). The transformation based on K 1 matrix, whose elements are proportional to (x − l i ) −1 and x 1 , obviously aims at removing these terms. As a matter of fact the new matrix, after the first transformation is applied, has, at order ε 0 , non-zero elements only in the rows 69 and 74, again proportional to (x − l i ) −2 and x 0 , with the exception of the row 74. Finally, after the second transformation, the whole matrix at ε 0 -order has only non-zero elements in the row 74 proportional to (x − l i ) −2 and x 0 , so by the application of K 3 transformation, the system is brought into the form (3.5). It should be noticed that the series of the above transformations do not correspond to the one described by the Moser algorithm [53][54][55][56].
Moreover, the final form does not only correspond to a Fuchsian form, but also factors out completely the ε dependence. Equation (3.5) can be straightforwardly solved and the result is given as 0 +b , the so-called alphabet, being rational functions of the independent variables { x}. Although this form looks very similar to our equation (3.5), with the identification α k (x) := x − l k , a comment is in order: usually in the canonical d-log form, DE are constructed with respect to all kinematical variables, and therefore the corresponding matrices A k should not depend on the kinematics. Contrary in our case, since the DE are derived with respect to only one variable, namely x, the matrices M a still depend on the (rest of the) kinematics {s ij }, defined in (2.6). As a consequence, in the multi-variable DE, the results are expressed in terms of GPs with coefficients that are rational numbers -commonly referred as universally transcendental (UT) -whereas in our case coefficients of GPs may still depend on {s ij }. In fact most of the elements of G in our solution do appear in a form that corresponds to the commonly defined UT form, with coefficients of GPs that are rational numbers independent of {s ij }, though others exhibit coefficients of GPs depending on {s ij }. In fact, working with DE in one independent variable, one can easily accommodate the dependence on the kinematics, as shown explicitly in Eq.(3.6), although it is an open question if the DE can be cast into a form so that all residue matrices M a are independent of the kinematics. In this context it is also interesting to mention that all the residue matrices M a , with a non-trivial dependence on the kinematics, do have eigenvalues that are negative integer numbers independent of the kinematics, related to the behaviour (x − l a ) −naε (n a positive integers) of the integrals at the corresponding limits x → l a .
The limit x = 1 represents the solution for all planar pentabox on-shell Feynman integrals. The limit can easily be obtained by properly resumming the log k (1 − x) terms. Interestingly enough we found a very simple formula for this limit given by with M 2 the residue matrix at x = 1 and G trunc derived from Eq.(3.6), by properly removing all divergencies proportional to log k (1 − x) and setting x = 1. Details can be found in Appendix A.
For the majority of the MI in the original basis (2.9), their boundary behaviour is captured by the DE itself, as it was also the case for the doublebox families [41]. This is achieved in a bottom-up approach. For several MI, such as the one-scale MI, including the genuine two-loop MI with three denominators, as well as for those MI that are expressed as product of one-loop integrals, their dependence on x is known in a closed form in terms proportional to x i+j . The rest of the MI, with a number m of denominators (m > 3), satisfy inhomogeneous DE, where the inhomogeneous part is completely fixed by MI with a number m of denominators, with m < m. Then, at each step of the iteration, we first calculate the part of the corresponding MI with m denominators, which we call the resummed part 5 , G (m) res , by introducing the following parametrization where i 0 and j are dictated by the DE itself in the limit x → 0. By putting the above expression (3.8) for the integrals in the DE and equating the coefficients of the terms x i+j with the same exponents on both sides of the DE in the limit x → 0, linear equations are obtained for the coefficients c i and d i . We note here that for the large majority of integrals we did not need to solve for the coefficients d i that correspond to x-suppressed terms. Their calculation were only required for those integrals whose DE had singularities of the form x −2+j at the boundary x = 0 (such singularities were also encountered for the DE of the one-loop pentagon discussed in [43]). Once the resummed terms in equation (3.8) were calculated, the resulting DE for G (m) res have no singularities at x = 0 and therefore can be safely expanded in ε and the result directly expressed in terms of Goncharov polylogarithms of the form G (l a , l b , . . . ; x).
For the pentabox family P 1 specifically, the majority of the coefficients (3.8) are fixed by the above-mentioned linear equations, while some others are not. We found in practice that for most of the integrals, the coefficients which are not fixed by the linear equations are zero and we confirmed this by the method of expansion by regions [48,49]. However, for some integrals we found that the method of expansion by regions predicts that some coefficients that are not determined by the linear equations, are in fact non-zero and require an explicit calculation. As described in the Appendix B, for those integrals we used other methods to calculate the unknown and nonvanishing coefficients. Once all boundary conditions were found for the integrals in the original basis (2.9), the boundary conditions for the canonical basis followed directly from the relation between the two bases that is described above.
It should be noticed that in the language of expansion by regions [48,49], MI in the limit x → 0, can be written, in general, as an expansion in terms proportional to x i+j log k (x), with k ≥ 0, which goes beyond our parametrization in (3.8). Nevertheless, we have explicitly checked, either by using the method of expansion by regions [57] or by comparing with previously known results, that in our case, for the pentabox family P 1 , no terms with k ≥ 1 were present, at the leading order in the limit x → 0.
The complete expressions for all MI are available in the ancillary files [44]. The solution for all 74 MI contains O(3, 000) GPs which is approximately six times more than the corresponding double-box with two off-shell legs planar MI. We have performed several numerical checks of all our calculations. The numerical results, also included in the ancillary files [44], have been obtained with the GiNaC library [58] and compared with those provided by the numerical code SecDec [59][60][61][62][63] in the Euclidean region for all MI and in the physical region whenever possible (due to CPU time limitations in using SecDec) and perfect agreement was found. For the physical region we are using the analytic continuation as described in [41]. At the present stage we are not setting a fully-fledged numerical implementation, which will be done when all families will be computed. Our experience with double-box computations show that using for instance HyperInt [64] to bring all GPs in their range of convergence, before evaluating them numerically by GiNaC, increases efficiency by two orders of magnitude. Moreover expressing GPs in terms of classical polylogarithms and Li 2,2 , could also reduce substantially the CPU time [42]. Based on the above we estimate that a target of O 10 2 − 10 3 milliseconds can be achieved.

Conclusions
In this paper we calculated, for the first time, one of the topologies of planar Master Integrals related to massless five-point amplitudes with one off-shell leg as well as the full set of massless planar Master Integrals for on-shell kinematics. We have demonstrated that based on the Simplified Differential Equations approach [43] these MI can be expressed in terms of Goncharov polylogarithms. The complexity of the resulting expressions is certainly promising that the project of computing all MI relevant to massless QCD, namely all eightdenominator MI with arbitrary configuration of the external momenta, is feasible. Having such a complete library of two-loop MI, the analog of A 0 , B 0 , C 0 , D 0 scalar integrals at one loop, the reduction of an arbitrary two-loop amplitude à la OPP can pave the road for a NNLO automation in the near future.
As experience shows, there are several issues that will need to attract our attention in order to accomplish our goal. First of all in order to systematize the whole procedure of reducing an arbitrary Feynman Integral in terms of MI in an efficient way, a deepening of our current understanding of IBP identities [65,66] is necessary [67]. Secondly, further standardising the procedure to obtain a canonical form of DE [29], which drastically simplifies the expression of MI in terms of GPs, is certainly a very desirable feature. Thirdly, the inclusion of MI with massive internal propagators, at a first stage with one mass scale corresponding to the heavy top quark, will provide the complete basis for NNLO QCD automated computations. Moreover, the calculation of boundary terms for the DE can benefit from further developments and exploitations of expansion-by-regions techniques, in conjunction with Mellin-Barnes representation of the resulting integrals. Finally, on the numerical side, a more efficient computation of polylogarithms is also necessary.
By equating the powers of log(1 − x) in the left-hand sides of (A.1) and (A.2) we obtain towers of equations for the coefficients X (n) and Y (n) . Consistency of these equations require the following identities to hold The above relation can be easily proved by mathematical induction. The first nontrivial relation corresponds to n = 3, namely −M 2 3 = 3M 2 2 + 2M 2 which is a direct consequence of the minimal polynomial x(x + 1)(x + 2) of the matrix M 2 . Equation (3.7) can now be easily derived from (A.2) by taking the limit x = 1 with G trunc ≡ G reg (x = 1). In other words, for the above integrals there is some behaviour at x → 0 which corresponds to coefficients c i and d i in equation (3.8) that are not all determined by the DE itself. For the above integrals we used various methods to compute the undetermined coefficients, which we explain further below.

B Methods of calculating the boundary conditions
Expansion by regions For the following integrals the method of expansion by regions was used to fix the non-zero coefficients in the boundary behaviour With the method of expansion by regions [48,49], the coefficients in equation (3.8) are expressed as Feynman integrals. For most of the integrals in (B.2), these integrals corresponding to the undetermined (by the DE) coefficients are reducible to known single scale integrals. However for a few, the integrals for the coefficients are non-trivial and require further calculation. For example let's consider the following MI and its DE expanded around x = 0: where B(s, ) is independent of the integral parameter x and the dots do not contain any further poles at x = 0. From expansion by regions it follows that at x = 0 the integral behaves as G ∼ c 1 x −1−4 + c 2 x −1−2 which corresponds to a soft and collinear region. The coefficient c 2 is completely determined by the DE (B.3) as explained in section 3, while c 1 is not. The coefficient c 1 is given by the method of regions as a Feynman integral:

(B.4)
We can calculate the above integral by the method of DE. As we did for the x-parametrization of the P 1 family, we introduce a parameter y in the denominators: (B.5) By solving the DE ∂ y I(y) and afterwards setting y → 1 we could calculate c 1 = I(1). We note that the boundary condition for the integral I(y) is completely determined by its DE. If this was not the case, we would use the method of regions again to express the missing coefficient in terms of a Feynman integral and repeat the above step, though this was never necessary.
Shifted boundary point For some other integrals we considered their limiting behaviour around another boundary point x 0 instead of at x 0 = 0: In practice, choosing the boundary at x 0 = ∞ corresponds to performing an inverse transformation x → 1/x at the level of the DE and afterwards considering x 0 = 0. This case was already discussed in [41] and therefore here we discuss the case of the MI G := G 01000001011 . Its DE is of the form: (1 − x ) and then afterwards integrating from x = 0 one can compute G as a function of x . Upon tranforming back to x, the limiting behaviour of G at x = 0 and therefore its corresponding boundary condition is found.
Extraction from known integrals The last three integrals in (B.1) correspond to taking the x → 1 limit of other known integrals (that lie in the same family P 1 ) and afterwards redefining the invariants as follows: The three integrals on the right hand side of equation (B.8) are MI that were previously calculated in [41]. From the exact result in x for the three integrals on the left hand side in (B.8) we could then compute their corresponding non-zero coefficients in (3.8) that are not determined 6 by the DE.