Integration-by-parts identities and differential equations for parametrised Feynman integrals

Integration-by-parts (IBP) identities and differential equations are the primary modern tools for the evaluation of high-order Feynman integrals. They are commonly derived and implemented in the momentum-space representation. We provide a different viewpoint on these important tools by working in Feynman-parameter space, and using its projective geometry. Our work is based upon little-known results pre-dating the modern era of loop calculations: we adapt and generalise these results, deriving a very general expression for sets of IBP identities in parameter space, associated with a generic Feynman diagram, and valid to any loop order, relying on the characterisation of Feynman-parameter integrands as projective forms. We validate our method by deriving and solving systems of differential equations for several simple diagrams at one and two loops, providing a unified perspective on a number of existing results.


Introduction
The calculation of high-order Feynman integrals is the cornerstone of the precision physics program at present and future particle accelerators [1].The systematic development of modern methods to compute Feynman integrals, going beyond a direct evaluation of their parametric expression, began with the identification and explicit construction of Integration-by-Parts (IBP) identities in dimensional regularisation, in Refs.[2,3], and reached a further degree of sophistication with the development of the method of differential equations [4][5][6].These two sets of ideas can be combined into powerful algorithms [7], and the procedure further streamlined and optimised by the identification of the linear functional spaces where (classes of) Feynman integrals live [8][9][10], and by taking maximal advantage of dimensional regularisation [11].The combined use of these tools has dramatically extended the range of processes for which high-order calculations are available, and has broadened our understanding of the mathematics of Feynman integrals, as reviewed for example in [1,12].
It is an interesting historical fact that the idea of studying and eventually computing Feynman integrals by means of IBPs and differential equations pre-dates all the developments just discussed, and was originally proposed not in the momentum representation, but in Feynman-parameter space.Of course, it is well-known that studies of Feynman diagrams flourished in the S-matrix era, as illustrated in the classic textbook [13].In particular, the projective nature of Feynman parameter integrands, and the importance of the monodromy properties of Feynman integrals under analytic continuation around their singularities, were soon uncovered, and attracted the attention of mathematicians [14,15] and physicists [16].In this context, Tullio Regge and collaborators published a series of papers [17][18][19] studying the 'monodromy ring' of interesting classes of Feynman graphs: first the ones we would at present describe as 'multi-loop sunrise' graphs in Ref. [17], then generic one-particle irreducible n-point one-loop graphs in Ref. [18], and finally the natural combination of these two classes, in which each propagator of the one-loop n-point diagram is replaced by a k-loop sunrise [19].All of these papers employ the parameter representation as a starting point, and make heavy use of the projective nature of the integrand.
At the time, these studies by Regge and collaborators did not immediately yield computational methods, but it is interesting to notice that, at least at the level of conjectures, several deep insights that have emerged in greater detail in recent years were already present in the old literature.For example Regge, in Ref. [16], argues, on the basis of homology arguments, that all Feynman integrals must belong to a suitably generalised class of hypergeometric functions, an insight that was sharpened much more recently with the introduction of the Lee-Pomeransky representation [20] of Feynman integrals and the application of the GKZ theory of hypergeometric functions [21][22][23][24][25][26].Regge further argues that such functions obey sets of (possibly) high-order differential equations, which he describes as 'a slight generalisation of the well-known Picard-Fuchs equations', also a recurrent theme [27].
While general algorithms were not developed at the time, two of Regge's collaborators, Barucchi and Ponzano, were able to construct a concrete application of the general formalism for one-loop diagrams [28,29].In those papers, they show that for one-loop diagrams it is always possible to organise the relevant Feynman integrals into sets (that we would now call 'families'), and find a system of linear homogeneous differential equation in the Mandelstam invariants that closes on these sets, with the maximum required size of the system being 2 n − 1 for graphs with n propagators1 .These systems of differential equations were of interest to Barucchi and Ponzano because they effectively determine the singularity structure of the solutions, and thus the monodromy ring, in agreement with the general results of Regge's earlier work.From a modern viewpoint, it is perhaps just as interesting to use the system directly for the evaluation of the integrals, as done with the usual momentum-space approach: this is the direction that we will pursue in our exploratory study.
In the present paper, we start from the ideas of Refs.[16][17][18][19] and the concrete results of Barucchi and Ponzano [28,29] to propose a projective framework to derive IBP identities and systems of linear differential equations for Feynman integrals.In order to do so, we need to generalise the Barucchi-Ponzano results in several directions.First of all, those results predate the widespread use of dimensional regularisation, and do not in principle apply directly to infrared-divergent integrals.Fortunately, the projective framework naturally involves the (integer) powers of the propagators appearing in the diagram.These can be continued to complex values, providing a regularisation that is readily mapped to dimensional regularisation 2 .We are then able to show that the projective framework applies directly to IR divergent integrals, and we provide some examples.Next, we observe that the procedure to derive IBP identities in projective space generalises naturally to higher loops.Clearly, at two loops and beyond it would be of paramount interest to have a generalisation of the Barucchi-Ponzano theorem, guaranteeing the closure of a system of linear differential equations, providing an upper limit for its size, and giving a constructive procedure to build the system.This would require a much deeper understanding of the monodromy ring of higher-loop integrals.Lacking this knowledge (a gap which certainly points to promising avenues for future research), we can nonetheless apply the parameter-space IBP technique, and derive directly sets of differential equation on a case-by-case basis.Indeed, we show that the method can be successfully applied to two-loop integrals, and we provide examples, including the two-loop equal-mass sunrise, for which we recover the appropriate elliptic differential equation.Finally, we note that our application of the projective framework highlights the importance of boundary terms in IBP identities: contrary to the momentum-space approach in dimensional regularisation, boundary terms do not in general vanish in the projective framework: on the contrary, they may play an important role in linking complicated integrals to simpler ones, as we will see in concrete examples.
We note that the work presented in this paper is part of a recent revival of interest in the mathematical structure of Feynman integrals in parameter space, and presents interesting potential connections to several current research topics in this context, including intersection theory [33][34][35][36][37], the concept of parametric annihilator [30], the use of syzygy relations in reduction algorithms [38,39], the study of generalised hypergeometric systems [40], and the reduction of tensor integrals in parameter space [41][42][43].More generally, for the first time in several decades we are witnessing a rapid growth of our understanding of the mathematical properties of Feynman integrals, in particular with regards to analiticity and monodromy (see, for example, [32,[44][45][46], and the lectures in Ref. [47]), with potential applications to questions of phenomenological interest, such as the study of infrared singularities [48] and the development of efficient methods of numerical integration [49,50].
The structure of our paper is the following.In Section 2 we introduce our notation for the parameter representation of Feynman integrals and for Symanzik polynomials, briefly reviewing well-known material for the sake of completeness.In Section 3 we introduce projective forms, and we use their differentiation and integration to lay the groundwork for the construction of IBP identities for generic projective integrals.In Section 4 we specialise our discussion to Feynman integrals, and give a general procedure to construct IBP identities in this case.In Section 5 and in Section 6 we validate our results by discussing several concrete examples at one and two loops.Four Appendices give some further technical details on these examples.Finally, in Section 7 we present an assessment of our results and perspectives for future work.

Notations for parametrised Feynman integrals
In this section we summarise some well-known basic properties of parametrised Feynman integrals, which will be useful in what follows.We adopt the notations of Refs.[12,51,52].
Consider a connected Feynman graph G, with l loops, n internal lines carrying momenta q i and masses m i (i = 1, . . ., n), and m external lines carrying momenta p j (j = 1, . . ., m).At this stage we do not need to impose restrictions on external masses, so p 2 j is unconstrained.On the other hand, momentum is conserved at all vertices of G, so one can parametrise the graph assigning l independent loop momenta k r (r = 1, . . ., l) to suitable edges of the graph.The line momenta are then given by where the elements of the incidence matrices, α ir and β ij , take values in the set {−1, 0, 1}.Working in d dimensions, with d = 4 − 2ϵ, and allowing for the possibility of raising propagators to integer powers ν i (i = 1, . . ., n), one may associate to each graph G a family of (scalar) Feynman integrals where we defined ν ≡ n i=1 ν i , and the integration must be performed by circling the poles in the complex plane of the loop energy variables according to Feynman's prescription.
The integration over loop momenta in Eq. (2.2) can be performed in full generality by means of the Feynman parameter technique, using the identity By virtue of Eq. (2.1), the sum in the denominator of the integrand in Eq. (2.3) is a quadratic form in the loop momenta k r , and can be written as where M is an l × l matrix with dimensionless entries which are linear in the Feynman parameters z i , Q is an l-component vector whose entries are linear combinations of the external momenta p j , and J is a linear combination of the Mandelstam invariants p i • p j and the squared masses m 2 j .Translational invariance of d-dimensional loop integrals allows to complete the square in Eq. (2.4): the integral over loop momenta can then be performed, leading to where the functions are called graph polynomials or Symanzik polynomials.References [12,51,52] discuss in detail the properties of graph polynomials: here we only note that both polynomials are homogeneous in the set of Feynman parameters, z i , with U being of degree l and F of degree l + 1; furthermore, both polynomials are linear in each Feynman parameter, with the possible exception of terms proportional to squared masses in F. These homogeneity properties set the stage for employing the tools of projective geometry, as discussed below in Section 3. Remarkably, Symanzyk polynomials can be constructed directly from the connectivity properties of the underlying Feynman graph.To do so, let us denote by I G the set of the internal lines of G, each endowed with a Feynman parameter z i .A co-tree T G ⊂ I G is a set of internal lines of G such that that the lines in its complement T G ⊂ I G form a spanning tree, i.e. a graph with no closed loops which contains all the vertices of G.The first Symanzik polynomial for the graph G is then given by (2.7) Note that, in the case of an l-loop graph, one needs to omit precisely l lines in order obtain a spanning tree: the polynomial U is therefore homogeneous of degree l, as announced.Similarly, we can consider subsets C G ⊂ I G with the property that, upon omitting the lines of C G from G, the graph becomes a disjoint union of two connected subgraphs.Clearly, each subset C G defines a cut of graph G, and contains l + 1 lines.One may further associate with each cut the invariant mass s (C G ), obtained by squaring the sum of the momenta flowing in (or out) one of the two subgraphs -by momentum conservation, it does not matter which subgraph we choose.The second Symanzik polynomial is then defined by As expected, F is homogeneous of degree l + 1 in the Feynman parameters.
To illustrate these rules, consider the one-loop box diagram depicted in Fig. 1a.As for any one-loop diagram, it is immediate to see that the first Symanzik polynomial is simply the sum of the Feynman parameters associated with the loop propagators.In this case (2.9) The second Symanzik polynomial depends on kinematic data.If for example one picks massless on-shell external legs, all cuts involving two adjacent propagators vanish.One is then left with the Cutkosky cuts in the s and t channels.Defining s = (p 1 + p 4 ) 2 and t = (p 1 + p 2 ) 2 (with all momenta incoming), and assuming all internal masses to be the same, one finds (2.10) At two loops, one may consider the sunrise diagram in Fig. 1b.In this case, each internal line is a spanning tree (complementary to a co-tree).This implies that the first Symanzik polynomial is (2.11) The cut-dependent part of the second graph polynomial is similarly straightforward, since only one cut exists.Taking equal internal masses, and denoting by p 2 the invariant mass of the incoming momentum, the second Symanzik polynomial is thus In what follows, the crucial property of Eq. (2.5) is the projective nature of the integrand.Indeed, one easily verifies that a change of variables of the form z i → λz i , with λ > 0, leaves the integrand invariant, except for the argument of the δ function.Since a change of variables cannot affect the integral, we see that one should properly look at Eq. (2.5) as the integral of a projective form over the (n − 1)-dimensional space PR n−1 .This statement will be further substantiated in the next sections: in Section 3 we will show some technology concerning such integrals, which will then lead to a general integration-by-parts formula for Feynman parameter integrals in Section 4.

Projective forms
In this section, we present a brief introduction to projective forms and to their integration and differentiation.Since the section is somewhat formal, it is useful to keep in mind from the beginning the announced correspondence between projective forms and Feynman-parameter integrands, which we will try to highlight with explicit examples.

Preliminaries
Let us begin by considering the Grassman algebra of exterior forms in the differentials dz i , where i ∈ D ≡ {1, 2, ..., N }, for a positive integer N .Let A be a subset of D, of cardinality |A| = a, and let ω A be its ordered volume form with i j ∈ A, and The volume form ω A can be 'integrated' by defining where A − i denotes the set A with i omitted, and we defined the signature factor ϵ k,B , for any B ⊆ D, and for any k / ∈ B, by means of Using the properties of the boundary operator d, one easily verifies that the differential of η A is proportional to ω A .Indeed As an example, consider again A = {1, 2, 3}: the form η A is then given by and its differential in fact is equal to 3 dz 1 ∧ dz 2 ∧ dz 3 .Consider next affine q-forms, defined by where R A is a homogeneous rational function 3 of the variables z i with degree −|A| = −q.The name affine form is a reference to their invariance under dilatations of all variables.Eq. (3.6) is readily seen to imply that also the (q + 1)-form dψ q is affine.Anticipating Section 6.1, an example of an affine form with q = 2 and N = 3, which is therefore the sum of three elements, is given by the integrand of the two-loop sunrise diagram, (3.7) where, as before, ν = ν 1 + ν 2 + ν 3 .The parameter λ, which at this stage is taken to be integer, will acquire a linear dependence on the dimensional regularisation parameter ϵ in the case of Feynman integrals, as discussed below.
An affine form is defined to be projective if it can be identically re-written as a linear combination of the 'integrated' forms η A , defined in Eq. (3.2).Then where the homogeneous functions T B (z i ) are obtained by suitably combining the functions R A (z i ) in Eq. (3.6) with appropriate factors of z i arising from the definition of η A in Eq. (3.2).As an example, the form ψ 2 in Eq. (3.7) is clearly projective, since the differentials reconstruct the form η A for A = {1, 2, 3}, given in Eq. (3.5).Another example of a projective form that will appear in the following sections is the integrand of the one-loop massless box integral, which reads where in the concrete application one will have λ = 2ϵ and r = t/s.

A useful theorem
A useful result in what follows is the statement that the set of projective forms is closed under differentiation.In other words Theorem 1.The boundary of a projective form is itself projective.
Proof.Consider an operator p trasforming an affine q-form into a projective (and therefore also affine) (q − 1)-form, according to First, we note that the operator p is nilpotent, i.e. p 2 = 0.This can be easily shown for a single term in Eq. (3.10), R A (z i ) ω A , and the generalization is then straightforward.In fact An example can serve the purpose of illustrating the cancellation in the last step: The cancellation clearly works for any subset A ⊂ D, since there are always two terms in the sum that are proportional to z i z j , and they contribute with opposite sign.Considering for example i < j, when the factor of z j is generated by the first action of p, the sign of the term is given by the position of the indices i and j in the ordered list of the elements of A. On the other hand, when the factor z i is generated by the first action of p, what is relevant is the position of j in A − i.
The nilpotent operator p can be combined with exterior differentiation to map affine q-forms into affine q-forms.One can then show that when acting on any affine q-form ψ q .In order to prove Eq. (3.13), we note that By manipulating the indices and combining terms, the sum of Eq. (3.14) and Eq.(3.15) becomes as desired.As an example of this last step, consider which implies On the other hand, one easily verifies that as desired.Eq. (3.13), just established, is actually sufficient to conclude the proof of the theorem.In fact, it can be shown [16] that all projective q-forms can be constructed by acting with the operator p on (q + 1)-forms: in other words, they are p-exact, and any ψ q can be written as ψ q = p ξ q+1 .An example can clarify this statement.Consider a generic affine two-form where R ij are homogeneous rational functions of degree −2 in the variables z 1 , z 2 and z 3 .By imposing that p(ψ 2 ) = 0 it is immediate to obtain which is a projective form.The generalization to affine n-forms is straightforward, as the condition generates a system of linear equations that is enough to fix all the rational functions appearing in the form, but one (the overall coefficient function multiplying the projective volume form).Using this information on the l.h.s of Eq. (3.16), one sees that the first term vanishes by the nilpotency of p; the second term must then also vanish, which implies that dψ q is itself p-exact and thus projective, as desired.
In the context of Feynman parametric integration, the theorem is significant for the following reason: given that Feynman integrals in the parameter representation are integrals of projective forms on a simplex (as discussed below), applying the boundary operator d on the integrand generates relations among forms with the same properties, i.e. other Feynman integrands, or generalisations thereof.These relations take the form of linear difference equations, which in turn can be used to build closed systems of differential equations to ultimately compute the integrals, just as normally done in the momentum-space representation.

Feynman integrals as projective forms
This section presents the core results of our paper.We identify the integrands of Feynman integrals as projective forms of a specific kind, we examine their properties, and finally we use the fact that the differential of a projective form is still a projective form to write a set of generic relations among parametric integrands that include and generalise those appearing in Feynman integrals.These relations take the form of integration-by-parts (IBP) identities relating different (generalised) Feynman integrals, and can be used to build and simplify systems of differential equations in parameter space.
To begin with, consider the projective form where η n−1 is the complete projective volume form of the projective space PC n−1 , while Q({z i }) is a polynomial of degree (l + 1)P − n and D({z i }) a polynomial of degree (l + 1).We recognise that the integrand of Eq. (2.5) is a specific instance of such a form, with the polynomial D given by the second Symanzik polynomial of the graph, F, and with P = ν − ld/2.A first important property of projective forms such as Eq.(4.1), and thus in particular of Feynman integrands, is the following Proof.This follows immediately from the fact that i) α n−1 is a closed form; ii) η n−1 is null on each surface defined by z i = 0 Indeed, if we denote by ∆ the subset of C n given by the surface connecting points in the boundaries of O and O ′ that have a common image in the projective space, then O+∆−O ′ α n−1 = 0 because of statement i), while ∆ α n−1 = 0 because of statement ii).
We note that this theorem provides, in particular, a proof of the so-called Cheng-Wu theorem [78].The latter, in its original form, states that in the argument of the delta function in the Feynman-parametrised expression, Eq. (2.5) it is possible to restrict the sum to an arbitrary subset of Feynman parameters z i .In fact, consider the integration of Eq. (4.1) on the projection of the i=1 z i = t} are equivalent for any positive value of t.In the limit t → ∞ the integration domain becomes independent of z n , and becomes a semi-infinite (n−1)-dimensional surface based on the simplex {z i | n−1 i=1 z i = 1}.Figure 2 provides an example in three dimensions of the two mentioned surfaces.Given these preliminary considerations, we can proceed using the conventional choice of S n−1 as integration domain.In that case so that Sn−1 where we use the shorthand notation z for the set {z i }.Any consistent choice of the polynomials Q(z) and D(z), yielding a projective form, provides a natural generalisation of Eq. (2.5).We can now use the fact that the boundary of a projective (n − 2)-form is a projective (n − 1)-form, to construct integration-by-parts identities in Feynman parameter space in full generality.To this end, consider the projective (n − 2)-form with any suitable choice of the polynomials H i (z) and D(z).Technically, we need to restrict the choice of D(z) so that the singularities of ω n−2 lie in a general position with respect to the simplex integration domain S n−1 , and in particular they do not touch the sub-simplexes forming the boundaries of S n−1 .This case was labeled as case (A) in [16].When the singularities reach the integration domain, it is necessary to perform a blow-up of the singular points and treat the singular regions separately.Note that dimensionally regularized UV and IR divergent integrals can be treated without difficulties: these divergences are regulated by the parameter ϵ, which features in the exponents of the parameters, while the restriction on D(z) is related to external kinematics, which may force the denominator to vanish on boundary simplexes.Acting now with the boundary operator d on the form ω n−2 gives This is the sought-for integration-by-parts identity: the integration of any projective form of the kind introduced in Eq. (4.1), with the choice can be reduced to the integration of forms with smaller values of P , modulo a possible boundary term, which can be integrated via the Stokes theorem on sub-simplexes (this is the reason for the requirement that the singular surface of α n−1 should not intersect the boundary).It is important to stress that the possible presence of non-vanishing boundary terms represents a substantial difference with respect to the conventional IBP identities in momentum space: in Feynman parametrisation, boundary terms do not in general integrate to zero.In terms of Feynman diagrams, the integration over sub-simplexes represents the shrinking of a line of the diagram to a point.Eq. (4.5) will be the basis for all the applications in the following sections.We note that, when applied to Feynman integrals, Eq. (4.5) is valid for any number of loops and external legs, since the structure of the integrands in parameter space can always be written as was done in Eq. (4.1).In order to explore its applications, we begin by specialising to one-loop graphs.

