Feynman Integrals and Intersection Theory

We introduce the tools of intersection theory to the study of Feynman integrals, which allows for a new way of projecting integrals onto a basis. In order to illustrate this technique, we consider the Baikov representation of maximal cuts in arbitrary space-time dimension. We introduce a minimal basis of differential forms with logarithmic singularities on the boundaries of the corresponding integration cycles. We give an algorithm for computing a basis decomposition of an arbitrary maximal cut using so-called intersection numbers and describe two alternative ways of computing them. Furthermore, we show how to obtain Pfaffian systems of differential equations for the basis integrals using the same technique. All the steps are illustrated on the example of a two-loop non-planar triangle diagram with a massive loop.


I. INTRODUCTION
The opportunity for revealing weak deviations from the Standard Model predictions within particle collision experiments at high-accuracy, or studying the properties of newly discovered particles, as well as the possibility of understanding formal properties of quantum theories not directly deducible from the basic structure of Lagrangians, or disclosing deeper similarities among theories which apparently describe the interactions among particles of different species, but also, more generally, the investigation of physical cases admitting a field-theoretic perturbative approach might depend on our capability of computing scattering amplitudes at higher order, therefore of evaluating Feynman integrals.
The exploitation of existing relations between multiloop integrals is of fundamental importance to minimize the computational load required for the evaluation of scattering amplitudes that, according to the number of involved particles and to their masses, depend on several kinematic scales. Feynman integrals are known to obey contiguity relations dubbed integration-by-parts (IBP) identities [1], which play a crucial role in the evaluation of scattering amplitudes beyond the tree-level approximation. IBP identities yield the identification of an elementary set of integrals, the so-called master integrals (MIs), which can be used as a basis for the decomposition of multi-loop amplitudes. At the same time, IBP relations can be as well exploited to derive differential equations [2][3][4][5][6][7][8][9], finite difference equations [10,11], and dimensional recurrence relations [12,13] obeyed by MIs. The solutions of those equations are valuable methods for the evaluation of MIs, should the direct integration be prohibitive, as it is usually the case in practical applications.
IBP identities can be generated by considering integrals of total derivatives that vanish on the integration boundary, which results in a system of linear relations between Feynman integrals that differ by the powers of denominators and/or scalar products in the numerator, see, e.g., [14][15][16]. This is an algebraic operation that can be easily automatized, without needing the explicit form of the integration domain, by simply imposing the qualitative information about null surface terms. Nevertheless, finding a basis of MIs often requires solving a huge linear system of equations [10]. For processes in which multiple kinematic scales are involved, it may represent an insurmountable task. We have been experiencing a vivid progress in performing such computations, by the refinement of the system solving strategy [10], either due to novel algorithms [17][18][19][20][21] or to the developing improved software [22][23][24][25][26], ending up with calculations of multiparticle amplitudes which were considered inaccessible a few years ago [27][28][29].
At the same time, one may want to seek alternative methods for obtaining the decomposition in terms of MIs, eventually searching for mathematical methods that allow for a direct integral reduction, which bypass the need of solving a system of linear equations.
In this work, we elaborate on a new method for establishing relations among Feynman integrals in arbitrary space-time dimensions, and for projecting them onto a basis.
An archetype of such a basis reduction is the Gauss' contiguous relation, e.g., It can be regarded as a basis reduction of the hypergeometric function 2 F 1 (a, b; c+1; z) on the left-hand side in terms of two MIs on the right-hand side. Of course, such an identity can be derived by considering linear relations between 2 F 1 's with different arguments, however modern approach to this problem is more direct [30]. One starts by writing down the integral representation: Here the integration domain C := − −− → (0, 1) together with the information about the branch of the integrand is called the twisted cycle, while the single-valued differential form ϕ := c c−b dlog x is called the twisted cocycle. 1 The above integral is understood as a pairing of these two objects [30]. B(a, b) is the Euler beta function. Similarly, we can consider two other logarithmic forms: which upon integration give rise to 2 F 1 (a, b; c; z) and 2 F 1 (a+1, b; c+1; z) respectively. The problem reduces to projecting ϕ onto a basis of ϕ 1 and ϕ 2 . This can be achieved by computing certain topological invariants called intersection numbers of pairs of cocycles, which are rational functions in the coefficients a, b, c. They become the coefficients of the basis expansion on the right-hand side of the Gauss' contiguous relation. Inspired by this approach, we wish to apply the same techniques of intersection theory to the study of Feynman integrals.
Among different representations of Feynman integrals, the one most closely resembling Aomoto-Gel'fand hypergeometric functions is the so-called Baikov representation [31]. It uses independent scalar products between external and internal momenta as the integration variables. This change of variables introduces a Jacobian equal to the Gram determinant involving both types of momenta, referred to as the Baikov polynomial B, raised to a power D−D * 2 for an integer D * . This selects a critical space-time dimension D * , in which the integral possesses nice mathematical properties. We review the Baikov representation in Section II and give more details of its derivation in Appendix A.
The Baikov polynomial fully characterizes the space on which the integrals are defined. Early signs of this fact were found by Lee and Pomeransky, who applied Morse theory to relate the number of MIs to counting of the critical points of B D−D * 2 [32]. 2 In fact, if one is interested in MIs in the top sector, it is enough to consider maximal cuts of Feynman integrals, which retain a lot of properties of the original integrals. For recent studies of maximal cuts in the Baikov representation, see [19,[36][37][38][39][40][41]. In this work we focus on maximal cuts in order to illustrate the basis reduction algorithm in the simplest possible setting.
After determining the size of the basis, which we denote by |χ|, we construct bases of twisted cycles C i and cocycles ϕ j , whose pairings give rise to |χ| 2 integrals: They form a minimal, linearly-independent, basis in which any other integral of the same type can be expanded.
Here B(z) is obtained by evaluating B on the maximal cut surface. Twisted cycles C i are chosen as certain regions with boundaries on {B(z) = 0}, while twisted cocycles ϕ j are differential forms with logarithmic singularities along all the boundaries of the corresponding C j . We describe them in Section III. Note that our basis integrals are not the conventional MIs, as they do not necessarily have an interpretation as valid maximal cuts themselves. They are more granular objects.
In Section IV we apply the tools of intersection theory of the appropriate homology and cohomology groups [42,43] to the problem of basis reduction. It can be done separately in the space of twisted cycles and cocycles. We focus on the reduction in the space of twisted cocycles and show how to apply two different techniques of evaluating intersection numbers: the special case of logarithmic forms has been treated by one of us in [44], while the case of general one-forms was discussed by Cho and Matsumoto in [43]. To further support the connection between Feynman calculus and intersection theory, we provide references to the relevant literature for the interested reader.
It is known that maximal cut of MIs obey homogeneous difference and differential equations [45][46][47][48]. From several studies focusing on the solution of differential equations for Feynman integrals [38,40,41,[46][47][48][49][50][51], the MIs have been identified with the independent components of the integration domain. Owing to the complete characterization of the integrand and of the integration domain, explicit solutions for the maximal cuts can be found in Baikov representation. Starting from a set of MIs whose system of equations is linear in the space-time dimension [52], the solutions of the homogeneous equations can be used to write transformation matrices (namely the resolvent of the homogeneous system around specific critical dimensions) [46,47], for building special set of MIs that obey canonical systems of differential equations [8,53], where the dependence on the space-time dimension is factored out of the kinematics. In Section V we discuss Pfaffian systems of differential equations satisfied by the basis integrals, and show that they can be derived using our basis reduction algorithm.
Even though the techniques described here are generally applicable, throughout the paper we consider an example of the two-loop non-planar triangle diagram with a massive loop for illustration purposes. In Section VI we study this integral directly in D=4 using intersection theory methods.

II. BAIKOV REPRESENTATION
Consider scalar Feynman integrals with L loops, E+1 external momenta, and N propagators in a generic dimension D: ( We focus on Euclidean space in all-plus signature for simplicity of discussion. Baikov considered a change of integration variables into all independent scalar products between loop and external momenta [31], (here we assume that D ≥ E, so that external momenta do not satisfy additional relations [54,55]). There are M := LE + 1 2 L(L + 1) such kinematic invariants, i · q j . In order to perform the change of variables, one needs to introduce M −N extra inverse denominators D a known as irreducible scalar products (ISPs) with exponents ν a ≤ 0. The original integral (3) is recovered when ν N +1 = · · · = ν M = 0. After the dust settles, one finds [13,56]: with a constant Jacobian c, which we will drop from now on. The integrand involves a rescaled Baikov polynomial : which is a ratio of two Gram determinants. Recall that a Gram matrix of Lorentz vectors {q i } is defined as The integration domain Γ is given by imposing L conditions: so that Γ := Γ 1 ∩Γ 2 ∩· · ·∩Γ L . Notice that this implies B > 0 everywhere on the integration domain. For convenience, we review a derivation of the Baikov representation in Appendix A. The representation (5) is particularly friendly towards computing cuts. By a linear transformation we can make a further change of variables into z a := D a for all M inverse propagators. A single cut corresponds to taking a circular integration contour {|z a | = ε}, which sets D a on-shell. Repeating this procedure N times we obtain a maximal cut, which takes the general form [36]: where we ignored an overall constant Jacobian. The integrand depends on the ISPs Baikov polynomial on the maximal cut B(z) is given by setting z 1 = · · · = z N = 0 in (6). Similarly, ϕ(z) is a differential (M −N )-form obtained as a result of the residue computation. For example, if the original integral had all propagators undoubled, we have ϕ(z) = M a=N +1 dz a /z νa a . Finally, the integration domain C is an intersection of Γ with the cut surface {z 1 = · · · = z N = 0}. Whenever this intersection is empty, the maximal cut vanishes and the diagram is reducible.
For one-loop diagrams the maximal cut fully localizes the integral and hence the size of the basis is one. Let us consider an interesting example of a two-loop non-planar diagram with internal mass m and p 2 1 = s, p 2 2 = p 2 3 = 0: We also choose an ISP D 7 = 2(p 2 + 1 ) 2 − p 2 1 , for later convenience. In this case, E = L = 2, M −N = 1 and the maximal cut in Baikov representation (8) becomes: where we relabelled z 7 → z, ν 7 → ν for clarity. The rescaled Baikov polynomial reads: where ρ := s(s + 16m 2 ). In the kinematic regime s, m 2 > 0 its roots are ordered as −ρ < −s < s < ρ. The constraints (7) imply that the integration region is: Here − −− → (a, b) denotes an oriented interval between a and b.

III. MINIMAL BASIS FOR MAXIMAL CUTS
Integrals of the type (8) admit a beautiful interpretation in terms of Aomoto-Gel'fand hypergeometric functions [57,58], where they are understood as pairings of twisted cycles C and cocycles ϕ. In order to see this, let us consider the following one-form: It defines a flat connection ∇ ω := d + ω∧. Its significance is that for any (M −N −1)-form ξ and cycle C. This means that we can define equivalence classes ϕ| of (M −N )-forms ϕ up to terms ∇ ω ξ that integrate to zero: Similarly, we have equivalence classes |C] of cycles C up to terms integrating to zero. 3 The two classes ϕ| and |C] encode all IBP identities and contour deformations. Their pairing, denoted by ϕ|C], is defined to be equal to the maximal cut (8).
In general, Baikov polynomial on the maximal cut B(z) admits a decomposition into irreducible components: From now on we will consider only the cases in which each b i (z) has degree at most M −N , so that the corresponding variety {b i (z) = 0} does not have self-intersections. We also assume that the dimension D in the exponent of the Baikov polynomial is generic, and in particular not an integer.
is a Morse function whenever its critical points are non-degenerate. Then the number of independent twisted cycles and cocycles, |χ|, is given by |χ| = where C γ is the number of critical points with a Morse index γ [59]. Under the assumptions given in [30], all critical points have the maximal index γ = M −N and hence: Here we used that ω = 0 determined the critical points. 4 When all irreducible components b i (z) are linear, |χ| is also equal to the number of bounded regions in R M −N \ {B(z) = 0} [60]. Let us discuss how to construct bases of twisted cycles C i and cocycles ϕ i for i = 1, 2, . . . , |χ|. The real section of the integration domain, R M −N \ {B(z) = 0}, decomposes into multiple disconnected regions, called chambers. Each chamber is a valid choice of a basis element C i . For each C i we can construct the corresponding twisted cocycle ϕ i 3 To be precise, on the space where Ω k ( * D) is the sheaf of smooth holomorphic k-forms on X with poles along the singular divisor D of B(z). Locally finite twisted homology groups H lf k (X, Lω) |C] with coefficients in a rank-1 local system Lω, are isomorphic by H k (X, ∇ω) Hom C (H lf k (X, Lω), C), which induces the pairing ϕ|C] := (8) [30]. 4 The idea of determining the size of the basis using the Euler characteristic χ was first considered in Feynman integral literature by Lee and Pomeransky [32], see also [35]. In the present work, however, we do not make any claims beyond maximal cuts in generic dimension D.
with logarithmic singularities along the boundaries of C i . For the algorithmic constructions, see, e.g., [30,[61][62][63][64] and the more recent study [65]. It is currently known how to do it in the cases when C i is bounded by hyperplanes and at most one hypersurface {b i (z) = 0} with degree up to M −N [65]. This is enough for considering many interesting examples of maximal cuts [40]. For instance, whenever a given chamber C i is simplexlike, i.e., bounded by exactly M −N +1 hyperplanes {b j (z) = 0}, j = 1, 2, . . . , M −N +1, we have: In situations when {B(z) = 0} is non-normally crossing, i.e., more than M −N hypersurfaces intersect at a single point, one needs to consider a blowup in the neighbourhood of this point. The choice of logarithmic basis comes with multiple advantages. For instance, its singularity structure in the variable ε := (D−E−L−1)/2 is manifest: neighbourhood of each co-dimension k ≤ M −N boundary of C j gives contributions at order ε −k , e.g., hypersurfaces contribute at order ε −1 , while their intersections at order ε −2 : Hence the leading divergence comes from the highest codimension boundaries (points), around which the integral behaves as ε −(M −N ) . To be precise, a given co-dimension k boundary of the form contributes to a given integral ϕ i |C j ] at order ε −k if and only if it belongs to the boundary of a given twisted cycle ∂C j and at the same time the twisted cocycle ϕ i has a non-vanishing residue along (20). At this stage let us remark that it is always possible to shift D → D + 2n for n ∈ Z in the exponent of the Baikov polynomial [12], at a cost at redefining ϕ → B(z) −n ϕ. (Requirement of single-valuedness of ϕ imposes n ∈ Z.) This changes ε → ε+n and hence the singularity structure, but does not affect the choice of the basis itself.
We can organize all possible pairings of twisted cycles and cocycles into the twisted period matrix P with entries: Its |χ| 2 components provide a minimal linearlyindependent basis for any integral of the form (8), as was first observed in [47], and later also [38][39][40] in Baikov representation. Let us stress that elements of the basis might not be necessarily recognized as Feynman integrals on their own.
For example, evaluating (14) for the diagram (10) we have: We find |χ| = 3 solutions of the condition ω = 0: Their positions are irrelevant for the counting, but they will be used later on for other purposes. (Size of the basis is alternatively determined to be 3 from the dimension of the pure braid group associated to (21).) A natural choice for the basis of twisted cycles is: The corresponding twisted cocycles are: which are designed to have residues ±1 on the two endpoints of each C i . They are special cases of (18). The matrix P consists of nine linearly independent Appell functions F 1 that form the basis. Note a symmetry for the two independent kinematic invariants, s → −s, ρ → −ρ, under which while B(z) is invariant. Therefore the entries of the matrix P are related as P ij → P 4−i,4−j and we end up with a five-dimensional functionally-independent basis.
In the Feynman integral literature it is customary to talk about bases of MIs, which have to consist of Feynman integrals, and hence have a fixed integration domain C and different powers of ISPs. In this language, the diagram (10) in generic dimension D has 3 linearly independent MIs and 2 after imposing the above symmetry, in agreement with [47,50].

IV. BASIS REDUCTION WITH INTERSECTION NUMBERS
The goal of basis reduction is expressing an arbitrary integral of the form (8) in terms of the |χ| 2 basis functions in P. This can be done separately in the space of cycles and cocycles. In order to do so, we introduce the notion of a metric on these spaces. Assuming existence of dual spaces |ϕ and [C|, let us consider pairings between their basis elements: These pairings are called intersection numbers. Using simple linear algebra, we can decompose an arbitrary twisted cocycle ϕ| into a basis of ϕ i | as follows: and similarly for a twisted cycle |C] in a basis of |C l ]: Here ϕ|ϕ j C −1 ij and H −1 kl [C k |C] are coefficients of the expansions. Putting these two decompositions together, we find that the original integral ϕ|C] is expressed in terms of basis functions in P as follows: In fact, this statement is completely general and holds for any Feynman integral in arbitrary parametrization, as long as one is able to identify |ϕ and [C| and their pairings. The advantage of Baikov representation of maximal cuts is that such identifications can be made, which allows for explicit computations.
To be precise, we define the dual space as equivalence classes |ϕ : ϕ ∼ ϕ + ∇ −ω ξ with the connection ∇ −ω and similarly for dual twisted cycles [C|. With this choice, intersection numbers [C i |C j ] are trigonometric functions of the dimension D [42]. They can be computed straightforwardly by considering all the places where C i and C j intersect geometrically. (Additional care needs to be taken when boundaries of C i and C j are non-normally crossing.) There exist numerous ways of evaluating them: we refer the reader to, e.g., [30,42,[66][67][68][69][70][71][72][73][74][75][76][77][78][79][80]. In the example at hand, the original Baikov integration domain C from (13) already decomposes as: and hence no detailed computation is necessary. Examples for other maximal cuts will appear elsewhere.

A. Intersection Numbers of Logarithmic Forms
Similarly, intersection numbers ϕ i |ϕ j can be evaluated in multiple different ways, see, e.g., [30,43,44,67,[69][70][71][72][73][74][75][76][77][78][80][81][82]. They are rational functions in kinematic invariants and the dimension D. It was recently found that for logarithmic forms ϕ i and ϕ j there exists a formula localizing on the critical points given by ω = 0 [44]: Here ω a are components of ω = M a=N +1 ω a dz a , and ϕ denotes a differential-stripped cocycle ϕ =: ϕ M a=N +1 dz i . Let use apply it to the two-loop example (10). For simplicity, we are going to choose the same representatives (25) for cocycle bases of both ϕ i | and |ϕ j . The above formula (31) becomes: where the sum goes over the three critical points (23) and we have a Jacobian ∂ ω/∂z coming from evaluating the delta function. Performing this computation for every combination of ϕ i | and |ϕ j gives us the matrix C from (26): It is always possible to choose the dual basis |ϕ j to be orthonormal, i.e., such that C = 1, which simplifies the decomposition (27). 5

B. Intersection Numbers of Non-Logarithmic One-Forms
In order to complete basis reduction (27), we ought to compute the remaining intersection numbers ϕ|ϕ j . Let us consider the maximal cut (11) with no numerators, ν = 0, for which we have ϕ = dz. This form has a double pole at infinity, which means that we cannot use (31). In this case we employ an alternative formula for general oneforms due to Cho and Matsumoto [43] (see also [81,83]): For a review of the derivation of (34), see, e.g., Appendix A of [44]. It involves a sum over all boundary components {B(z) = 0} = {−ρ, −s, s, ρ, ∞}. For each of them, one needs to compute ψ p := ∇ −1 ω ϕ p , i.e., find a unique holomorphic function ψ p solving Since the connection ∇ ω decreases the order of the pole by one, it is enough to consider an ansatz 5 For instance, an orthonormal basis to (25) is given by: though we will not make use of it in the text.
where ord p ϕ denotes the order of the zero of ϕ around p, e.g., ord p dz/(z − p) = −1. (Analogous expansion in 1/z is done when p = ∞.) Plugging it into (35) one can solve for the first few coefficients ψ (α) p order-by-order in (z−p). Notice that ord p ψ p = ord p ϕ + 1. The residues in (34) are non-zero only if ψ p ϕ j has at least a simple pole, or in other words: In our case, basis elements ϕ j have at most simple poles around each p, while ϕ = dz has only a double pole at infinity. Using the condition (37) we conclude that p = ∞ is the only contributing point in (34). Hence we compute: using the procedure above, but expanding around infinity. Plugging the result into the formula for intersection numbers (34) we find: 6 This completes a decomposition of ϕ| into a basis, which after evaluating (27) reads: Finally, using the cycle decomposition (30), we find: ρ P 11 +P 21 +P 31 + s P 21 (40) Recall that the terms P •3 = ϕ • |C 3 ] are related to P •1 = ϕ • |C 1 ] by the symmetry s → −s, ρ → −ρ, so in practice only three basis integrals need to be computed. The factor of ε = (D−5)/2 in front guarantees that the integral ϕ|C] is finite when ε → 0. Notice that the formula (34) is specific to one-forms, and in order to use it in other examples, one needs to generalize it to cases with more ISPs. This can be achieved using the algorithm of Matsumoto [81] and will be discussed in a future publication.

V. PFAFFIAN SYSTEMS
Derivation of Pfaffian systems of differential equations follows the same algorithm. Let us denote with d a differential on the appropriately chosen space of kinematic invariants K. Acting with it on basis elements, we have: We want to bring it into the Pfaffian form: 7 Here Ω is a one-form on K. In order to find it, it is enough to project each Φ i on the right-rand side of (41) in the basis P back again using (27). Hence the entries of Ω are computed with intersection numbers: We expect that for the logarithmic bases described in this work the matrix Ω should be proportional to ε := (D−E−L−1)/2. (Recall that whenever E+L is odd, we can dimension shift the system into ε = (D−4)/2.) In addition, it can be shown that the leading order in ε of P coincides with the intersection numbers, i.e., if the bases of cocycles are logarithmic. Let us derive differential equations in the example at hand. We choose K = {x ∈ C | x = ρ/s}, which gives: where ε = (D−5)/2. Using the definition (34), we evaluate: Notice that singularities in x can only occur when two punctures out of {−sx, −s, s, sx, ∞} collide, i.e., where x is −1, 0, 1 or ∞. The matrix C was already computed in (33), which after plugging in (43) gives the one-form Ω: with , Note the symmetry Ω ij = Ω 4−i,4−j , which descents from symmetries of P. This is a Fuchsian system, which can be solved using standard techniques [8,9,52,84,85].
The linear system (42) can be converted into a higherorder differential equations for each of the i's separately. For fixed i, the solutions of such an equation can be expressed in terms of the basis of P ij for different j's, i.e., different solutions correspond to distinct choices of the integration cycles.

VI. FOUR DIMENSIONS
So far we have been working in a generic dimension D / ∈ Z. In this section we illustrate how to apply the same techniques directly in the strict D → 4 limit for the maximal cut of the diagram (10). Setting D=4 we have: which behaves as z −2 when z → ∞ (as opposed to z 2(D−5) in the generic case). Hence the integral has a trivial monodromy at infinity, which changes its topological properties. We will see that this means the size of the basis |χ| drops from 3 to 2.
In order to be able to consistently use the techniques of the previous sections and make the connection to elliptic curves known from the literature [47,50,86], let us redefine This leaves the combination B(z) −1/2 ϕ = B(z) −1/2 ϕ invariant, but makes the integral defined on with the point at infinity included (hence we require that ϕ does not have poles at z = ∞). Indeed, the doublecover of X is an elliptic curve branching at the four points {−ρ, −s, s, ρ}. For simplicity, we continue to work directly on the base space X.
That the size of the basis is 2 can be counted by solving ω := − 1 2 d log B(z) = 0, which has two solutions, z * = 0, ∞. Alternatively, it is easily seen that |C 1 ] and |C 2 ] from (24) define two independent cycles on the elliptic curve, while |C 3 ] is homologous to |C 1 ]. Finally, let us see the same fact from intersection numbers, by choosing the bases of twisted cocycles ϕ i | and |ϕ j as in (25). (This might not be the optimal choice for D=4, but we use it for consistency.) Using the definition (34) with ω, but summing over p ∈ {−ρ, −s, s, ρ} we find: where The rank of this matrix is 2, which signals that there are only two linearly-independent twisted cocycles. From here it is also seen that ϕ 1 | and ϕ 3 | are linearly dependent, and hence we choose the basis to consist of ϕ 1 |, ϕ 2 | and similarly for the dual cocycles. This means we can no longer exploit the symmetry properties under s → −s, ρ → −ρ.
As in the previous sections, we decompose ϕ| = dz/(z 2 −ρ 2 ) into the basis by computing (27): where we used: and we used the corresponding 2×2 minor of C from (53). Alternatively, we could have computed that ϕ 3 | = ϕ 1 | and obtained the same result (55) from (39). Let us derive differential equations for this basis. The forms Φ i from (41) are related to those in (45) by Φ i = Φ i /(z 2 −ρ 2 ) and ε = −1/2. Hence we have: with Plugging into the expression for Ω in (43) (upon replacing ϕ → ϕ), we find: This provides a linear system of two coupled differential equations, d P i• = Ω ik ∧ P k• . We can solve them to obtain second-order decoupled equations for P 1• and P 2• separately: In each case, basis for the two solutions is provided by different choices of twisted cycles.

