D -Module Techniques for Solving Diﬀerential Equations in the Context of Feynman Integrals

Feynman integrals are solutions to linear partial diﬀerential equations with polynomial coeﬃcients. Using a triangle integral with general exponents as a case in point, we compare D -module methods to dedicated methods developed for solving diﬀerential equations appearing in the context of Feynman integrals


Introduction
Differential equations are ubiquitous in physics.Their solutions are often expressed in terms of special functions.Well-known examples are the hypergeometric function, or the Hermite polynomials that appear in the solution of the quantum mechanical hydrogen atom.In perturbative quantum field theory, a crucial bottleneck is the evaluation of Feynman integrals.In particular, obtaining analytic or numerical results for higher-loop Feynman integrals is paramount for deriving more precise theory predictions for processes at the Large Hadron Collider.Because of this, substantial research efforts are dedicated to evaluating Feynman integrals.These integrals can be seen as solutions to linear partial differential equations (PDEs) with polynomial coefficients.Analyzing these PDEs and extracting useful information from them is an important problem in this field.In this paper, we analyze this problem from two perspectives.The first one consists of D-module techniques from algebraic analysis in mathematics, and the second one consists of dedicated tools developed for solving the specific differential equations satisfied by Feynman integrals.
In the first approach, one uses the fact that systems of linear PDEs generate ideals in the Weyl algebra, which is denoted by D. Powerful D-module techniques [55,56] allow us to extract precious information from the PDEs.For instance, the singular locus of a D-ideal I describes where the solutions to the system of PDEs encoded by I may have singularities.The notion of a D-ideal being regular singular is inherently linked to its solutions having moderate growth-which is expected for Feynman integrals.The holonomic rank of a holonomic D-ideal tells us about the number of linearly independently solutions.It corresponds to the number of master integrals in the physics terminology.Finally, Gröbner basis methods allow us to compute canonical series solutions.These are solutions of a very special form: they are polynomials in the variables x i 's and log(x i )'s, with coefficients in formal power series in the x i .In the regular holonomic case, the solution space is fully spanned by solutions of that particular form.We compute these solutions with the help of algebro-geometric properties and the Gröbner data of the D-ideal via an algorithm of Saito, Sturmfels, and Takayama (SST).Given a regular holonomic D-ideal I and a real weight vector w, this algorithm allows us to compute all terms up to a specified w-weight k in the canonical series solutions to I.These truncated series take the form where we use the multi-index notation x a = x a 1 1 • • • x an n , and similarly for the logarithms.The index set C * Z over which the sum in (1.1) runs can be read completely from the D-ideal; this will be made precise in the article.
The second approach builds on various techniques that have been developed for evaluating Feynman integrals.The latter satisfy Picard-Fuchs equations [45,52], which are higher-order linear PDEs for individual Feynman integrals, or, equivalently, systems of firstorder equations for master integrals, see [6,38] for reviews.The Fuchsian nature of the singularities-which is expected for Feynman integrals-allows one to derive asymptotic expansions around singular points using a method of Wasow [65].
There have been various interactions between the relevant mathematics and physics communities already, resulting in several interesting works, for example on determining the number of master integrals [1,9,47], on deriving Picard-Fuchs equations for Feynman integrals [45,52], and on relating Feynman integrals and GKZ systems [3,21,27,42,43,64].However, to our knowledge, a systematic synthesis and comparison of D-module techniques and methods employed in high energy physics has not yet been done.In this work, we initiate such a comparison, with the aim of providing insights for both communities.For this initial study, we focus on a particular system of physically motivated PDEs and show how one obtains canonical series solutions of the form (1.1), both using the SST algorithm and Wasow's method.In the latter, we capture the weight vector from the Gröbner techniques via an auxiliary variable.
Outline.Section 2 introduces the systems of PDEs we are investigating here, together with their origin in physics.Section 3 provides background on ideals in the Weyl algebra and their (multi-valued) holomorphic solutions that will be needed throughout this article.It also outlines how to compute canonical series solutions of regular holonomic D-ideals via the SST algorithm, which is based on Gröbner basis computations in the Weyl algebra.In Section 4, we compute canonical series solutions for the three-particle case using the SST algorithm.In Section 5, we recover these solutions via a method of Wasow for computing solutions of systems that are in Fuchsian form.We conclude with an outlook to future work in Section 6. Appendices A and B contain some background on conformal symmetry as well as on integrable connections.In Appendix C, we provide an application of this method to a four-loop Feynman integral.

Running example: a triangle integral and conformal differential equations
As a test case for this work, we study a particular class of Feynman integrals: the one-loop "triangle" Feynman integrals, which correspond to the Feynman graph shown in Figure 1.These integrals are relevant to describe the scattering of three particles with momenta p 1 , p 2 , p 3 ∈ R d (with p 3 = −p 1 −p 2 due to momentum conservation) in a d-dimensional Minkowski spacetime.We denote by |v| 2 the Minkowski norm of the R d -vector v, i.e., |v| 2 := v ⊤ • g • v, where g = diag(1, −1, . . ., −1) is the metric tensor of Minkowski spacetime.The triangle Feynman integrals depend on the momenta of the scattering particles only through three independent variables x = {x 1 , x 2 , x 3 }, This becomes apparent in the Feynman parameterization, Figure 1: The Feynman graph representing the one-loop triangle Feynman integrals defined in Equation (2.1).Due to momentum conservation, p 1 + p 2 + p 3 = 0. Next to the internal edges, we record the corresponding exponent ν i as well as the loop momentum.
We approach the computation of these integrals by exploiting symmetries, which are encoded via differential equations.The triangle integrals are in fact solutions to a system of linear second-order PDEs, originating from conformal symmetry, where Here, c 0 = d is the number of spacetime dimensions, and c 1 , c 2 , c 3 are the conformal weights.We give some background on conformal symmetry and derive these PDEs in Appendix A. For the triangle integrals J triangle d;ν 1 ,ν 2 ,ν 3 are solutions to the conformal PDEs (2.5).Let us review a few useful properties of the operators P i in (2.6).First of all, the operator P 3 implies that the solutions are homogeneous of degree (2c 0 − c 1 − c 2 − c 3 )/2.Secondly, the system (2.5) is symmetric under permutations of the variables x i (together with the corresponding conformal weight c i ).While the operator P 3 is manifestly symmetric, the permutations map P 1 and P 2 into Q-linear combinations of themselves.As a result, the solution space is symmetric as well.For computing the solutions, we will focus on the physically motivated case c 0 = 4, c 1 = c 2 = c 3 = 2, which corresponds to the (classically) conformal φ 4 -theory in four spacetime dimensions. 1 In this case, the operators from (2.6) become Later on, we will consider the D 3 -ideal generated by the operators in (2.8) and will denote it by I 3 .Through Equation (2.7), we see that the relevant triangle integral has d = 4 and unit exponents, ν 1 = ν 2 = ν 3 = 1.A convenient analytic expression for the latter is from [67], where the function Li 2 denotes the dilogarithm (cf.[36]), λ is the polynomial and The function in Equation (2.9) is a solution of I 3 .It is a multivalued function and its analytic continuations are solutions of I 3 as well.The C-vector space of the analytic continuations of J triangle 4;1,1,1 is spanned by the four linearly independent functions (2.12) The functions f 2 , f 3 , f 4 are called discontinuities of f 1 in physics, i.e., they are differences of f 1 and an analytic continuation of f 1 .This can be shown conveniently using the symbol method [33] (for a review, see e.g.[30] and the references therein).We will see in Section 4 that the holonomic rank of I 3 is 4 and that this implies that the four functions in (2.12) span the whole space of holomorphic solutions to the system of PDEs encoded by I 3 .All these functions have moderate growth when approaching points of the singular locus and we will argue later that the D-ideal I 3 is indeed regular singular.As anticipated, the solution space is symmetric under permutations of the variables x: the functions f 1 and f 4 are invariant, while f 2 and f 3 transform into Q-linear combinations of themselves.In Section 5.2, we will recover these solutions by solving the Pfaffian system associated with I 3 .

D-ideals and canonical series solutions
We here recall basics about ideals in the Weyl algebra and their (multi-valued) holomorphic solutions that will be needed for the computation of canonical series solutions later on.This theory will also allow us to derive crucial information about our system of PDEs, such as the singular locus and number of holomorphic solutions, before even computing the solutions.
For further details about the theory of D-modules and holonomic functions, we refer our readers to [40,55,56] and the references therein.

