Decomposition of Feynman Integrals by Multivariate Intersection Numbers

We present a detailed description of the recent idea for a direct decomposition of Feynman integrals onto a basis of master integrals by projections, as well as a direct derivation of the differential equations satisfied by the master integrals, employing multivariate intersection numbers. We discuss a recursive algorithm for the computation of multivariate intersection numbers and provide three different approaches for a direct decomposition of Feynman integrals, which we dub the straight decomposition, the bottom-up decomposition, and the top-down decomposition. These algorithms exploit the unitarity structure of Feynman integrals by computing intersection numbers supported on cuts, in various orders, thus showing the synthesis of the intersection-theory concepts with unitarity-based methods and integrand decomposition. We perform explicit computations to exemplify all of these approaches applied to Feynman integrals, paving a way towards potential applications to generic multi-loop integrals.


Introduction
Feynman integrals in dimensional regularization admit parametric integral representations which expose their nature as Aomoto-Gel'fand integrals, thereby enabling a novel form of investigation of their algebraic structure by means of intersection theory of twisted de Rham (co)homology for general hypergeometric functions [1][2][3]. Accordingly, intersection numbers of differential forms [4] can be employed to define a scalar product on a vector space of Feynman integrals [1], such that projecting any multi-loop integral onto a basis of master integrals (MIs) becomes conceptually identical to decomposing a generic vector into a basis of a vector space. Univariate intersection numbers, as shown in the the original studies [1,2], were sufficient to validate a novel method based on intersection theory for deriving integral relations, which was used for the direct derivation of contiguity relations for Lauricella F D functions, as well as for Feynman integrals on maximal cuts, i.e. with on-shell internal lines, that admit a one-fold integral representations. As proposed in [2], applications of this novel method to the decomposition of full Feynman integrals in terms of a complete set of MIs, including the ones corresponding to subdiagrams, as well as deriving contiguity relations for special functions admitting multi-fold integral representation, required the use of multivariate intersection numbers [5][6][7][8][9][10][11][12][13].
A recursive algorithm for computing multivariate intersection numbers was proposed in [14] and later refined and applied to a few paradigmatic cases of Feynman integral decomposition [3]. This recursive algorithm was developed in order to compute intersection numbers for twisted cohomologies associated to n-forms, which in the general case may contain poles that are not necessarily simple. In the case of logarithmic (dlog) differential forms, owing to the presence of simple poles only, the computation of the intersection numbers is known to be simpler [6,12].
Recent complementary work [15,16] shows that intersection numbers play a fundamental role in the definition of a diagrammatic coaction for MIs, which combined with the master integral decomposition studied in this paper, as well as in Refs. [1][2][3] paves a way towards comprehensive computations of scattering amplitudes using the tools of intersection theory.
The intersection theory-based decomposition has also been recently applied to the study of Feynman integrals in d = 4 ± 2 space-time dimensions, from which an unexpected relation between the behaviors around → 0 and → ∞ emerged [17] and was used to investigate the properties of canonical systems of differential equations [18]. A further interesting step for the construction of canonical integrals with intersection theory has been reported in [19]. Moreover, it was observed that using recursion relations for computing intersection numbers can be further refined by relating them to dlog forms at each step of the recursive algorithm [20]. Other recent intersection-theory approaches include [21][22][23].
In this work, we elaborate on the vector space structure of Feynman integrals, whose complete framework was presented in [3], by providing the mathematical details that brought us to its formulation. We show how intersection numbers can be used to establish linear and quadratic relations for Feynman integrals, and, more generally, for Aomoto-Gel'fand generalized hypergeometric functions.
In particular, we focus on different ways of using intersection theory in order to derive linear relations for Feynman integrals, as well as the systems of differential equations and the finite difference equations they obey. Moreover, we present here for the first time, a novel algorithm for Feynman integral decomposition, which we will refer to as top-down decomposition, showing that the coefficients of MIs can be suitably extracted by projections via intersection numbers within an iterative strategy, starting from the integrals that correspond to graphs with the highest number of internal lines, and ending with those corresponding to graphs with the lowest possible number of internal lines (given, in the general case, by the product of as many tadpoles as the number of loops). Following [2,3], we also make use of two other algorithms: the bottom-up decomposition and straight decomposition to similar aim. All these strategies combine the advantages of the integrand decomposition techniques [24][25][26][27][28][29][30][31], the unitarity-based methods [32][33][34][35][36][37][38][39][40][41][42][43][44][45][46], and the intersections theory-based decomposition.
As originally defined [4], intersection theory for twisted cohomologies and the evaluation of multivariate intersection numbers are applicable to differential forms obeying certain genericity conditions, whose purpose is to regulate boundaries of integration and ensure that they integrate to analytic functions. In the physics language, this corresponds to the analytic regularization of Feynman integrals [47]. To simplify computations, we employ this regularization whenever necessary. It has the additional benefit of resolving the ambiguities that arise when there is a non-trivial overlap between critical points and singularities [2,20,48]. Recent mathematical developments, employing the notion of intersection numbers for the relative twisted cohomology [49], seem to offer the possibility of studying the vector space properties of Aomoto-Gel'fand hypergeometric integrals in absence of analytic regulators. This creates a natural path for further investigations of the connections between intersection theory and Feynman integrals, which are left for the future.
The paper is organized as follows: In Sec. 2 we begin by recalling the basics of the Feynman integrals in terms of twisted de Rham (co)homologies and their intersection theory. We show the representation of both the integral and its dual, together with the master decomposition formula needed for their direct decomposition. We discuss different ways to compute the dimension of the cohomology group. The differential equations satisfied by the forms and the dual forms are also provided. Then follows Sec. 3 in which we discuss multivariate intersection numbers. We start with an explicit construction of the 2-variable intersection numbers, which is expressed in terms of the univariate ones recursively. This procedure is generalized, resulting in the final formula for the n-variable intersection numbers. We also present an explicit example showing all the steps of the computation of a specific 2-variable intersection number, and discuss a few properties satisfied by the intersection numbers, as well as the simplified formula valid in the case of dlog forms. In Sec. 4, we discuss strategies for the decomposition of an arbitrary Feynman integral. Specifically, we show three different approaches, namely the straight decomposition, the bottom-up decomposition, and the top-down decomposition. Sec. 5 is dedicated to examples. We first consider the one-loop massless box and perform the decomposition with all these three approaches to show the steps involved explicitly. Moreover, we show the decomposition for the QED triangle as well as the differential equation for the QED sunrise. After that, we provide a few tables with all the key ingredients necessary for the computation of the multivariate intersection numbers needed to obtain the direct decomposition, as well as their differential equations, for the cases of the 1-loop box with 4 different masses, the 2-loop sunrise with 3 different masses, the 2-loop planar and non-planar massless triangle-boxes, as well as 2-loop massless double-box on a triple cut. Finally, Sec. 6 contains our conclusions and discussion. The paper ends with Appendix A containing the explicit forms of the multivariate intersection numbers used for the 1-loop massless box, the QED triangle, and the QED sunrise.

