Complete integration-by-parts reductions of the non-planar hexagon-box via module intersections

We present the powerful module-intersection integration-by-parts (IBP) method, suitable for multi-loop and multi-scale Feynman integral reduction. Utilizing modern computational algebraic geometry techniques, this new method successfully trims traditional IBP systems dramatically to much simpler integral-relation systems on unitarity cuts. We demonstrate the power of this method by explicitly carrying out the complete analytic reduction of two-loop five-point non-planar hexagon-box integrals, with degree-four numerators, to a basis of $73$ master integrals.


Introduction
High precision in the theoretical predictions for cross sections is necessary for the analysis of the Large Hadron Collider (LHC) Run II experiments. For many processes, the nextto-next-to-leading-order (NNLO) contributions are needed for precision-level predictions. NNLO in general requires the computation of a large number of two-loop (or, in some cases, higher-loop) Feynman integrals, which is often a bottleneck problem for particle physics.
Integration-by-parts (IBP) reduction is a key tool to reduce the large number of multiloop Feynman integrals to a small integral basis of so-called master integrals. Schematically, for an L-loop integral in dimensional regularization with the inverse propagators D 1 , . . . , D m , we have where the v µ j are vectors constructed from the external momenta and the loop momenta. Expanding the action of the derivative in eq. (1.1) leads to an IBP identity, which is a linear relation between multi-loop Feynman integrals. After collecting a sufficient set of IBP identities, one can carry out reduced row reduction in a specific integral ordering to express a target set of integrals as a linear combination of the master integrals.
The row reduction of IBP identities can be achieved by the Laporta algorithm [1,2]. There are several publicly available implementations of IBP reductions: AIR [3], FIRE [4,5], Reduze [6,7], LiteRed [8], Kira [9], as well as various private implementations. In the standard Laporta algorithm and many of its variants, integrals with doubled propagators appear in the intermediate steps and also in the master integral basis. To trim the linear system of IBP identities for row reduction, one may impose the condition that no integrals with doubled propagators appear in the IBP identities [10]. This condition can be solved by the syzygy computation [10] or linear algebra [11]. In refs. [12,13], methods for deriving IBP reductions on generalized-unitarity cuts was introduced. IBP reduction can also be sped up by the finite-field sampling method [14]. For several multiple-scale Feynman integrals, inspired choices of IBP generating vectors v µ j , leading to simple IBP identities, can be obtained from dual conformal symmetry [15]. In a recent development, IBP identities which target integrals with arbitrary-degree numerators [16] were introduced. An initial (and simpler) step to carrying out the IBP reductions is that of determining the integral basis itself. This can be done by the packages Mint [17] or Azurite [18]. The dimension of the integral basis can be efficiently determined by the computation of an Euler characteristic which naturally arises in the D-module theory of polynomial annihilators [19].
For multi-loop, high-multiplicity or multi-scale amplitude computations, the IBP reductions frequently become the bottleneck which requires extensive computing resources and time. Schematically, the complexity of IBP reductions mainly originate from the following facts, 1. Multi-loop IBP identities involve many contributing integrals, and the complexity of row reduction grows quickly with the number of integrals.
2. There are many external invariants (i.e., Mandelstam variables and mass parameters), and the algebraic computation of polynomials or rational functions in these parameters requires a substantial amount of CPU time and RAM.
In this paper, utilizing our new powerful module-intersection IBP reduction method, we present solutions to the above problems: 1. With the module-intersection computation [48], we dramatically reduce the number of relevant IBP identities and integrals involved by imposing the no-doubled-propagator condition and applying unitarity cuts.
2. With the classical technique of treating parameters as new variables and imposing a special monomial ordering, originating from primary decomposition algorithms [49] in commutative algebra, we can efficiently compute analytic IBP reductions with many parameters.
In this paper, we demonstrate our method by treating the cutting-edge example of the analytic IBP reduction of two-loop five-point non-planar hexagon-box integrals, with all degree-1, 2, 3, 4 numerators, to the basis of 73 master integrals. This computation was carried out using private implementations in Singular [50] and Mathematica of the ideas presented in this paper. Our result of fully analytic IBP reductions can be downloaded from, https://github.com/yzhphy/hexagonbox_reduction/releases/download/1.
This paper is organized as follows. In section 2 we present our module-intersection method for trimming IBP systems and describe a simple way of obtaining the individual modules needed for our method. In section 3, the central part of our paper, we present a highly efficient new method for computing these intersections. We then explain the technical details of performing row reduction for the trimmed IBP systems from module intersections in section 4. In section 5, we present our example, the module-intersection IBP reduction of hexagon-box integrals, in details. The conclusions and outlook are given in section 6.