One-loop parameter-based IBP
Consider Eq. (2.5) for a one-loop diagram with n internal propagators.In this case one can write where we introduced the matrix s ij (i, j = 1, . . ., n + 1), defined by as well as the auxiliary quantities The s ij represent then the invariant squared masses of the combinations of external momenta flowing in or out the diagram between line i and line j.We can now consider Eq. (4.4), and choose for H i simply the numerator of Eq. (4.7).Furthermore, we can consider separately the n forms obtained by omitting from the projective volume the variable z i , (i = 1, . . ., n), in turn.This amounts to setting for some h ∈ {1, ..., n}.Supposing ν h > 1, this leads to where s n+1,n+1 = 0.In order to clarify the structure of this expression, we introduce an index notation, following Ref.[28].We first define Then we write to denote the same function as in Eq. (4.12), where however the indices ν i ∈ {I, J , K} have been respectively decreased by one, left untouched, and increased by one.Note that we consider sets such that {I} ∪ {J } ∪ {K} = R. Furthermore, the raising and lowering operations are defined in order to preserve the character of f as a projective form, so the exponent of the denominator is re-determined after raising and lowering the indices.With this notation, Eq. (4.11) can be written as This is the desired 'integration by parts identity' at one loop, which at this stage is kept at integrand level to emphasise the fact that the boundary integral is not a priori vanishing.Notice that at the one-loop level one also has the constraint immediately following from the definition of f and the one-loop Symanzik polynomial.
Ref. [28] shows that, when the boundary term is zero, Eq. (4.14) and Eq.(4.15) allow for the systematic construction of a closed system of first-order differential equation.The proof proceeds by considering a set of parametric integrals containing all integrals obtained by raising an even number of parameter exponents by one unit (including the one of the first Symanzik polynomial).The derivatives of these integrals with respect to s ij are then included in a linear system of equations, as in Eq. (4.11) and Eq.(4.15), which is then solved in terms of the original set of integrals.This procedure is constructive and algorithmic, but one notices empirically that the number of integrals in the system is often higher than the number of actually independent master integrals, in concrete cases with specific mass assignments.In Section 5, we will use a similar construction, tryings to minimize the over-completeness of the resulting bases.