Feynman integrals and differential forms
We consider Aomoto-Gel'fand generalized hypergeometric integrals of the form In the context of the Feynman integrals addressed in this manuscript B is the Baikov (graph) polynomial, which has the property that it vanishes on the boundary of the integration domain in (2.1) while γ depends on the space-time dimensionality d, and on the number of loops and external legs. We assume γ to not be an integer, γ / ∈ Z, which follows from dimensional regularization. On the other hand, ϕ(z) is a single valued differential form whereφ L (z) denotes its differential-stripped version, f (z) is a rational function and a i are integer exponents, a i ∈ Z.
One of the key assumptions is that all the poles present in ϕ L must be regulated by u(z).
In genuine Feynman integrals this assumption is often violated; in this work we present two different strategies for overcoming this apparent obstacle. It is possible to identify equivalence classes of differential n-forms entering the integral (2.1). Forms in the same class are those that differ by a covariant derivative and give the same result upon integration, as will be explained below.

The cohomology group and its dual
Consider an (n−1)-differential form ξ L . In the absence of boundary terms due to (2.2) we have: where Thus we can write The forms ϕ L and ϕ L + ∇ ω ξ L , which give the same result upon integration, are in the same equivalence class Differential n-forms modulo the equivalence relation (2.7) belong to a vector space, the twisted cohomology group H n ω , and elements in this vector space are denoted by ϕ L |. In a similar way one can define an equivalence relation among integration contours which give the same result upon integration. Integration contours modulo the equivalence relation, are denoted by |C R ] and belong to the vector space H ω n , referred to as the twisted homology group. The integral of eq. (2.1) can be regarded as a paring between ϕ L | and the function u(z), integrated over the contour |C R ] (2.8) Given this terminology, we may now define a dual integral, given bỹ and consider the covariant derivative In analogy to (2.7) we can derive the equivalence relation such that differential n-forms modulo the equivalence relation eq. (2.11) belong to the dual vector space (H n ω ) * = H n −ω ; the elements of this space are denoted by |ϕ R . As done above, one can also consider an equivalence relation among integration contours, which leads to the vector space (H ω n ) * = H −ω n whose elements are denoted by [C L |. The dual integral of eq. (2.9) is interpreted as paring between |ϕ R and the function u(z) −1 , integrated over the contour [C L | (2.12)

Intersection numbers for twisted (co)homology classes
To summarize, within twisted de Rham theory, ϕ L | and |ϕ R are elements of the twisted cohomology class H n ω and the dual cohomology class H n −ω respectively. Because of a duality between twisted cycles and co-cycles [50], [C L | and |C R ] can be considered as elements of the homology class H ω n and the dual homology class H −ω n . Beside the two type of pairings that defined the integrals and the dual integrals above, we can consider • intersection numbers of twisted cycles [C L |C R ], as introduced in [51]; • intersection numbers of twisted co-cycles ϕ L |ϕ R , which were first considered in [4].
For the complete mathematical definitions of these objects, we refer the reader to the mathematical literature cited above.

Linear and quadratic relations
The reduction of a given integral, I = ϕ L |C R ], in terms of a set of ν MIs, can be interpreted in terms of differential forms, as since the integration cycle is the same for all the integrals of eq. (2.13). Likewise, the decomposition of a dual integralĨ = [C L |ϕ R in terms of a set of ν dual The coefficients c i , andc i in eqs. (2.14), (2.16) are determined by the master decomposition formulas where we introduced the (inverse of the) metric matrix By substituting eq. (2.17) in eq. (2.14) (or eq. (2.18) in eq. (2.16)), we obtain a representation of the identity operator in the cohomology space ν i,j=1 Similarly, in the homology space, the resolution of the identity is where H ij = [C L,i |C R,j ] is the metric matrix for the twisted cycles. The operators I c and I h can be inserted either in the bilinear pairing between the twisted cocyles or the twisted cycles, to obtain the quadratic identities which are known as Twisted Riemann's Period Relations [4]. Let us emphasize that while the choice of the dual basis |h j does not matter for the final result, it might greatly affect the complexity of intermediate steps of the computation. In particular, since the matrix C needs to be inverted, it is especially beneficial to choose |h j maximally orthogonal to e i | (in the sense of the intersection pairing), in order to make C as diagonal as possible. It may not always be clear how to choose such |h j and in fact for the purposes of this paper we will ignore this issue and set e i = h i throughout. This may give large expressions for intermediate intersection numbers, which tend to simplify after plugging in (2.17).

Dimension of twisted cohomology groups
In ref. [48], the number of MIs within the IBP-decomposition was related to the number of independent "contours" of integration, generating no surface terms. In particular, using as correspondence between the basis cycles and the critical points of the graph-polynomial of the considered integral parametrization, the number of MIs was related to the rank of the homology groups H ±ω n . In refs. [1][2][3], we considered a dual, equivalent description of the same problem, in terms of independent "differential forms". Accordingly, we define ν as the dimension of the twisted cohomology group, respectively, H n ±ω , here considered as a vector space, The complex Morse (Picard-Lefschetz) theory allows us to determine ν as the number of critical points of the function log u(z) [48]. We define 2.25) and the number of critical points is given by the number of solutions of the (zero dimensional) systemω i ≡ ∂ z i log u(z) = 0 , i = 1, . . . , n. (2.26) The number of solutions of (2.26) can be determined without computing explicitly its zeros [48]. In our applications the function u(z) always takes the form u(z) = j B γ j j (z), which gives the equations:ω In the absence of critical points at infinity, the number of solutions of (2.26) equals to the dimension of the quotient space for the ideal 1 In the special case when u(z) = B γ (z), it becomes simply [48] Considering a Gröbner basis G generating I, the Shape Lemma (see, e.g. [62], and [29] for an application to physics) ensures that the number ν of zeros of I, and hence the number of the solutions of the system (2.26), is the dimension of the quotient ring, where C[z] is the set of all polynomials that vanish on the zeroes of I (they identify a discrete variety, V ⊂ C ν ). In particular, the lemma ensures that the degree of the remainder of the polynomial division modulo G is ν + 1. In ref. [3], we recalled that ν can be computed using one of the many ways of evaluating the topological invariant Euler characteristics χ(X): X = CP n −P ω , where P ω ≡ {set of poles of ω} in projective space, the above relation can be written as where we used χ(CP n ) = n+1 together with the inclusion-exclusion principle for Euler characteristics. In other words, to compute ν, it is sufficient to evaluate χ(P ω ) of the projective variety P ω (see also refs. [63][64][65]).
In the following, we will compute the dimension of the cohomology groups to determine the size of the basis of differential forms for different choices of H n ±ω , each characterized by ω, or correspondingly by u.

Differential equation for forms and dual forms
Following [3,14], we provide the algorithm for the direct derivation of the systems of differential equations using multi-variate intersection numbers. Let us consider the system of differential equations in an external variable x for the basis e i | and the dual basis |h i , where σ = ∂ x log u. Here Φ i | and |Φ i can be decomposed in terms of e j |, and |h j respectively, by means of intersection numbers using eq. (2.17), and similarly where summation over indices j, k is implied. Using the above ingredients one can relate the matrices Ω and Ω through the identity 38) or in the matrix notation In particular, if the bases were orthonormal such that C = I then Ω = Ω.

