Feynman integrals as A-hypergeometric functions

We show that the Lee-Pomeransky parametric representation of Feynman integrals can be understood as a solution of a certain Gel’fandKapranov-Zelevinsky (GKZ) system. In order to define such GKZ system, we consider the polynomial obtained from the Symanzik polynomials g = U + F as having indeterminate coefficients. Noncompact integration cycles can be determined from the coamoeba—the argument mapping—of the algebraic variety associated with g. In general, we add a deformation to g in order to define integrals of generic graphs as linear combinations of their canonical series. We evaluate several Feynman integrals with arbitrary non-integer powers in the propagators using the canonical series algorithm.


Introduction
The analytic evaluation of Feynman integrals in dimensional regularization is still one of the main challenges to compute higher order corrections to observables in collider experiments. Methods for evaluating Feynman integrals involve a good understanding of their analytic properties. These have been important from the very beginning in order to develop techniques to evaluate them. Long time ago, some of these properties led to the recognition that Feynman integrals satisfy systems of differential equations analogous to those of hypergeometric functions 1 . The modern method to evaluate Feynman integrals is indeed based on differential equations. It is the combination of Integration by Parts (IBP) identities and differential equations [3][4][5][6][7][8]. The result is particularly simple [9] if the system of differential equations evaluates to combinations of Multiple Polylogarithms [10][11][12][13].
In favorable cases, Feynman integrals can be evaluated in terms of classical hypergeometric functions and their generalizations such as Appell, Horn, and Lauricella hypergeometric functions. Typically, these functions appear as infinite sums through the Mellin-Barnes representation of Feynman integrals (see e.g., [14][15][16]). Once evaluated in terms of hypergeometric functions, they can be utilized to recover the -expansion [17][18][19][20], although it is a nontrivial task when hypergeometric functions of many variables appear (see e.g., [21,22]). More recently, the Mellin-Barnes representation has been used to find systems of differential equations without resorting to IBP identities [23], through the differential reduction method [17,[23][24][25][26][27][28][29][30], which is based on Refs. [31,32]. Using the Mellin-Barnes representation it has been shown that a large class of Feynman integrals can be expressed in terms of Horn-type functions [33]. The differential reduction approach has led to the application of techniques coming from D-module theory as has been recently shown in Ref. [34]. This point of view is close to the one we will adopt in this paper.
A-hypergeometric functions were introduced by Gel'fand-Kapranov-Zelevinsky (GKZ) in 1990 [35] as a generalization of the well-known Appell, Lauricella, and Horn series. One important aspect of the theory is the study of polynomials with indeterminate coefficients associated with an integer matrix A [36]. This matrix and a vector of complex parameters furnish a system of partial differential equations (PDEs) known as a GKZ system. Solutions of these systems of PDEs are called A-hypergeometric functions. They can be represented as Euler-type integrals and hence these integrals define A-hypergeometric functions [35,[37][38][39]. On the other hand, series solutions can be computed by a generalization of the Frobenius method known as the canonical series algorithm due to Saito, Sturmfels, and Takayama [40]. GKZ systems and Feynman integrals have been object of recent interest to mathematicians, who have studied the relation among GKZ systems, Feynman integrals, and their regularization [41,42]. Maximal cuts in this language have been studied in Ref. [43].
In this paper, we will consider the parametric representation of Feynman integrals employed by Lee and Pomeransky to relate critical points of the sum of the Symanzik polynomials and the number of master integrals arising from IBPs [44]. This representation is based on the polynomial g = U + F, where U and F are the first and second Symanzik polynomials, respectively. We will show that the polynomial g defines a GKZ system, which is constructed by considering its coefficients to be indeterminate. We will use this information to obtain series expansions of the Euler integrals solutions using the Saito-Sturmfels-Takayama canonical series algorithm. Our definition of a GKZ system based on g allows the evaluation of integrals with arbitrary noninteger powers in the propagators using computational algebra. The polynomial g may lead to a matrix of codimension 0. In such cases, we introduce a deformation of g to ensure a canonical series representation. Once the canonical series are computed and integration constants are obtained, Feynman integrals can be recovered at the end of the computation by taking the limit of the deformation going to zero and setting its coefficients to their kinematic values.
This paper is organized as follows: In Section 2 we review the GKZ framework, present their solutions, and introduce the canonical series algorithm. We end this section with examples. In Section 3 we introduce scalar Feynman integrals and furnish a GKZ system based on them suited for canonical series. We give examples at the end of Section 3. We give our conclusions in Section 4.

A-hypergeometric functions and their representations
In this Section, we review basic aspects of the GKZ approach to hypergeometric functions. The main references are the book [40] and the lectures [45]. Reviews of the main concepts can also be found in Refs. [46,47].
This review contains four main parts. We first introduce polynomials with indeterminate coefficients (toric polynomials), which is one of the main ideas behind the GKZ approach. Polynomials with indeterminate coefficients lead to generalizations of discriminants, resultants, and determinants through the study of toric varieties 2 [36].
In the second part, we associate a system of partial differential equations(PDEs) to polynomials with indeterminate coefficients. We introduce Euler integral representations of solutions of GKZ systems. The system of PDEs may be formally defined as a holonomic ideal in a Weyl algebra D and it is the proper language for computational algebra purposes [40].
In the third part, we introduce series representations of solutions of GKZ systems and review the Saito-Sturmfels-Takayama algorithm to compute canonical series [40]. The fourth part contains examples of the methods introduced. A short example using the computer algebra system Macaulay2 can be found in Appendix B.
The connection with Feynman integrals will be made in Section 3. The Saito-Sturmfels-Takayama algorithm will be our main tool to evaluate Feynman integrals in Sec.3.

Notation
Throughout this paper we will employ multi-index notation, i.e., z α :=z α 1 1 · · · z α N N , c γ := c γ 1 1 · · · c γn n , where α ∈ K N , γ ∈ K n . For polynomials b 1 (z), . . . , b M (z) in the variables z = (z 1 , . . . , z N ), the multi-index notation for products of polynomials reads where β ∈ K M . Typically we will have K = C. For β = (1, . . . , 1), we simply write b(z) (1,...,1) = b(z). In the case M = 1, we write b 1 (z) β 1 = b(z) β . We call polynomials with indeterminate coefficients toric polynomials 3 and emphasize their dependence on their coefficients c by writing b(c, z). The multi-index notation for differential operators reads ∂ α = ∂ α 1 1 , · · · ∂ αn n , where ∂ i = ∂/∂c i . Euler operators are defined by Integrals where toric polynomials are involved will be identified by a subscript, e.g., I b indicates that the integral under consideration has a toric polynomial b(c, z) on its integrand. Finally, the Pochhammer symbol is defined by Useful identities following from this definition are summarized in Appendix A.

Polynomials, varieties, and their coamoebas
Let us consider q Laurent polynomials in N variables of the form where C * = C\{0}, α ij ∈ Z N , and n i is the length of the set of exponent vectors A i = {α i1 , · · · , α ik , · · · , α in i } associated with the polynomial b i (z). By abuse of notation, we denote by A i the N × n i configuration matrix of the exponent vectors of the i-th polynomial as where each (column) vector α ik is associated with a monomial term c ik z α ik in b i (z). Therefore |A i | = n i is the total number of monomials (columns in A i ) 4 . Let us give an example. Suppose we have a single polynomial in N variables with n monomial terms b ex (z) = n j=1 c j z α j , where we have labeled the coefficients simply by c j , j = 1, . . . , n. The N × n configuration matrix of b ex (z) reads The monomial term related with the first column reads With the above identifications of the columns and rows, the arrangement of the rows and columns (terms) is irrelevant as they define the same polynomial.
Let us now consider the product where each polynomial and its configuration A i are taken independently, in other words we do not expand b(z). Expanding the polynomials would lead to a single polynomial and hence it is a special case of the above. Let n := n 1 + · · · + n q be the total number of monomials and let us define the (N + q) × n matrix associated with b(z). Here 0 = (0, . . . , 0) and 1 = (1, . . . , 1) are row vectors of length |A i |. We define the codimension of A as The matrix A is the main object of study of the GKZ approach to hypergeometric functions. This definition allows us to consider integrals and PDEs later. This matrix is interpreted as representing a toric variety 5 . 4 Here we are adopting a rather unusual route by defining first Laurent polynomials and then the matrices Ai. The reason behind this is that in Feynman integrals we first consider polynomials with determinate coefficients and then consider the indeterminate case. In the mathematics literature, typically one first considers a configuration matrix and associates a polynomial to it. 5 See Ref. [48] for an introduction to toric ideals.
For later purpose we introduce the Newton polytopes of b i (z) following Ref. [38]. For each b i (z), the Newton polytope is the convex hull ∆ b i = conv(α i1 , . . . , α in i ) in R N . Like any other polytope, we can represent ∆ b i as the intersection of a finite number of halfspaces 6 : where µ i j ∈ Z N are primitive integer vectors in the inward normal direction of the facets of ∆ b i , and the ν i j ∈ Z are integers (See Fig.1). Here the product X · Y stands for the standard scalar product in R N . The polytope of b(z) is the Minkowski sum 7 (1, 1) Written as intersection of half planes, we have Let us now review the concept of a coamoeba of a variety. Consider the ideal I generated by some polynomials f 1 , . . . , f j ⊂ C[z 1 , . . . , z N ]. The zero set of the ideal I = f 1 , . . . , f j defines the algebraic variety Let us consider the case of a single polynomial f . The amoeba A f of the algebraic variety V(f ) is the image of V(f ) under the log mapping where Log(z) = (log |z 1 |, . . . , log |z N |). Amoebas were introduced by GKZ in Ref. [36]. Similarly, the coamoeba of V(f ) is the image of V(f ) under the argument mapping 6 Any polytope has a vertex representation and a half space representation. Going from one representation to another involves conversion algorithms that are beyond the scope of the present work. The interested reader can consult e.g., Ref. [49]. 7 For two polytopes P and Q in R N , their Minkowski sum P + Q is the set of all vectors p + q, such that p ∈ P , q ∈ Q [49].
where Arg(z) = (arg(z 1 ), . . . , arg(z N )). Here arg(z) is defined as usual. For example, the real positive line is given by R + = Arg −1 (0). Unlike the log mapping, the argument mapping is multivalued and we can think of it as a multiple periodic subset of R N or equivalently it can be viewed as the N -dimensional algebraic torus T N = (R/2πZ) N [38]. Coamoebas were introduced by Passare in 2004 8 . Our motivation to introduce coamoebas is that they define noncompact integration cycles for Euler-Mellin integral representations of A-hypergeometric functions. These cycles are particularly useful in generic cases where polynomials may vanish on the integration region. In many cases, coamoebas are difficult to study analytically and one has to consider a rough version, which is called the lopsided coamoeba [52][53][54][55]. The relation between A-hypergeometric functions and coamoebas has been studied in Refs. [37,38,55,56]. See Refs. [46,52] for an introduction to coamoebas.

GKZ systems and A-hypergeometric functions
We denote by H A (κ) the GKZ system associated with a matrix A and a vector of parameters κ. It is defined by the following data: 1. A (N + q) × n matrix A such that the vector (1, . . . , 1) lies in its row span. Hence, there is a vector ξ ∈ Z N +q , such that ξA = (1, · · · , 1). By definition this matrix is obtained from 2. A system of partial differential equations(PDEs) associated with A. Let u, v ∈ N n and consider where a ij denotes the components of A. Recall that θ j = c j ∂/∂c j and ∂ u = ∂ u 1 1 · · · ∂ un n .
A holomorphic function F (c) or formal series is called A-hypergeometric if it satisfies the above system of PDEs. GKZ systems can be rigorously defined in the language of holonomic ideals in the ring of differential operators with polynomial coefficients-the so called Weyl algebra D = K c 1 , . . . , c n , ∂ 1 , . . . , ∂ n modulo commutation rules. In this sense, the toric ideal associated with A is defined by where K[∂ 1 , . . . , ∂ n ] is a commutative polynomial ring. In addition, we construct the ideal generated by the column vectors κ T and θ = (θ 1 , . . . , θ n ) T . This ideal is given by where each generator of the ideal has the form n j=1 a ij θ j −κ i . The GKZ system H A (κ) denotes the left ideal on the Weyl algebra D generated by I A and Aθ − κ T . In this language a holomorphic function F (c) or formal series is called A-hypergeometric of degree κ if H A (κ) • F (c) = 0, where • denotes the action of the Weyl algebra on polynomials [40]. The language of holonomic ideals and D-modules is the appropriate one to treat the problem using computational algebra.
Let us denote by vol(A) the normalized volume-w.r.t the volume of the standard simplex which is equal to 1-of the convex hull of A 9 . Then, for generic parameters 10 κ, the rank of the system satisfies the inequality (Theorem 3.5.1 in [40]) which corresponds to the dimension of the solution space. This number can also be computed from the degree of the toric ideal I A .
For generic parameters rank(H A (κ)) = vol(A). In cases where the parameters have certain special values (nongeneric) the rank of the system can jump and we have rank(H A (κ)) > vol(A) (Example 4.2.7 in Ref. [40].).

Euler-type integral solutions
Solutions of GKZ systems have representations as Euler-type integrals. Accordingly, we call them A-hypergeometric functions. They can be constructed by taking the vector of parameters and then write the integral associated with A as follows: where the integration cycle is such that Ω ⊂ (C * ) N \V(b). It is usually assumed that the cycles are compact [35]. In Ref. [38], Berkesh, Forsgård, and Passare (BFP) constructed explicit noncompact cycles for these type integrals. In order to reach the above type of integral, BFP proceed in three steps. First, we consider that b(z) does not vanish in the positive orthant and consider the Euler-Mellin integral where in this case polynomial coefficients are fixed. BFP showed that if the polynomials b(z) are nonvanishing 11 , then Eq. (20) converges and defines and analytic funcion with parameters κ = (−β, −α) on the tube domain where int(τ ∆ b ) is the interior of of the weighted Minkowski sum of the Newton polytopes of b j weighted by τ , i.e., τ ∆ b = q j=1 τ j ∆ b j . The second step is to consider the less favorable case where the polynomials b(z) vanish on the positive orthant. Here, we can take a connected component Θ of R N \A b , where A b denotes the closure of the coamoeba of b and consider the integral where θ ∈ Θ is a representative of the connected component of complement of coamoeba of b(z). This is essentially a change of variables and a slight perturbation of θ does not impact the result of the integral. We write an analogue of Theorem 2.4 in Ref. [38], where the proof can be found.
The third step is the transition to A-hypergeometric functions by promoting the coefficients of b(z) in Eq. (22) to indeterminate, therefore we consider them as variables. The integral 12 is a representation of an A-hypergeometric function (Theorem 4.2 in [38]). For generic parameters κ, Eq.(24) provides a basis of solutions of H A (κ), where each integral is evaluated on a representative of Θ for each connected component of R N \A b (See Fig.2 13 ).
generate the series solutions from it. The role of the indicial equation in the Frobenius method is played by an indicial ideal ind w (I A ) with respect to w, together with the ideal Aθ − κ T . As in the Frobenius method, the problem consist of finding the roots γ of those ideals and generate the coefficients of the series. Series thus obtained belong to the Nilsson ring, i.e, series of the form In this paper, we will be chiefly interested in the case where the resulting series are logarithm-free, i.e., where β = 0 in the above equation. Let us turn our attention to these cases. Let be a lattice of rank m and let κ T = Aγ T 14 . For u ∈ L we can write u = u + − u − , where u ± ∈ N n have disjoint support. Now, for γ ∈ C n we define the following quantities expressed as falling factorials where (a) x are Pochhammer symbols. Then for γ ∈ C n , such that no element in γ is a nonnegative integer, the series 14 Notice that we will not be interested in general vectors γ ∈ C n but only those which are the roots of an ideal of H A (κ) with respect to some weight vector w. See Algorithm below.
is a solution of H A (κ) (See proposition 3.4.1 and Theorem 3.4.2 in Ref. [40]). When the vectors γ contain negative integers, we define the negative support of γ as In those cases, we consider the lattice and perform the sum over N γ in Eq. (29). In this paper, we will assume that the vectors κ are generic and thus none of the roots γ are negative integers 15 . This implies that supp(γ) = ∅ and therefore the sum runs over u ∈ L.
The roots γ can be obtained by finding the roots of an ideal known as fake indicial ideal. Accordingly, the roots are called fake exponents and we give the algorithm to compute them below. Series thus obtained are called canonical series. The series obtained from this algorithm have a common domain of convergence U w ∈ C n , which is characterized by a weight vector w. For generic roots γ, canonical series provide a basis of holomorphic solutions of H A (κ) in U w (For a description of this domain see Theorem 2.5.16 in Ref. [40] and Theorem 3.19 in Ref. [45]).
The relation with Euler integrals goes as follows. Suppose we find r different roots. The general solution of the integral (19) is given by a linear combination of its canonical series [40], i.e., where the canonical series φ r can be computed independently of the cycle Ω. The integration cycle plays a role in the computation of the integration constants K i . We discuss how to obtain integrations constants for the particular case of Feynman integrals in the next section. We proceed with the canonical series algorithm to compute fake exponents γ. A comprehensive review of this algorithm can be found in Ref. [45]. We start with some definitions.

Definitions
Toric ideal. Let D be a Weyl algebra over K, which is a free (noncommutative) associative K-algebra generated by K c 1 , . . . , c n , ∂ 1 , . . . , ∂ n modulo the commutation rules Let A be as in Eq.(5), then generates the toric ideal associated with A. Toric ideals can be computed using reduced Gröbner bases [57]. Initial ideal. Let w ∈ R n be a weight vector and let I A be a toric ideal. For w generic, the ideal in w (I A ) is a monomial ideal generated by the leading terms of I A with respect to the partial ordering w .
Standard pairs. Let R = K[∂ 1 , . . . , ∂ n ] and let I be a monomial ideal in R. Furthermore, let ∂ α be a monomial and F ⊆ {1, . . . , n}, where α ∈ N n . A standard pair of a monomial ideal I is a pair (∂ α , F ) satisfying three conditions: Let us denoted by S(I) the set of all standard pairs of I. The decomposition of I into irreducible monomial ideals can be obtained from the identity.

Algorithm (Saito-Sturmfels-Takayama)
In this algorithm we will set K = C.
Input: Matrix A, weight vector w, and complex parameters κ. Output: Roots of the fake indicial ideal fin w (H A (κ)).

Compute the toric ideal associated with
Notice that this is an ideal in the commutative polynomial ring C[∂ 1 , . . . , ∂ n ]. This ideal is the input to compute the combinatorial object of standard pairs S and the initial ideal with respect to w.
2. Let w ∈ R n be a generic weight vector. Compute the initial ideal in w (I A ) with respect to w and obtain its standard pairs S(in w (I A )). Standard pairs are combinatorial objects which tells us the types of solutions.
3. Use the standard pairs to construct the indicial ideal where 5. The fake indicial ideal with respect to w is given by 6. Compute the roots of fin w (H A (κ)). These are called fake exponents and we denote them by γ.
The canonical series are then given by Eq. (29). In order to write solutions, we compute the kernel of A to obtain the generating lattice Finally, we set [γ] u − = 0 whenever w.u < 0 (see Lemma 3.17 in Ref. [45]). Clearly, this property allows us to choose certain weights w that simplify the sums. We can use Macaulay2 to perform the above operations. An example is provided in Appendix B.

Examples
In order to illustrate the methods, we will take two representations of the Gauss hypergeometric function.

Double integral Gauss hypergeometric function
The first part of this example follows [37], where the reader can find details. Suppose we are interested in the following integral From the point of view we are adopting, we should think on the more general situation where the polynomial in the denominator has indeterminate coefficients, i.e., Its associated configuration matrix reads Then we consider the integral where Ω is a suitable cycle of integration. Following [37] the Newton polytope of b(c, z) can be represented by the inequalities(See Eq. (7)) and hence, we can read off the vectors µ i and numbers ν i Integration cycles can be obtained taking θ = (arg(c 1 /c 2 ), arg(c 1 /c 3 )). The BFP Theorem states that where Θ ∈ T 2 \A b . Taking the representative θ ∈ Θ, the entire function reads Taking the point c = (1, 1, 1, c), we have the coamoeba of b(z) shown in Fig.3, where we see that It is instructive to recover this result from the canonical series algorithm. Let w = (0, 1, 1, 1) 16 and κ = (−β, −α 1 , −α 2 ). The toric ideal associated with A reads We have in w (I A ) = ∂ 2 ∂ 3 , which we use to obtain its standard pairs 16 Choices of w are tied with the property that for w · u < 0 we have Thus, we obtain Computing its roots leads to In addition, computing ker( In order to start the sum from n = 0, we set Let us work out the case of γ 1 . Eq. (29) gives hence where we have used the series representation of the Gauss hypergeometric function (191). Similarly, for γ 2 we have Therefore, the solution as linear combination of canonical series reads We proceed to compute the integration constants. We take the noncompact cycle Ω = R 2 + . From the roots, we observe that if c 3 = 0, then φ 2 vanishes and if c 2 = 0, φ 1 vanishes. Hence we make which leads to Similarly, we obtain Let us now recover the original integrals by setting c 1 = c 2 = c 3 = 1, c 4 = c in Eq. (58). We have which equals Eq.(46) after using the identity (197) with z = 1 − c.

Single integral Gauss hypergeometric function
Let us now study the univariate integral of the Gauss hypergeometric function. We have In order to evaluate the integral we consider the (toric) polynomial and its associated GKZ system. The integral under consideration reads where κ = (−β 1 , −β 2 , −α). We have L = Z(1, −1, −1, 1). Taking the weight vector w = (0, 1, 1, 1), we obtain Hence, we have the two series solutions Specializing to our case, we have c 1 = c 2 = c 3 = 1 and c 4 = c, therefore the solution of the integral is given by where the constants can be obtained by setting c 2 and c 3 to zero in Eq. (65) choosing Ω = R + . The final result reads which evaluates to (63) after using Eq.(197) with z = 1 − c.

Feynman integrals as A-hypergeometric functions
In this Section we will interpret Feynman integrals as particular points of A-hypergeometric functions. First, we will introduce the parametric representation based on g = U + F used in Ref. [44] by Lee and Pomeransky. Later, we will give our proposal for defining GKZ systems based on g. Examples will be presented at the end of this Section.

Lee-Pomeransky representation of Feynman integrals
Let us consider a typical Feynman integral in Euclidean space in dimensional regularization. A L-loop integral with N propagators and E independent external momenta may be written as follows where the inverse propagators are of the form The matrices, M i , Q i , and J i have dimensions L×L, L×E, and 1×1, respectively. Integration over loop momenta through Feynman parameters and a Mellin transform leads to the Lee-Pomeransky parametric representation [44] I F (α) = ξ Γα where we have used the multi-index notation in the second equality and R + = (0, ∞). The overall factor and the polynomial g(z) are defined as follows: where U and F are the Symanzik polynomials. In order to compute them, we construct the matrices We then have where F is appropriately scaled in order to make it dimensionless. The polynomial U is a homogeneous polynomial of degree L and the polynomial F is homogeneous of degree L + 1.

Feynman integrals and canonical series
A suitable definition for computational algebra purposes would be to consider Feynman integrals as a linear combinations of their canonical series, whenever we can compute them. As they stand, Feynman integrals will not satisfy the system of PDEs associated with a GKZ system (Eqs.(13)- (14)) because the polynomial g(z) has fixed coefficients. Accordingly, the first step in our construction will be to take g(z) and consider its associated toric polynomial g(c, z). We will demand that this polynomial furnish GKZ system in the sense of Eq. (29). In particular, this means that g(c, z) must lead to matrices of co(A) > 0 in order to compute a generating lattice and follow the canonical series algorithm. Since we have a single polynomial, the codimension of the matrix A associated with g(c, z) is given by co(A) = n − N − 1, where n is the number of monomial terms in g(c, z) with common exponent vectors. Therefore, we will have co(A) = 0 when the number of terms in g(c, z) equals the number of rows in A, or equivalently when n = N + 1. A typical example of this situation is the polynomial of the L-loop massless cantaloupe graph (Fig.6). For instance, let us consider the case L = 1, i.e., the massless bubble graph (Fig.4). The g(z) polynomial of this graph is given by g(z) = z 1 + z 2 + sz 1 z 2 , therefore The codimension of A bubble is zero and hence ker A = ∅. Therefore, we cannot use this matrix to represent a solution of a GKZ system as canonical series in the sense of Eq. (29). However, in this simple example we know that the result of the massless bubble integral can be recovered from the one-mass bubble integral by taking the limit of the mass going to zero. Furthermore, as we will see in the examples, the matrix associated with the one-mass bubble integral is given by which has co(A one-mass ) = 1. Since the one-mass bubble evaluates to a Gauss hypergeometric function, we may consider the massless bubble as a special case of a 2 F 1 (a, b, c; p 2 /m 2 ) function, for some a, b, c depending on the powers of the propagators and the dimension. An equivalent alternative provided by GKZ systems is to deform the polynomial by adding an arbitrary constant r(z) = c 1 to g bubble (c, z) and consider instead where the constant c 1 can be set to zero at the end of the computation. This matrix corresponds to the GKZ system of a Gauss hypergeometric function. We will work out explicitly this example in Section 3.3.1. There, we will take the limit c 1 → 0 of a Gauss hypergeometric function as we would do for the mass in the case of a one-mass bubble. The choice of a deformation can be made systematically as we will see for the L-loop cantaloupe graph. Therefore, in cases where g(c, z) leads to a configuration matrix of co(A) = 0, we will consider a deformation of g(c, z) demanding that it defines a solution of a GKZ system in the sense of Eq.(29), thus allowing us to construct canonical series solutions. Recall that deg (U) = L and deg (F) = L + 1. We may choose any polynomial r(z) of deg(r) < L and define with the requirement that the associated matrix has co(A) > 0. Here r(c, z) denotes the toric polynomial associated with r(z). Similarly, U(c) and F(c) denote the Symanzik polynomials where the coefficients now are considered as variables. We have checked for the integral of a massless L-loop cantaloupe graph, up to 5-loop, that this definition leads to well behaved GKZ systems and that we can recover the original integral using canonical series. The second step in our construction will be to consider Euler-type solutions and series solutions of the GKZ systems furnished with g r (c, z). Let us first define the class of Euler integrals associated with the Feynman integral (74). We will drop the overall factor ξ Γα , which is a nonzero constant independent of the kinematics and instead consider which is a toric generalization of Eq.(74), i.e., where the deformation has been introduced and the coefficients of g r (z) are considered as variables. Notice that at this point, the noncompact cycle Ω will not be in general R N + and it will be determined by the coamoeba of g r (c, z) (Sec.2.3.1). Let us first show that the class of integrals obtained from g r (c, z) is A-hypergeometric 17 . We have the following theorem, which is a specialization of theorem 2.7 in Ref. [35] for the case of a single polynomial g r (c, z). That this theorem still holds for noncompact cycles is discussed in Ref. [38,39].

Theorem.
Let g r (c, z) be the deformed polynomial in N variables obtained from g(c, z) = U(c) + F(c), where F(c) and U(c) are obtained by considering the coefficients appearing in the Symanzik polynomials as variables. g r (c, z) is obtained by introducing a deformation r(c, z) demanding that its matrix satisfies co(A) > 0. Let A = (a 1 a 2 · · · a n ) be the configuration matrix associated with g r (c, z) and consider the polynomial with indeterminate generic coefficients Let A be its associated (N + 1) × n matrix A = 1 1 . . . 1 a 1 a 2 . . . a n .
The Euler-Mellin integral is a solution of the A-hypergeometric system H A (κ) of degree κ = (−d/2, −α). Noncompact cycles Ω can be obtained by taking the coamoeba of g r (c, z) and choosing representatives θ of connected components Θ ∈ R N \A gr (See Sec.(2.3.1).) Proof. The proof goes along the lines of Ref. [39] specializing to the nonhomogeneous case and the case of a single polynomial. Let us consider first Eq.(13). We have 17 A related construction-which does not introduce deformations-has been considered in Ref. [41] for the special case where the powers of the propagators are given by α = (1, . . . , 1). Here we consider generic powers in the propagators as dictated by canonical series solutions.
where (ρ) (x) denotes the falling factorial. Similarly and from Au = Av, the result (∂ u − ∂ v )I gr (κ) = 0 follows. Let us now focus on the second set of differential equations (14). Consider the first row of A in Eq.(85). We have where we have used integration by parts in the last equality. This completes the proof. In this way, we have generated a class of integrals related with the Lee-Pomeransky representation of Feynman integrals. Feynman integrals will correspond to special cases (points) of A-hypergeometric functions whenever the cycle Ω can be taken as R N + and the coefficients c can be taken as functions of the kinematic invariants. The behavior of the integral as the coefficients c vary can be studied through Eq. (22). This representation is also useful as it provides noncompact cycles at the cost of imposing conditions on κ. However, this can be sorted out by analytic continuation and it can be shown that, as a function of the variables c, Eq. (22) is A-hypergeometric everywhere. Here c is taken in C n \Σ A , where Σ A denotes the singular locus of all A-hypergeometric functions(Theorem 4.2 in [38]). The main difficulty in this approach is to choose a representative of some connected component Θ in R N \A gr . In other words, we have to select a point such that the cycle is nonvanishing on the set Arg −1 (θ), where θ ∈ Θ. The integration region of the Feynman integral, namely R N + , suggest taking any θ = (arg(f 1 (c)), . . . , arg(f N (c))), f i (c) > 0 provided θ / ∈ A gr (See e.g. Fig.3). This gives one connected component of vol(A) many ones and a possible integration cycle 18 .
On the other hand, logarithm-free canonical series allows us to study the behavior of the solution space under variations of the parameters κ at nonsingular points [57,61]. More important, once we have identified the above integrals as solutions of a GKZ system, we can use the statement in Eq.(32) relating Euler-type integrals and canonical series. For our case, this statement reads where φ 1 , . . . , φ M are the canonical series associated with A. Fake exponents γ for each φ 1 , . . . , φ M can be obtained from the SST algorithm in Sec.2.4.1. This equality is fundamental for our purposes as it allows both to take the limit of the deformation to zero and computing the integrations constants. Let us discuss these limits. The form of the series solution for some fake exponent γ reads These are characterized by a weight vector w ∈ R n , which selects a common domain of convergence U w (See Theorem 2.5.16 in [40]). In addition, we have the restriction [γ] u − = 0 for u.w < 0. Taking the limit of the deformation to zero amounts to take some of the coefficients in c = (c 1 , . . . , c n ) to zero in Eq.(88) and similarly for the remaining canonical series. On the LHS of Eq.(87), taking this limit amounts to recover the undeformed integral we started with. Let us clarify this point. The integral in the LHS of (87) is a well defined A-hypergeometric function provided c are generic and the cycle is taken in some θ / ∈ A gr . The RHS of (87) is an asymptotic expansion of such integral. Therefore, the limits on the LHS correspond to special values of the A-hypergeometric functions in the RHS. A judicious choice of the deformation ensures that this limit can be taken systematically as we will see in the examples.
In general, this limit will not be smooth as it will require the evaluation of A-hypergeometric functions on their singular points 19 and therefore we may require analytic continuation of the corresponding A-hypergeometric function. This can be seen as follows. A generic canonical series solutions can be understood as a function in co(A) variables. Setting one of these variables to zero may require a transformation of a variable to, say, its inverse and hence we require analytic continuation. A judicious choice of the weight vector w can simplify taking this limit, since it provides the condition [γ] u − = 0 for u.w < 0, thus determining the arrangement of the powers of c u in Eq.(88). As we will see for the L-loop cantaloupe graph, we can systematically choose a deformation and a weight vector to take this limit. However, we will leave as a conjecture that this this limit can always be taken.
For integration constants, we will use the information available on the fake exponents following Ref. [40]. The initial series in Eq.(88) are given by c γ . Therefore, the positions of the 0's on each γ will tell us which elements in c we have to set to zero in order to compute the integration constants.
Let us give some final comments. Notice that in order to recover the original Feynman integralin agreement with our definition of the polynomials associated with the GKZ system and of the integral (86)-we must set coefficients c to their kinematic values at the end of the computation. The integral (86) can be understood as a holomorphic function on U w . This observation gives us a nice way of deducing the long known fact that convergent Feynman integrals are functions of a Nilsson class [2], which is clear from their canonical series.
We end with our prescription to compute integration constants.

Integration constants
Suppose that we obtain M fake exponents from the SST algorithm. Let us isolate an exponent, say, γ = (γ 1 , . . . , γ n ) and suppose it contains k < n zeros in positions σ 1 , . . . , σ k . Take the coefficients associated with those positions to zero and consider where we take Ω = R N + as the cycle in the LHS. This procedure will compute a single coefficient. We repeat the process until we have computed all of them.

Examples of co(A) = 0
The purpose of co(A) = 0 examples is to show how to deal with the appearance of a deformation.

Structure of the examples
In the following examples we will omit the overall gamma factors ξ Γα and hence consider I(α) := I F (α)/ξ Γα along with its toric version I gr (κ). The vector κ has the form where N is the number of propagators. We assume that the powers of the propagators have generic noninteger complex values. This simplifies the discussion since the sum runs over L = ker Z A. The kernel of A leads to a rank co(A) lattice. For a single generator, we write L := Z(a 1 , . . . , a n ), u = n(a 1 , . . . , a n ), a i ∈ Z.
We set In order to indicate the i-th component of a root vector γ r we write γ i r .

Massless bubble
The simplest example is the massless bubble ( Fig.4) with inverse propagators where s = −p 2 . We have This polynomial leads to a matrix of codimension co(A) = 0, hence we introduce a deformation r(z) = c 1 and consider instead g r (c, z) = c 1 + c 2 z 1 + c 3 z 2 + c 4 z 1 z 2 ⇐⇒ A = Thus, we have for θ = (arg(c 1 /c 3 ), arg(c 1 /c 2 )). We have studied this integral in Sec.2.5, where we established that In order to recover our Feynman integral we have to take the limit c 1 → 0. The limit has to be taken carefully as c 1 appears both as a factor and in the argument of the hypergeometric function.
Using the identity (195) we have Therefore (96) Finally, taking the limit c 1 = 0, setting c 2 = c 3 = 1, and c 4 = s we recover the desired result where we have used identity (194). Notice that in order to take this limit, we have evaluated Eq.(94) on one of the singular points of the Gauss hypergeometric function, thus requiring an Euler transformation.

The single-scale massless triangle graph
Let us now turn our attention to those limits through canonical series. Let us consider the triangle graph in Fig.5. The inverse propagators read where −p 2 1 = −p 2 2 = 0. Computing the relevant polynomial and taking the deformation r(z) = c 1 leads to where, at the end of the computation, we will make the identifications c 2 = c 3 = c 4 = 1 and Computing ker(A), we have Choosing w = (1, 0, 0, 0, 0), we have u.w = n, hence [γ i ] u − = 0 for n < 0. Setting A = α 1 +α 2 +α 3 , the fake indicial ideal and its roots read Inserting the roots in Eq. (29) gives [γ i ] (0,n,n,0,0) [γ i + (n, −n, −n, 0, n)] (n,0,0,0,n) which leads to the series Integration constants can be easily computed by setting c 1 = 0 and c 5 = 0 in Eq.(100) with Ω = R 3 + . We write them collectively as Taking the limit c 1 → 0 and setting c 2 = c 3 = c 4 = 1, c 5 = s, we obtain Let us remark that in this example we have set c 1 to zero in order to compute one of the integration constants-as can be seen from the roots-therefore reaching a tautology. We can fix this by choosing an appropriate weight vector such that none of the roots contain zero in position 1. This leads to a more complicated version of the canonical series which will require analytic continuation in order to take the limit c 1 → 0. Since we want to interpret Feynman integrals with codimension zero matrix as a certain limiting cases of A-hypergeometric functions, we do not worry about this situation. In fact, as we will see in the next example, we can take advantage of this by choosing a weight vector such that all massless cantaloupe graphs are defined by the coefficients of a linear combination of two Gauss hypergeometric functions.

The massless L-loop cantaloupe graph
Let us consider the case of the L-loop cantaloupe graph (a.k.a. banana graph) shown in Fig.6.
We parametrize the inverse propagators as follows . . .  The polynomial g(z) of this graph can be written in general as where s = −p 2 . The integral to be computed reads It is easy to check that the g(z) polynomials of this graph lead to codimension zero matrices and hence they fall under the class of problems where we must introduce a deformation to define a GKZ system in the sense of Eq. (29). In order to perform such deformation systematically, let us introduce some notation. Let 1 i denote a sequence of 1's of length i and similarly for 0 j . We have the relation i + j = L + 1. Furthermore, let v : At each loop, we set a deformation monomial hence we have where c L+3 = s. Let us give an example. For L = 3, v = (1, 1, 0, 0) and r(z) = c 1 z 1 z 2 , then we have the deformed toric polynomial g r (c, z) = c 1 z 1 z 2 + c 2 z 1 z 2 z 3 + c 3 z 1 z 2 z 4 + c 4 z 1 z 3 z 4 + c 5 z 2 z 3 z 4 + c 6 z 1 z 2 z 3 z 4 .
After introducing the deformation, the (L + 2) × (L + 3) matrix associated with the L-loop cantaloupe graph can be written in the general form Notice that each row contains L + 3 elements. There are L rows with have a single zero and 2 rows with two zeros. The integral under consideration is given by where κ = (−β, −α 1 , . . . , −α L+1 ). Computing the kernel of the above matrix leads to where by definition 0 0 := ∅. We choose w = (1, 0 L+2 ), thus obtaining The roots can be written as which lead to the canonical series where The relevant integration constant reads which after a change of variables corresponds to the formula of the Mellin transform of a linear function. Setting c 1 = 0, c 2 = · · · = c L+2 = 1, and c L+3 = s we arrive at The formulas for the fake indicial ideal and its roots have been checked up to 5-loop. In this example, we have seen that deforming g(z) leads to a GKZ system of co(A) = 1 and taking the limit of the deformation to zero at the end of the computation allows us to interpret Feynman integrals as the limit of a linear combination of their canonical series.

One-mass bubble
Let us work now with the bubble integral with one massive internal line (Fig.7). The inverse propagators of this integral read Omitting the overall Γ factors we have the Lee-Pomeransky representation Now, let us consider the related problem from the GKZ point of view. We are interested in the more general situation Taking w = (0, 1, 1, 1), we obtain the following fake indicial ideal Then, the roots of this ideal read In addition, we have L = Z(−1, 1, 1, −1), and hence u = n(−1, 1, 1, −1). We also have u.w = n, which implies [γ i ] u − = 0 for n < 0. The canonical series simplify to the functions The result of the integral is a linear combination of φ 1,2 . Constants of integration can be obtained by integrating (129) taking c 3 = 0 and Ω = R 2 + , and similarly for c 2 . This leads to the multiplying factors Setting c 3 = (s + m 2 ) and c 4 = m 2 , the resulting integral reads The form of the result as a sum of two hypergeometric functions is reminiscent of the negative dimension approach [62]. Using the Eq.(196), we obtain
We can now write the canonical series (29) as The constants of integration read In this case the canonical series evaluate to hypergeometric functions 3 F 2 (a, b; c; d, e; x) (see definition in Eq.(192)). Let A = α 1 + α 2 + α 3 + α 4 and x = (c 2 c 4 c 5 )/(c 1 c 3 c 6 ). We can write the result in the condensed form 1 2 3 Figure 11: Three-scale triangle: Setting the constant to the values of the original polynomial, i.e., c 1 = · · · = c 4 = 1, c 5 = s and c 6 = t, we obtain Hence, we have recovered the result obtained by the differential reduction method using Gröbner bases given in Ref. [63] based on [64,65]. Interestingly, the computation of integrations constants in Ref. [63] follows a similar prescription as ours. However, in our case, because of our choice of w, we only set c 5 = 0 to compute K 1 but we do not use c 6 = 0 as in [63]. This results also matches nicely the results obtained in Ref. [66] using the negative dimension approach.

A co(A) = 2 example
Finally, let us consider the triangle integral shown in Fig.11 where all external momenta are nonvanishing. The inverse propagators read where we have −p 2 1 = s 1 , −p 2 2 = s 2 , and −(p 1 + p 2 ) 2 = s 3 . Using momentum conservation, we have The associated matrix reads This matrix has co(A) = 2 and hence the resulting A-hypergeometric functions will be functions of two variables. The integral under consideration then reads Computing ker(A), we have which means that in this case the lattice is generated by u = m(−1, 0, 1, 1, 0, −1) + n(−1, 1, 0, 0, 1, −1), m, n ∈ Z.
We choose w = (0, 0, 1, 0, 0, 0), which has the advantage of restricting the sum over m. We have u.w = m, hence [γ i ] u − = 0 for m < 0. The fake indicial ideal reads with roots Inserting the roots in Eq. (29) gives four solutions where we have defined x = (c 3 c 4 )/(c 1 c 6 ) and y = (c 2 c 5 )/(c 1 c 6 ). Due to our choice of weight vector, the sum has been only restricted in m. However, it is easy to see that the sum over negative integers n vanish due to the presence of (1) n . Hence the above sums have the form which are sum representations of the Appell hypergeometric function F 4 (See Eq.(193)). In order to compute the integration constants we can follow our prescription and set c 4 = c 5 = 0 in Eq.(177) and take Ω = R 3 . This computes the first coefficient in the expansion K 1 . We follow the same procedure for the remaining integration constants. They can can collectively be written as Setting c 1 = c 2 = c 3 = 1 and c 4 = s 3 , c 5 = s 1 , and c 6 = s 2 , we can write the Feynman integral as which agrees with the results obtained via the Mellin-Barnes integral representations [67,68] and the negative dimension approach [62].

Conclusions and outlook
In this paper have studied the relation between the Lee-Pomeransky representation of Feynman integrals and GKZ systems. We have shown that in generic cases we can associate a matrix A of co(A) > 0 to a deformed polynomial g r (c, z) = r(c, z) + U(c) + F(c), where r(c, z) is introduced to ensure a canonical series representation. U(c) and F(c) are toric polynomials associated with the Symanzik polynomials. Under these restrictions, we can interpret a large class of Feynman integrals as furnishing a solution of a GKZ system based on A. The canonical series algorithm then allows us to evaluate integrals with arbitrary powers in the propagators as linear combinations of A-hypergeometric functions. Feynman integrals are recovered at the end of the computation by identifying the coefficients of the toric polynomials with their kinematic values. Using the canonical series method, we have evaluated several integrals for arbitrary noninteger powers in the propagators. A particularly nontrivial example is the on-shell massless box with arbitrary powers in the propagators. With this method the result was obtained as a particular case of an A-hypergeometric integral and it matches the result based on the recurrence relations method based on Gröbner bases [63] and the negative dimension approach [66]. Another nontrivial example is the co(A) = 2 example of the three-scale massless triangle, which is in agreement with the negative dimension approach as well. It would be interesting to study the relation between those methods and the canonical series algorithm.
Computing canonical series is a straightforward computational algebra problem. However, these series heavily depend on the choice of a weight vector w, which sets the initial ideal and effectively chooses a domain of convergence. This choice is tied with the available information about the integration cycle and ultimately with our ability to compute integration constants. Our recipe of setting Ω = R N + is the obvious choice and was motivated by the cycles, which one obtains by studying the coamoeba of g r (z) in Euclidean kinematics. A full characterization of the coamoeba of g r (z) might be necessary when non-Euclidean kinematics is considered and going to higher codimensions.
The rank of the system arising from the toric version of g r (z) is bounded by vol(A) which, in general, is greater than the number of master integrals arising from IBP identities or from the Euler characteristic [69]. The canonical series algorithm produces vol(A) hypergeometric series in co(A) variables. They collapse to simpler expressions once we set the coefficients c to their kinematic values. This typically involves setting one or more of these variables to unity, which amounts to evaluate hypergeometric functions at singular points. At co(A) = 1 the functions appearing at those limits are Γ-functions. At co(A) = 2, those limits lead to hypergeometric functions of one variable, which then collapse to simpler expressions. It would be interesting to study the mechanism which relates the number of master integrals and the number of canonical series solutions.
We believe that the application of GKZ systems and canonical series to Feynman integrals is not limited to the Lee-Pomeransky representation. Indeed, it would be interesting to apply these ideas in representations where an algebraic definition of the integration cycles is available. For instance, this is the case of the representation due to Baikov [70], which has recently been studied specially in the context of maximal cuts [71][72][73][74][75]. This approach is closely related to Ref. [76], where bases of Pfaffian systems for GKZ systems are constructed using twisted cohomology groups.
We leave these explorations for future work.  (198)