Simplification of IBP systems by module-intersection method
In this section we present the details on trimming IBP systems by requiring the absence of integrals with doubled propagators [10], and by furthermore applying unitarity cuts. To demonstrate our method clearly, we first explain how to efficiently impose the nodoubled-propagator condition without applying unitarity cuts. Then we show that the same approach can be applied on unitarity cuts. In both cases, the algebraic constraints for trimming IBP systems are reformulated as the problem of computing the intersection of two modules over a polynomial ring [48]. At the end of this section, the algorithm for computing each individual module is introduced, while the intersection algorithm, the heart of this method, will be introduced in section 3.

Module intersection method without unitarity cuts
The objects under consideration are the multi-loop Feynman integrals, Here ν i ∈ Z, l 1 , . . . , l L denote the loop momenta, and p 1 , . . . , p E the independent external momenta. By the procedure of integrand reduction, we can set m = L(L + 1)/2 + LE. The D i denote inverse propagators.
To study IBP reduction of integrals of the form (2.1), we focus on the family of Feynman integrals associated with a particular Feynman diagram with k propagators (where k ≤ m) and all its daughter diagrams, obtained by pinching propagators. Without loss of generality, the inverse denominators of this Feynman diagram can be labeled as D 1 , . . . , D k . Therefore, the family of integrals is, Traditionally, IBP reduction is carried out within such a family. However, since the majority of the Feynman integrals that contribute to a quantity in perturbative QFT are integrals without doubled propagators, it is natural to consider the subfamily and the IBP relations for integrals in this subfamily [10]. We find that it is convenient to use the Baikov representation [51] of Feynman integrals for their integration-by-parts (IBP) reduction for several reasons: 1) the integrand reduction is manifest in this representation, 2) it is easy to apply unitarity cuts, 3) most importantly, it is surprisingly simple to trim IBP systems analytically in this representation. Here we briefly review the Baikov representation.
We collect the external and internal momenta as, The Gram matrix S of these vectors is, where the elements are defined as The upper-left E × E block is the Gram matrix of the external momenta, which is denoted as G. Defining z i ≡ D i and integrating out solid-angle directions, the Feynman integral (2.1) takes the following form in Baikov representation, where P ≡ det S, U ≡ det G, and the factor C L E originates from the solid-angle integration and the Jacobian for the transformation x i,j → z. For the derivation of IBP identities, U and C L E are irrelevant, so we may ignore them in the following discussion. Note that in this representation, the inverse denominators D i become free variables, and so the integrand reduction can be done automatically.
An IBP identity in this representation reads, Here the a i denote polynomials in the ring A = Q(s)[z 1 , . . . , z m ]. (We use s to represent the independent Mandelstam variables and mass parameters collectively.) We remark that in the Baikov representation, the Baikov polynomial P vanishes on the boundary of the integration domain, and hence there is no surface term in the IBP identity. Note that the terms with the pole 1/P appearing inside the parenthesis in eq. (2.7) correspond to dimension-shifted integrals, and as such are not favorable in deriving simple IBP relations.
To avoid such poles, we may impose the following constraints on the a i [12,13,52], where b is also required to be a polynomial in A. This constraint is known in computational commutative algebra as a "syzygy" equation [10]. In the following discussion, we only focus on the polynomials a i , since once they are known it is straightforward to recover the polynomial b. The solutions of eq. (2.8), taking the form, (a 1 , . . . , a m ) (2.9) form a sub-module of the polynomial module A m , which we denote M 1 in the following. Furthermore, to trim the IBP system, it is possible to work with integrals in F ndp , i.e. integrals without doubled propagators [10]. Note from the second line of eq. (2.7) that, even if the integral inside the differential operator has ν i ≤ 1, the derivative will produce integrals with doubled propagators. This can also be prevented by a suitable choice of the a i . For example, if we consider the integral of the parent diagram, ν 1 = · · · = ν k = 1, ν k+1 ≤ 0, . . . , ν m ≤ 0 inside eq. (2.7), we can avoid doubled propagators by requiring that a i is divisible by z i , Such (a 1 , . . . , a m ) also form a sub-module of A m , which we denote M 2 . We need to solve eqs. (2.8) and (2.10) simultaneously to find IBP relations which involve neither dimension shifts nor doubled propagators in Baikov representation. The strategy in ref. [13] to solve these conditions is to replace a i by b i in eq. (2.8) and then solve for (b 1 , . . . , b k , a k+1 , . . . , a m , b) as a single syzygy equation, employing Schreyer's theorem. However, this approach, although it works well for simple two-loop four-point integrals, becomes less practical for more complicated kinematics. It was suggested in the reference [48] that a better strategy is to determine the generators of M 1 and M 2 individually, and then calculate the module intersection, In subsection 2.3, we will see that it takes no effort to find the generators of M 1 and M 2 . In section 3 we present a highly efficient algorithm for computing the intersection Once eqs. (2.8) and (2.10) are solved, we obtain the simplified IBP identities without doubled propagators. They take the following form, (2.13) Once the generators of M 1 ∩ M 2 are obtained, we multiply them by monomials in the z i in order to get a spanning set of IBP identities for the reduction targets. Alternatively, it is possible to apply D-module theory to get the IBP reductions directly from the generators of M 1 ∩ M 2 . We leave this direction for future research.

