Feynman integral reduction using Gr\"obner bases

We investigate the reduction of Feynman integrals to master integrals using Gr\"obner bases in a rational double-shift algebra Y in which the integration-by-parts (IBP) relations form a left ideal. The problem of reducing a given family of integrals to master integrals can then be solved once and for all by computing the Gr\"obner basis of the left ideal formed by the IBP relations. We demonstrate this explicitly for several examples. We introduce so-called first-order normal-form IBP relations which we obtain by reducing the shift operators in Y modulo the Gr\"obner basis of the left ideal of IBP relations. For more complicated cases, where the Gr\"obner basis is computationally expensive, we develop an ansatz based on linear algebra over a function field to obtain the normal-form IBP relations.

family of integrals. The integers z i are called the indices of the integral. A given physical process is usually expressed in terms of integrals of several different families. Note that we have suppressed the integral's dependence on the kinematic invariants, since we are here mostly concerned with the dependence on the indices. For a modern account on the calculation of loop integrals see for example [2].
In a typical calculation, one has to evaluate thousands of loop integrals belonging to several integral families. An indispensable tool in the calculation of multiloop integrals is therefore the integration-by-parts (IBP) method [3,4], which provides relations between loop integrals with different propagator and numerator powers. These recurrence relations are called IBP relations and can be used to express all integrals of a given family in terms of a small number of so-called master integrals. Nowadays this is usually applied in the form of Laporta's algorithm [5], which solves the system of linear equations generated by plugging in numerical values for the indices z 1 , . . . , z n .
In order to generate a linear system of equations for Laporta's algorithm one has to specialize a set of IBP relations to a range of indices z 1 , . . . , z n . Typically this linear system contains a large number of integrals (as unknowns) that are oftentimes not directly needed but must be included to ensure the full reduction of the desired integrals. This problem can be (partially) avoided by starting with a set of so-called unitarity-compatible IBP relations based on syzygies [17][18][19][20][21][22] (over a polynomial ring). We refer to these IBP relations in Section 6 as special IBP relations. They reduce the size of the linear system and improve the performance. Additional new ideas towards a more direct reduction procedure have been developed: They rely on algebraic geometry [23][24][25] and intersection theory [26][27][28][29][30][31][32][33][34].
One limitation of Laporta's algorithm is that the reduction is only found for a given list of integrals. Thus, when additional integrals are needed at a later point, the program has to be run again. This can be overcome by deriving a full solution to the system of IBP recurrence relations. Since a parametric solution by hand is clearly not feasible for multiscale problems, it would be desirable to have an algorithmic way of solving the IBP relations once and for all. LiteRed [35] is a publicly available program that performs this task using tailored heuristics that is able to reduce several complicated integral families. In our work, we investigate the application of Gröbner bases to the solution of IBP recurrence relations.
Previous application of Gröbner bases in this context can be found in [36][37][38][39][40][41]. In [36], the IBP relations are first transformed into a system of partial differential equations for which the Gröbner basis is then computed. This method requires that all propagators have different, nonzero masses and that no external momentum squared is equal to one of the masses squared. Thus, it is not always possible to apply the result to cases with zero or equal masses or onshell momenta, since such limits can be singular. Reference [38] uses a modified version of Buchberger's algorithm to obtain so-called sector bases [39]. Finally we would like to emphasize that -in contrast to previous noncommutative Gröbner basis approaches -we work in the rational double-shift algebra defined in Section 2.
This paper is organized as follows. In Section 2 we define the (rational) double-shift algebra in which the IBP relations form a left ideal. In Section 3 we introduce the first-order normalform IBP relations and highlight in which sense they differ from other well-known sets of IBP relations. Furthermore we introduce the notion of a first-order family, i.e., an integral family for which the left ideal of IBP relations is generated by the first-order normal-form IBP relations. In Section 4 we define the notion of formally scaleless monomials and relate them to scaleless sectors. In Section 5 we demonstrate on selected examples of first-order families the computation of Gröbner bases, normal-form IBP relations, and the detection of scaleless sectors using Gröbner basis reductions. Section 6 recalls the construction of special IBP relations using syzygies in a polynomial ring. The special IBP relations turn out to provide a more efficient set of generators as a starting point for the Linear Algebra Ansatz in Section 7 to compute the normal-form IBP relations without precomputing a Gröbner basis. Finally we conclude in Section 8.
and define the Lorentz invariant quadratic expression For p = ℓ the vectors ℓ 1 , . . . , ℓ L will refer to the L loop momenta. For p = k the vectors k 1 , . . . , k E will refer to the E external momenta.
Consider the polynomial algebra Q[d, m 2 i ] with coefficients in the field Q of rational numbers. Its elements are polynomial expressions in the dimension symbol d and the symbols of squared masses m 2 i . Define the field of rational functions where the numerators Z and nonzero denominators N are polynomials in Q[d, m 2 i ]. Consider the Lorentz invariant expressions that are polynomial expressions in the scalar products of the L+E momenta ℓ 1 , . . . , ℓ L , k 1 , . . . , k E with coefficients in F. Each such expression can be written as a polynomial in the n = L(L+1) 2 + LE propagators P 1 , . . . , P n and socalled extra Lorentz invariants S 1 , . . . , S q with coefficients in F. The extra Lorentz invariants are constructed from the external momenta in such a way that the set {P 1 , . . . , P n , S 1 , . . . , S q } is algebraically independent over F. This means that P 1 , . . . , P n , S 1 , . . . , S q generate a polynomial algebra over F, which we denote by Since from some point on we do not need the special form of the P i 's and S j 's we replace them by symbols D 1 , . . . , D n and s 1 , . . . , s q , respectively. Likewise we replace the polynomial algebra T by the isomorphic polynomial algebra For more mathematical details on the construction of the polynomial algebras T and R see Appendix A.
The IBP relations are obtained from the fact that the operator ∂ ∂ℓ µ i v µ i turns the loop integrand into a divergence, i.e., annihilates the loop integral in dimensional regularization. More precisely: The standard IBP relations are obtained by C running through the standard basis of T L(L+E)×1 (see (2.26) below).
In the following we will rewrite the expression ∂v µ i ∂ℓ µ i in (2.7) and the differential operator v µ i ∂ ∂ℓ µ i in (2.8) in terms of the ring R. To this end we define the IBP-generating matrix 1 as the product matrix where J := ∂Pc ∂ℓ µ i ∈ T n×Ld ′ is the Jacobian matrix of the propagators, and where T is defined in Appendix A. Like the propagators, and unlike the Jacobian matrix, the entries of the IBPgenerating matrix belong to the subring T and can therefore be effectively rewritten as matrices over R ∼ = T using the subalgebra membership algorithm. The latter can be replaced by simple linear algebra due to the affine nature of P i as expressions in p i · p j . The dimensions of E are already independent of d ′ . However, its entries as expressions in the generators of the subring T formally still depend on d ′ . But once E is rewritten as a matrix over R, the initial dependency of E ∈ R n×L(L+E) on the dimension d ′ disappears 2 .
For the coefficients vector C ∈ R L(L+E)×1 consider the Jacobian (2.12) and the square matrix The divergence summand in (2.7) becomes 1 The name is motivated by equation (2.24). 2 Physically, d ′ should be thought of as the symbolic regularizing dimension d rather than an integer. Furthermore, the second summand (2.8) To determine the action of the differential operation v µ The next section introduces the shift algebra which contains the IBP relations as shift operators. with the relations (no summation over repeated indices) The prefix "double" refers to the simultaneous occurrence of both the lowering operators D i and the raising operators D − i . The action is partial since D − i cannot be applied to a scaleless integral. 3 Our choice of the right action will be justified in Remark 2.3 and the definition of scaleless integrals is deferred to Section 4.

