Differential equations from unitarity cuts: nonplanar hexa-box integrals

We compute ϵ-factorized differential equations for all dimensionally-regularized integrals of the nonplanar hexa-box topology, which contribute for instance to 2-loop 5-point QCD amplitudes. A full set of pure integrals is presented. For 5-point planar topologies, Gram determinants which vanish in 4 dimensions are used to build compact expressions for pure integrals. Using unitarity cuts and computational algebraic geometry, we obtain a compact IBP system which can be solved in 8 hours on a single CPU core, overcoming a major bottleneck for deriving the differential equations. Alternatively, assuming prior knowledge of the alphabet of the nonplanar hexa-box, we reconstruct analytic differential equations from 30 numerical phase-space points, making the computation almost trivial with current techniques. We solve the differential equations to obtain the values of the master integrals at the symbol level. Full results for the differential equations and solutions are included as supplementary material.


Introduction
In the coming years, Run 2 of the LHC and the future high-luminosity LHC will accumulate large datasets and deliver high-precision experimental measurements. These must be matched by high-precision theoretical predictions for the relevant cross sections, at NNLO for a number of observables. A fundamental ingredient in obtaining these predictions is the evaluation of two-loop Feynman integrals with several external legs. Feynman integrals are also crucial for studying scattering amplitudes in supersymmetric theories, where there has been great progress in constructing loop integrands but much less is known about integration, especially beyond the planar limit. For instance, while the nonplanar four-particle amplitude is known at three-loop order [1], at five-points only the integrand is known at two-loops [2,3].