Multivariate intersection numbers
In this section we provide the details of the recursive algorithm employed for the evaluation of intersection numbers of multivariate differential forms introduced in [14]; the algorithm was successfully applied in the context of Feynman integrals as well as hypergeometric functions in [3]. The recursive algorithm expresses the n-variable intersection number in terms of (n − 1)-variable intersection numbers and so on, where the last term of this sequence is the univariate intersection number discussed in refs. [1,2].
In particular, we consider integrals with n integration variables {z i 1 , . . . , z in }, which can be seen as iterative integrals, with a nested structure that follows from the chosen ordering {i 1 , . . . , i k } of the integers {1, . . . , n}. In order to compute multivariate intersection number for n-differential forms, we need to compute the dimension of the cohomology groups for all k-differential forms, from k = 1 to k = n. They can be obtained, for instance, by counting the number ν k of solutions to the equation system given by eq. (2.26), . . , n} with k distinct elements. In this way, one obtains a list of dimensions ν 1 , ν 2 , . . ., ν n , respectively corresponding to the iterative integration in . . , z in } variables, and where we used the vector notation, to indicate the integration variables.
It is interesting to observe that, while ν n is trivially independent of the ordering of the integration variables, the dimensions of the subspaces ν k may indeed depend on which specific subsets k of {1, 2, . . . , n} are chosen and in which order. As a working principle, we choose the ordering that minimizes the sizes of ν k for all k-forms (k = 1, . . . , n).
Before delving into the n-variable intersection number, we start with the example of 2-variable intersection numbers, written recursively in terms of univariate intersection numbers; an approach which will later be generalised to the n-variate case.