Module intersection method with cuts
In practice, for Feynman diagrams with high multiplicity or high loop order, instead of working with the Feynman integrals directly, it is more convenient to apply unitarity cuts [12,13]. In this section we show how to apply the module intersection method in combination with unitarity cuts, thus simplifying the construction and subsequent Gaussian elimination of the IBP identities. We follow the notation of ref. [54].
Consider the c-fold cut of eq. (5.2) with c ≤ k. Let S cut , S uncut and S ISP denote the sets of indices of cut propagators, uncut propagators and irreducible scalar products (ISP) respectively. Explicitly, (2.14) In Baikov represention, the c-fold cut of integrals with ν 1 = · · · = ν k = 1 takes a simple form, where,P = P | z ζ 1 →0,...,z ζc →0 . (2.16) We can derive IBP identities for the integrals on the cut S cut , Once again, to derive simple IBP relations without dimension shifts or doubled propagators [12,13], we may impose the conditions, Again, the module intersection method provides an efficient tool to solve these constraints simultaneously. However, for the purpose of automation, it is useful to make use of the following slight reformulation. Definẽ (2.22) By comparing these equations and their counterparts without applied cuts, we observe the following shortcut. Recall that the modules M 1 and M 2 defined in the previous subsection are the solution sets for eqs. (2.8) and (2.10) respectively. We definẽ whereby the c-fold cut on the elements of M 1 and M 2 has been imposed. Then clearly, solve the equations (2.21) and (2.22) simultaneously. Note that any element in (q 1 , . . . , q m ) ∈ M 2 has the property After imposing the cut we have, which is consistent with the requirement (2.20). This automates the computation: to obtaiñ M 1 andM 2 , we first compute M 1 and M 2 , and then simply apply the rules for the various cuts. An algorithm for computation of the intersection will be introduced in section 3.
Once the constraints (2.21) and (2.22) are solved, the IBP identities on the cut take the simple form, As for the uncut case, once the generators ofM 1 ∩M 2 are obtained, we multiply them by monomials in the z r i , i = 1, . . . , m − c to get a spanning set of IBP identities. The cuts necessary for reconstructing the complete IBP identities can be determined from a list of master integrals [13,18]: they are the maximal cuts of "uncollapsible" master integrals (master integrals which cannot be obtained from other integrals in the basis by adding propagators). In practice, the list of master integrals can be quickly determined by the packages Mint [55], Azurite [18]. The total number of master integrals can also be determined by the D-module theory method of ref. [19].
Upon applying Gaussian elimination to the IBP identities evaluated on each cut, we can subsequently merge the obtained reductions to find the complete IBP reductions without applied cuts.

