Macaulay Matrix for Feynman Integrals: Linear Relations and Intersection Numbers

We elaborate on the connection between Gel'fand-Kapranov-Zelevinsky systems, de Rham theory for twisted cohomology groups, and Pfaffian equations for Feynman integrals. We propose a novel, more efficient algorithm to compute Macaulay matrices, which are used to derive Pfaffian systems of differential equations. The Pfaffian matrices are then employed to obtain linear relations for ${\cal A}$-hypergeometric (Euler) integrals and Feynman integrals, through recurrence relations and through projections by intersection numbers.


Background
Within the perturbative approach to quantum as well as classical field theory, the evaluation of multiloop Feynman integrals is necessary for the determination of scattering amplitudes and related quantities (see [1] for a recent review). The by-now standard evaluation techniques of Feynman integrals (in momentum-space representation) exploit loop-momentum shift invariance to establish integrationby-parts (IBP) relations [2,3] among integrals whose integrands are built from products of the same set of denominators (and scalar products), but raised to different powers. IBP identities play a crucial role in the calculation of multi-loop integrals because they can be used to identify a minimal set of elements, dubbed master integrals (MIs), that can be used as a basis for the decomposition of multi-loop amplitudes. At the same time, IBP relations can be exploited to build systems of equations solved by the MIs: differential equations [4][5][6][7][8][9][10][11][12][13], dimensional recurrence relations [14,15], and finite difference equations [3,16] are linear relations emerging from the IBP reduction of integrals which spring from the action of special (polynomial and differential) operators on the integrands of the MIs. Solving such equations amounts to the actual determination of the MIs themselves, as an alternative to direct integration.
The derivation of IBP-decomposition formulas requires the solution of a system of linear relations, generated by imposing that integrals of a total differential (w.r.t. integration variables) vanish on the integration boundary [3,17]. For multi-loop, multi-scale scattering amplitudes, solving the system of IBP relations may, however, represent a formidable task. This problem has motivated the development of important techniques based on reconstruction of rational functions using finite fields [18][19][20][21].
Already during the early developments of S-matrix theory, it was recognized that topology and cohomological methods offer a connection between analytic properties of Feynman integrals and the geometry of their singularity structure, which emerge from the varieties associated to their graph polynomials. In more recent studies, intersection numbers for twisted de Rham cohomology groups [22][23][24][25][26][27][28][29][30][31][32][33] have been exploited to uncover the vector space structure of Feynman integrals and to derive novel algorithms for the direct projection of multi-loop integral onto MIs [34][35][36][37][38] (see also [39][40][41][42]). Important developments have been carried out in [43][44][45][46][47][48][49] and example of applications to non-linear relations can be found in [50,51]. See [1,42,52,53], for recent reviews. Within this approach, the number of MIs corresponds naturally to the dimension of the vector space of Feynman integrals, and can be related to topological quantities such as the dimension of the cohomology groups, the number of certain critical points [34,54], Euler characteristics [37,[55][56][57][58], as well as to the dimension of quotient rings of polynomials for zero dimensional ideals (Shape lemma) [38]. Linear relations (for integrals with shifted indices), Pfaffian systems of differential equations, finite difference equations, as well as quadratic relations, such as Riemann twisted period relations, can be derived by means of intersection numbers for differential (twisted) forms.
The pivotal role that twisted de Rham theory seems to have in controlling the algebra of Feynman integrals, and more generally of Euler integrals, has been stimulating us to elaborate on the isomorphism between cohomology groups, whose elements are differential forms, and D-modules, whose elements are partial differential operators in a Weyl algebra.
Within momentum-space representation, differential operators acting on multi-loop integrands and integrals are used to establish IBP relations, system of partial differential equations (SPDE) for MIs, and Lorentz invariance identities [9]. In the former type of relations, differentiation is carried out w.r.t. the integration momenta, while for the latter two ones, w.r.t. the momenta of external particles and masses. Also, within parametric representations of Feynman integrals, it is possible to study the action of differential operators acting on multi-fold integrands and integrals, distinguishing between the cases for which the differentiation is carried out w.r.t. integration or external variables.
In [57,58], it was shown that linear relations between Feynman integrals having shifted indices, similar to IBPs, can arise from the action of parametric annihilators of the integrand, within the Lee-Pomeransky (LP) representation. The algebraic properties of these special partial differential operators, w.r.t. the integration variables, are derived by studying the ideals of parametric annihilators within the language of D-modules. Remarkably, the number of master integrals, known to be finite [59], was found to be identical to the Euler characteristic of the complement of the hypersurface determined by the LP polynomial [57,58].
While the previous study focused on D-module theory for differential operators in the internal variables, in this work, we investigate the properties of differential operators in the external variables. In particular, we exploit the properties of GKZ-hypergeometric systems, introduced by Gel'fand, Kapranov, Zelevinsky [60], also known as A-hypergeometric systems, to derive Pfaffian equations for the generators of the corresponding D-module. As established in [61], GKZ systems provide a dictionary between Euler integrals and differential operators, which was later made fully algorithmic in [62]. Feynman integrals can be considered as special cases of Euler integrals, therefore the Pfaffian matrices derived in the context of D-module theory correspond to the matrices of the SPDE satisfied by the MIs [63].
The connection between Feynman integrals and the generalized hypergeometric functions was first proposed by Regge [64] and it was better understood with the knowledge of the Feynman integrals satisfying a holonomic differential equation [65], where the singularities of the differential equation were governed by the Landau singularities. The relation between Feynman integrals and hypergeometric system were studied in [66][67][68][69][70][71][72] (and reference therein). The isomorphism between the GKZ system and the Feynman integral was established in [63] and later realized within the LP representation in [73,74].
Within the proposed approach, many problems, including the derivation of differential equations, are translated into a ring theoretic computation where we may utilize various notions of computational ring theory such as Gröbner bases. As Gröbner basis calculations may be computationally expensive, we propose a different approach based on the Macaulay matrix, which is of central interest in this work. Usually, relations among integrals are employed to derive systems of partial differential equations. Instead, in the current work, we reverse the perspective, and show how Pfaffian matrices, built from Maculay matrices, can be used to derive linear relations for integrals. Moreover, through the secondary equation introduced in [33], we exploit the role of Pfaffian matrices for building cohomology intersection matrices to be used in the master decomposition formula presented in [34-37, 39, 43], for the direct integral decomposition.

Macaulay matrix for GKZ system
GKZ systems can be rewritten into systems of differential equations dubbed Pfaffian systems. By extension, the same holds true for Feynman integrals. A Pfaffian system for a holonomic function f in variables z = (z 1 , . . . , z N ) is of the form ∂ i F = P i · F, F = (s 1 f, . . . , s r f ), i = 1, . . . , N , (1.1) where ∂ i := ∂/∂z i , s i are differential operators acting on f , and P i are r × r matrices with entries being rational functions of z. The P i are called Pfaffian matrices (or simply Pfaffians) in this paper.
Our presentation is organized as follows. In Section 2, we review basic notions of the GKZ hypergeometric systems and their Euler integral representation. In Section 3, we discuss Pfaffian systems of differential equations, which are intimately related to GKZ systems. We present the Macaulay matrix algorithm, based only on linear algebra, to compute Pfaffian matrices in Section 4. We show its application to examples of differential equations for Feynman integrals in Section 5. In Section 6, we show how Pfaffians can be used to derive linear relations for GKZ systems, similar to IBP identities for Feynman integrals. Finally, in Section 7, we present the integral decomposition via intersection numbers, using Pfaffians to compute the required intersection matrices.
All algorithms in this paper are implemented in the computer algebra system Risa/Asir [85], Maple [86] and Mathematica [87] with FiniteFlow [88], while the calculations involving Feynman integrals are checked with LiteRed [89,90]. Programs used in this paper and machine readable data can be obtainable from [91].

GKZ hypergeometric systems
In this section, we briefly review some basic properties of the GKZ-hypergeometric systems to fix our notation. Section 2.1 introduces a particular integral representation related to the GKZ systems we work with, and Section 2.2 covers its relation to the algebraic de Rham cohomology groups. In Section 2.3, we describe how to represent a cohomology class by an element of Weyl algebra. Finally, in Section 2.4, we discuss the homogeneity property of GKZ systems, which allows us to reduce the number of independent variables.

Integral representation of GKZ-hypergeometric system
In this work, we consider Euler integrals of the form Here Γ is a twisted cycle 2 , β = (β 0 , . . . , β n ) ∈ C n+1 are complex parameters, and g(z; x) is a Laurent polynomial in x The monomials above are written in multivariate exponent notation: given an integer vector α i ∈ Z n we set where α i,j stands for the j-th component of the vector α i . Crucially, in (2.2) we regard each coefficient z i as an independent variable of f Γ (z). Let us construct the (n + 1) × N matrix whose columns a i are built from the monomial exponents α i as a i := (1, α i ), with the assumption that Span{a 1 , . . . , a N } = Z n+1 . Moreover, we introduce the (left) kernel of A, defined as, Then, by using A and β as input, we build the following set of differential operators: The function f Γ (z), defined in (2.1), satisfies the system of partial differential equations (PDE) therefore it is dubbed an A-hypergeometric function [60].

GKZ D-modules and de Rham cohomology
The operators in (2.6)-(2.7) can be regarded as elements of a Weyl algebra (2.10) In multivariate exponent notation, the elements of D N take the form k∈K h k (z)∂ k for some finite are polynomials in z with complex coefficients. The symbol ∂ i is an alias of ∂ ∂zi . We introduce the GKZ system as the left D N -module D N /H A (β), where H A (β) is the left ideal generated by E j and u , (2.11) Further details on D-modules theory can be found in the Appendix B. Let us list a few important properties of GKZ systems and their relations to de Rham cohomology groups, which are expressed through the following theorems and propositions.
First, we recall a theorem on the number of solutions to GKZ systems. Let ∆ A denote the convex polytope spanned by the columns of A. We say β is non-resonant when it does not belong to any set of the form span For example, if we take a 2 × 2 matrix The holonomic rank equals the number of independent solutions to the system of PDEs (2.6)-(2.7) at a generic point z ∈ C N . The first statement ensures that the rank is finite, while the second statement gives an exact formula for computing it in terms of combinatorial data.
Next, letting G m (resp. A) stand for the complex torus (resp. complex Affine line) equipped with the Zariski topology 5 and we denote by π : X → Y the natural projection from the space of GKZ and integration variables to the space of GKZ variables only.
. . x −βn n denotes the left D N -module O(X) endowed with this action. Formally, we have the identities Let us make this isomorphism explicit [95]. We let represent the space of relative k-forms 6 onto which we act with the covariant derivative in integration variables We then obtain a chain complex (2.21) The k-th relative de Rham cohomology group is defined as follows: is isomorphic to the n-th relative de Rham cohomology group H n , for which reason the latter is a left D N -module by Theorem 2.2. In fact, Theorem 2.2 can be rephrased as Proposition 2.3. Suppose that β is non-resonant. Then there is a unique isomorphism of D N -modules 3, which will be essential for our application of D N -module theory to Feynman integrals, is the following: given a cohomology class [ω(z)] ∈ H n , there exists a differential operator P ∈ D N , which is unique modulo H A (β), such that (2.24) The partial differential operators ∂ i in P act on a cohomology class [ω(z)] ∈ H n via The action (2.25) comes from differentiation under the integral sign: Since a Feynman integral can be represented by a cohomology class [34], we may equally well consider the operator P as representing that integral. An algorithm for computing P was developed in [62] and will be outlined in the following section. Moreover, in the view of the relation of GKZ-systems and Feynman integrals (to be elaborated on in Section 5), we observe that the finiteness of the rank, established by the first statement of Theorem 2.1, can be related to the finiteness of the number of master integrals [59]. The formula for its evaluation, given in the second statement of Theorem 2.1, offers an alternative way of determining dim(H n ), which is ordinarily computed in terms of Betti numbers, by counting the number of certain critical points, or by Euler characteristics -all of which are related to the number of master integrals. 6 X/Y does not denote a quotient space. Rather, it is a symbol for relative differential forms, by which we mean that we do not consider differential forms dz J where z ∈ Y .

Representing integrals inside D N -modules
Let us define the following family of differential forms: Pairing ω q with a cycle Γ yields the integral (2.28) We observe that differentiating the integral w.r.t z i shifts the q vector by +a i (defined after (2.4)): It is also possible to construct an operator which does the opposite, i.e. shifting q by −a i . In [96] it was shown that there exists a so-called step-up (or creation) operator U i and polynomial b i (β) such that (2.30) The operator U i then acts on the integral (2.28) via The polynomial b i (β) is called the b-function [96]. Its computation is algorithmic and is implemented in computer algebra packages such as mt gkz.rr [62]. Now, assume an integral of the form (2.28) is given for some choice of the integer vector q. Since the a i 's span Z n+1 by assumption, we may write (2.32) The coefficients r i ∈ Z allow us to balance the shifts by ±a i in the q of ω q Γ via equations (2.29) and (2.31). We may therefore construct the operator with the property that The two functions B(β) and B (β) are built from the prefactors on the RHSs of (2.29) and (2.31), and their explicit expressions can be found in [62]. The operator P (q) will be used extensively in our applications to Feynman integrals, in which case the vector q will be related to the propagator powers and space-time dimension of a given integral. Let us observe that the vector (r 1 , . . . , r N ) in (2.32) is not unique, because we can shift it by adding a term like i c i u i for some choice of c i ∈ Z, and u i ∈ Z N span Ker(A). This freedom can be exploited to simplify the expression of the operators P (q) appearing in (2.34), and for practical purposes, we find it convenient to choose r i > 0 for all i, in order to avoid the appearance of U i in (2.34), because it may contain monomials in ∂ i of large degree.

Homogeneity and integrand rescaling
Under rescaling of the z i variables, Euler integrals scale as where we used the multivariate exponent notation from (2.3) on the LHS. Indeed, (2.35) is equivalent to homogeneity equations E j f Γ = 0 , j = 1, . . . , n + 1 .
This property can be exploited to work on a simpler set of integrals obtained from the original definition by freezing (n + 1) out of N variables z i to e.g. 1. The residual functional dependence is then given by N − (n + 1) ratios of z i -variables. More details can be found in Appendix A.
Example 2.4. Consider the following Euler integral: which, according to eq. (2.1), corresponds to a 3 × 4 matrix: The homogeneity property from (2.35) involves three rescaling parameters {t 0 , t 1 , t 2 }: To derive this relation we need to rescale integration variables as The homogeneity property allows us to factorize the dependence on z 1 , z 2 , and z 3 by the following choice of the rescaling parameters t i : The original Euler integral from (2.37) then becomes and so we have effectively rescaled away (n + 1) of the z i variables. Finally, we show the differential equation obeyed by the rescaled Euler integral f Γ (1, 1, 1, z). The generator (2.7) of the original GKZ system with the A-matrix shown in (2.38) reads as Upon substitution of the representation (2.41) on the RHS, we can determine what the two terms in u map to when changing variables to w: These expressions can also be derived from Proposition A.1. Thus the original equation (2.42) becomes which we recognize as the differential equation for the Gauß hypergeometric function.
This concludes our short overview of the basic properties of GKZ systems and Euler integrals. In the following section, we describe how to obtain a Pfaffian system of first-order PDEs given a basis of Euler integrals.

Pfaffian systems
In this section, we introduce notation related to Pfaffian systems to be used in the rest of this work. In Section 3.1, we recall the passage from a holonomic D N -module to a Pfaffian system, while in Section 3.2 and Section 3.3 we present two methods for basis change.
For any positive integer N , we write C(z) = C(z 1 , . . . , z N ) for the field of rational functions in variables z 1 , . . . , z N . Let R := C(z 1 , . . . , z N ) ⊗ C[z1,...,z N ] D N be the rational Weyl algebra 7 ; note that R also has the structure of a non-commutative ring. Any element Q ∈ R can be uniquely written in the so-called normally ordered form with all the partial derivatives commuted to the right: (3.1) As we review in Appendix B, given a Gröbner basis G of a left ideal I in R, we denote by 8 Std the standard monomials w.r.t G and a term order ≺. For our purposes it is enough to consider a zero-dimensional left ideal I generated by a finitely many elements: (3. 2) The ideal I is zero-dimensional iff the corresponding set of standard monomials Std is finite (see [97, §6]).

From D N -module to DEQ system
We consider a left R-module M which is finite dimensional over C(z). Such a left R-module arises as an extension R ⊗ D N M for some other holonomic D N -module M . When {u 1 , . . . , u r } ⊂ M is a C(z)-basis, one can find an r × r matrix P j (z) with entries in C(z) such that We call P j (z) the Pfaffian matrix in direction z j . Let us now consider the GKZ D N -module. Namely, we consider a left R-module M = R ⊗ D N (D N /H A (β)), where H A (β) is as in (2.11). In view of Proposition 2.3, each u i corresponds to a cohomology class [ω i ] ∈ H n . Integrating (3.3) over a cycle Γ, we obtain an identity which is precisely the system of differential equations that the master Euler integrals are subject to. A standard method of computing the Pfaffian matrix P j (z) is to use a Gröbner basis of R ([97, Chapter 6]). Unfortunately, this approach is far from effective. We develop a different approach based on the Macaulay matrix in Section 4. Our construction of Pfaffian matrices by the Macaulay matrix method, to be given in Section 4, will assume that a set of standard monomials S for GKZ ideal RH A (β) is given. Note that the set S gives rise to a C(z)-basis of M . Although one usually computes a Gröbner basis of the ideal RH A (β) to find the set S, we claim that the set S can be found by computing a Gröbner basis of an ideal in a commutative subring of R. This reduces the cost of computation in a significant way. 7 The subscript of the tensor product ⊗ denotes the ring over which the product is defined. For example, it implies the identity: r ⊗ p Q = r p ⊗ Q, for r ∈ C(z), p ∈ C[z], and Q ∈ D N . 8 When considering Std as a basis, we instead write e (Std) .
For this purpose, we fix a term order ≺ on the commutative ring C[∂ 1 , . . . , ∂ N ]. It naturally induces a term order ≺ on the ring R because a term order on the ring R is a total order on the set of monomials ∂ k (cf. [97, §6.1]). We denote by in ≺ (I ) the ideal of C[∂ 1 , . . . , ∂ N ] generated by initial terms of I = u by the order ≺, where the box operators are given in (2.7). Any generator of in ≺ (I ) can be written as i . Let us write θ i for the Euler operator z i ∂ i and consider a subring of R generated by θ 1 , . . . , θ N over C. We will denote this ring by C[θ], which is clearly a commutative ring. By the correspondence θ k ↔ ∂ k , we equip a term order ≺ on the ring C[θ] induced from ≺. We then define the distraction of the element ∂ k ( [93, p.68] We consider the ideal I of the ring C[θ] generated by the distractions of the generators of in ≺ (I ) and the E 1 , . . . , E n+1 from (2.6). The following theorem is essential whose proof is based on a technique called Gröbner deformation as we will see in Appendix B.
Theorem 3.1. The set of standard monomials of I coincides with that of GKZ ideal RH A (β) when β is generic. More concretely, if {θ a1 , . . . , θ ar } is a set of standard monomials of I , the set {∂ a1 , . . . , ∂ ar } is a set of standard monomials of RH A (β).
The ideal I is known as toric ideal. Computation of a Gröbner basis of such an ideal is effective as it is generated by binomials. Thus, the computation of a set of generators of the ideal I is efficient. Running Buchberger's algorithm for the commutative ideal I , we obtain a set of standard monomials of GKZ ideal RH A (β) by Theorem 3.1. We emphasize that the process explained above does not involve any computation in the non-commutative ring R. Finally, we remark that the process above is implemented as a function cbase by euler in the Risa/Asir package mt gkz.rr.

Example 3.2. Let us consider a matrix
Let ≺ be the lexicographic order such that ∂ 1 · · · ∂ 4 . The ideal in ≺ (I ) is generated by monomials The distractions of them are We compute the Gröbner basis of an ideal generated by these distractions and Thus, the set of the standard monomials of GKZ ideal RH A (β) is given by

Basis change
Assume that we know Pfaffians P (Std) i in a basis of standard monomials e (Std) . In this section we show how to gauge transform the system into an arbitrary basis e.
Consider an ideal I over a rational Weyl algebra R. Let a basis e = (e 1 , . . . , e r ) T ⊂ R/I be given. In the following, the symbol • denotes multiplication inside of R, i.e. we write ∂ i • h when h ∈ R. We introduce • to distinguish composition in the ring from the action ∂ i H when, say, H is a matrix with entries in the field C(z).
We seek the Pfaffian matrices P i ∈ C(z) r×r satisfying the differential equation for the basis e: where e (Std) = {∂ α } ⊂ R/I are the standard monomials. If we could find a matrix G ∈ C(z) r×r relating the two bases, then we could relate the Pfaffian matrices in different bases via a gauge transformation: Note that G is invertible by assumption, since both e (Std) and e are bases. By (3.9), our problem is reduced to finding the transformation matrix G. To this end, we begin by expanding e in terms of e (Std) : Here we relied on the fact that any monomial of e is divisible by an entry of e (Std) . Next we transform each ∂ k • e (Std) into the form G (2) k · e (Std) via repeated application of (3.7). Let us illustrate the last step with an example.
We hence have that G We conclude that the gauge transformation matrix G decomposes as follows: where each G (2) k is built from the matrix products between Pfaffians and derivatives thereof. Thus, equations (3.9, 3.14) constitute one way to translate the Pfaffian matrices between different bases.

Basis change without derivatives
The matrix G (2) k appearing in the gauge transformation matrix (3.14) involves derivatives of Pfaffian matrices. If we were to work over a finite field, there would be no z i variables with respect to which we could take the derivative. This situation is amended, in the GKZ case, by the following proposition.
Proposition 3.4. The derivative of a Pfaffian matrix can be computed solely via addition and matrix multiplication from be a solution to a GKZ Pfaffian system given some monomial basis The following relations hold true: where we omitted the argument β for clarity. On the other hand, differentiating (3.18) we get where we applied the identity P k (β)F (β) = F (β − a k ) twice in the last step. The proposition follows upon equating (3.20) and (3.23) and isolating ∂ i P j (β).
Pfaffian systems introduced above are systems of linear partial differential equations (SPDE), satisfied by the solutions of a given GKZ system. As we will see later on, these equations are extremely useful in physical applications. Next we present an efficient way to calculate the Pfaffian systems, essentially via linear algebra.

Constructing Pfaffian systems from Macaulay matrices
In this section we describe a method for building the Pfaffian systems defined in eq. (3.3). The method amounts to first building an auxiliary matrix M called the Macaulay matrix, and then solving a special system of linear equations. In Section 4.1 we derive the Macaulay matrix (4.5) and the linear system (4.12, 4.13) that it satisfies. In Section 4.2 we then present Algorithm 1 for calculation of Pfaffian systems. In Section 4.3 we give several remarks about the algorithm and its efficiency. We close this section with several examples in Section 4.4, showcasing the steps and runtime statistics of the algorithm in practice.

From Pfaffian to Macaulay matrix
We present how the Macaulay matrix arises from a Pfaffian system in the basis of standard monomials. Since we will focus our discussion on the case of GKZ systems later on, based on the comments at the end of Section 3.1 we may safely assume that a set of standard monomials Std := {∂ k } is given, and that its size equals the holonomic rank |Std| = r, defined in (2.12). We remind that ∂ k denotes a monomial in derivatives, while ∂ i denotes a single derivative w.r.t. z i . Now, the Pfaffian matrix P i in the direction z i , defined in (3.3), is specified by the following relation: where • denotes composition in a ring R, α ∈ {1, . . . , |Std|} and I is a general holonomic ideal for now. Using the expression (3.2) for I, we carry over the relation (4.1) from the quotient R/I to the whole ring R as where each Q αj is an element of R and j ∈ {1, . . . , d} labels the ideal generators h j from (3.2). Now we bring the operators Q αj into normal ordered form as in (3.1): where q αjk ∈ C(z) |Std|×d×|Der| is a newly introduced rank 3 tensor 9 , and we denote by Der the set of all partial derivatives ∂ k appearing on the RHS. Now we are ready to introduce the main player of this section: the Macaulay matrix. It arises due to the normal ordered form of the first sum in (4.2). Substitution of the Q αj operators from (4.3) back into (4.2) results in the action of ∂ k ∈ Der on the ideal generators h j , which we bring to normal ordered form as follows: Here we introduced a rank 3 tensor M kjL ∈ C(z) |Der|×d×|Mons| and a set Mons, which collects the partial derivatives ∂ L emerging on the RHS. Upon concatenation 10 of the first two indices k and j, the rank 3 tensor M kjL turns into the Macaulay matrix where the combined index R := (kj) runs over a product of sets R ∈ Der × {1, . . . , d}. Plugging equations (4.3, 4.4, 4.5) back into the main relation (4.2) for the Pfaffian, we observe that the first term can be written as a product of matrices: where C ∈ C(z) |Std|×(|Der|·d) stands for the coefficient matrix C αR := q α(jk) of (4.3), and we regard the set Mons as a column vector. The set Mons, containing all monomials in ∂ i , can be partitioned into two disjoint sets consisting of standard and, what we call, exterior monomials: This partition naturally induces a column partition of the Macaulay matrix: with the two column blocks M Ext ∈ C(z) (|Der|·d)×|Ext| and M Std ∈ C(z) (|Der|·d)×|Std| . Similarly, the LHS of (4.2) can be decomposed as 11 where we introduced the coefficient matrix C ∈ C(z) |Std|×|Mons| in addition to its column blocks C Ext ∈ C(z) |Std|×|Ext| and C Std ∈ C(z) |Std|×|Std| , such that C = (C Ext |C Std ) in accordance with (4.8).
Finally substituting equations (4.6, 4.8, 4.10) into the main relation (4.2), we arrive at the following matrix relation: (4.11) Due to the linear independence of Ext and Std, this is equivalent to a system This system relates the Macaulay matrix M to the Pfaffian P i in the direction z i , and is central to our computational method presented next.

From Macaulay matrix to Pfaffian
Now we reverse the logic in the derivation above: although we started from (4.2), we would like to regard it as the goal. We view (4.2) and its matrix version (4.11) as a system of linear equations for the unknown matrix of coefficients C. Given the Macaulay matrix M , the idea is to solve first the system (4.12), and then substitute the solution into (4.13) in order to obtain the desired Pfaffian matrix P i . Let us now turn this idea into a practical algorithm. Say we would like to compute P i in direction z i given a zero-dimensional ideal I. First, given a D ∈ Z ≥0 we construct the set of all derivatives whose degree is bounded by D: (4.14) Next, similarly to (4.4), we write Mons D for the set of monomials in derivatives appearing in the normal ordered form of 11 Note that C , C Ext and C Std depend on the direction z i , since these matrices are born from the expression ∂ i • Std.
where we call M D the Macaulay matrix of degree D (compare to (4.5)). We adjust the maximal degree D and the sets Der D and Mons D until the linear system (4.12) can be resolved w.r.t. the unknown matrix C. Finally, the desired Pfaffian follows from using C in (4.13). We summarize this procedure in Algorithm 1. Mons ← Mons D ; Proposition 4.1. Algorithm 1 outputs the Pfaffian matrix P i in a finite number of steps. [97,Chap 6]). This means that the LHS can be expressed in terms of the generators {h 1 , . . . , h d } inside the ring R, as in (4.2). Let N 0 be the maximal degree of the monomials ∂ k ∈ Der from (4.3). Then Algorithm 1 outputs the Pfaffian matrix P i until the degree parameter D reaches N 0 .

Discussion
Some comments about Algorithm 1 are in order.
Other directions. The condition on line 4 can be extended by checking ∂ i • Std for all i = 1, . . . , N to generate a Macaulay matrix M D large enough for calculating Pfaffian matrices in all directions Rank of M Ext . The Macaulay matrix M Ext produced by the algorithm may not be of full rank. This makes the matrices M Ext , M Std smaller and therefore the whole algorithm more efficient than the algorithm of [78].
Degree of freedom. The unknown coefficient matrix C, which solves the system (4.12), can be chosen up to the left nullspace of the Macaulay matrix M . Indeed, suppose that both C = C (1) and C = C (2) solve the system (4.12). Clearly, we have (C (1) − C (2) ) · M Ext = 0. Uniqueness of the Pfaffian matrix alongside (4.13) imply (C (1) − C (2) ) · M Std = 0 as well. Therefore, we deduce according to the partition (4.8) of the Macaulay matrix.
Solvability condition. The system (4.12) has a solution when the rows of the C 1 matrix lie in the row space of the M Ext block of the Macaulay matrix. This can be formulated as a condition on the rank of the following matrix: Probabilistic method. Checking the solvability condition (4.18) can be computationally demanding. In order to reduce the computation time and memory usage we may employ the so-called probabilistic method. In the GKZ case, it amounts to setting the all the variables z i and parameters β (see (2.1)) to numbers in Q or a finite field F p , which greatly simplifies the row reduction required for the solvability check.
Smaller Macaulay matrix. Once the condition (4.18) is proved (almost surely, if one employs the probabilistic method), we further simplify the system (4.12) by selecting a subset of rows of the Macaulay matrix necessary for solving the system. In practice, we perform such a selection using the probabilistic method together with a determination of the left nullspace of the matrix on the LHS of (4.18).

Finite fields.
For fixed values of the variables z i and the exponents β, numerical values of the Pfaffian matrix P i can be efficiently obtained via the Macaulay matrix method. Note that one only needs numerical values of P i to obtain the values of normalizing constants for a certain class of statistical distributions, see [78]. In our applications, however, we would like to (re)construct the full rational dependence of P i on all the variables z i and β. When it becomes too demanding to solve the system (4.12) analytically, the probabilistic approach can be used in with tandem with rational function reconstruction over finite fields F p . Indeed, using modular arithmetic (over prime numbers) and the Chinese remainder theorem, the rational entries of P i ∈ C(z) r×r can be efficiently obtained by means of functional reconstruction algorithms for multivariate, rational functions, which were developed in the context of scattering amplitudes in [19]. This algorithm reconstructs analytic expressions for polynomials and rational functions via repeated numerical evaluation over finite fields by means of iterative application of Newton's polynomial representation and Thiele's interpolation formula. The reconstruction algorithm is implemented in the public code FiniteFlow [21], which uses so-called dataflow graphs for building numerical algorithms over finite fields.
Tests. Once we reconstruct the Pfaffian matrices from their numerical values over the finite fields, we check correctness as follows. Recall that Pfaffian equations (4.1) can be regarded as a Gröbner basis in R (see the discussion below eq. (1.2)). Therefore, if the solution rank of the left ideal (3.2) is known to be r, the correct r × r Pfaffian matrices should satisfy the integrability condition (1.2), and should also reduce the generators h i of the ideal (3.2) to 0 via the Pfaffian equations (4.1).
Further improvements. Let us comment on potential ways to improve the efficiency of the Algorithm 1. By virtue of (4.16), to generate the Macaulay matrix of degree D we are instructed to act with Der D , the set of all derivatives of degree ≤ D (see (4.14)), on the generators of ideal. The set Der D might be larger than necessary, so one way to improve the Algorithm 1 is to find a smaller set of derivatives needed for the construction of the Macaulay matrices. Another important direction of study is the prediction of singularities for the Pfaffian matrix in question, since such knowledge would greatly simplify the solution of (4.12) using rational reconstruction.
Setup Consider the following integral representation: and ∆ 2 is the two-dimensional simplex The ideal of the corresponding GKZ system is generated by three independent differential operators I = h 1 , h 2 , h 3 . However, to compute the Pfaffian matrices with generic parameters {α, β, β , γ} and {z 1 , z 2 }, it is enough to consider only the following pair: where the abbreviated coefficients are

Basis
The GKZ system in question has rank 3. The column vector of standard monomials reads Macaulay matrix For brevity, let us only consider the Pfaffian matrix P 1 in the direction z 1 . According to (4.1), our task is to express as a linear combination of Std over the rational function field Q(α, β, β , γ)(z 1 , z 2 ). Following the steps of Algorithm 1, we build the Macaulay matrix of degree D (4.16) as follows. First we act with Der = {1, ∂ 1 , ∂ 2 }, the list of all possible derivatives of degree D ≤ 1, on the two generators {h 1 , h 2 } from (4.22) Then we bring the ensuing expressions into normally ordered form (see (3.1)) and rearrange the result into a matrix. The rows of the matrix are labeled by the list of operators {h 1 , ∂ 1 h 1 , ∂ 2 h 2 , h 2 , ∂ 1 h 2 , ∂ 2 h 2 }, while the columns correspond to the monomials Mons appearing in the normally ordered form of these operators. Explicitly, we get the Macaulay matrix of degree 13 D = 1 where p i,j = ∂ zj p i and q i,j = ∂ zj q i . As dictated by (4.8), we decompose the Macaulay matrix into two blocks: the second block M Std has columns labeled by Std, and the first one M Ext by the rest of the derivatives Explicitly, we obtain the blocks where z − 2 := 1 − z 2 , and Pfaffian matrix We now turn to the linear system (4.12, 4.13). The coefficient matrices C Ext and C Std , defined in (4.10), are obtained by expressing (4.25) in terms of the sets (4.24) and (4.27): (4.30) 13 Had we used all three generators of the ideal, it would have been enough to consider the degree D = 0 Macaulay matrix. However, the shortened setup (4.22) considered here provides a better showcase for the steps of Algorithm 1.

Example 4.4. Consider the matrix
, whose solution rank is r = 33. Our choice of standard basis is Using the homogeneity propery discussed in Section 2.4, we may fix the following variables: z 1 = z 7 = z 8 = z 9 = z 10 = z 12 = z 13 = 1. The block M Ext of the Macaulay matrix of degree D = 2 turns out to be a sparse 945 × 958 matrix, whose rank is 534 by the probabilistic method. Therefore, there are many rows that are not needed for solving the system (4.12). Again, using the probabilistic method over the finite field Z/65537Z, we determine a maximal set of independent rows of the matrix M Ext . We find that a 534 × 958 matrix N Ext is enough to solve a smaller version of the system (4.12), namely C Ext = C · N Ext . As was mentioned in the Smaller Macaulay Matrix paragraph of Section 4.3, we can further reduce this new system as follows: since there are exactly 27 independent row vectors in C Ext , we may choose only a subset of the row vectors in N Ext , whose span includes the independent row vectors of C Ext . The Macaulay matrix is then obtained in 15.549 sec. It took 0.66 sec to calculate the Pfaffian matrix in the z 6 direction for fixed numerical values of variables z and parameters β. We leave the problem of full functional reconstruction of the Pfaffian matrices of this example for future work.
We have proposed an efficient Macaulay matrix method to construct Pfaffian systems. It can be applied to relatively large systems which are not feasible by Gröbner basis methods in the ring of differential operators. The latter methods are usually feasible for systems of up to, approximately, rank 10.

Macaulay matrix method and generalized Feynman integrals
The relation between Feynman integrals and GKZ-hypergeometric functions has been studied in [63,68,83,98,99]. Here, we show applications of GKZ systems combined with Macaulay matrices to Feynman integrals.
Let 0 < , δ 1, d 0 ∈ 2 · N, L ∈ N and ν := (ν 1 , . . . , ν n ) ∈ Z n . Moreover, fix the exponents of the Euler integral (2.1) to β = ( , − δ, . . . , − δ) − (d 0 /2, ν 1 , . . . , ν n ) . (5.1) We define a generalized Feynman integral as , |ν| := ν 1 + . . . + ν n . (5.4) The polynomial G takes the form and we set a i := (1, α i ) ∈ N n+1 0 as usual. The introduction of the δ parameter in (5.1) follows from the analytic regularization of [100,101] Using the notation of equations (2.27) and (2.28), we can equally well define the generalized Feynman integral as a pairing between a twisted cycle Γ and a twisted form In particular, we have The generalized Feynman integral reduces to the Lee-Pomeransky (LP) representation [54] of an L-loop integral in d = d 0 − 2 dimensions, with propagator powers ν i when: and where m i and p i stand for masses and external momenta, so that the polynomial G becomes the LP polynomial built from Symanzik (or graph) polynomials U, F, as G = U + F. 15 Let us observe that within the LP polynomial G(z; x), each z i may not necessarily be independent from each other. This is different from the polynomial g(z; x) appearing in the Euler integral f Γ (z), for which each z i is considered independent. Their independence ensures a non-degenerate correspondence between the monomials in x i variables and the partial differential operators in the z i variables -a crucial property to establish the isomorphism between twisted de Rham cohomology group and D-modules. This observation implies that Feynman integrals are restrictions of GKZ integrals, obtained from the latter by choosing suitable values of the variables z i [70,74,[80][81][82][83]. We note that the z i limits in (5.9) are generically not smooth at the level of the Pfaffian system for a GFI. In a forthcoming publication [105], we show that it is nevertheless possible to construct the Pfaffian system with all z i variables identified with their proper kinematic counterparts.

Massless one-loop diagrams
Let us consider a massless one-loop diagram G with n external legs and n internal edges. We use the symbol p i for each external momentum and x i for the Schwinger parameter of the i-th edge. The first and the second Symanzik polynomials read (1,2),...,(n−1,n), (1,n) σ ij x i x j , (5.11) 15 The limit δ → 0 of a generalized Feynman integral may yield to ill-defined expressions of the form In these cases, we adopt the replacement , with δ D (•) being the Dirac delta function, as proposed in [102][103][104].
where we set σ ij := −(p i + · · · + p j−1 ) 2 . Since the Lee-Pomeransky polynomial G = U G + F G consists of n 2 terms, it gives rise to an (n + 1) × n 2 A-matrix. In view of the discussion of Section 2.4, the number of reduced variables is which coincides with the number of independent Mandelstam variables s ij (as was also noted in [70]). In sum, we conclude that the master integrals for G are precisely subject to the GKZ system of PDEs. Moreover, the discussion above can be easily extended to the cases when external masses m 2 k = p 2 k are non-zero but treated as independent variables. For example, in the case of the one-mass n-gon diagram, all external masses are zero but one, say p 2 k = 0, therefore a new term −p 2 k x k x k+1 is added to the second Symanzik polynomial (5.11). The number of reduced variables, in this case, reads n(n−3) 2 +1 which is equal to {σ ij } ∪ {p 2 k }.

Example: One-loop box
In the following two sections, we illustrate the Macaulay matrix method for obtaining Pfaffian systems for Feynman integrals. We have chosen examples where the number of GKZ variables correspond to number of independent kinematic variables. The general case will be investigated in [105].
Setup Let us first derive the Pfaffian system for the one-loop massless box diagram. This example was studied in Lee-Pomeransky representation using twisted cohomology in [43]. The kinematics are The inverse propagators in momentum space are given by 16 (5.14) The corresponding generalized Feynman integral reads The s-dependent prefactor in (5.15) comes from rescaling integration variables as Within the LP-representation, the coefficients z i in (5.16) take the explicit form z 1 = · · · = z 5 = 1 , z 6 = t s =: z , (5.17) 16 Here 2 = 2 0 − 2 1 − · · · − 2 d .
leaving just one external kinematic variable z. This counting is compatible with the number of GKZ variables after using the homogeneity property discussed in Section 2.4. Indeed, out of the N = 6 variables z i appearing in the generalized Feynman integral, n + 1 = 4 + 1 = 5 of them can be rescaled, thereby yielding z 1 = · · · = z 5 = 1 and z 6 = z. Basis where Λ is a diagonal matrix containing prefactors: Applying the correspondence in (2.34) (or rather (A.10), which employs homogeneity) to the basis e (dR) , we obtain the D-module basis where ∂ := ∂ z .
Macaulay and Pfaffian matrices Using the Macaulay matrix method, we construct the Pfaffian system in the basis of standard monomials, Afterwards, we perform a gauge transformation to the basis e (D) , so that P can be read off the equation ∂ • e (D) = P · e (D) .
The Macaulay data for the system (5.22) is found to be M Ext = z 2 (z + 1) , (5.23a) The Pfaffian matrix is then = − (3δ + + 1) 2 + z (9δ + 2) 2 δ + 6δ + + 1 z 2 (z + 1) , Finally, according to the algorithm of Section 3.2, it is possible to build a suitable gauge transformation matrix G such that where the limit δ → 0 has been taken. The matrix (5.29) is canonical and in agreement with LiteRed.

Example: One-loop pentagon
Setup We consider a one-loop massless pentagon integral with one massive leg. The kinematics are Note that the identity (p 1 + p 2 + p 3 + p 4 ) 2 = (−p 5 ) 2 = p 2 imposes a relation on the Mandelstam variables s ij . The propagators are given by We can write the generalized Feynman integral as 33) and the corresponding A matrix A 5 is given in Example 4.3. The s 12 -dependent prefactor in (5.32) comes from rescaling integration variables by In the Lee-Pomeransky representation, the monomial coefficients are given by z 1 = . . . = z 6 = 1 , z 7 = 1 + y 2 + y 4 , z 8 = y 1 , z 9 = y 4 , (5.34) where the y i are ratios of kinematic variables: The Lee-Pomeransky polynomial hence contains 5 monomial coefficients different from unity. Similar to the previous example, we can obtain equally many non-unity GKZ variables via homogeneity. In particular, we start with N = 11 generic variables z i , after which we rescale n + 1 = 5 + 1 = 6 GKZ variables as z 1 = . . . z 6 = 1, leaving us with 5 variables z 7 , . . . , z 11 .
The dimensions d 0 are chosen so as to have non-negative r-vectors, as per the ending remark of Section 2.3. Using (2.34) (or formula (A.10) employing homogeneity), we find the D-module basis corresponding to e. We write e (D) = Λ · e (D) for a diagonal matrix Λ ij containing prefactors, , .
Macaulay and Pfaffian matrices The Macaulay matrix M Ext has dimensions 189 × 113, so in this case there are 113 exterior monomials Ext. The basis of 13 standard monomials is given by e (Std) = ∂ 9 ∂ 2 11 , ∂ 2 9 , ∂ 2 10 , ∂ 8 ∂ 11 , ∂ 9 ∂ 11 , ∂ 10 ∂ 11 , ∂ 2 11 , ∂ 7 , ∂ 8 , ∂ 9 , ∂ 10 , ∂ 11 , 1 T (5.40) We identify 133 independent rows of M Ext by row reducing it numerically. We therefore construct an unknown matrix C of dimensions 13×113 which must satisfy C Ext,i −C ·N Ext = 0, where N Ext contains the 133 independent rows of M Ext . This linear system can be solved in reasonable time on a laptop using FiniteFlow [88]. The Pfaffian matrix in direction i then follows from C Std,i −C ·N Std = P (Std) i , where N Std is built from M Std by taking the same rows as in N Ext .
Using the basis change algorithm outlined in Section 3.2, we perform a gauge transformation to relate P (Std) to P (e) i , written in the basis e. At this stage, we may safely take the limit δ → 0. We have verified that the resulting Pfaffian matrices are in agreement with LiteRed [89,90].
A similar analysis can be extended to diagrams with either more legs or more loops. For the purposes of this first work on this subject, we limited the application to one-loop integrals up to six external legs. In fact, Example 4.4 refers to the one-loop massless hexagon. Application to higher-loop cases will be discussed in future works.

Linear relations for generalized Feynman integrals
Let e = (I 1 , . . . , I r ) be a basis of master integrals, where each I i is a generalized Feynman integral of the form (5.2), i.e depending on generic variables z i . The vector e will be regarded as the column vector in the sequel; when the distinction between row and column vectors is clear from context, we will omit the transpose symbol • T .
Denote by e = (I 1 , . . . , I r ) a set of integrals in the same family as e. In both e and e , we replace the integrand factors x ν+ δ by x ν1+ δ1 1 · · · x νn+ δn n for a set of new parameters δ i . Theorem 6.1.
1. There exists a matrix U ∈ Q r ×r (z) such that e = U · e . (6.1) 2. Let L(e ) = ( 1 I 1 , . . . , r I r ), where the i are differential operators w.r.t. z with rational function coefficients. There exists a matrix V ∈ Q r ×r (z) such that Moreover, there are construction algorithms for U and V .
Remark 6.2. This theorem can be regarded as an analogue of IBP relations for generalized Feynman integrals. Recall the rank r for generalized Feynman integrals is possibly larger than for conventional Feynman integrals. However, once we set the GKZ variables z i equal to their physically relevant values in e.g. e = U · e, some of the master integrals in e might vanish or become equal, in which case we arrive at a conventional IBP relation.
Although the theorem can be proved by a method analogous to Algorithm 1 of [62], here we present a different approach based on Pfaffian matrices and the matrix factorial [84]. We begin by constructing recurrence relations in the general framework of GKZ hypergeometic systems, after which we specialize to the case of generalized Feynman integrals.

Recurrence relations for GKZ systems
For convenience, in this section we will work with a rescaled version of the Euler integral (2.1): where β is assumed to be non-resonant. We have the relation  (β), . . . , s r f (β) is then a solution to a Pfaffian system. Now, because ∂ i s j = s j ∂ i and ∂ i f (β) = f (β − a i ) we have that In other words, P i (β) yields a difference equation for F (β). Moreover, since β is non-resonant, the matrix P i (β + a i ) is invertible, in which case we also have To make the analogy clear between equations (6.5) and (6.6), let us introduce the operator ∂ −1 i acting as We define the falling matrix factorial as 8) and the rising matrix factorial as These matrix products can be swiftly computed using rational reconstruction over finite fields. When the subscript j is replaced by an integer vector κ ∈ N N 0 , we extend the definition of matrix factorials 18 to: κ j a j κ .
(6.11) Lemma 6.3. For κ ∈ N N 0 we have the recurrence relations Proof. By induction, let us derive that 14) The case j = 1 is shown in equation (6.5). Suppose that (6.14) holds for j − 1. Then we have Applying (6.14) to each (κ) j , we obtain the relation (6.12) for the falling matrix factorial. The case of the rising matrix factorial (6.13) is proved in a similar fashion.
, we note that the commutation relation P i (β − a j ) P j (β) = P j (β − a i ) P i (β) holds. The explicit matrix factorial representation of (6.12) is hence not unique. The same statement holds for (6.13).
Observing that F (β) i = s i f (β), we can write F (β) as where the vectors k i ∈ N N 0 are fixed according to s i = ∂ ki . Suppose that we want to obtain a recurrence relation for f (β − A · k 0 ) given some choice of k 0 ∈ Z n . Notice that f (β − A · k 0 ) is the first element of F (β − A · (k 0 − k 1 )). We propose to obtain the recurrence relation by Algorithm 2, whose correctness follows from Lemma 6.3. Note that the algorithm does not perform differentiation with respect to z i , meaning that it can still produce recurrence relations when the z i are fixed to numbers in the Pfaffian matrices.
Output the first element of step 3.