2.3.
Generating the left ideal of IBP relations. For an arbitrary coefficients vector C ∈ R L(L+E)×1 we get, using the above summation convention, the IBP (shift) operator Note that due to the Jacobian expression entering E C the map is not R-linear, but merely linear over the subalgebra F[s 1 , . . . s q ] < R.
The L(L + E) standard IBP relations are obtained by C running through the standard basis {e 1 , . . . , e L(L+E) } of R L(L+E)×1 , i.e., The left ideal of IBP relations annihilates all loop integrals of the given family, i.e., IBP whenever the partial action is defined. Since the map r in (2.25) is not R-linear one needs to formally prove that the left ideals I pol IBP and I IBP are finitely generated, more precisely: Proposition 2.2. The left ideals I pol IBP and I IBP are generated by the standard IBP relations: The right hand side of (2.30) results in the IBP operator Dictated by the loop diagram, some of the propagators play the role of numerators, i.e., only their nonpositive exponents z i are considered. These u many propagators are called irreducible numerators and are conventionally grouped at the end: However, annihilators of partial right (left) actions are merely left (right) ideals. But since the software we use for computing noncommutative Gröbner bases only supports left ideals, we had to opt for partial right actions.

GRÖBNER BASES IN THE NONCOMMUTATIVE DOUBLE-SHIFT ALGEBRAS
3.1. Gröbner bases, standard monomials, and master integrals. Below we need the notion of a Gröbner basis of the left ideals I pol IBP and I IBP in the respective noncommutative algebras In both cases we use generalizations of Buchberger's algorithm [42] to the context of GR-algebras and Ore algebras, respectively.
Replacing a set of generators of a left ideal by a Gröbner basis (with respect to a monomial order) might introduce redundant generators. However, these are necessary for the reduction procedure to produce unique remainders, independent of possible choices of the reduction steps.
Let G pol and G denote the Gröbner bases of the left ideals I pol IBP ¢ Y pol and I IBP ¢ Y , respectively. As customary we denote by NF G pol (f ) ∈ Y pol and NF G (g) ∈ Y the normal forms of f ∈ Y pol and g ∈ Y with respect to the Gröbner bases G pol and G, respectively.
One would generally expect the number of elements in G to be smaller or equal to that of G pol . This might fail in trivial cases as in Example 3.3. More involved cases like Example 5.2 show that the difference of cardinalities can be significant. For Gröbner basis and normal form computations in the polynomial double-shift algebra Y pol we use SINGULAR's subsystem PLURAL [43] and for the rational double-shift algebra Y we use Chyzak's Maple package Ore_algebra [44]. For the technical implementation we developed the package LoopIntegrals [45]. LoopIntegrals is currently written in GAP [46] and relies on the homalg-project packages [47] which offer a unified interface to SINGU-LAR [48] and Maple. A MATHEMATICA package that can perform the required Gröbner basis computations over the rational double-shift algebra is HolonomicFunctions [49,50]. The homalg-project does not offer an interface to MATHEMATICA yet.
The L(L + E) = 1 standard IBP relation is expressed as element of the polynomial double-shift algebra One can verify that the cyclic generator r 1 is already the reduced Gröbner basis G pol = {r 1 } of the left ideal I pol IBP ¢ Y pol . Switching to the rational double-shift algebra Y and computing the reduced Gröbner basis G of I IBP ¢ Y we get the two generators