D-ideals and their solutions
The (n-th) Weyl algebra, denoted D n or just D, is the free C-algebra generated by x 1 , . . ., x n , ∂ 1 , . . ., ∂ n modulo the following relations: all generators are assumed to commute, except ∂ i and x i .Their commutator is Elements of D are linear differential operators with polynomial coefficients: The rational Weyl algebra R n = C(x 1 , . . ., x n ) ∂ 1 , . . ., ∂ n is the ring of linear differential operators with coefficients in the field of rational functions C(x 1 , . . ., x n ) = {p/q | p, q ∈ C[x 1 , . . ., x n ], q = 0}.We will denote the action of a differential operator P on a function f (x) by the symbol •.For instance, ∂ • f denotes ∂f /∂x.
By ord (u,v) (P ), we denote the largest (u, v)-weight among the monomials appearing in P .The characteristic ideal of a D-ideal I is the ideal in (0,1) (I) ⊂ C[x 1 , . . ., x n ][ξ 1 , . . ., x n ] generated by the parts of highest (0, 1)-weight of all differential operators in I, where (0, 1) denotes the vector (0, . . ., 0, 1, . . ., 1) ∈ R 2n , i.e., one puts weight 0 to all x i and weight 1 to all ∂ i .The characteristic ideal lives in the polynomial ring, where one replaces ∂ i by ξ i to stress that they are now commuting variables.This is due to the fact that the commutator of two operators P, Q has (0, 1)-weight strictly smaller than deg(P ) + deg(Q).For n = 1, the initial in (0,1) (P ) is precisely the principal symbol of the differential operator P .The characteristic variety of I is Here, C 2n is the affine 2n-space C n x × C n ξ with coordinates x 1 , . . ., x n , ξ 1 , . . ., ξ n .A D-ideal is holonomic if dim(Char(I)) = n.By Bernstein's inequality, all components of Char(I) have dimension at least n; hence, a D-ideal I is holonomic if the dimension of its characteristic variety is smallest possible.The Zariski closure of the projection of Char(I) to the x-coordinates is the singular locus of I and is denoted by Sing(I) ⊂ C n x .Algebraically, Sing(I) is computed as the elimination ideal where (I : J (∞) ) denotes the saturation of an ideal I by an ideal J, the definition of which we recall now.Let I, J be ideals in a polynomial ring C[x 1 , . . ., x n ].For k ∈ N, the ideal quotient is (I : Then, the saturated ideal (I : J ∞ ) with respect to J is the C[x 1 , . . ., x n ]-ideal If I is holonomic, it follows that rank(I) < ∞.The reverse implication is not true.The holonomic rank can be computed as the number of standard monomials for a Gröbner basis of I in the rational Weyl algebra, see Section 3.2 for more details.Theorem 3.3 (Cauchy-Kovalevskaya-Kashiwara).Let I be a holonomic D n -ideal.On a simply connected domain U ⊂ C n \ Sing(I), the C-vector space of holomorphic solutions to I on U has dimension rank(I).Remark 3.4.In Theorem 3.3, finite holonomic rank is sufficient for the statement to hold.One arrives at a holonomic D-ideal by taking the Weyl closure (see Definition 3.5). ⋄ is a holonomic D-ideal.Numerous functions in the sciences are holonomic, e.g., hypergeometric functions, many trigonometric functions, some probability distributions, and many special functions such as Airy's or Bessel's functions, and polylogarithms.Zeilberger [66] was the first to study holonomic functions in an algorithmic way.If all solutions to a D 1 -ideal I have moderate growth (cf.[62, p. 146]) when approaching the singular locus of I, including at ∞ (and at boundary components of a compactification, respectively), the D 1 -ideal is called regular (singular), and irregular singular otherwise.For n > 1, it is more involved to read if a D n -ideal is regular singular.
Given an annihilating D-ideal of some functions, there is a straightforward way to slightly enlarge it, which is at the same time of theoretical interest: the Weyl closure.Hence, techniques that are valid for holonomic D-ideals can be applied to non-holonomic D-ideals of finite holonomic rank by passing to their Weyl closure.

Pfaffian systems
If I is a holonomic D-ideal, the D-module D/I gives rise to a vector bundle of rank rank(I), say rank(I) = m, with an integrable connection induced by the action of D on D/I.We refer the interested readers to Appendix B, where we explain this geometric interpretation.
In particular, the origin of the integrability conditions of Pfaffian systems becomes visible there.Here, we stick to what is needed for the sake of our computations in this paper.Let S = {s 1 , s 2 , . . ., s m } be the set of standard monomials of I in R n for a chosen term ordering ≺, i.e., those monomials ∂ b , b ∈ N n , that are not contained in the initial ideal of I; cf.[55, p. 33] for more details.Without loss of generality, we can assume s 1 = 1.Let f ∈ Sol(I) be a solution to I and let F = (f, s 2 • f, . . ., s m • f ).Since rank(I) = m, there exist unique matrices P 1 , . . ., P n ∈ C(x 1 , . . ., x n ) m×m such that The matrices P i can be computed via a Gröbner basis reduction (cf.[56, p. 23]).The system in (3.12) is the Pfaffian system of I for the chosen term ordering on the Weyl algebra.The matrices P i obey the integrability condition, which translates as On the right hand side of (3.13), entry-wise differentiation the matrices P i is meant.

Indicial ideals and the Nilsson ring
We first recall some results from [55, Sections 2.2-2.6] in Theorems 3.6-3.16below, in order to present the overall strategy to compute solutions of our D-ideals.Here and throughout the rest of the file, weight vectors for the Weyl algebra are allowed to be taken from Among others, W contains the set {(−w, w)|w ∈ R n }.For a weight vector (u, v) ∈ R 2n and P ∈ D n , the initial form in (u,v) (P ) of P denotes the part of P of maximal (u, v)-weight.
For weights of the form (−w, w), one denotes the initial form simply by in w (P ).Each weight vector (u, v) induces an increasing filtration F • (u,v) (D n ) of the Weyl algebra by the (u, v)-weight.For weights of the form (−w, w), the associated graded ring is isomorphic to the Weyl algebra itself, cf.[55, p. 4].For (u, v) = (0, 1) ∈ R 2n , the associated graded ring gr (0,1) (D n ) is the polynomial ring C[x 1 , . . ., x n ][∂ 1 , . . ., ∂ n ], in which case one typically replaces ∂ i by ξ i to highlight the commutativity of the variables.The initial ideal in (u,v) (I) of a D n -ideal I is the ideal generated by in (u,v) (P ) of all P ∈ I.It naturally lives in the graded ring gr (u,v) (D n ).The initial ideal for the weight (0, 1) ∈ R 2n is exactly the characteristic ideal of the D n -ideal, which was introduced in Section 3.
In what follows, we focus on weights of the form (−w, w).The following theorem summarizes Theorems 2.2.1 and 2.5.1 of [55] about the holonomic rank of the initial ideal.Theorem 3.6.Let I be a holonomic D n -ideal and w any weight vector in R n .Then also the initial ideal in (−w,w) (I) is holonomic and rank(in (−w,w) (I)) ≤ rank(I).If I is regular holonomic, then rank(I) = rank(in (−w,w) (I)).(1) I is torus-fixed if and only if in (−w,w) (I) = I for all w ∈ R n .
(2) If w is generic for I, then in (−w,w) (I) is torus-fixed.
For the precise definition of what "generic" means here, we refer to the later Section 4.1.
We are going to exploit that distractions of torus-fixed ideals (such as indicial ideals) take on a very special form.The indicial ideal of I with respect to w is the distraction of the initial ideal, i.e., the C[θ 1 , . . ., θ n ]-ideal where θ i = x i ∂ i denotes the ith Euler operator.
Definition 3.8.The zeros of the indicial ideal are called the exponents of I with respect to w.
A D n -ideal F which is generated by elements of C[θ 1 , . . ., θ n ] is called Frobenius ideal.Frobenius ideals hence are of the form F = D n J with J an ideal in C[θ 1 , . . ., θ n ].A Frobenius ideal F = D n J is holonomic if and only if J is Artinian, i.e., if C[θ]/J is a finite-dimensional C-vector space.Now assume J is Artinian.The primary decomposition of J is of the form where Q A is an Artinian ideal, and Q A (θ − A) denotes the ideal obtained from replacing each θ i by θ i − A i in Q a .Solutions of Frobenius ideals take on the very special form x A • g(log(x)), which can be read from the primary decomposition of J. Namely, A ∈ V (J) and g runs over the finite-dimensional C-vector space the orthogonal complement of the Artinian ideal Q A .Since we will make use of that strategy regularly later on, we summarize what was said above in the following proposition.
Proposition 3.9.Let F = D n J, where J ⊂ C[θ], be a holonomic Frobenius ideal.The solution space of F is spanned by the functions x A • g(log(x)), where A runs over the points of the variety V (J), and g runs over the orthogonal complement Proposition 3.10.Let I be a holonomic ideal and w ∈ R n generic for I. Then ind w (I) is a holonomic Frobenius ideal whose rank equals rank(in (−w,w) (I)).
By N, we denote the ring of functions of the Nilsson class, i.e., those functions which can be represented by an element of for suitable vectors u 1 , . . ., u n , β 1 , . . ., β n ∈ C n (see [55, (2.31)]).The coefficients lie in the ring C x u 1 , . . ., x u n of formal power series in the x u i 's.
Definition 3.11.The w-weight of a monomial x A log(x) B is the real part ℜ(w • A) of w • A.
The initial series of a function f = A,B c AB x A log(x) B ∈ N, denoted in w (f ), is defined to be the finite subsum of all terms of minimal w-weight.Note that a weight vector w ∈ R n induces a partial order on functions of the Nilsson class as follows: (3.21) Since the w-weight does not give a monomial ordering, one needs a monomial order ≺ as a tie breaker and denotes the resulting monomial order by ≺ w .We will take ≺ to be the lexicographical ordering on N obtained as restriction of the lexicographic ordering on C n ⊕N n .The set of starting monomials of I with respect to ≺ w is where in ≺w (f ) = x A log(x) B for some A ∈ C n and B ∈ N n .Here, (3.23) , then A is an exponent of I with respect to w.For each exponent A, the number of starting monomials of the form x A log(x) B is the multiplicity of A as a root of the indicial ideal ind w (I).
Theorem 3.13 ([55, Theorem 2.5.12]).For each starting monomial , there exists a unique f ∈ N with the following properties: (1) f is annihilated by I, i.e., f ∈ Sol(I). ( (3) The monomial x A log(x) B is the only starting monomial that appears in f with non-zero coefficient.
Solutions to I as in the theorem above are called canonical (series) solutions of I with respect to ≺ w .The solution functions f in Theorem 3.13 can be shown to actually live in the Nilsson ring N w ⊂ N, which is the content of the next proposition.
and [55]).The elements of C C w (I) * Z are power series in x whose exponent vectors lie in C w (I) * Z .More precisely, the canonical solutions to I with respect to ≺ w have the form x A • g, where A is an exponent of I and g is an element of C C w (I) * Z [log(x 1 ), . . ., log(x n )], such that the degree of each log(x i ) in g is at most rank(I) − 1 (see [55,Theorem 2.5.14]).Definition 3.15.Let w ∈ R n .A finite subset G of a D-ideal I is called a Gröbner basis with respect to w if I is generated by G and in (−w,w) (I) ⊂ D is generated by the initial forms in (−w,w) (g), where g runs over the set G.
Moreover, if G is a Gröbner basis of I with respect to ≺ w , where ≺ denotes any term order on D (see [55, p.5]), G is also a Gröbner basis of I with respect to w.