2-variable intersection number
We start by considering an integral with two integration variables {z 1 , z 2 }, written as follows, L is a differential 2-form in variables z 1 and z 2 , while C R is a twodimensional integration domain embedded in some ambient space X with complex dimension 2. We assume that X admits a fibration into one-dimensional spaces X 2 z 2 and X 1 z 1 2 , and correspondingly ϕ R can be decomposed in a similar manner. Similarly, we can consider a dual integral, given bỹ with all the variables defined analogously to the ones above.
As before, we have and employing eq. (2.26), we can count the number of MIs in the X 1 space, which we label as Then the goal is to determine the 2-variable intersection number ϕ in terms of the univariate intersection numbers on the X 1 space, which are calculable with the univariate methods discussed in [1,2] and assumed to be already computed. We start by decomposing the differential forms as into an arbitrary basis forms e (1) i | and their duals |h In the above expressions ϕ (2) L,i | and |ϕ (2) R,j are one-forms in the variables z 2 , and they are treated as coefficients of the basis expansion. They can be obtained by a projection similar to eq. (2.17), using only univariate intersection, namely (sum over repeated indices is understood) with the metric matrix, which is also a univariate intersection matrix From refs. [1,2] we know that the univariate intersection number is given as is the local solution of the differential equation around every pole p of ω 1 , denoted by the set P ω 1 . Here the connection ω 1 is just the dz 1 component of ω, and ∇ ω 1 = (d + ω 1 ∧). In [14] it was shown that putting these ingredients together one can write the 2-variable intersection number as is the local solution of the differential equation around each point q from the set of poles of Ω (2) denoted by P Ω . In eq. (3.13), the 2variable intersection number ϕ has been expressed in terms of the known univariate intersection numbers and a new connection matrix Ω (2) . To determine Ω (2) , we will follow the same trick as adopted in the case of single variable, namely starting from the integral and defining the equivalence class of the single valued differential form. We want to find an analogue of the fact that with du = ωu, but for two-fold integrals. Let us consider the original integral I from eq. (3.3) and apply the decomposition we used in eq. (3.6): where we defined Now, there could exist many forms ϕ (2) L,i that integrate to give the same result. Let us consider a total derivative of u i times any function (0-form) ξ i (z 2 ) with poles correctly regulated, where d z 2 denotes the differential acting only on z 2 , i.e. d z 2 = dz 2 ∂ z 2 . Let us notice that u i (z 2 ) satisfies the following differential equation in z 2 following Sec. 2.5: where Ω (2) is a ν 1 × ν 1 matrix. Inserting this into eq. (3.18), we obtain: where the final equation defines our new connection ∇ Ω (2) .
The Ω (2) can be obtained directly from computing the z 2 -differential of u i (z 2 ), The final line can be further simplified by using the master decomposition formula in eq. (2.17) in the z 1 -variable, such that Using eq. (3.19), we can identify Ω (2) through Let us discuss an alternative recursive formula for intersection numbers, which uses the dual connection matrix Ω (2) instead of Ω (2) . This amounts to repeating the same steps, but now using the decomposition of the differential forms as described in eq. (3.6). Following [14] the 2-variable intersection number can be written as In the above equation, the 2-variable intersection number ϕ has been expressed in terms of the known univariate intersection numbers and a new connection Ω (2) . To determine Ω (2) one follows steps similar to those describe above.
Let us consider the dual integral with two variables as follows: where we use the decomposition of ϕ R from eq. (3.6) in the first step and defined We then consider a total derivative of u ∨ i times a function ξ i (z 2 ) Using the results from Sec. 2.5, the vector u ∨ i (z 2 ) satisfies the following differential equation in z 2 : where Ω (2) is a ν 1 × ν 1 matrix. Inserting this into eq. (3.28), we obtain: Finally, the matrix Ω (2) can be obtained directly by computing the z 2 -differential of u ∨ i (z 2 ), which shows that its components are given by

n-variable intersection number
Following the above discussion, we can generalize the 2-variable intersection number to the n-variable case, where we start by considering an integral with n integration variables (z 1 , z 2 , . . . , z n ), written as with the notation n = {1, . . . , n}. The ϕ (n) L is an n-variable differential form on some space X. Similarly, one can can define a dual form ϕ (n) R . We assume that the n-complexdimensional space with coordinates (z 1 , . . . , z n ) admits a fibration into a (n−1)-dimensional subspace parametrized by (z 1 , . . . , z n−1 ), denoted by n−1, which we call the inner space, and a one-dimensional subspace with z n , which we refer to as the outer space. We have and employing eq. (2.26), we can count the number of MIs on the inner space, which we define as ν n−1 . The aim is to express the n variable intersection number ϕ in terms of intersection numbers in (n−1)-variables on the inner space, which are assumed to be known at this stage, following the recursive nature of the algorithm. The choice of the variables (and their ordering) parametrizing the inner and outer spaces is arbitrary: as before, we use the generic notation k ≡ {i 1 , i 2 , . . . , i k } to denote the variables taking part in a specific computation.
Thus, the original n-forms can be decomposed according to where ν n−1 is the number of master integrals on the inner space with arbitrary bases e R,i are one-forms in the variable z n , and they treated as coefficients of the basis expansion. They can be obtained by a projection similar to eq. (2.17), giving We stress again that the (n−1)-variable intersection numbers are assumed to be known at this stage. The recursive formula for the intersection number reads [14]: where the functions ψ (n) i are the solution of the system of differential equations andφ L,i are obtained through eq. (3.36). Here,Ω (n) is a ν n−1 × ν n−1 matrix, whose entries are given byΩ and finally P n is the set of poles ofΩ (n) defined as the union of the poles of its entries (including a possible pole at infinity). We observe that the solution of eq. (3.40) around z n =p can be formally written in terms of a path-ordered matrix exponential for a vector ψ (n) with entries ψ (n) i . Nevertheless for its use in eq. (3.39), it is sufficient to know only a few leading orders of ψ (n) around each p ∈ P n . Therefore, it is easier to find the solution of the system eq. (3.40) by a holomorphic Laurent series expansion, using an ansatz for each component ψ i , see [1,2]. Such a solution exists if the matrix Res zn=p Ω (n) does not have any non-negative integer eigenvalues, which we assume from now on (when this is not the case one can employ a regularization discussed in Sec. 4.1). Moreover, the number of critical points of the determinant of the Ω (n) provides the dimension of that cohomology group, i.e. the number of the corresponding MIs, see also [20].
The recursion terminates when n=1, in which case the inner space is trivial: ν 0 = e (0) 1 | = |h (0) 1 = 1, and we impose the initial conditionŝ (3.43) In this case eq. (3.39) reduces to a computation of an univariate intersection number [4,6] previously studied in refs. [1,2]. Let us notice also that combining eqs. (3.39) and (3.37) gives which is suitable for practical calculation purposes. Using the above identity recursively, the intersection number can be expressed as