Recurrence relations for generalized Feynman integrals
We can employ Algorithm 2 to find relations among generalized Feynman integrals by specializing β to (6.4). The matrices Q i = P −1 i exist for this generic choice of and δ i 's. In certain cases, we can even take a limit δ i → 0 to obtain conventional IBP relations as we observe in the example below.

Example: One-loop bubble
Setup We use Algorithm 2 to obtain a linear relation for the one-loop bubble integral with one massive propagator (see also Example 1.1 of [74]), whose denominators are defined as (6.17) From the Lee-Pomeransky representation we read off We pick master integrals e = I(4, 1, 1), I(4, 2, 0) and wish to decompose e = I (4, 1, 2) . For this diagram, the choice δ 1 = δ 2 := δ is allowed. In order to apply Algorithm 2, we require k i corresponding to a relevant set of standard monomials. We use auxiliary vectors k i and α such that k i correspond to master integrals e. Their definitions are given by (6.28) and (6.29) in the proof of Theorem 6.1.
Step 0: Input The master integrals correspond to the choice 20) Given that these k i only have positive entries, we have α = (0, 0, 0, 0) using (6.29). Then the A · α terms drop, and we have using notation defined in (6.30). The vector k 0 represents the integral we want to reduce, I(4, 1, 2).
Step 3 Next, we are instructed to compute the matrix factorial corresponding to the operator ∂ +1 3 ∂ −1 1 , where the ∂ 3 stems from κ + and the ∂ 1 stems from κ − . Explicitly, we calculate P 3 (β + a 1 )Q 1 (β) . (6.25) Step 4: Output The recurrence relations follow from the 1st element of the matrix factorial above, apart from an adjustment of Γ-prefactors. This adjustment is due to the constants c and c (i) in equations (6.32) and (6.33c). Given all of the above, Algorithm 2 produces the result in agreement with LiteRed.
Critical case. We observe that our method may not give linear relations for some A and non-generic choices of and δ i 's. Put δ 1 = . . . = δ n := δ and take  [96]. In fact, we have Note that ∂ 3 and ∂ 4 stand for P i . We have b ( , − δ, − δ) = 0, for which reason Q i , i = 3, 4 acquires a 0 in the denominator. Since Q i acts similarly to ∂ −1 i , the recurrence relations involving Q i do not exist. Moreover, b(β) = 0 has the effect of breaking several isomphism theorems for D-modules (see, e.g., [96], [93, §4.5] and [106] for recent references on isomorphism theorems of GKZ systems).
We close this section with proving Theorem 6.1.
Proof. (Theorem 6.1.) We begin by noting that the matrices Q i exist since P i are invertible for generic parameters δ i .
Without loss of generality, suppose that r = 1. Then e = (I 1 ) is given by an integral I 1 = I(d 0 , ν) in the notation of Section 5. We week to reduce I 1 in terms of master integrals e = (I 1 , . . . , I r ), where we use the notation I i = I(d (i) 0 , ν (i) ). Take vectors k i ∈ Z N such that (6.28) We decompose k i into a sum of positive and negative parts: Then we may define the vector α = max and take k 0 ∈ Z N such that A·k 0 = A·α+(d 0 /2, ν). We apply Algorithm 2 using input k 0 , β = β +A·α for β = ( , − δ 1 , . . . , − δ n ) and s i = ∂ ki . The output of the algorithm takes the form Let us inspect the expressions for f on both sides of (6.30). On the LHS, 19 α is constructed in order to have a common shift for all k i such that their entries become non-negative.
So we obtain the generalized Feynman integral I 1 apart from constant prefactors c = c(d 0 , ν) defined in (5.4). On the RHS, In other words, we obtain I i apart from Γ-prefactors The coefficients u i (β) multiplied by Γ-prefactors give the matrix U . Thus, we conclude Statement 1. of Theorem 6.1.
Statement 2 can be proven by noting that ∂ i induces the parameter shift and applying Statement 1.