The Saito-Sturmfels-Takayama algorithm
We now state the theorem which allows us to compute canonical solutions starting from a Gröbner basis of I with respect to a generic weight vector w ∈ R n .This is [55, Theorem 2.6.1], and we here recall it together with its proof, which contains the algorithm to lift solutions of the indicial ideal to canonical series solutions of I up to arbitrary w-weight.We are going to refer to this algorithm as SST algorithm.
Theorem 3.16.Let I be a regular holonomic ideal in Q[x 1 , . . ., x n ] ∂ 1 , . . ., ∂ n and let w ∈ R n be a generic weight vector for I. Let I be given by a Gröbner basis G w.r.t.w.There exists an algorithm which computes all terms up to specified w-weight in the canonical series solutions to I with respect to ≺ w .
For the convenience of our readers, we summarize the main steps of the SST algorithm as a procedure before turning to the proof.We also give an example of the algorithm running on a one-variable hypergeometric system.We here already mention Gröbner fans; the definition is contained in the later Section 4.1.
Procedure 3.17 (Computing canonical series solutions of a D-ideal up to a chosen order).Input: A regular holonomic D n -ideal I, its small Gröbner fan Σ in R n , a weight vector w ∈ R n that is generic for I, and the desired order k ∈ N.
Step 1. Compute the indicial ideal ind w (I) and its rank(I) many solutions.They are the form x A log(x) B with A ∈ V (ind w (I)), and will be the starting monomials of the canonical series solutions.For each starting monomial, carry out Steps 2-5.
Step 2. Determine a Gröbner basis G of I with respect to the weight w.
Step 3. Create the recurrences.By a recurrence, we mean a way of writing each element .., h r as matrices acting on the vector spaces L p , L p−β(1) , . . ., L p−β(r) , which contain the coefficient vectors c p of x p log(x) b , 0 ≤ b i < rank I.
Step 4. Apply the recurrences.Assuming one has the coefficient vectors c p−β(1) , ..., c p−β(r) , this amounts to solving the system of matrix equations to obtain c p .
Step 5.If the matrix of f is singular, then one may need to use the condition that for i = j, the series expansion s i must have coefficient 0 for the starting monomial of s j .
Output: The canonical series solutions of I with respect to w, truncated at w-weight k.
Step 1.The order of x with respect to w = 1 is −1.Thus in (−w,w) (I) = θ(θ −3) .The initial ideal is already torus-fixed, so ind w (I) = in (−w,w) (I).Solutions to the indicial ideal are x 0 = 1 and x 3 .We select the starting monomial Step 2. The ideal is principal, hence its generator is a Gröbner basis of the ideal.
Step 3. We write P = f − h, where f = θ(θ − 3) and h = x(θ + a)(θ + b).It suffices to compute the action of θ on each element of L p and extend it k[θ]-linearly.We have Thus, the matrix of the operator θ in the basis {x p+3 , x p+3 log(x)} is p + 3 1 0 p + 3 .
Step 4. Let c p,1 and c p,2 be the coefficients of x p+3 and x p+3 log(x) in the power series expansion.Then we can write our operators as matrices, and our recurrence as with initial values c 0 = 1, d 0 = 0. Solving the recurrence yields the explicit formulae where (a) p = a(a + 1) is the Pochhammer symbol.
Step 5.If we choose the starting monomial x 0 = 1 instead, the matrix of f is singular for p = 3.We leave it to the reader to find the series expansion, or to see [55, pp. 98-99].

⋄
We now turn to the proof of Theorem 3.16, making the steps of Procedure 3.17 precise.Proof and algorithm.Let w ∈ R n be generic for I. Compute the roots of ind w (I) and extend the field of coefficients by them.Denote the resulting, computable field extension of Q by K. To compute the canonical solution of I whose starting monomial is x A log(x) B , one proceeds as follows.For p ∈ Z n , denote by L p the K-vector space The monomials of L p are a K-basis of it.They are ordered by the term order ≺ w on the Nilsson ring, starting with the smallest.The matrix of f in this basis is an upper triangular square matrix.Let L ′ p denote the set of monomials in L p that are not contained in Start ≺w (I).Now let {f 1 , . . ., f d } be any generating set of ind w (I) and restrict f i : L p → L p to L ′ p ; this corresponds to deleting some of the columns in the associated matrix.Denote the resulting matrix by F i .Then the map is injective and is represented by the matrix obtained as vertical concatenation of F 1 , . . ., F d .Now let G = {g 1 , ..., g d } be a Gröbner basis of I with respect to w; its Gröbner cone in R n is denoted by C w .For each g ∈ G, choose a Laurent monomial x α such that where Here, ord (−w,w) (h) denotes the largest w-weight of a monomial appearing in h.Then the set of operators f 1 , . . ., f d obtained this way generate ind w (I) and, as maps as in (3.27), they are injective.The Laurent monomial x α for a Gröbner basis element g is obtained by taking a highest-weight term x a ∂ b m in g, where m is a monomial in the θ i 's and for all k, at least one of {a k , b k } is zero.Intuitively, this corresponds to pulling out as many θ's as possible into m.Then x α = x b−a .The h i may have terms of different weights; in this case, we get a recurrence that involves L p for more than two different p's.The coefficients of the canonical series solution are now computed by induction on the w-weight k.We start from a canonical series solution x A log(x) B + • • • and assume the coefficients c pb are already known for 0 and let M k be the space of terms with w-weight greater than k, i.e.M k = p•w>k,p∈C * Z L p .Then, by definition, we have that Assuming we know F k (x) for some k, we are going to construct a recursion which allows us to determine the additional terms which are needed to lift F k (x) to F k+1 (x).The starting point of that recursion will be the starting monomials.We hence look for an element To achieve this, observe that ord (−w,w) (h i ) < 0 implies that h i • x ℓ has higher w-weight than x ℓ .We can show this on monomials as follows.Suppose that h i = x q for some q ∈ N n with q • (−w) < 0. Then q • w > 0, so h i • x ℓ = x ℓ+q has higher w-weight than x ℓ .Similarly, suppose that h i = ∂ r where r • w < 0. Then (−r) • w > 0. Thus h i • x ℓ = Cx ℓ−r for some constant C, and x ℓ−r has higher w-weight than x ℓ .Together with (3.28) and (3.30), this implies that , which gives the desired recursion relation for By the injectivity of the map F from (3.27) and the existence of a canonical series solution, there exists a unique solution E k+1 (x) to (3.32), and this lifts F k (x) to F k+1 (x).
Remark 3.19.The choice of α in (3.28) is in general not unique, and sometimes there might not be an α at all which satisfies all of the conditions.For example, consider w = (−1, 0) and the operator 4 Computing solutions of I 3 with Gröbner techniques In this section, we compute the canonical series solutions of the D-ideal I 3 generated by the operators in (2.8) using the SST algorithm from the proof of Theorem 3.16.We begin in Section 4.1 by showing that I 3 is holonomic, and by computing some preliminary data, such as its holonomic rank, its singular locus, and its Gröbner fan.In Section 4.2, we discuss thoroughly the computation of the canonical series solutions for one of the cones of the Gröbner fan.Thanks to the symmetry of the system, the computations for the remaining two cones work analogously; we give more details for that at the end of the section.