An explicit example in two variables
Let us consider intersection numbers based on the following function: which givesω We will focus on the steps required for the computation of the intersection number given by These forms only have poles at infinities. Letting 1 = {1} define the inner space, we find corresponding to the fact thatω 1 = 0 has one solution. We choose the inner basis for the left and right forms, denoted by e (1) | and |h (1) respectively, aŝ Given two arbitrary forms ϕ L | and |ϕ (2) R the following decompositions hold: where ϕ with In the recursive approach we assume the one-variable intersection numbers w.r.t. z 1 , to be computed in the previous step. They are given by: while the new 1 × 1 connection matrixΩ (2) is given by: and we see that the poles ofΩ (2) are located at Next, we consider the differential equation: The full analytic solution of (3.62) is not required, but rather a power series around each p ∈ P 2 is sufficient. Denoting by y the local coordinate around the pole, the solutions of (3.62) to leading orders in y read: • Solution around p = 0 (y = z 2 ): (3.63) • Solution around p = 1 (y = z 2 − 1): • Solution around p = ∞ (y = 1/z 2 ): , . (3.66) Finally we may evaluate the bi-variate intersection number as a sum of univariate residues, as given by eq. (3.13): giving the final result for the intersection number: . (3.68) We notice that, in the case at hand, only the residue at p = ∞ gives a non-zero contribution to the intersection number.

Properties of intersection numbers
In this subsection we briefly review some relevant properties of intersection numbers. Their rigorous definition is given by [4] where in order to make the integral well-defined, one needs to regularize at least one of the forms ϕ L or ϕ R by imposing compact support near the boundaries of X, ϕ L/R → ϕ c L/R . Performing this regularization and manipulating the result leads to the concrete expressions in terms of residues given in the previous subsections, see, e.g., [66,Sec. 3.2] for the intermediate steps. This definition can be used to directly prove the following properties.
• Intersection numbers are invariant under a change of differential forms within the same equivalence classes, namely 72) and the covariant derivatives ∇ ±ω defined in eqs. (2.5) and (2.10), explicitly read while ξ L and ξ R are arbitrary (n−1)-forms with poles regulated by u.
• Intersection numbers obey the symmetry relation which follows directly from the definition and the fact that commuting ϕ L with ϕ R yields a sign change of (−1) n . We stress that the right-hand side is evaluated with respect to −ω rather than ω.

Intersection numbers of logarithmic forms
Intersection numbers for multivariate logarithmic forms were first considered in [6]. Alternative formulas for more direct calculations were later presented in [12,14]. In particular, if ϕ L and ϕ R are both dlog, we have where the sum goes over all the ν critical points given by the solutions of the system of equationsω as in eq. (2.26). When at least one of the forms is non-logarithmic, the formula (3.75) is only valid asymptotically in the limit γ → ∞. In those cases one can still calculate intersection numbers as a series expansion in 1/γ, which was successfully applied to the computation of differential equations for certain Feynman integrals in [17]. The recursive algorithm for the computation of the multivariate intersection numbers presented in Sec. 3 is applicable for any rational form. However, at each step of the recursive algorithm, the coefficientsφ Thus, under the assumption that the connection matrices Ω (n) and Ω (n) contain only simple poles, its possible to replace the coefficientsφ (n) L,R containing higher-degree poles, with a suitably chosenφ (n) L,R belonging to the same equivalence class, but containing simple poles only. One may exploit this fact to compute intersection numbers in one variable as a univariate global residue, without introducing any algebraic extensions as observed in [20].

Feynman Integral Decomposition
As proposed in refs. [1-3, 17, 20], the use of multivariate intersection numbers yields a direct decomposition of a given Feynman integral I in terms of an a priori chosen set of MIs J i , with i = 1, . . . , ν. The decomposition given by eq. (2.13) is on the form where the determination of the coefficients c i is the goal of this section. We identify three possible strategies which can be adopted in order to achieve this task. They all employ the master projection formula from eq. (2.17), which is applied to differential forms constructed differently in the three cases. We name them the straight decomposition, the bottom-up decomposition, and the top-down decomposition. All the approaches have the first step in common: finding the number of MIs which appear in the decomposition and choosing them accordingly. We introduce the following definitions: • Σ denotes the set of integers used to label the full set of denominators; • σ denotes a set of integers that label a subset of denominators, σ ⊆ Σ; • sector is the set of integrals for which only the subset of propagators specified by σ appear in the denominator (thus, a sector is unambiguously identified by σ).
There is a one-to-one correspondence between sectors and (generalized unitarity) cuts. On the level of the function u, this correspondence is manifested by setting all z j 's belonging to σ to zero in the original u(z), where we work in Baikov representation. Given u σ , the number of MIs in the corresponding sector, ν σ , can be determined through the criteria given in Sec. 2.4. The total number of MIs (without taking into account any symmetry relations) is then given by where the sum is over all sectors. Finally we can choose the forms e i | associated to the (arbitrarily chosen) MIs J i , through the identification (4.4)