Algorithm for computing individual modules
The modules M 1 and M 2 defined in subsection 2.1, and henceM 1 andM 2 defined in the previous subsection, can all be determined without effort.
The condition (2.8) for M 1 is a syzygy equation for the Baikov polynomial P and its derivatives. Schreyer's theorem [56] guarantees that solutions for syzygy equations can be obtained from Gröbner basis computations. However, for the Baikov polynomial P this is not necessary, owing to the special structure of P . A convenient way to find the solution, or equivalently the tangent vectors for the hypersurface P = 0, is to use the basic canonical IBP identities [12]. Here alternatively, we use the Laplace expansion method [54] 1 of the Gram determinant S (2.5) to determine M 1 , since the results are manifestly expressed in the z variables.
Laplace expansion of P = det S yields, These L(E + L) relations are syzygy relations between P and its derivatives in the x i,j variables. It is straightforward to convert them to solutions of eq.
where a i,j and b are obtained by applying the chain rule to the expressions in eq. (2.27). Explicitly, they are given by, It is proven in ref. [54] via Józefiak complexes [57] that the L(E + L) tuples of polynomials (a i,j ) α form a complete generating set of M 1 . We remark that, 1. The generating set (2.29) is at most linear in the z i .
2. The generating set (2.29) is homogenous in the z i and the Mandelstam variables/mass parameters, as can be inferred from dimensional analysis.
The second property is crucial for our highly efficient algorithm for computing intersections of modules, to be explained in section 3. The generating set for M 2 is trivial. There are m generators, Here e j is the j-th m-dimensional unit vector. The cut cases,M 1 andM 2 defined in eq. (2.23), can then be obtained from eqs. (2.29) and (2.30) by simply setting z ζ i → 0, i = 1, . . . , c.

Determining the module intersection
In this section, we introduce a highly efficient algorithm for computing generators of the module intersections and an algorithm to trim the generating sets.
Given the two submodules M 1 , M 2 ⊂ A t over the multivariate polynomial ring A = F [z 1 , . . . , z n ] with coefficients in the multivariate rational function field F = Q(c 1 , . . . , c r ), our goal is to obtain a generating system of the module M 1 ∩M 2 (which is finitely generated since A is Noetherian). We apply Gröbner basis techniques to address this problem. Our algorithms are implemented in the computer algebra system Singular, which focuses on polynomial computations with applications in commutative algebra and algebraic geometry [50].
We first recall some terminology: Using the notation z α = z α 1 1 ·. . .·z αn n for the monomials in A, we call z α e i with a unit basis vector e i ∈ A t a monomial of A t . An F -multiple of a monomial is called a term. A monomial ordering on A t is an total ordering > on the set of monomials in A t , which respects multiplication, that is, z α e i > z β e j implies z α z γ e i > z β z γ e j for all α, β, γ, i, j, and z α e i > z β e i ⇔ z α e j > z β e j for all α, β, i, j. Then > induces a monomial ordering on the monomials of A, which we again denote by >. In turn, any monomial ordering > on A induces two canonical monomial orderings on A t , position over term z α e i > z β e j :⇐⇒ i < j or (i = j and z α > z β ) (3.1) and analogously term over position. We call a monomial ordering on The monomials of A t come with a natural partial order, which we call divisibility For terms c 1 z α e i and c 2 z β e j with z α e i | z β e j we define their quotient as Moreover, we define the least common multiple lcm(z α e i , z β e j ) as zero if i = j, and as lcm(z α , z β ) otherwise. For a subset G ⊂ A t , the leading module L(G) is the module of all A-linear combinations of the lead monomials of the non-zero elements of G.
By iteratively canceling the lead term of f via multiples of lead terms of the divisors in G = {g 1 , . . . , g l } ⊂ A t with respect to a fixed global ordering, we obtain a notion of division with remainder yielding a division expression , and the lead monomial of f is not smaller than that of any a i g i .
Let U ⊂ A t be a submodule and > a global monomial ordering Theorem 1 (Buchberger). With notation as above, the following conditions are equivalent: is the so-called S-polynomial (or syzygy polynomial ) of f and g.
For a proof of this standard fact, see for example section 2.3 of ref. [58]. A generating set G of U can be extended to a Gröbner basis by means of Buchberger's algorithm, which according to the above criterion computes the remainder r = NF > (spoly > (g i , g j ), G) for all g i , g j in G, adds r to G if r = 0, and iterates this process with the updated G until all remainders vanish. This process terminates, since A t is Noetherian and, hence, any ascending chain of submodules becomes stationary (section 2.1 of ref. [58]). Along this process, we can determine all relations between the g i : with regard to the position over term order, h 1 , . . . , h m are the elements of H in t+l i=t+1 e i , and π : A t+l → A l is the projection onto the last l coordinates, then is generated as an A-module by π(h 1 ), . . . , π(h m ).
Algorithm 2 in conjunction with Lemma 3 can be used to determine a generating system of M 1 ∩ M 2 . For the module intersection problems arising from the non-planar hexagon-box diagram however, the performance of Buchberger's algorithm over F is not sufficient to yield a generating system in a reasonable time-frame. It turns out that a classical technique (which to our knowledge dates back to ref. [49]) to simulate computations over rational function fields via polynomial computations is much faster. We apply this technique to the Gröbner basis computation yielding the syzygy matrix used to determine the module intersection.
Let G ⊂ B t be a Gröbner basis of U with respect to a global block ordering (> 1 , > 2 ) with blocks z 1 , . . . , z n > c 1 , . . . , c r . Then G is also a Gröbner basis of U with respect to > 1 .
Proof. Denote the block ordering (> 1 , > 2 ) by >. Every f ∈ U can be written as By h · f ∈ U and G being a Gröbner basis of U , we know that NF(h · f, G) = 0. Hence, there is a g ∈ G with L > (g) | L > (h · f ) and is a unit (invertible) in A and > is a block ordering, we obtain that L > 1 (g) | L > 1 (f ). This argument shows that L > 1 (U ) = L > 1 (G), that is, G is a Gröbner basis of U with respect to > 1 .
Remark 6. The method of Lemma 5 turns out to be efficient since the input modules in our setting are homogeneous in the variables z 1 , . . . , z n , c 1 , . . . , c r , which allows for efficient sorting of the S-polynomials in Buchberger's algorithm by degree.
Remark 7. For our setting, modular techniques, which compute over finite fields, combine the results using the Chinese remainder theorem and apply rational reconstruction (as developed in a general setting in refs. [59] and [60]) seem not to be useful since very large constant coefficients occur. As a result, this approach would require considering a large number of primes for lifting.
The module intersection algorithm resulting from Lemma 3 and Lemma 5 produces generating sets which usually are not minimal in any sense. For a homogeneous module, a minimal generating system can be determined, however, the computation is very expensive. Another option is to determine the unique reduced Gröbner basis. A Gröbner basis g 1 , . . . , g l is called minimal if L(g i ) does not divide L(g j ) for all i = j. Such a minimal Gröbner basis is called reduced if none of the terms of the tails g i − LT(g i ) is divisible by some L(g j ). It does, however, also not make much sense to pass to a minimal or reduced Gröbner basis, since Gröbner bases can be much larger than generating systems and usually cannot be obtained in a reasonable time in our setting. We hence employ a randomized algorithm for trimming the generating systems to remove extraneous generators. This algorithm is based on determining reduced Gröbner bases after passing to a finite field and specific values of parameters c i : and Multiple runs of the algorithm with different p and p i reduce the chance of a bad prime or a bad parameter value. We apply Algorithm 8 iteratively to drop generators, starting with generators of large (byte) size.