The Gröbner fan of I 3
In order to employ the SST algorithm for computing the canonical series solutions, we first need to certify that the ideal is holonomic.Let I 3 (c) be the D-ideal generated by the operators in (2.6), with c = (c 0 , c 1 , c 2 , c 3 ).A computation in Singular [26,35] approves that I 3 (4, 2, 2, 2) is holonomic.The output of the Singular code below proves that I(c) is holonomic also for generic c.For that, we use the Weyl algebra with the degree reverse lexicographical order for The code uses the D-module libraries [5], and can be carried out by running the following lines, which also compute the holonomic rank, the singular locus, and the characteristic variety of I 3 (c).Since we compute over the field of rational functions in the parameters c, the results are valid for generic c.The characteristic ideal in (0,1) ]-ideal generated by the 8 operators and The cones C 1 , C 2 , and C 3 correspond to the following cones in R 3 : and ) ranges over any of its cones.The small Gröbner fan is the restriction of the Gröbner fan to Hence, the cones of the small Gröbner fan of I are in one-to-one correspondence with the initial ideals in (−w,w) (I), with w ∈ R n .
Definition 4.1.A weight vector is generic for I if it lies in an open cone of the small Gröbner fan of I.
To compute the small Gröbner fan of our system I 3 , our strategy is to first approximate the fan by computing the initial ideals for various lattice vectors, and finding where they change.Once we have a good enough approximation to guess our fan, we can verify our guess by computing initial ideals along the candidate rays.The system I 3 is homogeneous with respect to the vector (1, 1, 1).Thus, our small Gröbner fan really lives in R 3 /(1, 1, 1).More specifically: if our Gröbner fan is Σ, then Σ = Σ ′ × (1, 1, 1) where Σ ′ is a fan in R 2 .We choose the orthogonal basis v 1 = (1, 0, −1) and v 2 = (−1, 2, −1) of R 3 /(1, 1, 1).By explicitly computing the initial ideals of I 3 , we obtain Σ ′ as depicted in Figure 2. Note that this picture is not canonical; it depends on the choice of basis of R 3 /(1, 1, 1).
The D-ideal I 3 is moreover regular singular, as we will read later from the system (5.28).
Remark 4.2.The triangle integral-read as a function of the coefficients of the graph polynomial-is a solution of a GKZ system [27, Section 3.5], in which case the SST algorithm simplifies substantially.The conformal ideal I 3 (4, 2, 2, 2) is a restriction of that GKZ system to a subspace of special values of some of the coefficients.Our further investigations will not rely on that fact; we here provide a general implementation of the SST algorithm.⋄