simple computation reveals that the Gröbner basis reductions modulo
In particular, the normalized IBP relation r 1 2m 2 takes the special form: This special form of the IBP relation also reflects the functional dependence of the integral's closed-form result: Using and Γ(z + 1) = z Γ(z) we obtain It is obvious from (3.8) that the special form of the IBP relation -regardless of its characterization using Gröbner bases-is ideally suited for performing IBP reductions. Finally, Hence, the set of standard monomials with respect to G is merely is a K-linear combination of the standard monomials, the first-order normal-form IBPs.
The adjective "first-order" refers to the linear occurrences of D − i and D i in the monomials in (3.12) for which the normal forms are to be computed.
Examples (with u = 0) show that NF G pol (a i D − i ) with respect to the polynomial Gröbner basis G pol includes expressions in the D − j 's. In contrast, in the examples treated in Section 5, the normal forms NF G (a i D − i ) with respect to the rational Gröbner basis G do not involve any D − j . However, the K-coefficients of NF G (a i D − i ) will often be true fractions in Remark 3.5. We expect the numerators num(R i ) ∈ Y pol to lie in I pol IBP . This can be verified for each specific example by checking that all num(R i ) reduce to zero modulo the polynomial Gröbner basis G pol ⊂ Y pol . This proves a posteriori that the reduction modulo G used to compute the R i 's can be obtained without dividing by polynomials in Q[a 1 , . . . , a n ]. In particular, the normal-form IBP relations are valid for all indices (z 1 , . . . , z n ) ∈ Z n .
We call the matrix τ a certificate of polynomiality.
To verify this it suffices to compute the (normalized) reduced minimal Gröbner basis G ′ of . . .