Sparse row reduction
In this section we turn to discussing linear algebra techniques. Although the generation of IBP identities using eq. (2.26) is very fast, achieving the reduction of target integrals to linear combinations of master integrals is highly non-trivial. This is because the latter step requires computing the row reduced echelon form (RREF) of the IBP identities, which becomes computationally intensive in cases with multiple external invariants (i.e., Mandelstam variables and mass parameters), as these enter the IBP system as parameters. Therefore, analytic computation of the RREF requires sophisticated linear algebra techniques.

Selection of relevant and independent IBP identities
The IBP identities generated from eq. (2.26) usually contain linearly redundant identities, as well as identities which are irrelevant for reducing the target integrals. To speed up linear reduction in the subsequent step, we make use of the following methods to remove redundant and irrelevant identities.
1. Removal of redundant linear identities. This can be done with the standard linear algebra algorithm of picking up independent rows of a matrix. We construct the matrix of all requisite IBP identities, sort the rows by their density (i.e., number of non-vanishing entries) or their byte count, and then compute the RREF of the transposed matrix numerically. The pivot locations correspond to the linearly independent IBP identities, giving preference to sparser IBP identities, or IBP identities of smaller sizes.

2.
Removal of irrelevant linear identities. Furthermore, given a target integral set, we can single out the relevant IBP identities which will ultimately reduce them to master integrals. This can also be done with a numeric RREF. We carry out the reduction numerically and record the rows used for reducing the targets in the computation (by recording the left-multiplying matrix of the row reduced matrix).
Regarding the numeric RREF above, in practice we assign generic integer values to all external invariants (Mandelstam variables and mass parameters) and the spacetime dimension and work with finite fields. The computation is powered by the highly efficient sparse finite-field linear algebra package SpaSM [61].