One-loop massless box
Let us consider the integral in Eq. (4.7) for the one-loop massless box integral, where n = 4, and for simplicity we set the renormalisation scale as µ 2 = (p 1 + p 4 ) 2 ≡ s, while (p 1 + p 2 ) 2 ≡ t (all momenta incoming).In particular, we focus on the simple case where all ν i = 1, and we define where, as in Eq. (3.9), we defined r = t/s, and the notation for four-point integrals is from now on I(ν 1 , ν 2 , ν 3 , ν 4 ; ν 5 ).This notation is set up so that the arguments of the function correspond directly to the exponents of the propagators in the corresponding Feynman diagrams.Notice that in this framework, as is well known, the dimension of spacetime becomes simply a parameter related to the exponents of the first Symanzik polynomial, and dimensional shift identities are naturally encoded in the parameter based IBP equation, Eq. (4.14).The matrix s ij for I box reads The construction in Ref. [28] shows that the integrals that appear in the final differential equations system are I(1, 1, 1, 1; 2ϵ) and the ones obtained from it by raising by one the powers ν i for an even number of parameters.Based on this result, we expect a basis set of integrals to be given by In this simple case, we know that this basis is overcomplete: only three linearly independent master integrals are needed for the calculation of the one-loop massless box [53].Since however our goal here is just to establish the viability of the method, we proceed with the ansatz in Eq. ( 5.3).
To verify that this is indeed a basis and that we can close the system, consider first the derivative of I(1, 1, 1, 1; 2ϵ) with respect to r, which indeed contains only integrals belonging to the desired set.On the other hand In order to proceed, it is necessary to express the integral I(3, 1, 3, 1; 2ϵ) in terms of integrals belonging to the chosen set.Eq. (4.14) for The boundary term is which follows from the fact that the projective form η 234 vanishes on all boundary sub-simplexes, except the one defined by z 1 = 0, where however the integrand is zero.This property holds for all the identities used in this section, and the reasoning will not be repeated.Consider now Eq.(4.14) for ν 1 = ν 2 = ν 3 = 2, ν 4 = 1 and h = 2, as well as the sum rule in Eq. (4.15) for I(2, 1, 2, 1; −1 + 2ϵ).One finds where the symmetry of the integrand under the exchange of (z 1 , z 3 ) with (z 2 , z 4 ) has already been taken into account.This system and Eq.(5.6) allow to find a solution for I(3, 1, 3, 1; 2ϵ), given by involving only integrals allowed in the system.Furthermore, one easily sees that the integral I(2, 2, 2, 2; 2ϵ) is also involved in the equation , 2, 2; 2ϵ) . (5.10) The last derivative to be computed in terms of the chosen set of basis integrals is ∂ r I(2, 2, 2, 2; 2ϵ), which is proportional to I(3, 2, 3, 2; 2ϵ).Using the same procedure adopted so far, it is possible to get a linear system of equation, whose solution for the desired integral is The system of differential equations for our basis set of integral is now complete, and it reads (5.12) Given Eq. ( 5.12), one can proceed using standard methods.In particular, Eq. (5.12) is not in canonical form [11]. Several techniques are available to solve this problem [54][55][56][57].Here, we simply follow the method of Magnus exponentiation [58]: the necessary steps are presented in Appendix B.
Once the system is in canonical form, it can be solved iteratively as a power series in ϵ by standard methods 4 .
In the spirit of a proof-of-concept, we have not developed a systematic approach to the search for useful boundary conditions to determine the unique relevant solution of the system.In the case at hand, continuity in r = −1, uniform-weight arguments, and the known value of the residue of the double pole in ϵ can be used to recover the known solution.We find matching the result reported, for example, in Ref. [53].The one to one correspondence between the two results is found by setting the overall constant k(ϵ) = 4 − π 2 3 ϵ 2 − 40ζ(3) 3 ϵ 3 .