num(Rn)
over Y and check that G ′ = G. Alternatively one could verify that reduces to zero modulo G ′ .
Definition 3.9. A bookkeeping of this reduction yields a matrix η ∈ Y L(L+E)×n such that We call the matrix η a certificate of first-order generation.
Due to Remark 3.7 one cannot expect that τ and η are mutually inverse over Y , even when L = 1 and n = L(L + E). means that all integrals of a sector have the same factors P i in the denominator, possibly raised to different positive powers z i . In a given sector with subset of positive indices U ⊆ {1, . . . , n − u} the integral with z i = χ U (i) is called the corner integral of the sector, where χ U is the characteristic function of U as a subset of {1, . . . , n}. This means, the corner integral of a sector has no numerator and all factors in the denominator have power one Let I(z 1 , . . . , z n ) (with nonnegative z i 's) be an L-loop integral which does not factorize into a product of lower-loop integrals. Setting the speed of light to unity, this integral has the physical dimension of the mass to the power Ld − 2 z i . After the integration, only masses and kinematical invariants can carry this mass dimension. If no such quantities are available after integration, either because they are set to zero by the external kinematics or because the (evaluated) integral does not depend on them, the integral is considered to be scaleless. Scaleless integrals are set to zero in dimensional regularization. If the L-loop integral I(z 1 , . . . , z n ) factorizes in a product of lower-loop integrals, then I(z 1 , . . . , z n ) is considered scaleless if any of its factors is scaleless.
In an integral family scaleless integrals typically appear in sectors with few propagators (note that terms in the numerator cannot provide a mass scale). It is well-known that if the corner integral of the sector is scaleless, then all integrals of this sector are scaleless as well. In such a case we call the sector a scaleless sector. Note that all subsectors of a scaleless sector are also scaleless.
For an algorithmic identification of scaleless sectors by means of linear algebra the Symanzik polynomials of the loop diagram (often denoted as U and F) can be used. Our package LoopIntegrals [45] uses ideas of [51] based on [52], see also [53] for a similar algorithm.
Motivated by computations in the noncommutative rational double-shift algebra we suggest the following definition and the conjecture below: is scaleless.
It hence suffices to look for the formally scaleless monomials among those of the form D i 1 1 · · · D in n for i j ∈ {0, 1}. These correspond to the corner integrals of the respective scaleless sectors defined by the subset {j | i j = 1}. The first two conditions ensure that the transformation leaves the evaluated integral invariant (modulo a global factor given by the Jacobian | det M LL |). However, the integrand might change due to the transformation and with the last condition we ensure that the transformed integrand can be expressed as linear combination of integrands of the same integral family. Note that for the performance of Laporta's algorithm it is important to include the information on the scaleless sectors and the symmetries while generating the linear system since it reduces the number of variables and the size of the system. The n = 2 propagators are The following transformation is a symmetry of the loop integral: where the lines indicate the different block matrices in (4.3) and the dot denotes an entry that is equal to zero.
The L(L + E) = 2 standard IBP relations are , expressed as elements of the polynomial double-shift algebra The reduced Gröbner basis G pol for the one-loop bubble over the polynomial double-shift algebra Y pol has 4 elements and was computed in less than a second using PLURAL: The reduced Gröbner basis over the rational double-shift algebra was computed in less than a second using Ore_algebra. It also has 4 elements and reads We can now compute the normal forms of the operators a i D − i with respect to the Gröbner basis G of the left ideal I IBP := ⟨r 1 , r 2 ⟩ ¡ Y (5.9) generated by the above two standard IBP relations: The C 2 -symmetry of the problem is obvious from the above normal forms. It further manifests itself in the normal forms of the indeterminates D i : The set of standard monomials with respect to G is merely Computing normal forms with respect to the Gröbner basis G one can easily verify that the minimal scaleless monomials D 1 , D 2 are indeed formally scaleless with respect to I(1, 1) in the sense of Definition 4.1.
As discussed in Section 3 the normal form of a i D − i with respect to the polynomial Gröbner basis G pol ⊂ Y pol includes expressions in the D − j 's. Still, the polynomial reduction verifies that for i = 1, 2 the numerator num( These reductions yield as in Definition 3.6 the certificate of polynomiality matrix The following nontrivial row-syzygy 16) in (Y pol ) 1×2 of ( r 1 r 2 ) shows that (5.15) does not uniquely determine τ . The computation of the Gröbner basis of N in Y pol verifies that N does not generate I pol IBP . However, r 1 , r 2 reduce to zero modulo N in Y , yielding the certificate matrix of first-order generation from Definition 3.9, proving that ⟨R 1 , R 2 ⟩ Y = I IBP := ⟨r 1 , r 2 ⟩ Y As mentioned in Remark 3.10, the numerators of the normal-form IBP relations R 1 , R 2 cannot lie in the image of the map r, since it is obvious from (2.24) that the IBP relations in the image of r are affine expressions in the a i 's.
We can now verify that the normal forms in (5.10) and (5.11) reflect the functional dependence of the integral's closed-form result on the indices z 1 and z 2 , analogously to (3.8). Using (5.18) and Γ(z + 1) = z Γ(z) we have, for example, the contiguous function relation which corresponds to NF G (a 1 D − 1 ) in (5.10). Note that 1 is always a standard monomial, unless the left ideal is the unit left ideal. This happens for example for the scaleless one-loop bubble (with s = k 2 1 = 0). There G pol has five elements but G = {1}, i.e., I IBP = Y . This means that the set of standard monomials is empty (NF G (1) = 0), i.e., there are no master integrals as expected for a scaleless integral family (contrary to what is claimed in [54]).
Example 5.2 (One-loop box). The one-loop box is defined by the loop momentum ℓ 1 and linearly dependent external momenta k 1 , k 2 , k 3 , k 4 . Due to momentum conservation k 1 + k 2 + k 3 + k 4 = 0 we take k 1 , k 2 , k 4 to be the set of linearly independent external momenta (in particular, L = 1, E = 3). The external lines are on-shell and massless implying k 2 i = 0 for i = 1, 2, 3, 4. Internal lines are also massless. This results in the independent external kinematic invariants s 12 = 2k 1 · k 2 and s 14 = 2k 1 · k 4 .
The following two involutions generate a Kleinian symmetry group V 4 of the loop integral:

(5.22)
Their product yields For completeness, we mention that additional symmetry relations can be found when some of the indices are negative, e.g., for encodes the symmetry  The reduced Gröbner basis G pol for the one-loop box over the polynomial double-shift algebra Y pol has 28 elements and was computed in less than a second using PLURAL. The reduced Gröbner basis over the rational double-shift algebra , s 12 , s 14 )(a 1 , a 2 , a 3 , a 4 ) was computed in less than five seconds using Ore_algebra. It has 9 elements and reads (5.29) with the abbreviations a i 1 ...i k := k j=1 a i j . We provide both G and G pol in computer-readable form in the ancillary files.
We can now compute the normal forms of the operators a i D − i with respect to the Gröbner basis G of the left ideal generated by the above four standard IBP relations: The normal forms of the indeterminates D i reveal the V 4 -symmetry of the problem: The set of standard monomials with respect to G is Computing normal forms with respect to the Gröbner basis G one can easily verify that the minimal scaleless monomials D 1 D 2 , D 1 D 4 , D 2 D 3 , D 3 D 4 are indeed formally scaleless with respect to I(1, 1, 1, 1) in the sense of Definition 4.1. The above computations can be found in the notebook [55].
As discussed in Section 3 the normal form of a i D − i with respect to the polynomial Gröbner basis G pol ⊂ Y pol includes expressions in the D − j 's. Still, the polynomial reduction verifies that for i = 1, . . . , 4 the numerator num( shows that the matrix τ is not uniquely determined. Furthermore, we have verified that the (normalized) reduced minimal Gröbner bases of {r 1 , . . . , r 4 } and {num(R 1 ), . . . , num(R 4 )} over Y coincide, proving that the latter (or equivalently {R 1 , . . . , R 4 }) generates the left ideal The rational Gröbner basis G can be used to perform fast reductions. The ancillary files of the arXiv submission include a FORM [56] program and a MATHEMATICA program which compute normal forms modulo G. In order to express the integral I(z 1 , . . . , z 4 ) in terms of the above three master integrals (5.34) the program computes the normal form of the monomial For example, the provided FORM-program is able to express I(10, 10, 10, 10) in terms of the master integrals in about 5 seconds.
We also provide a MATHEMATICA program which uses the four normal-form IBPs . . , 4} ⊂ I IBP ¡ Y to perform the reduction. This program illustrates that the normal form IBPs can be used to construct an algorithm for the reduction to master integrals similarly to the standard IBPs. However, due to the much simpler structure of the normal form IBPs, it is much easier to design the reduction algorithm. In particular, they already have the correct form to reduce integrals in the top-level sector, where all indices z i with i ∈ {1, . . . , n − u} = {1, . . . , 4 − 0} are positive.
The symmetric group S 3 is a symmetry group with the following three transpositions as generators: The L(L + E) = 4 standard IBP relations are (5.40) expressed as elements of the polynomial double-shift algebra The reduced Gröbner basis G pol for the two-loop tadpole with three massive lines over the polynomial double-shift algebra Y pol has 47 elements and was computed in a few seconds using PLURAL. The reduced Gröbner basis over the rational double-shift algebra was also computed in a few seconds using Ore_algebra and has 9 elements. Both G and G pol are included in the ancillary files.
We can now compute the normal forms of the operators a i D − i with respect to the Gröbner basis G of the left ideal I IBP := ⟨r 1 , . . . , r 4 ⟩ ¡ Y (5.43) generated by the above four standard IBP relations: (5.44) The S 3 -symmetry of the problem is obvious from the above normal forms. (5.45) NF G (1) = 1 Y 1 is (trivially) a standard monomial, The set of standard monomials with respect to G is of which the last three are equal. It is interesting to note that the homogeneity of the integral with respect to m 2 manifests itself through the relation: 48) or equivalently the IBP operator Computing normal forms with respect to the Gröbner basis G one can easily verify that the minimal scaleless monomials D 1 D 2 , D 1 D 3 , D 2 D 3 are indeed formally scaleless with respect to I (1, 1, 1) in the sense of Definition 4.1.
As discussed in Section 3 the normal form of a i D − i with respect to the polynomial Gröbner basis G pol ⊂ Y pol includes expressions in the D − j 's. Still, the polynomial reduction verifies that for i = 1, 2, 3 the numerator num( IBP . These reductions yield as in Definition 3.6 the certificate of polynomiality matrix The following nontrivial row-syzygy shows that (5.52) does not uniquely determine τ .
The computation of the Gröbner basis of N in Y pol verifies that N does not generate I pol IBP . However, we have verified that the (normalized) reduced minimal Gröbner bases of {r 1 , . . . , r 4 } and {num(R 1 ), num(R 2 ), num(R 3 )} over Y coincide, proving that the latter (or equivalently {R 1 , R 2 , R 3 }) generates the left ideal I IBP ¡ Y .