Sparse REF and RREF strategies
We find that the IBP system that arises from eq. (2.26), after removing the linearly dependent and irrelevant IBP identities for targets is, in general, very sparse. To find the row echelon form (REF) and RREF efficiently, it is crucial to apply a sophisticated pivoting strategy to keep the linear system sparse in the intermediate steps.
First, we write the IBP identities in the form of a matrix, with the columns sorted according to some integral ordering, for example like that of Azurite. Then the RREF will eventually reduce the target integrals to the Azurite master integral basis. However, it is important to swap rows and columns during the REF computation, and find suitable pivots for row reduction, in order to keep the matrix sparse. This can be achieved by the heuristic Markowitz strategy [62], provided that all entries in the sparse matrix are of a similar size. However, in our cases, the entries are polynomials or rational functions in Mandelstam variables and mass parameters, and so a weighted pivot strategy, considering both the sparsity and the byte sizes, is used.
Note that we use a total pivoting strategy for which both row swaps and column swaps are used. The row swaps will not change the final result of the RREF, whereas the column swap will change the final result of the RREF. This means that the target integrals will typically not be reduced to the desired master integral basis, but a different integral basis. If we require the target integrals to be reduced to a specific pre-determined basis (say, the Azurite basis), a basis change must be carried out after RREF computation. We find that it is more efficient to allow both row and column swaps, and then compute the basis change, than to allow row swap only (partial pivoting strategy).
The REF and RREF algorithm is implemented in our primitive Mathematica code. A more efficient implementation in Singular is in preparation and will become available soon.

The non-planar hexagon-box diagram example
To demonstrate the power of our new method, we consider a cutting-edge integration-byparts reduction problem: the reduction of two-loop five-point non-planar massless hexagonbox integrals. Recently, this diagram has attracted a great deal of interest, and the hexagonbox integral with a chiral numerator has been analytically computed by use of the bootstrap method and by superconformal Ward identities [63,64]. Here we consider the analytic IBP reduction of hexagon-box integrals with arbitrary numerators with the degree up to four.
The hexagon-box diagram, and the necessary cuts for deriving the IBP identities, are illustrated below in figure 1.
The advantage of applying cuts is that the number of integrals in the IBP relations on each cut will be much less than when no cuts are applied.  Table 1. "Pre"-master integrals supported on each of the 10 cuts necessary to construct the complete IBP reductions.
Hence the set of 73 master integrals is, In this paper, we first reduce integrals to the linear combination of 75 "pre"-master integrals and then apply eq. (5.6) to further achieve the reduction to the linearly independent 73 master integrals.

Reduction of IBP identities on cuts
After computing the module intersections on the 10 cuts, we use eq. (2.26) to generate IBP identities without doubled propagators on each cut in turn.  We calculated the row reduced echelon form (RREF) of these linear systems using the total pivoting strategy, in order to retain the sparsity in intermediate steps. For the linear systems corresponding to some cuts, we are able to directly obtain the RREF analytically in s 12 , s 13 , s 14 , s 23 , s 24 and the spacetime dimension D. For linear systems corresponding to other cuts, we calculate the RREF with integer values of one or two s ij repeatedly. The fully analytic RREF is then readily obtained by our private heuristic interpolation algorithm. These computations are implemented in our primitive Mathematica code. A much more efficient RREF code in Singular is in preparation.
We also remark that for the computation of different cuts, a useful trick is to work with different choices of independent Mandelstam invariants. For a specific cut, a suitable choice of Mandelstam invariants can speed up the computation and also save the RAM usage. After the RREF is obtained, we use the program Fermat [65] to replace the new Mandelstam variables by the original choice s 12 , s 13 , s 14 , s 23 , s 24 .
The running time and resources required of our Mathematica code depends on the size of the linear systems and also the complexity of coefficients in the reduced IBPs. For the smallest linear system, the one corresponding to the quadruple cut {1, 4, 6, 7}, our Mathematica RREF code obtained the fully analytical RREF in 31 minutes with one core and 1.5 GB RAM usage on a laptop with 16 GB RAM. For the largest linear system, that corresponding to the triple cut {3, 6, 7}, we run our Mathematica code and assign integer values to two Mandelstam invariants. This finished in 2.5 hours and used 1.8 GB RAM. We parallelize the running with various integer values, evaluating 440 points on the IRIDIS High Performance Computing Facility. The semi-analytic results are then interpolated to the fully analytic RREF result by our heuristic multivariate interpolation algorithm, requiring a CPU time of 23 minutes with one core and 15 GB RAM usage.