Solutions to the indicial ideal
We will compute the initial and indicial ideal of the D-ideal I 3 from (3.16) with respect to the weight vector w = (−1, 0, 1) ∈ C 1 .In the basis v 1 = (1, 0, −1), v 2 = (−1, 2, −1) of Section 4.1, w is (−1, 0).From Figure 2, we read that w is contained in the relative interior of the cone C 1 in the Gröbner fan, hence w is generic for I 3 .The initial ideal in (−w,w) (I 3 ) is generated by the three operators We compute the indicial ideal Users of Singular may compute the initial ideal and a Gröbner basis with respect to the weight vector w with the following code.
LIB "dmod.lib";int c0 = 4; int c1 = 2; int c2 = 2; int c3 = 2; ring r = 0,( This returns the initial ideal as in (4.3).We deduce that the indicial ideal ind w (I 3 ) is generated by the three operators They are a Gröbner basis of in (−w,w) (I 3 ), and in this case in (−w,w) (I 3 ) = ind w (I 3 ).One reads that i.e., the only zero of the indicial ideal is the point (−1, 0, 0).For the next step of the algorithm, we need to compute solutions to ind w (I 3 ).These will become the starting monomials of the canonical series solutions to I. Recall that ind w (I 3 ) = DJ for J the C[θ 1 , θ 2 , θ 3 ]-ideal generated by the operators in (4.5).We need to compute the primary decomposition of J.In Singular, one obtains the primary decomposition with the library primdec by running the following code, where we encoded the Euler operators θ 1 , θ 2 , θ 3 as variables th1, th2, th3.
LIB "primdec.lib";ring r = 0,(th1,th2,th3),dp; ideal J = th1 + th2 + th3 + 1, th2^2, th3^2; list pr = primdecGTZ(J); pr; From the output, one reads the ideals coinciding with the Gröbner basis in (4.5).Recalling that V (J) = {A} with A = (−1, 0, 0), we read that Q A from (3.18) is the ideal generated by θ 1 + θ 2 + θ 3 , θ 2 2 , and θ 2 3 .The orthogonal complement Q ⊥ A from (3.19) hence is This is a finite-dimensional vector space spanned by the 4 polynomials Hence, by Proposition 3.9, the solution space to ind w (I 3 ) is spanned by the 4 functions In Macaulay2 [34], the generators of the indicial ideal in (4.5) and its solutions are obtained conveniently using the commands distraction and solvePDE as follows.Running this code outputs the following: As is explained in [2], the solutions to ind w (I 3 ) are then obtained by replacing dthi by log(x i ) and multiplying the resulting functions by x A = x −1 1 .We will take advantage of the homogeneity of the system to write it in fewer variables.Consider the change of variables Lemma 4.3.Every solution f (x 1 , x 2 , x 3 ) to I 3 can be written in the form where f (y 2 , y 3 ) denotes the function f (y 2 , y 3 ) = f (1, y 2 , y 3 ).
Therefore, we can write each canonical series in the form We will call an element of the form (log y 2 ) i (log y 3 ) j a monomial.We now encode our D-ideal in the new variables.Now D denotes the Weyl algebra in the y-variables, i.e., In the y-variables, we thus obtain the following basis of the solutions space of ind w (I 3 ): 1 log(y 2 ) log(y 3 ) . (4.15) After substituting back to the x i variables, these functions are-up to reordering and a sign-exactly the functions found in (4.10).Since each m i is entirely contained within a single L p , we have Start ≺w (m i ) = m i (cf.[55, Lemma 2.5.10]).

Lifting the initial monomials to solutions of I 3
We proceed using the algorithm described in Theorem 3.16.The space L p is 16-dimensional and defined as the span of elements of the form x −1 1 y p (log y 2 ) m (log y 3 ) n where p ∈ Z 2 is fixed and 0 ≤ m, n ≤ 3.In order to execute the algorithm, we must first find the matrix of the transformation from L ′ p to (L p ) d .As we will see below, this matrix consists of 3 vertically concatenated blocks, where each block is an upper triangular 16 × 16 matrix.First, we compute that a Gröbner basis G w of I 3 with respect to w is given by Switching to the y-variables, the operators in (4.16) turn into Replacing θ y 1 by −1 gives (4.18) Observe that This justifies to substitute θ y 1 by −1.We point out that the ideals generated by the operators in (4.17) and (4.18), respectively, are not the same ideals; but the quotients by them are isomorphic as D-modules on the algebraic torus, and hence their solution spaces are isomorphic.⋄ We now check that the order condition on the h i 's is satisfied.Indeed, for w = (−1, 0, 1), the variable Thus ord (−w,w) (h i ) = ord (−w,w) (y i ) < 0 for i = 2, 3, as desired.In other words, g 1 gives rise to a recurrence relation between L p and L p+(1,0) , while g 2 gives a recurrence relation between the spaces L p and L p+(0,1) .
We have implemented the SST algorithm from the proof of Theorem 3.16 in Sage for PDEs in two variables.Running it produces the canonical series solutions that are shown in Equation (4.19) at the very end of this section.
In summary, we achieved a basis of holomorphic solutions to I 3 purely from the Gröbner data of I 3 .We carried the computations out for a weight vector w from the cone C 1 of the Gröbner fan of I 3 .For our computations, we changed to y-variables.To refine the weight-ordering of w to a term ordering on the Weyl algebra, we used the degree reverse lexicographical ordering for In order to compute a basis of the solution space to I 3 , one could equally have started from cone C 2 or C 3 , resulting in different orderings.One would change the refining order to The Gröbner bases of I 3 w.r.t.weights from C 2 and C 3 are then obtained by a relabelling of the indices, and accordingly for the starting monomials and canonical series solutions.In Lemma 4.3, one would factor out x −1 2 for C 2 , and x −1 3 for C 3 .
We end this section by giving the canonical series solution of I 3 for w from C 1 obtained by our implementation of the SST algorithm.We display them from w-weight 0 to 4 on the next page.We recall that these four functions are a basis of the space of solutions to the D-ideal I 3 .In order to pick out a specific solution, such as the one corresponding to the triangle Feynman integral from which we started, one needs to prescribe rank(I 3 )-many initial conditions.This can be done in a number of ways.For example, we may numerically evaluate the triangle Feynman integral J triangle 4;1,1,1 at four arbitrary kinematic points. 3For each cone, the kinematic points x = (x 1 , x 2 , x 3 ) must be chosen so that consecutive truncations of the canonical series evaluated at x converge.In this way, the truncation of the canonical series at w-weight k provides an approximation of the solutions, with an error comparable to the size of the omitted terms of w-weight k + 1.For example, in the cone C 1 with weight vector w = (−1, 0, 1), we must have that |x 1 | ≫ |x 2 | ≫ |x 3 |.Comparing the numerical evaluations obtained using AMFlow [51] to a linear combination of the canonical series in (4.19) allows us to determine that We validated this against further numerical evaluations, as well as against the analytic solution in (2.9).

Series solutions using Wasow's method
In this section, we discuss how to construct series solutions of the ideal I 3 using the toolkit developed for computing Feynman integrals in high energy physics.We provide an alternative approach for constructing the canonical series solutions discussed in the previous sections without resorting to Gröbner basis techniques, and give a dictionary between the two methods.We also touch base with the problem of solving a first-order Pfaffian system of PDEs, namely a system of PDEs of the form where F (x) is a vector-valued function of x = (x 1 , . . ., x n ).This is a typical problem that arises in the evaluation of Feynman integrals.Often, one aims to approximate the solution F (x) for small values of some variable, say of x 1 .If the Pfaffian system (5.1) is in a manifestly Fuchsian form at x 1 = 0, i.e., if the matrices P i (x) have at most simple poles at x 1 = 0, there is an algorithm [65] to construct the asymptotic series of the solution F (x) around x 1 = 0; see e.g.[67] for a discussion in the context of Feynman integrals.They will be of the form up to an arbitrary power of x 1 , with m max a positive integer which depends on the specific Pfaffian system, and c k,m are Laurent series in x 2 , . . ., x n .The series in (5.2) is a particular case of (1.1): for the weight w = (1, 0, . . ., 0), the series in (1.1) is the truncation of (5.2).
A number of issues needs to be addressed in order to apply this method to D-ideals.First of all, we need to construct a Pfaffian system associated with the considered D-ideal.In [56, p. 23], it is explained how to achieve this by Gröbner basis computations.In Section 5.1, we present an alternative algorithm which relies on linear algebra only, and apply it to the D-ideal I 3 .This method is inspired by the problem of integration by parts (IBP) reduction for Feynman integrals [22,58].We spell this analogy out, and provide a dictionary of the relevant concepts.As a by-product, this method allows us to determine the holonomic rank and the singular locus of the D-ideal without computing Gröbner bases.The resulting Pfaffian system is in general not in a manifestly Fuchsian form; double or higher-order poles may appear in the DEs, preventing us from applying the algorithm from [65] to construct the asymptotic series of the solutions.Finding a gauge transformation which puts the Pfaffian system into Fuchsian form is a well known problem in Feynman integrals [37], and a number of strategies have been developed in that context, see e.g.[38,50].In Section 5.2, we argue that, if a D-ideal is regular holonomic, for every singular point x * it is possible to find a gauge transformation such that the associated Pfaffian system has at most simple poles at x = x * .In the case of I 3 , we show that we can even construct a Pfaffian system which is manifestly Fuchsian everywhere, namely which has at most simple poles at any singular point.The resulting system is particularly simple; we can solve it analytically, and by doing so, we will recover the solutions to I 3 given in (2.12).Finally, we need to fill an important conceptual gap.The canonical series computed by Gröbner bases are truncated by the w-weight of the terms.The asymptotic series expansions, on the other hand, rely on the notion of a "small parameter", e.g.x 1 in Equation (5.2), and are truncated at a given power of it.To link the two notions, in Section 5.3, we introduce an auxiliary variable t, which captures the w-weight of the monomials.We compute the asymptotic solutions of the Pfaffian system around t = 0, and verify that they coincide with the canonical series computed in Section 4.

Pfaffian systems
In this section, we present a method for constructing a Pfaffian system of PDEs associated with a given D-ideal using linear algebra.We are going to apply that method to the conformal D-ideal I 3 .A closely related method to efficiently construct Pfaffian matrices is discussed in [21].Our presentation builds on the analogy with IBP reduction, to the benefit of the readers who are familiar with Feynman integrals.
In Section 3.2, we have seen that the vector-valued function F (x), which satisfies a Pfaffian system (5.1)associated with a D n -ideal I with rank(I) = m is given by where f (x) is a general holomorphic solution of I, and p(1), . . ., p(m) ∈ N n such that {∂ p(1) , . . ., ∂ p(m) } are the standard monomials of the left R-ideal RI for a given monomial ordering.We may assume ∂ p(1) = 1, so that the first entry of F (x) is a holomorphic solution of I.The standard monomials can be determined by Gröbner bases in R, and the matrices P i (x) in the Pfaffian system (5.1) can be obtained through reduction in R, Our goal is to construct a Pfaffian system of I using linear algebra only.Let Q 1 , . . ., Q N be generators of I. First of all, the R-monomials ∂ p(i) defining F (x) as in (5.3) need not to be standard monomials for F (x) to satisfy a first-order Pfaffian system of PDEs.It suffices that they are a C(x)-basis of R/RI.One way to find such a basis is to write down the relations among the monomials ∂ a in R/RI, and solve them.There are infinitely many such relations in R/RI, obtained by multiplying the generators Q i of I by any ∂ a , We view this as a linear system of equations in the "variables" ∂ a .In Section 3.2, we have seen that, if I has finite holonomic rank m, the linear system (5.5) has m independent variables.In other words, any ∂ a can be expressed as a linear combination of m ∂ a 's in R/RI by solving the linear system (5.5).The issue is that there are infinitely many equations.We cannot solve them all, nor do we need to.We instead adopt an approach inspired by Laporta's algorithm [46] for solving integration-by-parts relations in the context of Feynman integrals.
We introduce an ordering > of the unknowns ∂ a graded by the total degree of the derivatives.The ordering among unknowns with the same degree is instead arbitrary.It will lead to a different basis of R/RI.For instance, we choose the graded lexicographic order.Explicitly, for n = 2, we have (5.6) Next, we multiply the generators Q i ∈ D by ∂'s from the left, up to a certain degree d max . 4his yields the equations where a ∈ N n runs over all vectors of natural numbers with We denote the resulting finite system of equations by S dmax .We dub this operation seeding, in analogy with the equivalent step in constructing IBP relations for Feynman integrals.Since we have truncated the system of equations, some unknowns cannot be solved for with the equations in S dmax .The key idea is to solve them only for a subset of the unknowns, in particular those with total degree up to some value d needed below d max (e.g.d needed = d max − 1), and to eliminate the other unknowns through Gaussian elimination.We denote by D the vector containing all unknowns ∂ a appearing in S dmax sorted according to the ordering defined above, and by D needed the needed unknowns, (5.8) By construction, D has the block form with D extra = D\D needed .We represent the equations in S dmax in matrix form as where M(x) is a [S dmax ] × [D] matrix whose entries depend on x.We row-reduce M(x) and drop the null rows.The resulting system of equations has the form Finally, we perform back substitution, but only for the needed unknowns.In other words, we solve the subset of equations (5.12) As a result, we obtain the expressions of the needed unknowns in terms of a subset of independent monomials ∂ a in R/RI.Let us make a few comments about this algorithm.
Construction of a Pfaffian system IBP reduction with Laporta's algorithm a set of master integrals Table 1: Analogies between the construction of a first-order Pfaffian system of DEs associated with a given D-ideal and IBP reduction of Feynman integrals.
(1) The value of d max must be higher than the maximal degree of the standard monomials.
In practice, we iterate the algorithm for various increasing values of d max until we see that the resulting independent monomials stay the same.These are a basis of R/RI.
(2) If the D-ideal I has infinite holonomic rank, the algorithm will not terminate.In practice, as we increase d max , we will obtain more and more independent monomials.Thus, this procedure is a linear algebra approach for computing the holonomic rank of a D-ideal.
(3) Let {∂ p(1) , . . ., ∂ p(m) } be a C(x)-basis of R/RI.If the needed unknowns include ∂ i ∂ p(j) for any i = 1, . . ., n and j = 1, . . ., m, the matrices P i (x) of the Pfaffian system (5.1) can be read off from the solution of the system of equations discussed above.
(4) This algorithm does not prove that the identified monomials ∂ a are a basis of R/RI.However, the resulting Pfaffian system of DEs must satisfy the integrability conditions (3.13).If the latter are not satisfied, either a higher value of d max must be used, or the ideal has infinite holonomic rank.If they are satisfied, the solution space has a C-dimension equal to the length of F (x) (see Section 3.2).Since the first component of F (x) is by construction a general solution of the original ideal I, we can conclude that the algorithm has terminated, and that a basis of R/RI has been identified.
(5) The common denominator of the entries of the matrices P i in the Pfaffian system (5.1)defines a hypersurface which contains the singular locus of I.
(6) There is a strong analogy with IBP reduction in Feynman integrals, which we summarized in Table 1.
We now apply this algorithm to the conformal ideal I 3 .It is convenient to use the variables y defined in (4.11), and to eliminate y 1 from the problem as follows.Lemma 4.3 guarantees that every solution of I 3 has the form f (y 2 , y 3 )/y 1 .Any function of such form is annihilated by the generator P3 in (2.8).The action of the other two generators can be expressed as follows: (5.13) Here, Q 1 and Q 2 are differential operators which depend on y 2 and y 3 only: (5.14) We denote by I y 3 the D-ideal generated by Q 1 and Q 2 .It is holonomic, it has holonomic rank 4, and its singular locus is where λ is obtained from the polynomial λ (2.10) as λ = λ| x 1 =1, x 2 =y 2 , x 3 =y 3 . (5.16) The general solution of I 3 is given by the general solution of I y 3 divided by y 1 .Hereinafter, we will thus focus on I y 3 .Since the latter depends on y 2 and y 3 only, we use the shorthand notations y = (y 2 , y 3 ), y a = y a 2 2 y a 3 3 , and ∂ a y = ∂ a 2 y 2 ∂ a 3 y 3 .We compute a basis of R/RI y 3 by running the algorithm described above.We use d needed = d max − 1.For d max = 1, we find only one basis element, namely 1. Starting from d max = 2 and onwards, we identify four basis elements.We define where f (y) is a holomorphic solution of I y 3 .It satisfies a system of first-order Pfaffian DEs and the 4 × 4 matrices Pi satisfy the integrability conditions.For example, P2 is the matrix From Equation (5.19), we see that the common denominator is y 2 y 2 3 λ, which is in agreement with the singular locus computed via Gröbner bases.Furthermore, some entries have double poles at y 3 = 0.The system is thus not in Fuchsian form.In the next section, we will perform a gauge transformation to obtain a Pfaffian system which is in a manifestly Fuchsian form.

Writing the system in Fuchsian form
If a D-ideal I is regular holonomic, all solutions have moderate growth when approaching the singular locus (cf.Section 3).This implies that the solutions cannot have essential singularities.This has implications for the form of the associated first-order Pfaffian system of DEs.A first-order ODE having poles of order higher than 1 implies that its solution has an essential singularity.The moderate growth condition on the solutions to regular holonomic D-ideals thus implies that the associated Pfaffian system has at most simple poles.The Pfaffian systems may however have "spurious" double poles.For example, the system exhibits a double pole at x = 0, yet the solution has no essential singularity at x = 0.Such spurious poles can be removed by a suitable gauge transformation, The goal is to construct the transformation matrix T (x) such that the new vector-valued function G(x) satisfies a first-order system of DEs which is manifestly Fuchsian at x = 0.For this simple rank-two example this can be achieved by which leads to The problem of finding a gauge transformation which puts a Pfaffian system into a manifestly Fuchsian form is central in the modern methodology for computing Feynman integrals [37].For our purposes, it suffices to put the Pfaffian system in a form which is manifestly Fuchsian at the singular point around which we compute the asymptotic expansion.If the entries of the required transformation matrix are in the field of rational functions, this problem is solved (see [38,50] and the references therein).For the Pfaffian system associated with the D-ideal I y 3 constructed in the previous subsection, we perform the gauge transformation with the following transformation matrix:5 As a result, G(y) satisfies a linear system of first-order DEs which is in Fuchsian form.It can be elegantly expressed as where In the Feynman integral literature, the arguments li of the logarithms in the Fuchsian Pfaffian system (5.27) are called symbol letters, and their ensemble { li } symbol alphabet [33].They encode the possible singular points of the solution.Indeed, the letters in Equation (5.29) vanish exactly on the singular locus of I y 3 given in (5.15).The system (5.27) is therefore in a manifestly Fuchsian form, i.e., it has at most simple poles at every singular point.In (5.29), the readers with some experience in Feynman integrals may recognize the alphabet of the three-mass triangle Feynman integrals [20], minus an additional letter equal to λ.The relation between this alphabet and conformal symmetry was already observed in [23], and is a promising hint that this alphabet may be valid at any loop order.The matrix in (5.28) exhibits the block triangular structure observed in [19] for the differential equations satisfied by finite Feynman integrals.We expect that differential equations corresponding to (5.27) can be obtained using the method of [19].
Before we move on to constructing the asymptotic series solutions in the next section, we point out that the form of the first-order Pfaffian system in (5.27) is particularly wellbehaved; its solutions can be written down explicitly in terms of logarithms and dilogarithms.The matrix B(y) (5.28) is upper block-triangular, which allows us to solve the system (5.27)iteratively.We obtain where f 1 (x), . . ., f 4 (x) are the solutions to I 3 from (2.12), and k 1 , . . ., k 4 ∈ C are arbitrary constants.Recalling the relation (5.25) between F (y) and G(y) with T (y) as in (5.26), and that the first component of F (y) divided by y 1 is a general solution to I 3 , we have proven that f 1 (x), . . ., f 4 (x) are indeed a basis of the space of solutions to I 3 .

Capturing the weight via an auxiliary variable
We now compute the asymptotic solutions to the first-order Pfaffian system (5.27).Since the latter is in Fuchsian form, we can use Wasow's algorithm [65] to expand around any singular point without resorting to Gröbner bases.In order to build a bridge to the canonical series solutions, which are defined with respect to a weight vector w ∈ R n , we introduce an auxiliary variable t as follows.Let f (x) be the general solution to the ideal under consideration.We define a new function f w , which depends both on the original variables x and on t as where t w x is a short-hand notation for (t w 1 x 1 , . . ., t wn x n ).
Remark 5.1.To motivate the construction in (5.31), consider the toy example where the function is a monomial, say f (x) = x a for some a ∈ N n .The exponent of t in the auxiliary function f w (t, x) = t w•a x a gives the w-weight of the monomial, namely w • a.In this sense, the exponent of the auxiliary variable t captures the notion of weight.⋄ We then compute the asymptotic expansion of f w (t, x) around t = 0, where m max is a natural number which depends on the specific ideal.By construction, the monomials in c k,m (x) have w-weight k.By definition, we have that f w | t=1 ≡ f.Hence, we expect that the asymptotic expansion (5.32) around t = 0, truncated at t k and evaluated at t = 1 equals the canonical series expansion truncated at w-weight k.
Example 5.2.Consider the simplest of the solutions to the ideal I 3 from Equation (2.12), where λ is the homogeneous polynomial of degree 2 defined in Equation (2.10).For the weight vector w = (−1, 0, 1) in cone C 1 , we introduce the auxiliary variable t via The asymptotic expansion around t = 0 is given by a Taylor expansion in t, We see that the monomials x a appearing as coefficients of t k have w-weight w • a = k.
Truncating the expansion (5.35) at order k and setting t = 1 gives the canonical series solution truncated at w-weight k, e.g., where the dots denote terms of w-weight 5 or higher.Indeed, this coincides with the series f1 (y 2 , y 3 ) in (4.19), upon substituting y in terms of x in the latter, and dividing it by x 1 to obtain a solution of I 3 through Lemma 4.

⋄
The derivatives of f w (t, x) with respect to x and t can be obtained from the derivatives of f (x) through the chain rule, so that we can straightforwardly construct a system of PDEs to which f w (t, x) is a solution, starting from that for f (x).In the next section, we will use this approach to compute the canonical series solutions to the Fuchsian system derived in Section 5.2, and verify that they reproduce those computed in Section 4.

Computations for the cone C 1
We now turn to the asymptotic solution of the Fuchsian system (5.27).We choose the weight vector w = (−1, 0, 1) from the cone C 1 of the small Gröbner fan.The corresponding weight for the y-variables is w y = (−1, 1, 2).We thus introduce the auxiliary variable t via where G(y 2 , y 3 ) is a holomorphic vector-valued function which satisfies the Fuchsian Pfaffian system (5.27).The auxiliary function Gw (t, y) satisfies a Pfaffian system in Fuchsian form; for i = y 2 , y 3 , t, The matrices Bi (t, y) are obtained from the derivatives of G(y) through the chain rule, In order to compute the asymptotic expansion of the solutions around t = 0, we need to study the behavior of the matrices Bi (t, y) around t = 0.The Laurent expansion of Bt (t, y) is (5.40) Since the system is in Fuchsian form, there can be no pole of higher order.The residue at t = 0, i.e., the constant matrix Bres t , has the eigenvalue 0 only.The matrices By 2 (t, y) and By 3 (t, y) are instead non-singular at t = 0.
The first step of the algorithm in [65] is to perform a gauge transformation Gw → G′ w , Gw (t, y) = Ũ (t, y) • G′ w (t, y) . (5.41) In the new basis, our vector-valued function satisfies the system for i = y 2 , y 3 , t, with The key idea is to choose the gauge transformation such that the DEs for G′ w (t, y) are simpler than those for Gw (t, y) in view of the asymptotic solution around t = 0.In particular, we wish to simplify B′ t (t, y) so that it has a simple pole only, i.e., is of the form B′ t (t, y) = Bres t t . (5.44) We can construct the transformation matrix Ũ (t, y) which achieves (5.44) as a series in t,6 Ũ (t, y) = k≥0 t k Ũk (y) . (5.45) We plug this expansion into Equation (5.43) and impose that the resulting B′ t (t, y) satisfies (5.44).The first term in the expansion, Ũ0 (y), must commute with Bres t .The simplest way to achieve this is to choose Ũ0 (y) = I 4 to be the 4 × 4 identity matrix.The terms of expansion (5.45) of higher order are determined recursively in terms of the lower ones through the solution of a linear system of equations, (5.46) As an example, we spell out the first few terms of the resulting transformation matrix: (5.47) As expected, each order in t involves only monomials whose w-weight matches the power of t.As a result, G′ w (t, y) satisfies the following simplified system of PDEs, with B′ y 2 (t, y) and B′ y 3 (t, y) determined by (5.43), and non-singular at t = 0. We can solve the simplified system (5.48) using the path-ordered exponential formalism (see e.g.[14] and the references therein).We integrate this system along a path in the (t, y)-space of the form 0 , y (0) −→ (0 , y) −→ (t, y) , ( for some arbitrary values y (0) of y which do not belong to the singular locus of the ideal.In other words, we pick (t = 0, y = y (0) ) as boundary point, 7 restore the dependence on y by integrating along the first piece of the path, and then on t by integrating along the second.
and h(y) is the vector-valued holomorphic function which solves the system (5.52) Explicitly, we have We now have all the ingredients to compute the asymptotic expansion around t = 0 of the solution of the first-order Pfaffian and Fuchsian system of DEs (5.27).It is given by Gw (t, y) = Ũ(t, y) • e Bres t log(t) • h(y) . (5.54) The gauge transformation matrix, Ũ (t, y), is given by a Taylor series around t = 0 starting with the identity matrix, and can be computed algorithmically up to any order in t.The matrix exponential, e Bres t log(t) , is irrelevant for our purposes, as we will eventually set t = 1, and e Bres t log(1) = I 4 .The matrix-valued function h(y) given in Equation (5.53) is therefore the only source of powers of log(y) in this approach.
We can now write down the asymptotic expansions of the solution to the ideal I y 3 .We recall that a solution to the latter, f (y), is by construction given by the first component of the vector-valued function F (y) (see (5.17)).The latter is related to the vector-valued function G(y) through the gauge transformation (5.25) with transformation matrix T (y).Putting everything together, we have where the subscript 1 denotes the first component of the vector.The series expansion of G(y) truncated at w-weight k is given by the asymptotic expansion of Gw (t, y) around t = 0 up to order k, given by (5.54), setting t = 1.Similarly, the series expansion of the transformation matrix is obtained by the Taylor expansion of Tw (t, y) = T (t w y) around t = 0, truncated at t k and evaluated at t = 1.Then, the four linearly independent series solutions to I y (5.56) In the above functions, the square brackets highlight terms of different weight, and the dots denote terms of w-weight 4 or higher.From these, we see that the starting monomials coincide with the solutions to the indicial ideal in (4.15).The series expansions in (5.56) match those computed by Gröbner bases in (4.19) (with hi (y) = fi (y)), and are related to the known solutions f i from (2.12) via (5.57) We thus find agreement among the known solutions, the canonical series solutions computed via the SST algorithm, and those computed by solving the Pfaffian system asymptotically.
The computations for the remaining cones of the Gröbner fan follow the same strategy.The results can be obtained by relabelling the indices as discussed at the end of Section 4.

Conclusion and outlook
In this paper, we studied algorithms for obtaining solutions of linear partial differential equations relevant for particle physics.In principle, this applies to all Feynman integrals, as the latter are known to satisfy systems of differential equations (see e.g.[52]), as well as to further special functions appearing as their solutions, such as multiple polylogarithms.
We implemented a method, originally proposed by Saito, Sturmfels, and Takayama [55], to compute canonical Nilsson series expansions.The SST algorithm had been in use already [27,57] for GKZ systems.To our knowledge, our paper is the first one to implement the algorithmic ideas of SST without requiring hypergeometricity.Using the particular example of conformal differential equations satisfied by a one-loop triangle Feynman integral with general propagator powers, we implemented the SST method and compared it to a method of Wasow [65] that is more commonly used in the context of differential equations for Feynman integrals [38,50].
Using the example of the triangle Feynman integral, we showed how methods [55,56] from the theory of D-modules can be used to determine key features of the differential equations, as well as canonical series expansions of the solutions.In particular, for holonomic systems, the holonomic rank corresponds to the number of master integrals in the physics literature.The singular locus determines the set of (possible) singularities of the corresponding Feynman integrals.The solutions to the indicial ideal are the starting terms of the series.Following SST, we use Gröbner basis methods to obtain canonical series expansions.Leveraging the techniques developed for computing Feynman diagrams, we present alternative methods which allow us to extract this precious information from the PDEs without relying on the computation of Gröbner bases.
It is interesting to ask about the scope of the method.Proxies for the complexity of the differential equations are: the number of variables; the order of the differential operators; the number of differential equations.A potential bottleneck in the SST approach could be the reliance on Gröbner bases.However, compared to other situations where Gröbner methods are used (e.g. for analyzing equations of high polynomial degree), one would expect the polynomials appearing in the PDEs of physical interest to have moderate degree.(As a proof of principle, we showed in Appendix C an application to a four-loop ladder integral.)Moreover, our comparison with Wasow's method shows that, in principle, these methods can be replaced by linear algebra operations, which could potentially scale better.Then again, Wasow's method relies on the Pfaffian system being in Fuchsian form at the regular singular point around which we wish to expand.While several strategies to achieve this have been developed especially in the context of Feynman integrals, they have limitations.This makes it even more interesting to have two complementary approaches.
There are various interesting practical and conceptual questions for follow-up work.On the practical side, the physics literature provides many important classes of Feynman integrals that are relevant for the phenomenology of elementary particles, in the context of gravitational wave physics, or in cosmology, for example.A prerequisite of our method is the knowledge of the relevant differential equations.Obtaining the latter is an interesting active area of research in itself, see e.g.[45] and the references therein.It would be interesting to study the scope of the SST algorithm for the Picard-Fuchs equations obtained in [45].
On the conceptual side, let us mention the following questions.Firstly, when dealing with equations in multiple variables, one may ask how one can restrict the differential equations to lower-dimensional submanifolds.Physical examples include on-shell limits, or restrictions to submanifolds corresponding to (spurious) singularities.Secondly, there are important situations where only a subset of the relevant differential equations is available (e.g. because they are easier to obtain than the complete equations, as e.g. in [29]).In other words, in these situations, the holonomic rank of the D-ideal is not finite, and hence there is functional freedom in the solutions.This happens for example for the differential equations that follow from conformal or Yangian symmetry beyond the three-particle case.It is then interesting to study how holonomic techniques can be used to obtain useful information, such as which additional constraints may be added to make the system holonomic.Finally, the triangle Feynman integral considered here has a finite value, but Feynman integrals typically exhibit divergences in the infrared and ultraviolet regions of the loop integration.In dimensional regularization, these divergences are regulated by analytically continuing to generic d ∈ C spacetime dimensions.The divergences are then manifested as poles in the Laurent expansion around d = d 0 for some integer d 0 (typically d 0 = 4).The coefficients of this Laurent expansion are solutions to regular holonomic D-ideals in the kinematic variables and, as such, can in principle be treated using the D-module techniques discussed here.It would therefore be important to find a way to set up the expansion around d = d 0 within the D-module approach.
Beyond these specific research questions, we believe that the exchange of methods promoted in this work will be fruitful for both the communities of mathematicians and theoretical physicists.Given the ubiquity of PDEs, we envision these methods will be beneficial for many further applications.

A Conformal symmetry
Conformal symmetry underlies many of the quantum field theories which are important to our understanding of fundamental physics: massless quantum electrodynamics, quantum chromodynamics, Yang-Mills theory, scalar φ 4 theory, and more, all have conformal symmetry at the classical level. 8Despite its ubiquity, little is known about its implications for on-shell scattering amplitudes.The latter are functions of the particles' momenta which play a central role in describing the interactions of fundamental particles.One of the obstacles which has hindered the study of this topic is that the constraints imposed by conformal symmetry have the form of second-order differential operators in momentum space.Determining the conformally-invariant functions relevant for the scattering of particles thus requires to solve systems of linear second-order PDEs.Given the large number of variables which describe the particles' scattering, this problem is challenging from a technical point of view.As such, it is an optimal testing ground for D-module techniques.
The conformal symmetry group is an extension of the Poincaré group SO(1, 3)⋉R, which is the semi-direct product of the Lorentz group and the Abelian group of translations.The Poincaré group plays a central role in fundamental physics, as it is the symmetry group of Einstein's theory of special relativity, and hence of all theoretical models which stem from it.The conformal group is obtained by extending the Poincaré group by two classes of transformations: dilatations and conformal boosts (alias "special conformal transformations").We refer the interested readers to [28] for a thorough discussion of conformal symmetry in field theory, and present here what is necessary to introduce the system of PDEs addressed in this work.We begin by introducing some notation.We denote by the vector of d-dimensional spacetime coordinates, where z 0 is the time coordinate, and z i for i = 1, . . ., d − 1 are the spatial coordinates.We define the scalar product of two coordinate vectors z 1 and z 2 as The conformal transformations are given by the composition of the coordinate transformations shown in Table 2. 9 They form a group with respect to composition.A function is conformally-invariant if it is annihilated by the generators of the conformal group.The latter are first-order differential operators in the spacetime coordinates which obey the commutation relations of the Lie algebra associated with the conformal group.The generator of the group of dilatations for n points of coordinates z 1 , . . ., z n is for instance given by where the c k ∈ R are parameters called conformal weights, and ∂ z k denotes the vector of partial derivatives, Scattering amplitudes are functions of the particles' momenta.Like the spacetime coordinates z k , the momenta p k are d-dimensional vectors.Their zeroth components, p 0 k , give the particles' energy, while the other components, p µ k with µ = 1, . . ., d − 1, are related to their velocity.Coordinate and momentum space are related by a Fourier transform.
The momentum-space realization of the generators of the Poincaré group are given by first-order operators also in the momenta p k .The constraints they impose are particularly simple.The invariance under translations implies the conservation of total momentum, namely that p 1 + • • • + p n = 0. We can use this constraint to eliminate one of the momenta, say p n .The invariance under Lorentz transformations instead implies that the invariant functions depend on the momenta only through the scalar products p k •p ℓ , called Mandelstam invariants in the physics literature.Therefore, a generic Poincaré invariant function of n momenta p 1 , . . ., p n depends only on n(n − 1)/2 variables, which can be chosen as where the integration domain is the Minkowski spacetime (see Section 2).The corresponding Feynman graph is shown in Figure 3.The normalization factor of (|p 1 + p 2 | 2 ) 2(6−d) |p 2 + p 3 | 2 is chosen to make the integral dimensionless and simplify the results.We take the momentum p 4 to be off-shell, i.e.Thanks to the chosen normalization, the ladder integral depends only on y 2 and y 3 .The state of the art for Feynman integrals with these kinematics is three loops, cf.[63,18,39].Any information about its analytic structure would therefore be of great interest.For this reason, we focus on its maximal cut in d = 4 dimensions.The maximal cut-which amounts to replacing all propagators, 1/|q| 2 for some momentum q, with delta functions δ(|q| 2 )-captures precious analytic information about the integral.In particular, it is central in the "method of differential equations" (see e.g.[38] for a review) to obtain the so-called "canonical form" [37], which then allows for a systematic solution in terms of special functions.
In our presentation here, we take a more general approach to derive differential operators which annihilate the integral.Instead of using conformal symmetry, we construct the Picard-Fuchs operators [45,52].We do this by following a path which is the reverse of what we discussed in Section 3.2.First, we construct a Pfaffian system of PDEs for the ladder Note that, for the given generators, the D-ideal P 1 , P 2 is already too complicated for Singular to compute a Gröbner basis within a couple of minutes.But, one can check that P 2 ∈ W (I) \ I. Indeed, I W (I) = W ( P 1 , P 2 ) as D-ideals.Their singular loci are Sing(I) = Sing(W (I)) = V (y 2 y 3 (y 2 + y 3 )) . (C.6) The operator G 1 encodes that the solution functions f (y 2 , y 3 ) of the D 2 -ideal I are homogeneous of degree 0, and therefore f (y 2 , y 3 ) = f (y 2 /y 3 , 1).This is not obvious by looking at the Feynman integral, and is therefore a powerful insight of the D-module approach.
Remark C.2.It is possible to deduce that the system can be rewritten in one variable purely from the geometry of the Gröbner fan.The Gröbner fan has rays generated by the vectors (1, 1) and (−1, −1), and is a one-dimensional fan in R 2 , perpendicular to (1, −1).The only cones in the dual fan are R ≥0 • (1, −1) and R ≥0 • (−1, 1).Thus, the SST algorithm tells us that there are two families of series expansions in y 2 and y 3 , with exponents only in N • (1, −1) and N • (−1, 1), respectively.From this, we see that we can rewrite these series expansions in the single variable y = y 2 y −1 3 .⋄ We hence reduce the system to a single variable, namely the ratio y := y 2 /y 3 .To do so, we change variables from (y 2 , y 3 ) to (y, z) via For solving the associated differential equation, we can ignore the pre-factor of z −1 and hence are in the ODE case, with a single variable y.
Since the system is in a single variable, we will only need to distinguish between w being positive or negative, and can hence restrict to the weight being w = ±1.The starting monomials for w = 1 can be obtained from the initial ideal in (−1,1) (G y 2 ) = θ 4 .They are 1, log(y), log(y) 2 , log(y) 3 . (C.10) As in Example 3.18, we encode the coefficients of our Puiseux series in a vector c p of length 4 specifying the coefficients of y p , y p log(y), y p log(y) 2 , and y p log(y) 3 .As before, we can encode the action of θ on the series as multiplication of the coefficient vector c p by the matrix One can easily verify that these are solutions to the operator G y 2 .Similarly, for w = −1 we obtain the exponents {0, −1} and starting monomials 1, log(y), y −1 , y −1 log(y) . (C.16) We use the recursion Here, there is a slight complication: with starting monomial log(y), the matrix F p has a singularity at a = 0 and p = −1.However, we use the definition of canonical series to impose that neither of the terms y −1 , y −1 log(y) appears in the series expansions of the remaining solutions.Explicitly, we set the vector c −1 equal to zero.As before, expressing the recurrences in closed form gives us the Laurent series

Example 3 . 1 .
The Airy equation f ′′ (x) − xf (x) = 0 (3.4) is encoded by the operator P = ∂ 2 − x ∈ D. The Airy functions Ai and Bi are solutions of P , i.e., P • Ai = P • Bi = 0.They span the 2-dimensional solution space of P .⋄ A (left) D-ideal is a subset I ⊂ D that is closed under addition and multiplication by elements of D from the left.They encode systems of linear PDEs.For a vector

Definition 3 . 5 .
The Weyl closure of a D n -ideal I, denoted W (I) is the D n -ideal W (I) := R n I ∩ D n .(3.10) Clearly, I ⊆ W (I). Hence, for the singular locus and the space of holomorphic solutions to the system of PDEs encoded by I, one has Sing(I) ⊇ Sing(W (I)) and Sol(I) ⊇ Sol(W (I)) .(3.11)Moreover, rank(I) = rank(W (I)), since R n I = R n W (I). Since every element Q of W (I) can be written as Q = r • P for some r ∈ C(x 1 , . . ., x n ) and P ∈ I, we also have the inclusion Sol(I) ⊆ Sol(W (I)).Hence, Sol(I) = Sol(W (I)).In summary, a D-ideal I and its Weyl closure W (I) have the same solution space, but W (I) might contain additional operators that annihilate all solutions of I.If rank(I) < ∞, it does not follow that I itself is holonomic; but there is a rescue: taking the Weyl closure turns I into a holonomic D-ideal.
Program (WASP) funded by the Knut and Alice Wallenberg Foundation.This project received funding from the European Union's Horizon 2020 research and innovation programs Novel structures in scattering amplitudes (grant agreement No. 725110) and High precision multi-jet dynamics at the LHC (grant agreement No. 772099), and from the European Union's Horizon Europe Research and Innovation Programme under the Marie Sk lodowska-Curie grant agreement No. 101105486.This research was supported by the Munich Institute for Astro-, Particle and BioPhysics (MIAPbP), which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy -EXC-2094 -390783311.