One-loop massless pentagon
We now turn to the natural next step, the one-loop massless pentagon.In dimensional regularisation, it is well-known that this integral can be expressed as a sum of one-loop boxes with one external massive leg (corresponding to the contraction to a point of one of the loop propagators), up to corrections vanishing in d = 4.In this section, we will recover this result, showing that, in this case, the method connects to the derivation of the pentagon integral first reported in Ref. [59].
From the point of view of projective forms, the analysis of this case is interesting because it involves a non-vanishing boundary terms, contrary to what happened in Section 5.1 and to the analysis of Ref. [28].Consider then Eq. (4.11) for a five-parameter integral with ν 1 = ν 2 = ν 3 = ν 4 = ν 5 = 1, and the exponent of the U polynomial equal to 2ϵ.Starting with the case h = 1, we obtain the equation with Using Stokes theorem, and considering the only subset of the boundary of the five-dimensional simplex where η {2,3,4,5} ̸ = 0, the boundary term of this equation becomes box is a one-loop box integral with one massive external leg, with a squared mass proportional to s 25 .Effectively, the propagator with index ν 1 has been contracted to a point.Note that, when applying Stokes theorem, the integration over boundary domains corresponds to the proper integration region, needed to obtain the lower-point Feynman integral, up to a sign arising from the orientation of the boundary.Reversing this orientation, when needed, produces a sign that, for example, cancels the minus sign in Eq. (5.15).
Considering, in a similar way, all the possible values for h, the following system of equations is obtained 0 0 s 13 s 14 0 0 0 0 s 24 s 25 s 13 0 0 0 s 35 s 14 s 24 0 0 0 0 s 25 s 35 0 0 where the integral I(1, 1, 1, 1, 1; −1 + 2ϵ) is proportional to the pentagon integral in d = 6 − 2ϵ.The solution to this system for the pentagon integral recovering the result of Ref. [59].The correspondence between the coefficients reported here and those of Ref. [59] can be derived using the definition c i = 5 j=1 S ij in their notation.A direct consequence of Eq. (5.18) is the well-known theorem stating that the one-loop massless pentagon can be expressed as a sum of one-loop boxes with an external massive leg, up to O(ϵ) corrections.This last statement is due to the infrared and ultraviolet convergence of the 6 − 2ϵ dimensional pentagon, which implies that the last line of Eq. (5.18) is O(ϵ).