THE SPECIAL IBP RELATIONS
In the previous section we computed the set of normal-form IBP relations {R i | i = 1, . . . , n} starting from the set of standard IBP relations {r i | i = 1, . . . , L(L + E)}. However, it turns out to be more efficient to start with a different set of IBP relations introduced in [18,22].
Recall that the standard IBP relations are obtained from (2.24) Remark 6.2. A column syzygies matrix S ∈ R c A ×c S of A ∈ R r×c A modulo B ∈ R r×c B can be computed by computing a column syzygies matrix S T ∈ R (c A +c B )×c S of the augmented matrix (A|B) ∈ R r×(c A +c B ) and then discarding the lower part T . Most computer algebra systems with a Gröbner basis engine provide procedures to compute S directly without (fully) computing T .  In all examples where we were able to compute the (reduced) Gröbner bases of the standard IBP relations and the special IBP relations in the noncommutative polynomial double-shift algebra Y pol they coincided, i.e., we were able to prove that This means, the set of special IBP relations generates the same left ideal in Y pol as the standard IBP relations. This is remarkable since ⟨S⟩ is a proper R-submodule of R L(L+E)×1 .
The reduced column syzygy matrix of in R 4×6 producing 6 special IBP relations in Y pol ⊂ Y . They are provided in an ancillary file attached to the arXiv submission. As mentioned above, a Gröbner basis computation in Y pol verifies that the 6 special IBP relations generate I pol IBP ¡ Y pol .