Merging of IBP reductions and final result
Having obtained the RREF of the 10 IBP systems on cuts, it is straightforward to merge the coefficients to obtain the complete IBP reductions without applied cuts [13]. For example, to determine the reduction coefficient of I 42 = I(1, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0), we search for I 42 in table 1 and find that it is supported on the cuts {1, 5, 7} and {3, 6, 7}. We find that, for every target integral, the reduction coefficient for I 42 on the cut {1, 5, 7} equals the corresponding coefficient on the cut {3, 6, 7}, as must be the case. Thus, for each target integral in turn, we obtain its reduction coefficient of each I j as the reduction coefficient of I j on each cut supporting I j . In this way we obtain the complete reduction of the target integral without applied cuts.
After merging the 10 RREFs we reduce the 32 target integrals in eq. (5.8) to the 75 "pre"-master integrals {I 1 , . . . , I 75 }. We then apply the global symmetry in eq. (5.6) to eliminate the redundant integrals I 63 and I 68 , and then finally obtain the complete analytic reduction to the 73 master integrals. The IBP reduction result, as replacement rules (in a compressed file of the size 280 MB), can be downloaded from the link provided in the introduction.

Comparison with other IBP solvers
We have checked our results with FIRE5 [5] in the C++ implementation with LiteRed [8], and KIRA (v 1.1) [9]. We note that it is not an easy task to perform the analytic IBP reduction of hexagon-box integrals with degree-four numerators. Due to the RAM limit we faced, we have not yet been able to obtain the analytic IBP reduction for degree-fournumerator hexagon-box integrals from FIRE5 or KIRA with the Rackham cluster, on a node with two 10 core Intel Xeon V4 CPU and 256 GB of RAM 2 . On the other hand, it is easy to obtain the numeric IBP reduction for degree-four-numerator hexagon-box integrals with FIRE5 and KIRA. For example, if all the five s ij are taken to be integers, FIRE5 is able to generate the purely numeric IBP reductions in about 6.0 hours. We ran FIRE5 purely numerically for many sets of integer values, and all the numeric IBP reduction results are consistent with our analytic IBP reductions, after a basis change between the Azurite integral basis and the FIRE basis.

Conclusion and outlook
In this paper we have presented a new and efficient method for computing integration-byparts (IBP) reductions based on the ideas developed in refs. [13,48]. We used a module intersection method [48] to trim IBP systems on unitarity cuts. The key idea is the efficient analytic computation of module intersections, which is achieved by the mathematical technique of treating the kinematical parameters as variables and using a block monomial ordering for which [variables] [parameters]. This trick could also be helpful for other types of multi-loop multi-scale amplitude computations. After solving the module intersection problems, we find linear IBP systems on cuts which are of very small sizes. This part is implemented in our highly efficient and automated Singular code. For example, the analytic IBP vectors for the hexagon-box integral with no doubled propagators on triple cuts can be computed in minutes.
Furthermore, we applied sophisticated sparse linear algebra techniques to compute the row reduced echelon form of the IBP identities on cuts. For example, we applied a weighted version of the Markowitz pivoting strategy to retain the sparsity in intermediate steps of the row reduction. We have implemented the sparse linear algebra part of our algorithm as a preliminary Mathematica code. A more efficient Singular implementation will become available in the near future.
In this paper, we have solved a cutting-edge IBP reduction problem fully analytically: that of the reduction of non-planar five-point hexagon-box integrals, with numerators of degree four, to the basis of 73 master integrals. Our result has been verified (numerically) with the state-of-art IBP reduction programs.
We are currently preparing an automated implementation based on open-source software, such as Singular, with which we expect to be able to solve yet more difficult IBP reduction problems. We expect that our method will boost the computation of NNLO cross sections for more complicated 2 → 3 scattering processes and higher-multiplicity cases.
Finally, we also expect that the ideas presented in this paper, such as the module intersection for trimming integral relations, the special ordering for variables and parameters, and sparse linear algebra techniques, can also be combined in various ways with other existing computational methods such as: the finite field sampling and reconstruction approach [14,66], the dual conformal symmetry construction of IBP vectors [15] and the newly developed D-module method [19], in order to determine tailored optimal reduction strategies for various classes of Feynman integrals.