VII. DISCUSSION
The mathematical structure of scattering amplitudes is richer than what is currently known. The study of geometric aspects related to the decomposition of multi-loop amplitudes in terms of a basis of independent integrals, and the recent development of ideas for the evaluation of the latter, seem to offer a new perspective on Feynman calculus.
In this work, we introduced the tools of intersection theory, borrowed from the theory of Aomoto-Gel'fand hypergeometric functions, to Feynman integrals. We identified a correspondence between the forms associated to generalized hypergeometric functions and Feynman integrals in Baikov representation, and demonstrated the applicability of intersection theory, by discussing choices of bases of twisted cycles and cocycles, basis reduction in both spaces, as well as derivation of the differential equations. We applied it to the special case of a two-loop non-planar three-point function with internal massive propagators, both in arbitrary D dimensions and in the four dimensional case, which, as studied in the recent literature, it is known to involve elliptic integrals. In particular, we analyzed maximal cuts of Feynman integrals, which have the property that the twisted cocycles have singularities only along the vanishing of Baikov polynomials, that defines the cut surface. For non-maximal cuts, by definition, there exist singularities of twisted cocycles that do not coincide with the cut surfaces. In these cases, one needs to consider relative twisted homology and cohomology groups. Definition of intersection numbers in such cases was discussed in the mathematics literature only very recently [87]. It is expected that it can also help with lifting the assumption of genericity of D, so that one can study integer dimensions directly in more general cases. We leave these questions for future research.
Even in the realm of maximal cuts, the variety defined by the vanishing of the Baikov polynomial, could become more involved at higher loops or with more kinematic scales, which complicates the determination of relevant bases of twisted cocycles. It would be interesting to study whether there exist consistent changes of variables that simplify the decomposition (16) into linear factors, perhaps at the cost of introducing different exponents, as is the case for, e.g., the Appell function F 4 [30].
Our analysis can be considered as a preliminary exploration of the possibility of developing a method for the direct decomposition of Feynman amplitudes in terms of a basis of independent integrals by means of projections, where the elements of intersection theory allow to define a metric within the space of Feynman integrals. In this appendix we review the derivation of the Baikov representation [13,31,56] for a scalar Feynman integral with L loops, E+1 external momenta, and N propagators:

ACKNOWLEDGMENTS
(A1) Out of the independent momenta (for simplicity we assume that D ≥ E, though it is not necessary) involved in the scattering process, we can construct the Gram matrix G {qi} := [q i · q j ], which has M := LE + 1 2 L(L + 1) independent entries depending on loop momenta, i · q j . We use them to express each of the N inverse propagators as: where the sum goes over M kinematic invariants i · q j and m a denotes the mass of an internal particle. Here A defines an N × M matrix, whose rows are labelled by propagators and columns by Lorentz products i · q j . For later convenience, we extend A into an M × M invertible matrix by introducing additional propagators D N +1 , D N +2 , . . . , D M called irreducible scalar products. We choose their powers to be ν a ≤ 0, so that they only appear in numerators. Baikov representation makes manifest the propagator structure of a given Feynman integral. In order to do so, we first change the integration variables from µ i to i · q j . This is done by decomposing each loop momentum i = i + ⊥ i , into the (E+L−i)-dimensional space spanned by { 1 , 2 , . . . , i−1 , p 1 , p 2 , . . . , p E } and the orthogonal complement [56]. (The order in which loop momenta are labelled does not matter for the final result.) The integration measure becomes: is a constant, the integration domain is Γ := Γ 1 ∩ Γ 2 ∩ · · · ∩ Γ L , and the rescaled Baikov polynomial reads B := det G { 1,..., L ,p1,...,p E } det G {p1,...,p E } , The constraints on Γ imply that B > 0 everywhere on the integration domain.
Since the inverse propagators D a are related to i ·q j by a linear transformation (A3), we can make an additional change of variables to express the Feynman integrals in terms of M variables z a := D a : where c = (det A) −1 is a constant Jacobian. We set c c = 1 for convenience. B and Γ are the same quantities as before, but expressed in terms of the new variables z a . The expression (A9) is called the Baikov representation.