THE LINEAR ALGEBRA ANSATZ
We now describe a method for producing the normal-form IBP relations without computing the Gröbner basis for I IBP ¡ Y . This method is based on linear algebra (LA-Ansatz).
We start with a generating set M = {g 1 , . . . , g t } of the left ideal I IBP ¡ Y , e.g., the set of standard IBP relations or the set of special IBP relations. For a fixed order o ≥ 0: • Multiply M by (D − 1 ) i 1 · · · (D − n ) in for i j ≥ 0 and i 1 + · · · + i n ≤ o from the left. . . , a n ] ⊂ K := F(s 1 , . . . , s q )(a 1 , . . . , a n ).   One can simulate the Gaussian elimination over K by running the algorithm over the subring A and dividing by the content of the resulting rows of reduction steps. This relies on fast multivariate gcd computations for which we use the Julia package HECKE. The complement of the subring Q[a 1 , . . . , a n ] < A is a multiplicative subset of A. Strictly speaking one must carry out the elimination in the localized subring (A \ Q[a 1 , . . . , a n ]) −1 A < K. We leave this for future work.
For the on-shell kite we have the IBP-generating matrix (cf. (2.11)) The reduced column syzygy matrix S of E ′ ∈ R 6×13 modulo diag(D 1 , D 2 , D 3 , D 4 , D 5 ) ∈ R 5×5 is producing 13 special IBP relations in Y pol ⊂ Y . The matrix S and the corresponding special IBP relations are provided in ancillary files attached to the arXiv submission. Since we were not able to compute G pol we were not able to verify that the 13 special IBP relations generate I pol IBP ¡ Y pol .
integral family. The latter makes this approach well-suited for parallelization. However, the bottleneck of this approach is the computation of the noncommutative Gröbner basis G of I IBP ¢ Y . Indeed, the existing implementations failed to produce the Gröbner basis for the on-shell kite. The bottleneck in computing the Gröbner basis G motivated the Linear Algebra Ansatz developed in Section 7 which enabled us to compute the first-order normal-form IBP relations even when we were not able to compute G. Here it turned out that the special IBP relations of Section 6 yield a more efficient set of generators of I IBP than the standard IBP relations. The special IBP relations are in this sense intermediate between the standard and the normal-form IBP relations.
Examples in Section 5 clearly demonstrate, why the normal-form IBP relations -regardless of their definition using Gröbner bases-are ideally suited for IBP reductions. We have provided as an ancillary file a proof of concept MATHEMATICA program which uses the normalform IBP relations to perform reductions for the integral family of the one-loop box.
Preliminary experiments show that the use of normal-form IBP relations in Laporta's algorithm improves the sparsity of the system of linear equations. We leave the systematic combination of both approaches for future work.
While we focused on first-order integral families in the work, we expect that integral families with more than one master integral per sector are "higher-order" generated. It would be interesting to have a simple characterization for first-order generation. We also leave this for future work.
Finally, we mention that the main ideas of Sections 2 and 3 have also been discussed in a somewhat more pedagogical form in [57]. are F-linearly independent elements of the (q + n)-dimensional F-linear subspace T 2 < T generated by the quadratic expressions p i · p j (for p ∈ {ℓ, k}), where q depends on the relations ρ o . This means, that (the representative in T of) each P c is an affine polynomial in the quadratic expressions p i · p j with constant term either a Q-multiple of m 2 i ∈ F or a rational number. Let S 1 , . . . , S q , P 1 , . . . , P n (A.5) be a basis of T 2 . The representatives of S 1 , . . . , S q in the (ambient) algebra F[ℓ 1 , . . . , ℓ L , k 1 , . . . , k E ] can be chosen homogeneous of degree 2 and are called the extra Lorentz invariants. Together with the masses they form the kinematic invariants of the process. Then T can equally be expressed as the subring having T as its image in T .