JHEP01(2019)006
A problem that has recently attracted much attention is the calculation of massless two-loop five-point amplitudes. The two-loop five-point amplitudes for pure Yang-Mills theory have been recently computed in the planar limit, first for the special case of all-plus external helicities [4][5][6] and subsequently for arbitrary helicities [7,8]. The calculation of loop amplitudes requires both the evaluation of the integrand, which for the five-point amplitudes mentioned above has often been performed in the framework of generalized unitarity [9][10][11][12], and the evaluation of Feynman integrals. An important property of dimensional regularization is that total derivatives with respect to the loop momenta of Feynman integrands in momentum space integrate to zero (boundary terms always vanish). This gives rise to powerful integration-by-parts (IBP) identities [13][14][15], which relate integrals with a specific set of propagators but with arbitrary numerators and different integer exponents of propagators to a basis of so-called master integrals. Using these relations, an amplitude can be reduced to a linear combination of master integrals that must then be evaluated. All master integrals required for planar massless five-point two-loop amplitudes have been evaluated analytically [16,17]. Beyond the planar limit, the amplitudes are as yet unknown and require in particular the evaluation of nonplanar five-point two-loop integrals. Some master integrals have been recently evaluated using bootstrap and (super-) conformal symmetry methods [18,19], with additional hidden symmetry properties explored in [20,21], but full results for the nonplanar hexa-box and double-pentagon topologies are still missing.
A very successful method for evaluating dimensionally-regularized master integrals is through the solution of differential equations with respect to kinematic invariants [22][23][24][25][26][27]. IBP relations also play a central role in this method. Indeed, the differential equations are obtained by differentiating Feynman integrals, and rewriting the resulting expression in terms of the basis of master integrals one is trying to compute. This last step is done with the help of IBP relations. The set of IBP relations that is needed for differential equations is in general simpler to obtain when compared to the one required for the evaluation of an amplitude because the integrals that must be reduced are simpler (indeed, only recently have IBP-reduction tables for planar five-point massless amplitudes been obtained [28,29]). Despite this, for nonplanar five-point two-loop integrals the IBP-reduction step is still the bottleneck in computing the full set of master integrals from the differential-equation approach.
In this paper, we present a method for constructing the differential equations which overcomes this issue. The main result is an approach to obtain a more compact IBP system than the one in used in more conventional methods [30][31][32][33][34] which is thus easier to solve. It builds on the recent observation that, in conjunction with computational algebraic geometry tools, unitarity-based methods can be applied to construct IBP relations [8,[35][36][37][38][39][40][41] and differential equations [42][43][44]. Our approach is based on the construction of unitarity-compatible IBP relations [8,[35][36][37][38][39][40][41], which allow us to control the powers of the propagators in the integrals appearing when generating the IBP relations. It is made even more efficient by working on specific unitarity cuts and then combining the results obtained on each of the cuts, in a procedure similar to the one used in [37,41]. As an example of our approach, we compute the differential equations for the nonplanar hexa-box topology. The differential equations for the maximal-cut integrals of this topology were previously derived by one of the authors using unitarity-based methods [43]. Here we present results for the differential equations of the uncut integrals.

JHEP01(2019)006
Although any basis of master integrals is in principle equivalent, it is well known that a basis of so-called pure master integrals is more convenient for solving the differential equations, as the dependence on the dimensional regularization parameter = (4 − d)/2 completely factorizes [27]. We thus construct such a basis and use it to write the differential equations. Explicit expressions were known in the literature for most of the pure master integrals we required, except for some 5-point examples. 1 Pure integrals are often constructed from the analysis of the integrand in strictly 4 dimensions. In this paper we observe that numerators constructed from specific Gram determinants that vanish in 4 dimensions can also be used. They lead to very compact expressions for the numerators, and we give explicit expressions for all planar five-point topologies, even beyond the ones required for the present calculation. The differential equations in our pure basis are given in the supplementary material diffEqMatrices.m.
We also propose an alternative and much more direct way of obtaining the differential equations of the nonplanar hexa-box, which relies on the knowledge of its so-called alphabet [45] and of a pure basis. In this approach, the analytic forms of the differential equations are highly constrained, leaving a small number of free parameters that are fixed by numerical fitting, with IBP reduction computed at numerical values of the kinematic invariants. This makes the computational resources required trivial.
To validate our results, we compute the solution of the differential equations at the symbol [45] level by constructing a generic solution, imposing the so-called first-entry condition [46] on the symbol, and finally computing a single trivial integral that allows to determine all initial conditions required at the symbol level.
The paper is structured as follows. In section 2, we describe our approach to generating the unitarity-compatible IBP relations required for computing differential equations, and the application to the case of the nonplanar hexa-box. In section 3, we discuss the use of the IBP reduction results to obtain the differential equations in dimensional regularization, the transformation into a pure basis, how we construct the pure basis, and our alternative method of constructing the differential equations. In section 4 we discuss our results for the differential equations of the nonplanar hexa-box and their solution at symbol level. In appendix A we describe some further details of our implementation of IBP reduction, namely the generation of IBP identities sector-by-sector. In appendix B we discuss an alternative set of variables leading to a rational alphabet. Finally, in appendix C we tabulate all the 73 pure integrals in our pure basis.
2 Unitarity-compatible integration-by-parts reduction A crucial step in constructing differential equations for a Feynman integral is obtaining the associated set of integration-by-parts (IBP) relations. These are required to rewrite the derivatives of the integrals in terms of a set of master integrals. This crucial step is currently the main bottleneck in calculating new Feynman integrals through the method of differential equations. Figure 1. The nonplanar hexa-box topology. We consider the case with all external and internal lines massless.

JHEP01(2019)006
In this paper we describe an approach to tackle this issue. The improvement we present is based on a more efficient construction of the IBP relations. We focus on the application of the approach to the construction of a system of differential equations for the nonplanar hexa-box topology shown in figure 1, which has so far not been achieved through more standard techniques. This serves as an illustration of the potential of our approach, which is completely generic and applicable beyond this example.

Loop momentum parametrization
Let us begin by defining our conventions to describe the nonplanar hexa-box topology. We choose the five independent kinematic invariants to be s 12 , s 23 , s 34 , s 45 , s 51 , We always set s 45 = 1 (2.3) in our calculations to eliminate one kinematic scale, which amounts to normalizing all variables to an overall scale. The analytic dependence on this overall scale can always be recovered from dimensional analysis. With the conventions of figure 1, the inverse propagators are written as There are three 'irreducible numerators' that cannot be written as a linear combination of inverse propagators. We choose to define them as ρ 9 = (l 2 + p 1 ) 2 , ρ 10 = (l 2 + p 1 + p 2 ) 2 , ρ 11 = (l 2 + p 1 + p 2 + p 3 ) 2 . (2.5)

JHEP01(2019)006
Our approach to the construction of the IBPs is inspired by generalized unitarity and thus heavily relies on putting some of the propagators 'on-shell'. In the context of this paper, putting a propagator 1/ρ on-shell simply means we set its inverse ρ to zero when it appears in the numerator. It is thus convenient to parametrize the loop momenta in terms of the inverse propagators and irreducible numerators [36,37,[47][48][49][50] ρ 1 , ρ 2 , . . . , ρ 11 . (2.6) We shall call these natural variables (in recent literature, they are also known as Baikov variables). There is an invertible linear map between the natural variables and the following eleven dot product, which we call the canonical variables, Therefore, any expression that is a polynomial in the canonical variables, such as the numerator of a loop integrand, is also a polynomial in the natural variables. It then becomes trivial to impose unitarity cuts on such an expression. For example, the maximal cut of a numerator of the nonplanar hexa-box topology is simply computed by setting ρ 1 = ρ 2 = · · · = ρ 8 = 0, leaving a polynomial in the three irreducible numerators.

IBP relations with controlled propagator powers
In dimensional regularization, integrals of total derivatives with respect to the loop momenta of a Feynman integral in the momentum-space representation do not have boundary terms. This leads to the so-called IBP relations, which relate Feynman integrals with different powers of propagators. For instance, for the nonplanar hexa-box integrals, IBP relations are obtained from for a generic vector v µ i , where there is an implicit summation over the loop momentum label i = 1, 2. We can expand v µ i as where the implicit summation ranges are j = 1, 2 and A = 1, 2, 3, 4. Therefore the v µ i have a total of 2 × (2 + 4) = 12 coefficients v ij and v iA . We require all 12 coefficients to be polynomials in the loop-momentum components such that the object we obtain after computing the derivative is still a Feynman integral. For generic vectors v µ i , the action of the derivative will naturally lead to propagators with increased power (squared propagators in the case of eq. (2.8)). Therefore, conventional implementations of integration-by-parts relations [30][31][32][33][34] involve 'auxiliary integrals' with raised propagators powers. Although one is a priori not interested in the reduction of these auxiliary integrals, in conventional approaches one is forced to also reduce them in order to have a complete IBP system. This leads to a proliferation of terms and can make solving the IBP system a very complicated task.

JHEP01(2019)006
A method that avoids the appearance of auxiliary integrals was proposed in ref. [35]. It is based on the observation that these will not be generated if one demands that the where f j are arbitrary polynomials in the natural variables. We will call vectors v µ i which solve the above equation IBP-generating vectors. We stress that the above condition applies to the inverse propagators but not to the irreducible numerators, i.e., in our example it applies to ρ 1 , ρ 2 , . . . , ρ 8 but not to ρ 9 , ρ 10 , ρ 11 . It is clear that, when using such an IBPgenerating vector in eq. (2.8), no terms with squared propagators will be generated on the right-hand side of the equation, and one thus obtains a relation between integrals with unit powers of the propagators.
Note that this procedure also allows us to deal with integrals that have squared propagators (or, in fact, a propagator raised to an arbitrary positive power). Indeed, consider the example If the vectors v µ i satisfy the condition in eq. (2.10), then the right-hand side will involve integrals with the propagator 1/ρ k raised to power two or less, and the remaining propagators raised to power one or less. Also in this case, we find that the IBP-generating vectors allow to control the proliferation of auxiliary integrals. We note that this case is important if one wants to use IBP relations for constructing differential equations, as integrals with squared propagators can be generated when differentiating an integral with respect to an external kinematic invariant.
In summary, using IBP-generating vectors, we are able to generate IBP relations in which we have control over the power of the propagators of the integrals. In particular, we never generate integrals with higher powers of propagators. IBP relations obtained in this approach are called unitarity-compatible IBP relations. These IBPs lead to a more compact IBP system that is thus easier to solve. 2 The final question we must address is how we construct the IBP-generating vectors. A full set of solutions to eq. (2.10) can be obtained from computational algebraic geometry [8,35,37,39,41,51], computational linear algebra [52], and in some cases, analytic formulas from loop-by-loop considerations [36] and dual conformal symmetry [53]. In this paper, we adopt the method of ref. [8], which directly solves eq. (2.10) after rewriting the lefthand side in terms of natural variables, using algorithms implemented in the computational algebraic geometry package SINGULAR [54]. In addition, we identify and remove redundant IBP-generating vectors found by SINGULAR, by reducing the vectors against each other at random numerical values of the kinematic invariants s ij , modulo an arbitrary large prime number. We only retain vectors that are not reduced to zero by other vectors. This ad hoc procedure, similar to the one in ref. [41] but without Groebner-basis computations, significantly reduces the number of IBP-generating vectors.

Simplification on unitarity cuts
In this subsection, we discuss IBP reduction on a spanning set of non-maximal cuts before full results are obtained from merging the results computed on specific cuts, in a procedure similar to the one used in refs. [37,41]. This approach further simplifies the construction of the IBP relations.
We first introduce a compact notation for integrals with generic powers of propagators and numerators (specialized to the nonplanar hexa-box topology), Here, the ν i denote the powers of propagators and numerators. We include a standard overall normalisation factor that depends on the dimensional regulator (we take d = 4 − 2 ) and the Euler-Mascheroni constant γ E . The indices ν 9 , ν 10 , ν 11 are required to be non-positive, since the three irreducible numerators ρ 9 , ρ 10 , ρ 11 never appear in the denominator with positive powers. On the other hand, the indices ν 1 , ν 2 , . . . , ν 8 are allowed to be any integers. Indeed, we include integrals of sub-topologies with a proper subset of the propagators where the missing inverse propagators can appear in the numerator with positive powers. Every IBP relation is a linear relation between the different F [ ν] integrals in eq. (2.12) with different sets of ν i indices.
In the following, we will consider the generalized unitarity cut of a subset ∆ of the propagators of a given For the purpose of this paper, it is sufficient to note that (2.14) In other words, the cutting operation sets to zero any integral that does not have all cut propagators raised to nonzero positive powers. As defined above, the cutting operation is a linear operation. If one has constructed an IBP relation for uncut integrals, one can then apply a cut to the relation and obtain the IBP relation between the cut integrals. Since some integrals are set to zero by the cutting operation, the cut relation will be simpler (i.e., involving fewer integrals) than its uncut counterpart. One thus expects that it should be simpler to compute IBP relations on a given cut, and this is indeed true. In particular, it is simpler to construct IBP-generating vectors by solving eq. (2.10) subject to cut conditions. One thus obtains IBP relations modulo integrals that do not depend on the cut propagators.
To choose which cuts to consider, we construct all subsets of propagators and check if, when removing any extra propagator, the corresponding Feynman integral becomes scaleless. If that is the case, we keep the cut associated with that subset. These cuts form a spanning set of cuts, and an explicit set for the nonplanar hexa-box will be given in the next section. In practice, working on a cut means we have to construct IBP-generating vectors that depend on fewer variables, which can greatly simplify their calculation using computational algebraic geometry.

JHEP01(2019)006
Having constructed an IBP relation on a given cut, there are two types of terms we need to worry about not capturing. The first type corresponds to integrals that only have a subset of the cut propagators. Given the way we chose the cuts, these are scaleless, vanishing integrals, and thus they do not modify the IBP relation. The second type of terms that we miss on a given cut are terms with a subset of the cut propagators and some other propagators. These terms will be captured by another cut in the spanning set. By merging the information on the spanning set of cuts, we construct the full IBP-reduction table which reduces any uncut integral to a set of master integrals. The coefficient of a particular master integral is computed from any cut which does not set this master integral to zero.
When we use eq. (2.11) to generate IBP relations with doubled propagators 1/ρ 2 k , care must be taken when we work with a cut that includes the k-th propagator. We temporarily treat the propagator 1/ρ k as uncut when constructing IBP-generating vectors from eq. (2.10). After the associated IBP relations have been generated, we re-impose the cut conditions and set any integral to zero if the two powers of 1/ρ 2 k are completely canceled in the integral. More explicitly, we set to zero any integral for which the k-th index in eq. (2.12) is strictly non-positive, rather than being 1 or 2.

Implementation for the 2-loop 5-point nonplanar hexa-box
In this section we give the details of the implementation of the general formalism of the previous subsections to the calculation of the IBP relations required for the differential equations of the nonplanar hexa-box. We need to generate two kinds of IBP relations: 1. IBP relations involving only integrals with no doubled propagators.
Furthermore, it is sufficient for our purposes to generate the IBP relations at numerical values of d. The reasons why this is sufficient will be made clear in the next section. We perform IBP reduction on 8 triple cuts and 2 quadruple cuts, 3 containing the following subsets of propagator lines as labeled by eq. (2.4), When generating IBP relations with a doubled propagator 1/ρ 2 2 using eq. (2.10), the propagator 1/ρ 2 must be uncut, i.e. deleted from the cut lines. The following double cuts are then also needed,∆

JHEP01(2019)006
The IBP relations are in fact generated sector-by-sector. This involves solving eq. (2.10) for every sub-topology supported on the cut, and details of this procedure are presented in appendix A. The two kinds of IBP relations are combined into a single linear system which we solve on each cut using a private code written in Wolfram Mathematica, calling Fermat [55] through A.V. Smirnov's FLink (part of the FIRE5 package [33]). We first set the kinematic variables s ij to numerical values, and perform Gaussian elimination with a partial (row) pivoting scheme that preserves the sparsity of the linear system. The pivotal rows correspond to a minimal set of independent IBP relations. We then perform Gaussian elimination with identical pivot choices but with full analytic dependence on the s ij , starting from the independent IBP relations identified during the previous numerical step. In agreement with ref. [41], we find that there are 73 master integrals (75 after IBP reduction, and 73 after accounting for symmetries).
The linear systems obtained for the quadruple cuts ∆ 9 and ∆ 10 in eq. (2.15) are relatively small and require only a few minutes each to be solved. For the linear systems of the remaining triple cuts, the computation times range from about half an hour to nearly two hours. The total time used to solve these linear systems is about 8 hours, using one CPU core. 4 In the end, the results for the 10 cuts are merged to produce the full IBP reduction table for all integrals with up to rank-3 numerators and up to one doubled propagator 1/ρ 2 2 . By applying a graph symmetry we obtain from it the reduction table for integrals with the doubled propagator 1/ρ 2 3 .
3 Differential equations and basis of pure integrals

Kinematic derivatives and momentum routing
As discussed above, it is sufficient to fix s 45 = 1 and compute the partial derivatives of the master integrals with respect to the other kinematic variables, s 12 , s 23 , s 34 , s 51 . These may be represented as partial derivatives with respect to the four independent external momenta, where k runs over 1, 2, 3, 4 in the implicit summation. The calculation is made simpler by a judicious choice of momentum routing. The only invariant depending on both p 4 and p 5 is s 45 , which we have fixed to be 1. Therefore, as Lorentz transformations leave all kinematic invariants unchanged, infinitesimal variations in the other 4 kinematic variables can always be represented by infinitesimal variations in p 1 , p 2 , p 3 while keeping p 4 and p 5 unchanged. The combination (p 1 + p 2 + p 3 ) = −(p 4 + p 5 ) is also unchanged under these transformations. We thus find that eq. (3.1) can be reduced to

JHEP01(2019)006
It is straightforward to find the β 1 µ ij and β 2 µ ij which give the correct kinematic derivative while preserving the on-shellness of p 1 , p 2 , and p 3 . 5 By inspection of eq. (2.4), we can see that eq. (3.2) annihilates all inverse propagators except ρ 2 and ρ 3 . Therefore, kinematic derivatives of any integral (without doubled propagators) can only have a doubled propagator 1/ρ 2 2 or 1/ρ 2 3 , while the other six propagators will always remain with single powers. As discussed in the previous section, it is thus sufficient for us to perform IBP reduction for integrals with a doubled propagator 1/ρ 2 2 . From those, we can then obtain IBPs for integrals with the doubled propagator 1/ρ 2 3 by applying a graph symmetry corresponding to flipping figure 1 upside down.
The outcome of the procedure are differential equations in the form where each M ij is a 73 × 73 matrix and the I a denote the 73 master integrals. Using differential forms, eq. (3.3) can be written in a more compact form, In the above equation, I is the column vector of the 73 master integrals, dI is its exterior derivative in the kinematic variables, and M is a matrix of one-forms A useful consistency check is the integrability condition, which follows from applying the exterior derivative to both sides of eq. (3.4), or more intuitively, from the fact that partial derivatives commute with each other when we differentiate the master integrals repeatedly.

Transformation to a pure basis
In this section we consider a transformation into a new basis where the differential equations take a particularly simple form, the so-called canonical form. Here we review this concept and introduce some notation. Under a general change of basis of master integrals, In component notation, the new matrix is As argued in ref. [27], for all integrals that evaluate to multiple polylogarithms [56,57], as is the case for the nonplanar hexa-box, we can choose a basis of 'pure' master integrals. These master integrals have the property that the coefficients in their Laurent expansion in have a uniform (transcendental) weight with no rational dependence on the kinematic variables (we refer the reader to e.g. refs. [58,59] for introductory discussions on these topics). The weight of each coefficient is related to the power in the expansion, such that if we assign weight −1 to , each term in the Laurent expansion has the same weight. If one finds a basisĨ where all master integrals are of this form, then the differential equations take a very special form, the so-called canonical form, with where the matrices M α only contain rational numbers, and all dependence on the kinematic variables appears in the dlog-forms d log r α . In component form, In such a basis it is sufficient to compute the differential equations for a numerical value of (or equivalently d) since the dependence on this variable is trivial. Each r α is a letter of the so-called symbol alphabet associated with the integralsĨ a [45]. In the case of the nonplanar hexa-box, the r α are algebraic functions of the s ij , involving in general the square-root of the 5-point Gram determinant. We note that the task of writing the solution to the differential equations in terms of multiple polylogarithms is greatly simplified if the alphabet is rational. This can be achieved for this particular differential equation using the so-called momentum-twistor variables discussed in appendix B. In particular, in these variables it becomes trivial to write a general solution in terms of multiple polylogarithms, leaving only the initial condition to be determined. Using known properties of the analytic structure of Feynman integrals, their symbol can be easily constructed from a differential equation in canonical form. While it is not the same as having a complete solution of the differential equation, the symbol encodes non-trivial information on the analytic structure of the solution of the differential equation.

Constructing a pure basis
The problem of finding a pure basis of master integrals is a complicated one, even though some strategies exist to attempt to find such a basis in an automatic way [60][61][62][63]. There exist also less automatic approaches based on constructing an integral with the appropriate analytic properties. For instance, one might construct an integrand that can be written in a dlog form, or an integral with 'unit leading singularities' which is conjectured to be pure (see [58] for a review of these different approaches with examples).
For the particular integrals that interest us in this paper, explicit expressions for most pure integrals were known in the literature, in particular in refs. [3,64]. We had to construct one pure master in the nonplanar hexa-box top sector, and all the pure master integrals for the five-point planar integrals. 6 For the latter, we observe that pure integrals can be obtained with numerators constructed from specific minors of the Gram matrix associated with the integrand. A basis of pure integrals can be found in appendix C and here we give details on how the missing pure master integrals were constructed.
For the nonplanar sectors, we had to construct a single pure integral for the nonplanar hexa-box topology. The missing pure integral, corresponding to the N 3 numerator in table 1 in appendix C, was computed by constructing an integral with unit leading singularities. We find that the numerator N 3 = s 12 s 23 (l 1 + p 4 ) 2 (l 1 + p 5 ) 2 − l 2 1 (l 1 + p 4 + p 5 ) 2 (3.14) on the nonplanar hexa-box has the required properties. For the planar topologies, there are two distinct five-point topologies, shown in figure 2, each with two master integrals, for which we had to construct pure integrals. These are relevant, for instance, for the sectors associated with (N 17 , N 18 ) and (N 28 , N 29 ) in table 1 in appendix C. In both cases the scalar integral is pure once it is normalized to have unit leading singularities (for the penta-bubble integral, there is also a simple d-dependent normalization).
The remaining pure master integrals are computed from insertions of loop-momentum components beyond 4 dimensions. As a motivation for this approach, consider the one-loop scalar pentagon in 6−2 dimensions. This integral is known to be pure (see e.g. refs. [65,66] for the explicit differential equations in canonical form), and can be written as a (4 −

JHEP01(2019)006
2 )-dimensional integral with a nontrivial numerator insertion through simple dimensionshifting identities. Up to normalization, this numerator is the Gram determinant of the linearly-independent momenta of the integrand [67], which is easily related to the square of the loop-momenta components beyond four dimensions (below, we show this explicitly at two loops). We apply the same principle at two loops and consider tensor insertions of minors of the Gram matrix of the linearly independent momenta of the integrand, denoted G = G(l 1 , l 2 , p 1 , p 2 , p 3 , p 4 ), and sometimes called the Baikov matrix. If we write the ddimensional loop-momenta as l i = (l [4] i , l and define where G [i,j] is the matrix obtained by deleting row i and column j of G. We have introduced the quantity tr 5 , which plays an important role in 5-point kinematics, and can also be defined as a trace of γ-matrices: The above equation will be taken as the definition of the branch of the square root of the 5-point Gram determinant, as well as the chirality of the external momentum configuration in the remaining of this paper. By inserting specific minors from the set listed in eq. (3.17) we are able to construct pure integrals for the two-loop 5-point integrals. As explicit examples, for a box-triangle and a penta-bubble we find that the same numerator insertion gives a suitable choice of master integral: where l 1 is the loop momentum of the box or the pentagon sub-loop respectively (we refer to appendix C for the details on the indexing of each master integral). We note that the factor of 1/ρ 6 means a doubled propagator 1/ρ 2 6 in the integral. Furthermore, we observe that this approach is more general. Indeed, we can find pure master integrals for all planar two-loop 5-point massless topologies using such numerators. For the double-box with five external legs, two pure master integrals are obtained in a simple generalization of the four-point massless case (see e.g. [27]). The third pure master can be chosen to be We note that these insertions are particularly suited to a calculation using the natural variables introduced in section 2.1.
To finish this section, we comment on how the change of basis is implemented in our calculation. We first construct the differential equations in an arbitrary basis of integrals, each of which is given by a single term F [ ν] in the notation of eq. (2.12) for some indices ν. Then we construct the matrix T that implements the transformation to the pure basis, see eq. (3.9). 7 Finally, the inverse matrix T −1 is computed using the software Fermat [55] and Mathematica, and we obtain the differential equations in the pure basis from eq. (3.10).

An alternative construction of differential equations in canonical form
We devised an alternative and much more direct way of obtaining the same system of differential equations, which we now describe. This alternative method only requires IBP reduction at numerical values of the s ij , making the computational-resource requirements almost trivial when combined with our unitarity-compatible IBP reduction. However, it crucially relies on the a priori knowledge of the alphabet and is thus not generally applicable. In the context of this paper, we have thus implemented it as a check of the completely generic approach discussed in previous sections.
The full alphabet for nonplanar 2-loop 5-point integrals has been conjectured to be composed of 31 letters given explicitly in ref. [18]. These are written in terms of the tr 5 defined in eq. (3.18) and the additional variables defined by with s 56 cyclically identified with s 51 . In term of these variables, the letters are with 1 ≤ i ≤ 5. These letters are a possible choice for the r α on the right-hand side of eq. (3.12). Imposing s 45 = 1 eliminates one letter. Using this alphabet and the pure basis of appendix C, we construct the differential equations at 30 randomly chosen rational numerical phase-space points s ij . The matrices JHEP01(2019)006 (M α ) ab on the right-hand side of eq. (3.12) can then be found by solving small linear systems. 8 The use of rational numbers removes the issue of rounding errors during the fitting procedure. Using purely numerical data we thus reconstruct the complete analytic result.
Whenever one is able to construct an ansatz for the alphabet of a given integral and has a pure basis, this method can be used to directly and efficiently construct the differential equations. In particular, it reduces the reliance on the efficiency of the IBP reduction. We note that the test of whether a basis is pure or not can also be done numerically, by verifying if the right-hand side of the differential equations are proportional to at a numerical phase-space point. Furthermore, the validity of the ansatz for the alphabet is automatically determined by whether the linear systems have solutions.

Results
We have presented a completely generic method for constructing differential equations of Feynman integrals using unitarity-compatible IBP relations. To show its potential, we applied it to the nonplanar hexa-box integrals and computed their differential equations. This highly non-trivial example illustrates the power of our approach as it is beyond what can currently be achieved with more standard approaches. In this section we give details on how the differential equations were obtained and the checks that were made. In particular, as a very strong check, we discuss the solution of the differential equations at the symbol level which we have compared against known results in the literature [18].

Differential equations
Following the steps described in the previous sections, we constructed the system of differential equations for the nonplanar hexa-box integrals in canonical form, using the basis of master integrals given in appendix C. The differential equations are obtained at the arbitrarily chosen value of d = 67/13, and working in the finite field of cardinality 7666667, also chosen arbitrarily. We checked that the integrability condition (3.7) indeed holds. We then manipulate the expressions to obtain the explicit decomposition in terms of the dlogforms d log r α and the matricesM α on the right-hand side of eqs. (3.12) and (3.13). We verified that the results remain unchanged when different numerical values of d are used, which confirms the factorization of in the differential equations. Dividing the results by = (4 − d)/2 and reconstructing rational numbers from their finite-field representations, using the extended Euclidean algorithm, we obtain fully analytic results for the differential equations. We note that the rational reconstruction from a single finite field is only possible once the system of differential equations is in canonical form. 9 Indeed, if one were to attempt to reconstruct the rational numbers in a non-pure basis, one would have to perform the evaluation in multiple finite fields and apply the Chinese remainder theorem.

JHEP01(2019)006
For the specific nonplanar hexa-box in figure 1, we observe that the differential equations actually depend on only 24 letters, instead of the 31 listed in eq. (3.23). The remaining 7 letters will contribute once permutations of the external legs are considered, as required e.g. when computing a physical scattering amplitude. For the particular permutation of external legs in figure 1 we find that the following letters are absent: r 8 , r 9 , r 10 , r 21 , r 22 , r 23 , r 24 . (4.1) The differential equations are derived in terms of the Mandelstam invariants s ij . In these variables we observe that, as expected, the r α of eq. (3.12) contain square roots (of the 5-point Gram determinant or, equivalently, tr 5 as defined in eq. (3.18)). As already discussed, it is known that if we transform the results into a new set of so-called momentumtwistor variables [68], under which all spinor components of external momenta have rational components, there should be no square roots. In other words, written in terms of the twistor variables, the r α should be rational functions. We give details of a possible choice of twistor variables in appendix B. As a consistency check of our calculation, we have changed to this new set of variables and observe that the alphabet is indeed rational. As the most stringent check of our result, we have verified that we reproduce exactly the same differential equations if we use the alternative method described in section 3.4. As mentioned there, this is a much more direct and computationally efficient method of obtaining the same result, but it relies on previous knowledge of the alphabet and of a basis of pure master integrals. The timings required to obtaining the differential equations through this approach are negligible, even compared to our approach based on unitarity-compatible IBP relations.
The differential equations are given in the supplementary material diffEqMatrices.m, which contains explicit expressions for the matrices (M α ) on the right-hand side of eq. (3.13), corresponding to the alphabet of eq. (3.23) and the basis of appendix C. The matrices are presented as two-dimensional arrays mat[n] in the Wolfram Mathematica format, with n = 1, . . . , 31 being the index that identifies the associated letter. 10 If a letter is absent, the corresponding matrix is the zero matrix. As an example, the matrix M 31 , corresponding to the letter r 31 = tr 5 , is extremely sparse, with only the following 15 nonzero entries among the 73 × 73 components:  The matrices related to the other letters, given in the supplementary material, also contain very simple rational numbers, whose numerators and denominators are never larger than 16. The amazing simplicity of the results is the reason why we are able to reconstruct these rational numbers from their image in a single finite field with a 7-digit prime modulus.

Symbols of the master integrals
As stated in section 3.2, once one has obtained the differential equation in canonical form for a Feynman integral one can easily construct its solution order by order in at symbol level [45]. Indeed, constructing a generic solution at symbol level from the differential equation is a simple exercise in combinatorics: using the different letters one constructs a set of tensors that are consistent with the differential equation (see e.g. refs. [58,59] for a pedagogical introduction). The weight-zero contributions must a priori be determined from the initial condition, as they are just rational numbers that vanish upon differentiation. However, given the iterative structure of the solution of the differential equation, where the solution at order n depends on the solution at order n−1 , see eq. (3.13), the weight-zero constants appear as rational coefficients of the higher-weight terms. We can then use arguments on the analytic structure of the answer, such as the position of the branch points, to find relations between the different weight-zero constants [27]. More precisely, we impose the so-called first-entry condition [46] which requires all symbol tensors to have physical branch points. In our example, the weight-zero constants appear at order 1/ 4 , and we should determine 73 such constants. The first entry condition means that the first entries of the tensors must correspond to the letters r 0+i or r 15+i in eq. (3.23), which completely fixes all the 73 weight-zero coefficients up to an undetermined overall factor. This factor can be fixed from the simplest master integral with three propagators, which is trivial to compute to all orders in . 11 Following these steps, we computed the symbol of the 73 master integrals from our differential equations in canonical form. As a sample result, the first integral in our pure basis, corresponding to the N 1 numerator in eq. (C.3), has no 1/ 4 or 1/ 3 poles. The coefficient of the 1/ 2 pole is given by weight-two symbol tensors built from the letters in eq. which satisfies the first-entry condition stated above.
To check the validity of our results, we re-computed the symbols of the two nonplanar five-point integrals studied in ref. [18] and found complete agreement with their results. The symbols of all the master integrals, including expansion terms from order 1/ 4 to order 0 , are given in the supplementary material symbols.m. Another supplementary material masters.m gives the pure master integral basis, as tabulated in appendix C.

Conclusions
In this paper we showed that unitarity-compatible IBP relations can be used to construct differential equations for Feynman integrals. These relations reduce the size of the IBP system by avoiding to raise propagator powers, which leads to the proliferation of auxiliary integrals in more standard techniques. Although the idea of using this type of IBP relation was introduced already some years ago [35], only recently have they been used as a powerful tool in fully-fledged practical calculations. In our approach to construct the differential equations, the efficiency gained by using such compact IBP relations was complemented by

JHEP01(2019)006
ideas from two-loop numerical unitarity [8,39]: we construct analytic IBP relations sectorby-sector, so that we do not have to obtain all IBP-generating vectors simultaneously. Using our refined implementation, we were able to solve the IBP system needed to construct differential equations for all the nonplanar hexa-box integrals in 8 hours on one CPU core. This required the reduction of tensor numerators of rank up to three and up to one doubled propagator. We believe this shows the potential of our approach: with very limited computational resources we are able to tackle a problem that is beyond the reach of current tools even when using much more computer power. We thus look forward to applying our method to new and more complicated problems. We present explicit results for the differential equations of the nonplanar hexa-box integrals and their solution at symbol-level. The differential equations are written for a pure basis of master integrals. Most of them were available in the literature, and we constructed the remaining ones ourselves. For the 5-point planar two-loop topologies, some pure integrals were obtained by inserting as numerators specific minors of the Gram matrix of the linearly-independent external and loop momenta, which vanish in 4 dimensions. This was shown to be a good approach even for the more complicated planar topologies, such as the double-box with five external legs and the penta-box, that are not required for the particular differential equations we study in this paper. Numerators obtained with these minors have a very compact representation. We believe it would be very interesting to better understand this type of numerator insertions, and why they are good candidates for pure integrals.
Furthermore, we formulated another method for constructing differential equations in a very efficient manner. This method relies on prior knowledge of the alphabet. In the present case, full analytic results are reconstructed from numerical evaluations at 30 numerical phase-space points, which makes the IBP reduction trivial (even with conventional Laportatype algorithms that do not use unitarity-compatible IBP relations). To make this method more generally applicable, one needs to be able to write a conjecture for the alphabet. We note nevertheless that the validity of the conjecture is verified by the success of the fitting procedure to reconstruct the differential equation. We expect that in the near future this method can be used in other 2-loop 5-point topologies for which the alphabet has been conjectured [18], but also that it will be applicable to more complicated diagrams once techniques to construct conjectures for the alphabet become more developed.
Our results constitute a significant step towards obtaining analytic expressions for all the master integrals relevant for 5-point massless amplitudes in the nonplanar sectors. Indeed, a convenient set of variables that rationalizes the alphabet is known and all that remains to obtain explicit expressions in terms of multiple polylogarithms is the calculation of boundary conditions for the differential equations. Combined with the recent progress in the reduction of 5-point massless amplitudes to master integrals, we believe the calculation of these amplitudes beyond the large-N c limit is now well within reach.

A Sector-by-sector generation of IBP relations
We generate IBP relations sector by sector. The sector of an integral refers to the set of propagators that are present in the integral, and the sector of an IBP relation is defined by the integral (in the relation) with a maximum number of propagators.
For the nonplanar hexa-box, there are 8 propagators, but the IBP relations contain integrals of lower sectors with only a subset of these propagators. Sectors will be arranged into 'levels' according to the number of propagators; those with more propagators will be called higher sectors, and those with fewer propagators will be called lower sectors.
As an example, consider IBP relations for the integral topology shown on the top right corner of table 1, where the propagator 1/ρ 1 of the nonplanar hexa-box has been canceled. Such IBP relations may be generated by eq.
Eq. (2.10) is now trivially satisfied for j = 1, so we only need to solve the modified equation 2) will appear as a proper subset. However, in practice we impose a degree bound when solving eq. (2.10) to speed up the computation. 12 As a result, we will miss the solutions of v µ i related to lower sectors with many canceled propagators, because many inverse propagators, not just ρ 1 , will have to appear in eq. (A.1). Fortunately, the limited solutions of v µ i still generate a full set of IBP relations on the maximal cut of the nonplanar hexa-box, i.e. ignoring integrals of lower sectors. Then we recursively descend into lower sectors and solve the corresponding versions of eq. (A.2), in order to generate a full set of IBP relations in the lower sectors. In the end, we combine all the relations into one linear system and solve them together.
In this paper, the above sector-by-sector procedure is carried out on cuts. For each cut, we only descend into those sectors where none of the cut propagators are canceled, and combining IBP relations from all sectors gives us a complete set of IBP relations on the particular cut.

JHEP01(2019)006 B Momentum-twistor variables
Although we use the kinematic invariants in eq. (2.1) to compute the differential equations, eventually one might be interested to change to momentum-twistor variables [68], under which all spinor components of external momenta have rational components. This ensures that the pure integrals, tabulated in appendix C, do not involve square roots in their alphabet. Having a rational alphabet is also an important step if one wants to find an explicit solution in terms of so-called Goncharov multiple polylogarithms [56,57].
The momentum-twistor parametrization of external kinematics is highly redundant due to a GL(1) rescaling symmetry for each twistor and a global SL(4) symmetry. We adopt a convenient choice of the independent variables x 1 , x 2 , . . . , x 5 given in [4], but with a cyclic shift s ij → s i+3 j+3 , The inverse map of eq. (B.1) is also found in ref. [4], but again slightly modified by us with the cyclic shift s ij → s i+3 j+3 , As always, we have set s 45 = 1 to eliminate a trivial overall scale from the problem.

C Table of basis integrals of uniform transcendentality
In this appendix, we explicitly write the numerators giving pure master integrals, omitting the propagators of the integrals. Some of our numerators include a propagator, which means that the propagator will be doubled in the corresponding master integral, see e.g. eq. (C.12). All nonplanar sectors are written explicitly. For planar sectors we write a single representative and list the ones that can be obtained from it. Our notation for the sectors and the
Many pure nonplanar 5-point integrals are given in ref. [3]. The numerators are often written in the spinor helicity formalism. 13 Once manipulated to become neutral under the little group, they can be converted to expressions in the invariants s ij and the tr 5 defined in eq. (3.18) using the following formulas, [ij] jk [kl] li = 2s il s jk − 2s ik s jl + 2s ij s kl + tr 5 2 sgn(i, j, k, l) .
In the second and third lines of the equation above, (i, j, k, l) are permutations of (1, 2, 3, 4), and the sgn symbol denotes the signature of the permutation. When dealing with expressions involving |5 or |5], we first rewrite the expressions using only the spinors of the legs 1, 2, 3, 4, and then apply the last two relations of eq. (C.1). After rewriting the numerators in terms of s ij and tr 5 , we define the even-parity and odd-parity parts of any numerator N as Recall that tr 5 is imaginary in Minkowskian signature, and is proportional to a contraction with the Levi-Civita ε-tensor, which is odd under parity. By projecting the numerators of ref. [3] onto even and odd parts to construct our basis, we make sure that the integrands have definite parity, and are either purely real or purely imaginary. Nontrivial complex values only arise after loop integration.

JHEP01(2019)006
Sector and diagram Numerators Sector and diagram Numerators  Table 2. Diagrammatic representation of each sector and associated numerators.

JHEP01(2019)006
Sector and diagram Numerators Sector and diagram Numerators  Table 3. Diagrammatic representation of each sector and associated numerators.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.