Decomposition via cohomology intersection numbers
Relations between Feynman integrals, equivalent to IBP identities, and, more generally, identities for Euler-Mellin integrals, equivalent to contiguity relations, can be derived by means of intersection theory for twisted de Rham cohomologies [34,35,37,38]. According to the mentioned algorithm, the decomposition of any given integrals in terms of an independent basis of MIs can be obtained from the projection of the twisted differential form appearing in the integrand of the integral to decompose into a basis of differential forms that generate a de Rham twisted cohomology group, via intersection numbers. For the case of generalized Feynman integrals (5.7), the covariant derivative (2.20) reads and we denote the associated n-th de Rham cohomology group as H n (see also (2.22)). We can also introduce a dual covariant derivative ∇ ∨ x = ∇ x →− and let H n∨ be the n-th (dual) de Rham cohomology group associated to it. The cohomology intersection number is a natural pairing between the elements of the two groups.
Let {e i } r i=1 be a basis for H n and {h i } r i=1 a basis for H n∨ ; the decomposition of any twisted form ϕ ∈ H n in terms of {e i } r i=1 can be obtained via chomology intersection numbers according to the master decomposition formula [34,35,37,38], This formula implies the decomposition of (generalized) Feynman integrals in terms of master integrals, upon the identification in (5.7).
Looking at (7.3) we infer that two distinct sets of intersection numbers are required, namely { ϕ, h i ch } r i=1 and {(I ch ) ij } r i,j=1 . Therefore, in order to apply the decomposition formula (7.3), it is required the determination of the matrix I ch and of the vector ϕ, h i ch .
According to the algorithm proposed in [62], given a basis , the cohomology intersection matrix (I ch ) ij = e i , h j ch can be obtained as a rational solution of the so called secondary equation, which is a (system of partial) differential equation(s) for the cohomology intersection matrix I ch , controlled by the Pfaffian matrices In particular, any non-zero rational solution of the secondary equation I is related to I ch as up to a constant κ 0 , corresponding to a boundary value, to be independently provided, see theorem 2.1 of [33].
Remarkably, also the vector ϕ, h i ch can be obtained by a second application of the same algorithm to a suitably chosen set of differential forms. For this purpose, it is sufficient to build the auxiliary basis 20 {ϕ}, as well as the auxiliary dual basis , and the corresponding Pfaffian matrices, say {P aux Then, the cohomology intersection matrix for the auxiliary bases, I aux ch , whose elements are defined as (I aux ch ) ij = e aux i , h aux j = e aux i , h j , can be obtained as rational solution of this matrix differential equation, The wanted vector of intersection numbers ϕ, h j ch , can be read off the r th -row of I aux ch , because (I aux ch ) rj = ϕ, h aux j ch = ϕ, h j ch , j = 1, . . . , r.
Actually the solution of the auxiliary secondary equation, I aux ch = κ aux 0 I aux (7.8) also requires the independent knowledge of the constant κ aux 0 . The proposed twofold procedure, yields the determination of all the cohomology intersection numbers needed in (7.3), up to the knowledge of the ratio of two boundary constants κ aux 0 /κ 0 . Let us observe that the master decomposition formula (7.3) requires the knowledge of the inverse matrix I −1 ch . By using the secondary equation (7.4), one can show that I −1 ch can be directly determined as rational solution of the following matrix differential equation, {ϕ} fails to be a basis when ϕ = 0. In the following discussion, we implicitly assume that the cohomology class ϕ is not zero and the set {ϕ} is a basis. This can be achieved when ϕ is represented by a monomial differential form ω d 0 /2,ν .
Determination of coefficients. Let us apply the master decomposition (7.3) to a vector formed by r differential forms e 1 , e 2 , . . . , e r−1 , ϕ, reading as, By construction, the product I aux ch · I −1 ch has the following form, (7.11) where I r−1 is the identity matrix in the (r − 1) × (r − 1) space, and the entries of the last row are the coefficients of the decomposition (7.3). On the other side, owing to the solution of the secondary equations, one has, This relation among matrices can be exploited to fix the value of the ratio κ aux 0 /κ 0 . In fact, any of the elements on the diagonal, say (I aux ch · I −1 ch ) kk , for any k ∈ {1, . . . , r − 1}, amounts to 1, therefore, yielding the relation κ 0 κ aux 0 = (I aux · I −1 ) kk , (7.14) namely fixing the ratio of the boundary constants from one of the elements of the product of two known matrices. Finally, (7.10) becomes, hence the coefficients of decomposition of the differential form ϕ in (7.3) reads as, for any chosen k ∈ {1, . . . , r − 1}. This result bypasses the determination of the complete cohomology intersection matrices, I ch and I aux ch , and relies only on the knowledge of the Pfaffian matrices, P and P aux , which control the secondary equations. The determination of their rational solutions, I and I aux , can be efficiently combined with the rational function reconstruction over finite fields, earlier discussed.
Let us remark that, if needed, the individual expressions of the normalization constants κ 0 , and κ aux 0 can be determined by means of [107, Theorem 8.1].

Conclusion
In this work, we presented a simplified algorithm for the determination of Pfaffian matrices from Macaulay matrix, making use of the theory of D-modules for Gel'fand-Kapranov-Zelevinsky systems. Then, we introduced two algorithms for the derivation of linear relations for GKZ-hypergeometric integrals making use of Pfaffians: i) the first one uses the properties of the matrix factorial, within the holonomic gradient method; ii) the second one uses the direct integral decomposition via cohomology intersection numbers and the rational solution of the secondary equation for cohomology intersection matrices. In the latter case, we derived a novel variant of the master formula for the projection of differential forms onto a set of independent forms.
Our investigation exploited the isomorphism between de Rham twisted cohomology groups, whose elements are differential forms, and D-modules, whose elements are partial differential operators in a Weyl algebra. Our results can be applied to GKZ-hypergeometric functions as well as to Feynman integrals, which can be considered restrictions of the former class of functions. Within our analyses the number of master integrals, corresponding to the dimension of the de Rham cohomology group, is related to the rank of the GKZ system, computed by polyhedra combinatorics in terms of volumes of polytopes. We showcased a few applications to simple mathematical functions and one-loop integrals.
Within the standard approaches, linear relations among integrals, such as contiguity relation for GKZ-hypergeometric functions and integration-by-parts identities for Feynman integrals, are employed to derive systems of partial differential equations for a chosen set of independent function. Instead, in the current work, we reversed the perspective, and showed how Pfaffian matrices for system of partial differential equations, built from Maculay matrices, can be used to derive linear relations for integrals.
We hope that our study could offer a novel, more complete view on integral relations which emerge from the application of differential operators acting both on internal and external variables, and could provide additional insights to the investigation of de Rham twisted co-homology theory in computational (quantum) field theory.
where A η θ w is a column vector, obtained from matrix multiplication of A η and a column vector θ w .
Let us close this discussion with one application of Proposition A.1, namely the action of a generic partial derivative operator on the cohomology class [dx/x]: (A.20)

B Holonomic D-modules in a nutshell
In this appendix, we summarize basics of the theory of holonomic systems utilized in our study of Feynman integrals. We cite introductory textbooks rather than the original papers or comprehensive textbooks.

B.1 Holonomic ideals
Let D be the Weyl algebra with polynomial coefficients: where the commutators are [z i , z j ] = 0, [∂ i , ∂ j ] = 1, and [∂ i , z j ] = δ ij . When we need to specify the number of variables, we also denote this algebra by D N . Any element L ∈ D can be written in the normally ordered form via the commutator relations. The principal symbol in (0,1) (L) of L is defined as the sum of the highest order differential operators in L: where |q| := q 1 + · · · + q N , and the commutative variable ξ denotes the derivative ∂ when a given expression is written in the normally ordered form. We see that the principal symbol is the sum of the highest weight terms with respect to the weight 0 for z and the weight 1 = (1, . . . , 1) for ∂. Example B.1. Take N = 2 and let I be a left ideal generated by z 1 and ∂ 2 . The ideal of principal symbols is then generated by z 1 and ξ 2 . The dimension of the zero set of z 1 = ξ 2 = 0 in C 4 is 2, and so we conclude that D/I is holonomic.

B.2 Standard monomials
Let R be the rational Weyl algebra When I is a holonomic ideal in D, the left ideal R I of R is a zero-dimensional ideal. In other words, the dimension of R/(R I) as a vector space over the rational function field C(z 1 , . . . , z N ) is finite. Similar to the polynomial case (B.2), any element L of the rational Weyl algebra R can be brought into the normally ordered form Let ≺ be a term order among the monomials in ∂ i 's. The largest monomial ∂ q in L (stripped of its rational function coefficient) is called the initial monomial in ≺ (L) of L. Monomials in derivatives naturally form the polynomial ring C[∂ 1 , . . . , ∂ N ]. A finite subset G of R is called a Gröbner basis of R I with respect to ≺ when the initial monomial ideal is generated by in ≺ (L), L ∈ G. A monomial ∂ α is called a standard monomial with respect to G when ∂ α does not belong to the initial monomial ideal (in ≺ (R I) in our case). Within the theory of Gröbner basis in R, the zero-dimensionality of R I is equivalent to the finitness of the set of standard monomials. In the following, an element ∂ α of in ≺ (R I) ⊂ C[∂ 1 , . . . , ∂ N ] is also denoted by ξ α to emphasize that it is an element of the polynomial ring. Let us now give a few examples of standard monomials for different ideals.
Example B.2. When N = 2 and I = z 1 , ∂ 2 ⊂ D, we have R I = 1 because 1 z1 z 1 = 1. Then the set of standard monomials is empty.
where β 0 and β 1 are complex numbers. These two operators annhilate the function (1 − z 1 z 2 ) −β0 z β1 2 . By a Gröbner basis computation in the Weyl algebra D, we find that the ideal of principal symbols is generated by the two elements ξ 1 − z 2 2 ξ 2 and z 1 ξ 1 − z 2 ξ 2 . Since it defines a 2-dimensional variety in C 4 , we conclude that the module D/ L 1 , L 2 is holonomic. The set {L 1 , L 2 } is a Gröbner basis in R with any term order ≺ and the initial ideal is generated by ξ 1 and ξ 2 . The set of standard monomials is simply {ξ 0 1 ξ 0 2 } ≡ {1}, corresponding to the point (0, 0) in the figure to the right, where the axes represent powers of ξ 1 , ξ 2 . Note that the second construction D ∩ J from J is a highly non-trivial. Algorithms realising such a construction are called Weyl closure algorithms.
We close this subsection with a note on the relation between Pfaffian equations and Gröbner bases. Let G be a reduced Gröbner basis of a zero dimensional ideal RI in the rational Weyl algebra R. Since it is zero-dimensional, the set of the standard monomials is a finite set. We denote it by {s 1 , . . . , s r }.
Computing the normal form of ∂ i s j by the Gröbner basis, we obtain The matrix P i = (P k ij ) j,k is the Pfaffian matrix for the operator ∂ i . Conversely, assume that we are given a basis S = {s 1 , . . . , s r } of R/RI as a vector space over C(z 1 , . . . , z N ). Here the s i 's are monomials. Let P i be the Pfaffian matrix for ∂ i . Then we have ∂ i s j = r k=1 (P i ) jk s k modulo RI. We assume that there exists a term order ≺ such that s ≺ t for all monomials s ∈ S and t ∈ S c where S c is the set of the monomials which do not belong to S. If such an order exists, it can be found by solving a system of linear inequalities [93, §2.1]. See also a note below. Then the set G = {g ij := ∂ i s j − r k=1 (P i ) jk s k | ∂ i s j ∈ S} is a Gröbner basis with repect to the term order ≺, and the set S is the set of the standard monomials for the Gröbner basis. Let us prove it. Since u 1 for any monomial u, we have ut t for any monomial t. Therefore, we have ut ∈ S c for t ∈ S c and any monomial u by the condition of the order. It implies that a Gröbner basis of the monomial ideal S c is {∂ i s j | ∂ i s j ∈ S}, and S is the set of standard monomials for that basis. It follows from the condition on the order that we have in ≺ (g ij ) = ∂ i s j ∈ S and in ≺ (G) ⊂ in ≺ (RI). Since dim C[∂ 1 , . . . , ∂ N ]/in ≺ (G) = r by in ≺ (G) = S c , in ≺ (G) = in ≺ (RI) holds. Thus, G is a Gröbner basis for the order ≺.
Note that it is necessary to pose the condition on the order and the basis S. For example, consider N = 1 and M = R/R(∂ 1 − z 1 ). S = {∂ 1 } is a vector space basis of M over C(z 1 ) and the Pfaffian system with respect to S is ∂ 1 − (z 1 + 1/z 1 ), because ∂ 1 (∂ 1 − z 1 ) = ∂ 2 1 − z 1 ∂ 1 − 1 and 1 z1 ∂ 1 − 1 belong to the ideal R(∂ 1 − z 1 ). However, the reduced Göbner basis of the ideal for any term order is {∂ 1 − z 1 }. In fact, the condition ∂ 1 ≺ 1 does not hold for any term order.
Finally, we explain a procedure to check if there exists a term order ≺ such that s ≺ t for all monomials s ∈ S and t ∈ S c . Let w ∈ R N ≥0 be a weight vector. We solve a system of linear inequlities with respect to w w · α ≤ w · β for all s = ∂ α ∈ S and t = ∂ β ∈ S where S is the set of the minimal generators of the monoid S c . If the solution space is an Ndimensional cone, then take w from the interior of the cone and define the order by ≺ w , where we may take any tie breaker ≺. Here, we mean ∂ a ≺ w ∂ b if and only if w · a < w · b or (w · a = w · b and ∂ a ≺ ∂ b ). If the cone is empty, there does not exist such an order. When the cone has dimension less than N , take w(1) ∈ R N ≥0 from the relative interior of the cone and also take the set P of (α, β)'s such that w(1) · α = w(1) · β. We solve a system of linear inequalities w · α ≤ w · β for all (α, β) ∈ P .
Take w(2) ∈ R N ≥0 from the relative interior of the solution cone. If this w(2) satisfies w(2) · α = w(2) · β for all (α, β) ∈ P , there exists no term order we want. If not, take ≺ w(2) as a tie breaker. We repeat this procedure. Since any term order can be expressed by a matrix of weight vectors [109], we can check if there exists a term order s ≺ t for all monomials s ∈ S and t ∈ S c and we can construct it if it exists.

B.3 Proof of Theorem 3.1
In this subsection, we prove Theorem 3.1. We follow the notation of Appendix 3.1. Let J be an ideal of C[θ]. It naturally lifts to a left ideal of R generated by the elements of J which we will denote by RJ. Recall that a term order ≺ on the ring C[θ] naturally corresponds to an order among monomials ∂ k by the correspondence θ k ↔ ∂ k . Thus, ≺ induces a term order on the ring R.
Lemma B.5. The set of ≺-standard monomials of the ideal J coincides with that of RJ through the correspondence θ k ↔ ∂ k .
Proof. Let us recall an identity z k i ∂ k i = θ i (θ i − 1) · · · (θ i − k + 1) (B.10) for any k ∈ N. For any k = (k 1 , . . . , k N ) ∈ N N , we set We claim that a ≺-Gröbner basis G ⊂ J of J is again a ≺-Gröbner basis of RJ. Let us take P, Q ∈ G and write them as P = aθ k + · · · , Q = bθ + · · · (B.12) where a, b ∈ C and underlines indicate the leading terms. Regarding P and Q as elements of the ring R, we obtain an expansion P = az k ∂ k + · · · , Q = bz ∂ + · · · . The second factor clearly belongs to the ideal J and the S-pair of P and Q is reduced to zero by G in the ring R. This shows that G is a ≺-Gröbner basis of RJ. Thus, the set of leading monomials of J is identical to that of RJ through the correspondence θ k ↔ ∂ k . This proves the lemma.
The following theorem is known as Gröbner deformation.
(B.15) Now, let us prove Theorem 3.1. Since z k is invertible in the ring R, Theorem B.6 shows that in ≺ (RH A (β)) is identical to the ideal RI . Clearly, the right-hand side is a lift of an ideal I in C[θ]. By Lemma B.5, we can conclude that the set of standard monomials of I is identical to that of in ≺ (RH A (β)) through the correspondence θ k ↔ ∂ k .

B.4 Integration and restriction
Now we discuss two important properties of the holonomic D-modules: the so-called integration and restriction constructions.
Let D M +N be the ring of differential operators in M + N variables: Finally, let us note that the integration and the restriction constructions commute: Restrictions of generalized Feynman integrals will be discussed in the forthcomming paper [105].