Straight decomposition
We consider the following decomposition Hereφ andê i correspond simply to the integrands of the integral I to decompose and of the chosen master integrals, J i , respectively. In order to evaluate the intersection numbers, all the poles present in the differential forms must be regulated in u. If this assumption is violated, we can introduce a regulated u, denoted by u ρ , which contains a monomial z ρ k k for each (non-regulated) pole present in the differential forms, that is and correspondingly where we emphasized the action of regulators. By analogy, we also introduce a regularized version ofΩ (n) , whenever Res zn=pΩ (n) has any non-negative integer eigenvalue. The Thus, we obtain a new system of differential equations, analogous to eq. (3.40), which is, in this case, controlled byΩ (n) Λ . We assume that the solution of the latter around a pole p, denoted by ψ (n) Λ,p , reproduces in the limit Λ → 0, a solution for the original system (around the pole p). The intersection numbers are computed through ω ρ , and lead to a set of coefficients, denoted by c ρ,i , which depend on the set of regulators, collectively indicated by ρ. The coefficients c i , which appear in the original decomposition eq. (4.5), are recovered in the limit ρ → 0 3 This approach requires the evaluation of intersection numbers, for which all the integration variables are present simultaneously. For ease of notation, whenever the regulated u is introduced, in the following we will omit the subscript ρ from the individual intersection numbers ϕ|h j ρ and e i |h j ρ .

Bottom-up decomposition
In this approach, proposed in [3], the decomposition is applied to the spanning set of cuts, defined as the minimal set of cuts such that each MIs appears at least once [3,67] (a cut behave like a pass-high filter, therefore MIs with a number of internal lines smaller than the number of cut variables do not contribute to the decomposition on the cut). We denote a given spanning cut (i.e. an element in the spanning set of cuts) by τ ; moreover S τ is the set of sectors which survive on that spanning cut Finally, the number of MIs which survive on the spanning cut τ , denoted by ν Sτ is On the spanning cut τ , we define and we consider the following decomposition with As expected,φ τ andê i,τ are inferred from the cut-integrals. As in any unitarity-based approach [68][69][70], the coefficients c i determined from a cut decomposition are identical to those appearing in the original decomposition -the coefficients are invariant under cuts. Therefore, the complete decomposition for the (uncut) integral I can be obtained by combining the coefficients determined from the individual spanning cuts. As described in Subsec. 4.1, all the poles present in the differential forms must be regulated in u τ . If this is not the case, we can introduce the regularized u τ , denoted by u ρ,τ which leads to used in the evaluation of the intersection number. We also use a regularized version ofΩ (n) , whenever Res zn=pΩ (n) has any non-negative integer eigenvalue, as explained above. Now, the coefficients of the decomposition, c ρ,i depend on the set of regulators ρ. The coefficients of the original decomposition (4.14) are recovered in the ρ → 0 limit: This procedure requires the evaluation of the intersection numbers only for the uncut variables, therefore it can be significantly less demanding than the previous case. As before, whenever the regulated u is introduced, we will omit the subscript ρ from the individual intersection numbers.

Top-down decomposition
This approach is new and combines the advantages of the decomposition by intersection numbers with the top-down subtraction algorithm traditionally used in methods of integrand decomposition [24,28,30,31]. In particular, as for the integrand decomposition, one can determine the coefficients of the MIs systematically, beginning from the ones with the highest number of internal lines (the top sector) and moving downward, ending with the sector with a minimal number of lines equal to the number of the loops (built from product of tadpoles). At any step, the determination of the coefficients of a given MI, say J i , is obtained on the corresponding cut, after subtracting off the known contributions coming from higher sectors, as the latter are written as a linear combination of the MIs with a higher number of internal lines (whose graph contain the one corresponding to J i as subdiagram), coming from the earlier steps of the decomposition. In particular, let us reconsider the complete decomposition, 19) and assume that, within the top-down approach, after at most n-steps, the coefficients c i , with i = 1, . . . , n have been determined, and can be considered as known. We can write, (4.20) which, in terms of pairings, reads, where φ n |, defined as, is a known differential form. By applying a properly chosen maximal cut, identified by τ , we can then determine the coefficients c i of a number ν τ MIs J i , whose graph contains exactly those lines that are cut. In fact, on the maximal cut τ , we can define and and the decomposition simplifies and becomes, Two important observations are in order. First, we notice that the subtraction in eq. (4.20), is similar in spirit to the subtraction performed in an integrand decomposition, although the known coefficients depend also on d, and not only on the kinematical variables. Second, after the subtraction of the known terms, the differential form φ n,τ may contain spurious poles, which are not regulated by u τ . These poles can be eliminated by redefining φ n,τ , φ n,τ → φ n,τ = φ n,τ + ∇ ωτ ξ L,τ , (4.27) using a suitable ξ L,τ , which can be systematically built. Thus, in this approach, the regulators are not introduced. At this point the determination of the coefficients via intersection numbers can proceed iteratively, top-down, until all sectors have had their c i coefficients determined.

Examples
In this section we illustrate the previously-discussed decomposition algorithms on a few examples. As the first example we will discuss the one-loop massless box. This diagram was discussed in the context of intersection theory already in ref. [3], but we will here add further details, and go through the reduction with each of the three methods presented in Sec. 4.

The one-loop massless box
The kinematics is such that