Two-loop examples
The first Symanzik polynomial for l-loop Feynman integrals, with l > 1, displays a much more varied and intricate structure compared to the one-loop case, corresponding to the factorially growing variety of graph topologies that can be constructed.Some classes of diagrams can still be described to all orders: a natural example is given by the so-called l-loop sunrise graphs, depicted in Fig. 3, contributing to two-point functions and involving (l + 1) propagators.The monodromy ring for these graphs was identified in Ref. [17], but this result was not (at the time) translated into a systematic method to construct differential equations.The simplest non-trivial graph of this kind corresponds to l = 2, and we will discuss it below, in Section 6.1, in the case in which the masses associated with the three propagators are all equal.We will then consider the other non-trivial topology contributing to two-point functions at two loops, the five-edge diagram depicted in Fig. 4.

Two-loop equal-mass sunrise integral
Sunrise graphs at l loops are characterised by the first Symanzik polynomial where ẑi is excluded from the product.Graphs of this class have generated a lot of interest in recent years.The two-loop sunrise graph with massive propagators is the simplest Feynman integral involving elliptic curves, and has been extensively studied both in the equal-mass case and with different internal masses [60][61][62][63][64][65][66][67][68][69][70]; furthermore, sunrise diagrams with massive propagators at higher loops provide early examples of integrals involving higher-dimensional varieties, notably Calabi-Yau manifolds [71][72][73][74][75].
In our present context, we would simply like to show how the projective framework that we are developing leads to the Picard-Fuchs differential equation obeyed by the (equal-mass) two-loop  [61].To this end, consider again equation Eq. (3.7), which gives the relevant integral.In our present notation where here r = p 2 m 2 , and p µ is the external momentum.For simplicity, we will work in d = 2, where the integral is finite both in the ultraviolet and in the infrared.In this case, the first Symanzik polynomial drops out, and the integrad is simply the inverse of the second graph polynomial.It is important to note that for this diagram both Symanzik polynomials vanish when approaching the boundary of the simplex S {1,2,3} , at the points z i → z j → 0 and z k → 1.In principle, this configuration invalidates the application of Stokes theorem, as discussed in Section 4, and one needs to introduce a regularisation, for example by deforming the boundaries of the simplex near the corners [61,64,69].In the equal-mass case, the domain deformation can be avoided, since the corresponding corrections cancel: we will therefore proceed with the general method, applying directly Eq. (4.5).For an explicit discussion of the differences between the two cases, see Ref. [76].
Continuing with the strategy adopted at one loop, we use the numerator of Eq. (6.2) (at this stage still for generic d) to define which gives Furthermore, denoting as before the square bracket in denominator of the integrand in Eq. (6.2) by D(r), we find Inserting Eq. (6.4) and Eq.(6.5) into Eq.(4.5), and picking the appropriate value of P to ensure projective invariance, we arrive at the IBP equations where in the second step we used the notation for raising and lowering operators in the function f as discussed in Section 4.1.The functions f are also related by the identity Using the sum rule in Eq. (6.7), and Eq. ( 6.6), we can build a linear system of equations involving the integrals I(0, 0, 0, 3ϵ), I(1, 0, 0, 1 + 3ϵ), and a non-vanishing boundary contribution B, arising from the IBP relation for I(1, 1, 0, 1 + 3ϵ) when taking h = 3 (at this point, it should be clear that boundary terms only survive when ν h = 1).The linear system is presented in Appendix D, and the boundary term contributes to the equation where Note that the minus sign in the denominator and the factor of (−1) 2ϵ come from the convention of including the masses with a minus sign in the second Symanzik polynomial.As stated above, we now set d = 2, so that the boundary term simply becomes B = 1 2 .The linear system given in Appendix D is sufficient to yield the following non-homogeneous differential equations, involving two master integrals (the third master integral appears here as the non-vanishing boundary term): We note that the differential equation system in Eq. (6.10) is the same reported in [77], up to a a different normalisation of the non-homogeneous term, which is solely due to our different normalisation of Feynman integrals.This system can be transformed into a single second-order differential equation of Picard-Fuchs type by using the OreSys package for Mathematica: the result is r 3 corresponding to the elliptic second order differential equation discussed in [60,77], up to our different normalisation.

Two-loop five-edge diagram
As a last example, we consider the two-loop, five-edge diagram represented in Fig. 4, with all internal edges taken to be massless.In this way, the only kinematic parameter is the squared momentum p 2 carried by the external legs.This diagram has been extensively studied, starting with the seminal discussion in Ref. [3].In this section, the result of [3] is re-derived by using the parameter-space method presented in this article.
The two polynomials (and the corresponding Feynman integral) are symmetric under the re-labeling a property which reflects the symmetries of the graph, and which can be used to simplify the expressions resulting from the integration by parts identities in Eq. (4.5).Indeed, it is a virtue of all approaches based on parameter space that such symmetries under permutations of the graph propagators are manifest from the beginning, and one does not have to deal with the degeneracy of possible graph parametrisations associated with different loop-momentum assignements, as is the case in momentum space.
To illustrate the use of these symmetries, consider the integral which is proportional to the Feynman integral associated with Fig. 4. Eq. (6.12) implies that and the use of the symmetry properties of the graph polynomials reduces this equation to the much simpler form Eq. (6.16) is the first step necessary for reducing integral in Eq. (6.14) to a linear combination of simpler integrals.Consider now the integration by parts identity in Eq. (3.4), with h = 1, and with (6.17) One obtains then the integration by part identity Upon integration over the simplex S 5 , this yields The integral Ω 1 can be calculated by means of Stokes' theorem, with the result where the sign in the definition of ω 1 , Eq. (6.17), is absorbed by the boundary ∂S 5 = −S {2,3,4,5} , since the integrand vanishes on the other sub-simplexes comprising ∂S 5 .The boundary term Ω 1 is proportional to the Feynman integral obtained from the diagram in Fig. 4 when the edge labelled 1 shrinks to a point, and the propagator corresponding to edge 3 is raised to the power of 2. This integral can be evaluated straightforwardly, yielding a product of Gamma functions (see for example Ref. [3]).
A similar strategy can be applied to find the three other equations that are necessary to reduce the two-loop five-edge integral to simpler integrals.The resulting equations are integral, as we did, but it would be interesting to explore variations on this theme with an eye to optimisation.On the other hand, in contrast to momentum-pace algorithms, we observe that the parameter-space approach bypasses the ambiguity due to the choice of loop-momentum routing, which can be non-negligible for complicated diagrams; similarly, the issue of irreducible numerators is implicitly dealt with at the momentum integration stage.These two aspects are among the consequences of the fact that parameter space offers a minimal representation of Feynman integrals, transparently related to the symmetries of the original Feynman graph.An especially promising aspect of the projective framework is its close connection to the most significant algebraic structures associated with Feynman integrals.The Barucchi-Ponzano analysis can indeed be seen as an application of the results of Ref. [18], and it is notable that it succeeds not only in constructing a system of differential equations for n-point one-loop integrals, but also in setting a bound on the size of the system, guaranteeing its closure, and providing an algorithmic construction.This is to be contrasted with the very large size of the systems of IBP identities that emerge in the intermediate stages of calculations in standard algorithms.It is clearly a goal of future research to extend these techniques and the analysis of Regge and collaborators to more complicated two-and higher-point integrals.In particular, studies on three-loop two-point functions and on two-loop three-point functions are currently ongoing, and steps towards the automation of the generation of IBPs in the projective framework are under way, with the goal of reaching stateof-the-art topologies such as two-loop penta-and hexa-boxes and three-loop four-point functions.When complex multi-scale examples of this kind become available, a more thorough comparison of the two approaches, including computational aspects, will become possible.After this last step, the system is finally reduced to its canonical form, which can be solved iteratively.We write These are seven equations involving nine independent integrals, two of which are the chosen basis integrals.
The system is of course easily solved with elementary methods.

Figure 1 :
Figure 1: a) One-loop box diagram b) Two-loop sunrise diagram

Figure 2 :
Figure 2: Example of two domains in R 3 that are equivalent under projective transformations.

5 ) 0 Db
which reduces the system to ∂ r b ′ (r) = The diagonal part D(r) of the ϵ-independent term is removed by the matrix C d (r) = exp r ′′ (r) .(B.7)One may now directly apply Magnus' theorem, with the final change of basis given by C red (r) =