Performing the Baikov parametrization yields
and performing the sector-by-sector analysis described in the beginning of Sec. 4 yields ν σ = 1 for the sectors and ν σ = 0 for the remaining sectors, corresponding to the well-known set of master integrals: the box and the s-and the t-channel bubble: The corresponding differential forms read In the following we will decompose the example which can be expressed in terms of the chosen master integrals as We will determine these coefficients with the three methods presented in Sec. 4.

Straight decomposition
As prescribed in Sec. 4.1 we may construct the regulated u as where in this case we pick the regulators to be all equal. From this definition we may construct the corresponding ω as Choosing the variable ordering to be, from the innermost to the outermost, z 4 , z 3 , z 2 , z 1 , we can compute the dimensions of the twisted cohomology groups corresponding to the individual layers of the fibration. The result is Corresponding to the order of variables given above, we pick the basis for each level to bê We choose the dual bases to beĥ i =ê i . In the following, we will decomposê The required intersection numbers are C ij = e i |h j , 1 ≤ i, j ≤ 3, (5.14) and The individual intersection numbers, up to the leading order in ρ, are presented in App. A.1.
Combining the intersection numbers as dictated by eq. (4.10), we obtain, after taking the limit ρ → 0, the coefficients These results are in agreement with the values obtained with FIRE [71].

Bottom-up decomposition
The first step of a bottom-up decomposition is to identify a spanning set of cuts τ . That set is easily seen to be the cuts corresponding the two bubbles Let us first consider the τ = {1, 3} cut. On this cut, the decomposition reads: We have where and ω ρ,τ =ω 2 dz 2 +ω 4 dz 4 witĥ The variable ordering, from the innermost to the outermost, is chosen as z 2 , z 4 . The dimensions of the cohomology groups read: The basis elements, on the cut, are: The dual basis elements are chosen asĥ i,τ =ê i,τ . We will show the decomposition, on the cut, of: This requires the intersection numbers and Expressions for the individual intersection numbers are presented in Appendix A.1. Combining them as prescribed by eq. (4.18), and considering the limit ρ → 0, we obtain the coefficients c 1 and c 2 in agreement with eq. (5.16).
• Cut τ = {2, 4}. Performing instead the decomposition on the second of the spanning cuts, τ = {2, 4} will allow us to reconstruct c 1 and c 3 in eq. (5.16), which means that in total all of the master integral coefficients c i have been extracted.

Top-down decomposition
The first step in the top-down decomposition is the extraction of the box-coefficient.
The coefficient c 1 can be computed as ϕ/e 1 on the maximal cut: Here we have We also get from which we can extract ν τ = 1 corresponding to the s-channel bubble. We know that has to be reducible to the s-channel bubble. This property is not apparent asφ contains poles in z 2 and z 4 that distinguishes the box and the bubble sectors. However, we know that φ is in the same equivalence class as a φ without these poles. Writing we may make the following ansatz for ξ, Fitting the free coefficients κ with the requirement that all poles of φ in z 2 or z 4 vanish, gives a solution The corresponding φ is of the formφ where P is a polynomial, so we see explicitly that the z 2 and z 4 poles are gone, and that no poles are present in φ that are not poles of ω. With this we may perform the bi-variate intersections, and we get in agreement with eqs. (5.16). The expressions for the two intersection numbers are listed in App. A.1, and please note that they are much simpler than for the other two approaches due to the absence of the regulator. For the t-channel cut one may proceed likewise, and extract the coefficient of the t-channel bubble, again in agreement with eqs. (5.16).
Let us note that one could use the subtraction in eq. (5.33), where κ 1 is a free coefficient. Then, the fitting of the unknown coefficients of eq. (5.35) generates a system whose solution does require the value κ 1 = c 1 . In other words, κ 1 , which in this case corresponds to the coefficient of a master integral in the higher sector (the box function) may be fixed together with the remaining κ-parameters 4 .

Discussion of the example
Considering the three intersection-based reduction methods of the one-loop massless box, we see that the straight decomposition required the computation of 12 4-variate intersection numbers, the bottom-up decomposition required 12 2-variate intersection numbers, and the top-down decomposition required 4 2-variate intersection numbers. Due to the recursive nature of the multivariate algorithm of Sec. 3, the computation of a 2-variate intersection number is much easier than a 4-variate intersection number, thereby showing the efficiency of the bottom-up algorithm compared to the straight decomposition. On the other hand, in the top-down decomposition, we compute fewer intersection numbers than in the other two approaches, but there is a trade-off since in this approach we may have to perform a fit of the extra κ-coefficients in the subtraction terms as done in eq. (5.34) for the one-loop massless box, which may become computationally expensive in a generic case.

The QED triangle
In this subsection we discuss the one-loop QED triangle [2]. The denominators are and the kinematics is such that p 2 1 = p 2 2 = m 2 , (p 1 + p 2 ) 2 = s. The Baikov parametrization yields: Performing the sector-by-sector analysis described in the beginning of Sec. 4 we obtain ν σ = 1 for the sectors and ν σ = 0 for the remaining ones.
The master integrals are chosen as: 44) and the corresponding differential forms read In the following we will decompose: which can be expressed in terms of the chosen master integrals as

Straight decomposition
We introduce a regularized u given by and then We consider the ordering of the variables, from the innermost to the outermost layer, as z 3 , z 1 , z 2 and the dimension of the twisted cohomology groups are Given the order of variables considered above, we chose the basis elements to bê while the dual basis elements are chosen asĥ i =ê i . The required intersection numbers are and ϕ|h k , 1 ≤ k ≤ 3. Combining the intersection numbers, and taking the ρ → 0 limit as in eq. (4.10), we obtain . (5.54) These coefficients are in agreement with the result obtained from FIRE [71] (before applying any symmetry relations). Here, we consider 2-loop QED sunrise diagram as shown in Fig. 3. The denominators are:

The QED sunrise
while the ISPs are chosen as: The Baikov parametrization gives: We choose the invariant p 2 1 = s and normalise it by the squared internal mass effect m 2 , effectively setting m 2 = 1, and the m 2 dependence can be recovered later by power counting. We perform the sector-by-sector analysis for each of the 7(= 2 3 − 1) sectors as described in Sec. 4, and obtain zero MIs in all sectors except for where for the sector {1, 2, 3} we obtain 3 MIs and for the sector {1, 3} 1 MI, amounting to a total of 4 MIs. The MIs are chosen as the following: 61) and the corresponding differential forms read Here, we will build the differential equation for the set of master integrals We will now determine Ω using the bottom-up decomposition as described in Sec. 4

Bottom-up decomposition
First we identify a spanning set of cuts τ . That set is easily seen to only contain the cut corresponding to the double tadpole: On this specific cut, we use: and ω ρ,τ =ω 2 dz 2 +ω 4 dz 4 +ω 5 dz 5 witĥ We consider the ordering of the variables, from the innermost to the outermost, as z 4 , z 2 , z 5 and the corresponding numbers of independent forms read: On the cut we have where J ρ,τ = C u ρ,τ e τ and the differential equation reads The twisted cocycles corresponding to the individual MIs on the cut arê Following eq. (2.35) we define σ = ∂ s log u ρ,τ and the corresponding twisted cocycles for the decomposition of eq. (5.70) read: For the inner spaces, we choose the basis elements as: Using eq. (3.39) we may then get the individual entries of the differential equation matrix The individual multivariate intersection numbers are provided in App. A.3. Using these intersection numbers, we obtain after taking the limit ρ → 0 which is in agreement with the result obtained from LiteRed [72].

Further Examples
In the following, we present the key information useful to perform the reduction by means of intersection theory, in a set of cases all corresponding to physically relevant Feynman integrals. In particular, for each case, we provide a table containing: the definition of the integral family; the spanning cuts (τ ); the dimensions of the vector spaces at each step of the recursive algorithm (ν) and the corresponding bases (e), for the evaluation of multivariate intersection numbers; a pictorial decomposition of a generic integral, whose coefficients can be determined by means of our master decomposition formula eq. (2.17). In all these cases, the reduction and/or the differential equations were computed successfully, in agreement with the results of public IBP codes [71][72][73][74].

Conclusions
In this work we elaborated on the vector space structure of Feynman integrals, and on the existence of what amounts to a scalar product among them first presented in Ref. [3], showing a detailed description of their systematic decomposition in terms of Master Integrals. In particular, we described the evaluation of multivariate intersection numbers for twisted cocycles, which are the key ingredient of the master decomposition formula eq. (2.17), in terms of a recursive algorithm boiling the computations down to univariate intersection numbers. We applied the master decomposition formula to derive integral relations and differential equations for a number of Feynman integrals. As shown in previous works [1,2], they can also be used for deriving dimensional recurrence relations (finite-difference equations) for Feynman integrals. We discussed algebraic properties of integrals and dual integrals as well as systems of differential equations they obey. We provided three different strategies for Feynman integral reduction, which we dubbed the straight decomposition, the bottom-up decomposition, and the top-down decomposition, which show possible combinations of the intersection-theory concepts together with unitaritybased methods and integrand decomposition.
The recursive computation of multivariate intersection numbers requires regulated integrals, not plagued by spurious irregular behavior which might emerge at the intermediate steps of the evaluation. For this purpose, we employed the analytic regularization procedure. On the other hand, using the richer mathematical structure of the relative twisted cohomology, the use of regulators might be avoided, thereby offering a very interesting new direction for future studies and applications to physics.
Let us conclude by observing that the decomposition formula, or better the corresponding formula for the identity resolution, in terms of multivariate intersection numbers, is applicable to generic parametric representations of Feynman integrals, including those not considered here. More generally, it can be used to derive linear and quadratic relations for Aomoto-Gel'fand type of integrals (and their duals), which are of broad interest and have applications in different contexts in physics as well as mathematics.  [76]. CloudVeneto is acknowledged for the use of computing and storage facilities.

A Intersection numbers for the three examples
In this appendix we provide the explicit form of intersection numbers needed for the Feynman integral decompositions performed in Sec. 5. Since we work in analytic regularization with a parameter ρ that is taken to zero at the end of the computation, it suffices to know only the leading ρ-orders of intersection numbers. While our algorithm computes them exactly in ρ, in order to save space in this appendix we list only the leading term for each intersection number individually. One can check that the orders given here are sufficient for reconstructing the coefficients c i to order O(ρ 0 ) and that their limit as ρ → 0 is in fact smooth.

A.1.1 Straight decomposition
Here we provide the intersection numbers, up to the leading order in ρ required for the decomposition presented in Subsec. 5.1.1: e 3 |h 1 = e 2 |h 1 , (A.8) and (A.14)

A.1.3 Top-down decomposition
For consistency with the straight decomposition and the bottom-up decomposition, we also provide here the intersection numbers needed for the top-down decomposition of Subsec. 5.1.3, on the τ = {1, 3} cut. They are . (A.23)

A.2 The QED triangle
Here we provide the intersection numbers, up to the leading order in ρ, required for the system of differential equations presented in Subsec. 5.2:

A.3 The QED sunrise
Here we provide the intersection numbers, up to the leading order in ρ, required for the system of differential equations presented in Subsec. 5.3: