Vector Space of Feynman Integrals and Multivariate Intersection Numbers

Feynman integrals obey linear relations governed by intersection numbers, which act as scalar products between vector spaces. We present a general algorithm for constructing multivariate intersection numbers relevant to Feynman integrals, and show for the first time how they can be used to solve the problem of integral reduction to a basis of master integrals by projections, and to directly derive functional equations fulfilled by the latter. We apply it to the derivation of contiguity relations for special functions admitting multi-fold integral representations, and to the decomposition of a few Feynman integrals at one- and two-loops, as first steps towards potential applications to generic multi-loop integrals.


INTRODUCTION
Scattering amplitudes encode crucial information about collision phenomena in our universe, from the smallest to the largest scales. Within the perturbative fieldtheoretical approach, the evaluation of multi-loop Feynman integrals is a mandatory operation for the determination of scattering amplitudes and related quantities.
Linear relations among Feynman integrals can be exploited to simplify the evaluation of scattering amplitudes: they can be used both for decomposing scattering amplitudes in terms of a basis of functions, referred to as master integrals (MIs), and for the evaluation of the latter. The standard procedure used to derive relations among Feynman integrals in dimensional regularization makes use of integration-by-parts identities (IBPs) [1], which are found by solving linear systems of equations [2] (see [3,4] and references therein for reviews). Algebraic manipulations in these cases are very demanding, and efficient algorithms for solving large-size systems of linear equations have been recently devised, by making use of finite field arithmetic and rational functions reconstruction [5][6][7].
In [8], it was shown that intersection numbers [9] of differential forms can be employed to define (what amounts to) a scalar product on a vector space of Feynman integrals in a given family. Using this approach, projecting any multi-loop integral onto a basis of MIs is conceptually no different from decomposing a generic vector into a basis of a vector space. Within this new approach, relations among Feynman integrals can be derived avoiding the generation of intermediate, auxiliary expressions which are needed when applying Gauss elimination, as in the standard IBP-based approaches.
In the initial studies, [8,10], this novel decomposition method was applied to the realm of special mathematical functions falling in the class of Lauricella functions, as well as to Feynman integrals on maximal cuts, i.e. with on-shell internal lines, mostly admitting a one-fold inte-gral representation. Those results concerned a partial construction of Feynman integral relations, mainly limited to the determination of the coefficients of the MIs with the same number of denominators as the integral to decompose, which was achieved by means of intersection numbers for univariate forms.
In this paper, we make an important step further, and address the complete integral reduction, for the determination of all coefficients, including those associated to MIs corresponding to sub-graphs. In the current work, we discuss the one-loop massless four-point integral as a paradigmatic case, although the algorithm has been successfully applied to several other cases at one-and two-loop.
Generic Feynman integrals admit multi-fold integral representations. Their complete decomposition requires the evaluation of intersection numbers for multivariate rational differential forms. Intersection numbers of multivariate forms have been previously studied in [11][12][13][14][15][16][17][18][19]. Recently, a new recursive algorithm was introduced in [20]. In this letter, we present its refined implementation and application to Feynman integrals, which provide a major step towards large-scale applicability of our strategy for the reduction to MIs. The results of this work show potential for further applications ranging from particle physics, through condensed matter and statistical mechanics, to gravitational-wave physics, while making new connections to mathematics.

INTEGRALS AND DIFFERENTIAL FORMS
In this work, we focus on integrals of the hypergeometric type, where: z = (z 1 , . . . , z n ) are integration variables; C is the integration domain; u is a multi-valued function of the form u = i B i (z) γi with γ i / ∈ Z, such that i B i vanishes on the integration boundary ∂C; and ϕ is a single-valued differential n-form, withφ being a rational function with all poles regulated by u(z). Then employing Stokes' theorem we find equivalence classes of n-forms, for any (n−1)-form ξ and where ∇ ω ≡ d + ω∧ is a covariant derivative with a one-form ω ≡ d log u. The space of n-forms modulo the relation eq. (3) forms a vector space called a twisted cohomology group 1 H n ω . We denote its elements by ϕ| ∈ H n ω . Within this framework, the integral I from eq. (1) can be interpreted as a pairing of ϕ| with the integration contour |C], Since in our applications |C] will always stay constant, the vector space of such integrals is the same as that of ϕ|.
Consider a set of ν MIs, say J i , defined as in terms of any independent set of differential forms e i |. Then, the decomposition of a generic integral I in terms of the MIs J i , can be interpreted as coming from the more fundamental decomposition of the differential form ϕ| in terms of the basis forms e i | , namely with the coefficients determined by the master decomposition formula [8,10], 1 We refer the interested reader to [20][21][22] for reviews of twisted (co)homologies and their intersection theory, as well as [8,10,18,20,[23][24][25][26][27][28] and [29][30][31][32] for some recent applications of these ideas to physics.
Using eqs. (6,8), our algorithm for expressing any integral of the type of eq. (1) as linear combinations of MIs proceeds along three steps: 1. Determination of the number ν of MIs.
2. Choice of the bases of forms e i | and |h i .
3. Evaluation of the intersection numbers for multi-variate forms, appearing in the entries of the C-matrix, and in ϕ|h j .
Number of Master Integrals. Under some assumptions one can show that all other vector spaces H k =n ±ω are trivial, which means that ϕ can only be n-forms [33]. In those cases the dimension of these vector spaces, i.e. the number ν of MIs, can be determined topologically 3 , in terms of the Euler characteristic χ(P ω ) of the projective variety P ω defined as the set of poles of ω. This connection allows us to use complex Morse (Picard-Lefschetz) theory to determine ν as the number of critical points of the function log u(z). Let us define then the number of critical points is given by the number of solutions of the system of equationŝ with the short-hand notation ∂ zi ≡ ∂/∂z i , provided that the set of solutions is finite. Additional details are provided in the App. A.

INTERSECTION NUMBERS
In this section we review a recursive algorithm for the evaluation of intersection numbers of multivariate differential forms introduced in [20].
We start by decomposing the n-dimensional space with coordinates (z 1 , . . . , z n ) into a (n−1)-dimensional subspace parametrized by (z 1 , . . . , z n−1 ), which we call inner space, and a one-dimensional subspace with z n , dubbed outer space. The aim is to express the original intersection number n ϕ (n) L |ϕ (n) R in terms of one-dimensional residues on the outer space and intersection numbers n−1 . . . | . . . on the inner space, which are assumed to be known at this stage. The choice of the variables (and their ordering) parametrizing the inner and outer spaces is arbitrary: in the following, we use the generic notation m ≡ (12 . . . m) to denote the variables taking part in a specific computation.
Thus, the original n-forms can be decomposed according to ϕ (n) where ν n−1 is the number of master integrals on the inner space with arbitrary bases e (n−1) i |, |h (n−1) j and the metric matrix In the above expressions ϕ R,j are dz n -forms treated as coefficients of the basis expansion. They can be obtained by a projection similar to eq. (8), for example: where from now on we use the implicit sum notation for repeated indices. The recursive formula for the intersection number reads where functions ψ (n) i are solutions of the system of differential equations where ϕ R,i dz n from eq. (15). The ν n−1 ×ν n−1 matrixΩ (n) given bŷ and finally P n is the set of poles ofΩ (n) given by the union of the poles of its entries (including possible poles at infinity). We observe that the solution of eq. (17) 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. (16), 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. (17) by a holomorphic Laurent series expansion, using an ansatz for each component ψ (n) i , see [8,10]. 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.
The recursion terminates when n=1, in which case the inner space is trivial: In this case eqs. (16,17) reduce to a computation of univariate intersection numbers [9,12] previously studied in [8,10]. Plugging everything together, eq. (16) can be expressed as where the ranges of summations are i m = 1, . . . , ν m and each ψ im−1im dz m coming from the projection: which is known a priori, because the bases of all inner spaces are arbitrarily chosen. The matricesΩ (m) needed in eq. (22) are computed analogously to eq. (18). Notice that all ψ (m) entering eq. (21) need to be computed only once for a given family of integrals. The multivariate intersection number given in eqs. (16,21) is the key formula used in this letter. Paired with the master decomposition formula eq. (8), the above recursion for intersection numbers yields an expansion of multi-fold integrals of the form in eqs. (1,4) in terms of MIs.

HYPERGEOMETRIC FUNCTION 3F2
In order to illustrate application of the above algorithm we start with a more familiar case of contiguity relation for the hypergeometric function 3 F 2 . Consider the function H, defined as, where β(a, b) = Γ(a)Γ(b)/Γ(a+b) is the Euler betafunction, In this case, ω is defined through eq. (10) witĥ where z i+1 should be understood as z 1 if i = 2. The systemω 1 =ω 2 = 0 has three solutions, corresponding to ν (12) = 3 MIs. We choose three master forms, e which correspond to the following set of MIs, At the same time, we define the dual basis, |h yields the decomposition of the function defined in eq. (24) in terms of those in eq. (28), which amounts to a contiguity relation for 3 F 2 functions. The coefficients c i are determined by means of eq. (8), requiring the computation of 12 intersection numbers for two-forms, that is 9 elements of the matrix (C (12) ) ij = (12) e for i, j = 1, 2, 3. To apply formula eq. (16), we consider the z 2 -subspace as the inner space. In turn, the number of MIs for the inner space is determined by counting the number of The individual intersection numbers are too large to be printed here. Yet, the final result is rather simple, and, in terms of 3 F 2 -functions, it reads, where the X i are the arguments of the functions H in eq. (28) and This relation has been (numerically) verified with Mathematica.

FEYNMAN INTEGRAL DECOMPOSITION
Within the Baikov representation (BR), Feynman integrals can be cast in the form 4 (1), with either u(z) = B γ (z) [39,40] or u(z) = i B i (z) γi [41]. The factors B and B i are graph (Baikov) polynomials, and their exponents depend on the dimensional parameter d, and on the number of loops and external momenta of the corresponding diagram. The number n of integration variables corresponds to the number of scalar products formed by the external and the loop momenta. In fact, within BR, the propagators of the diagrams, supplemented by a set of auxiliary propagators (related to the irreducible scalar products), are the integration variables. Therefore, ϕ in (2) can be generically written as where a i ∈ Z, and where f is a rational function of z.
Multiple-cut integrals, identified by the on-shell conditions z i1 = . . . =z i k =0, are also of the form (1), but their integrands depend on fewer integration variables (and their integration contour is modified), see [8,10]. Bottom-up Decomposition. For a given integral, any set of its denominators identifies a sector. Therefore, one maximal-cut (when all denominators are cut) corresponds to each sector. The number of MIs in each sector can be determined by counting the number of critical points of the corresponding maximal-cut 5 , using eq. (11). After determining the number of MIs, the decomposition of Feynman integrals can be obtained by means of eq. (8). This is done by redefining u by multiplying it with a regulating factor z ρi i for each uncut denominator, where the exponents ρ i are regulators to be put to zero at the end of the calculation. This is done in order to allow the theory to access sectors where the regulated variables appear as propagators. The determination of coefficients can be performed on unitarity cuts, where the integrands are simpler, and the evaluation of the multivariate intersections requires fewer iterations. A minimal set of spanning cuts will be sufficient to retrieve the information on the complete decomposition [42], and then, using the regulated u, the master decomposition formula (8) yields the coefficients of those MIs that survive on the cut. As in the case of IBP-based approaches, additional relations may be obtained from the symmetries of the diagrams, in order to minimize the number of independent integrals.
As discussed in refs. [8,10] also differential equations in kinematic variables, e.g. ∂ s J i = j a ij J j , can be obtained with the above techniques.
so that any integral I of the form of eq. (1), with u given in eq. (34), and ϕ defined in eq. (33) (with n = 4), can be decomposed as, 6 If the Baikov polynomial B is a non-zero constant on the maximal cut, the integral is fully localized by the cut-conditions. In this case, the condition ω = 0 is always satisfied, and there is ν = 1 master integral.

Differential Equation.
Let us consider the differential equation: where we restore the s-dependent prefactor: On the Cut 1,3 we obtain: ∂s . On this cut we have: On the Cut 2,4 we have: withφ 2,4 = f z1z3 . On this specific cut we obtain: where a 1 is in agreement with eq. (52) and Let us finally remark that we have successfully applied the aforementioned algorithm to the complete decomposition of a few one-and two-loop integrals associated to the diagrams shown in Fig. 2, involving the evaluation of up to six-variable intersection numbers, and that the resulting expressions are in agreement with the IBP relations [43][44][45][46]. Further examples are provided in the App. B.

CONCLUSIONS
Elaborating on the original proposal of [8] and on the wider studies [10,20], we have shown that Feynman integrals can be expressed in terms of a complete basis of integrals, by making use of intersection numbers, which act as scalar products for the vector space of integrals, through the pairing of differential forms appearing in their integrands. Let us notice that the final result of the recursion eq. (21) should not depend on the parametrization of the inner and outer space. Nevertheless, we observed that suitably chosen variable orderings may simplify and fasten the recursive procedure. This is a feature of the proposed algorithm that requires a dedicated study, which goes beyond the goal of the present work. Within Baikov representation, one-loop and multi-loop integrands have a similar structure, and therefore we expect that our decomposition algorithm can be applied to the case of integrals associated to more complex diagrams than the ones considered here, which we plan to investigate in the near future.
Scattering amplitudes are analytic functions, determined by their singularities. Intersection numbers, and their relation to Stokes' and Cauchy's residue theorems, embed what we believe is a clean role of analyticity in the amplitudes decomposition. We investigated the geometric origin of master integrals within the formalism of twisted (co)homology, where it was possible to relate them to the number of critical points and Euler characteristics in the connection to Morse/Picard-Lefschetz theory (as a special case of the Poincaré-Hopf index theorem). Applications to Feynman integrals in representations other than Baikov will also constitute topics of future works. The present study can be broadly applied in the context of theoretical particle physics, condensed matter and statistical mechanics, gravitational-wave physics, as well as mathematics. comments. H. Let us consider a single-valued k-form ϕ k and a multivalued function u(z) integrated over a k-real-dimensional submanifold C k ⊂ X inside of some space X of complex dimension n, If u(z) regulates all boundaries of C k then by Stokes' theorem: where ∇ ω ≡ d + ω∧ is a covariant derivative with a oneform ω ≡ d log u(z). Thus adding terms of the form ∇ ω ϕ k−1 to ϕ k does not change the value of the integral of eq. (56). Similarly, we can impose that integrals over boundary terms of the form ∂C k+1 vanish: which corresponds to ∇ ω ϕ k = 0. These two requirements define a set of natural vector spaces for k = 0, 1, . . . , 2n: called twisted cohomology groups [21]. Under some assumptions amounting to the fact that u(z) regulates all boundaries of X, one can show that in fact H n ω is the only non-trivial space and all other H k =n ω vanish [33]. From now on we consider only such cases, even though Feynman integrals are known to sometimes violate these assumptions [10,35].
One can also construct a dual vector space (H n ω ) * = H n −ω , with the same properties, given by a replacement ω → −ω in the above definition eq. (59). In this work we consider ϕ L | ∈ H n ω and |ϕ R ∈ H n −ω and a scalar product ϕ L |ϕ R called the intersection number [9] 7 .
The Euler characteristic χ(X) of the space X can be computed as an alternating sum of dimensions of H n ω , Since all H k =n ω vanish, we find that the dimension of H n ω , and hence also the number ν of MIs is given by Thus ν can be computed using one of the many ways of evaluating the topological invariant χ(X). We review a few of them below. Since X = CP n −P ω , where P ω ≡ {set of poles of ω}, we can simplify the above relation to where we used the fact that χ(CP n ) = n+1 and the inclusion-exclusion principle for Euler characteristics. The computation thus amounts to evaluating the Euler characteristic χ(P ω ) of the projective variety P ω , see [36][37][38] for related approaches. Let us introduce a simple function u(z) that will serve as an instructive example in the remainder of this appendix: which arises physically from the maximal cut of a two-loop non-planar triangle diagram [8] and gives rise to Appell F 1 functions with some constants s, ρ, γ. Computing ω = d log u(z) gives straightforwardly P ω = {±ρ, ±s, ∞}, and hence X = CP 1 −P ω is a one-dimensional space parametrized by an inhomogeneous coordinate z. The point at infinity is removed from X since Res z=∞ (ω) = 0.
Since the Euler characteristic of 5 distinct points is simply χ(P ω ) = 5, using eq. (62) we find: which is the correct number of MIs in this case [8]. 7 Similarly, eq. (56) is a scalar product ϕ k |C k ] between H k ω ϕ k | and the twisted homology group H ω k |C k ], which is non-zero only for k=n. Since |Cn] is always constant in Feynman integral computations, H n ω can be also regarded as the vector space of Feynman integrals in a given family with the same ω.
Let us now consider a real-valued function h(z) ≡ Re(log u(z)), called a Morse function, which assigns a "height" to every point z ∈ X. Special role in this construction is played by critical points z * of h(z) defined by dh(z * ) = 0. Using Cauchy-Riemann equations it is straightforward to show that this condition is the same as ω = n i=1ω i dz i = 0 and thus the critical point equations readω We assume that all critical points are isolated and nondegenerate. To each of them the Morse function assigns a pair of flows, labelled by a sign ± and parametrized by an auxiliary "time" variable τ , (66) In the − case we have dh(z)/dτ < 0 and hence it corresponds to a downward flow from the α-th critical point z * (α) , which defines a submanifold of X called a Lefschetz thimble (or a path of steepest descent) J α with some real dimension λ α . Similarly, the + case defines an upward flow, which generates a path of steepest ascent K α through the critical point z * (α) , with real dimension 2n−λ α . Here λ α is the number of unique negative directions extending from the α-th critical point, called its Morse index.
One of the key results in complex Morse theory (often called Picard-Lefschetz theory) is that the Euler characteristic can be expressed as [49]: where M λ is the number of critical points with the Morse index equal to λ. Since u(z) is a holomorphic function, near each z * (α) we can pick local coordinates w (α) (with the critical point at w (α) =0) such that the Morse function admits an expansion: Treating X as a real manifold with coordinates w (α) = x (α) + iy (α) we find and hence every critical point has a shape of a saddle with exactly n upward and n downward directions, or the Morse index λ α = n. This means that only M n is non-vanishing and hence using eqs. (61) and (67) we find [33]: In the context of Feynman integrals these arguments were first given in [35]. The critical points can be also used to compute asymptotic behavior of intersection numbers [18]. Let us mention that Lefschetz thimbles are integration contours along which eq. (56) converges the most rapidly for k=n, and thus the set {J α } n α=1 can be used a basis of integration cycles 8 . For explicit examples of projecting cycles onto such bases using homological intersection numbers see App. A of [20].
In the example at hand, eq. (63) gives ν=3 solutions of the critical point equations, in agreement with eq. (64). The form of Lefschetz thimbles depends on the values of s, ρ, γ and here we choose ρ>s>0 and γ>0 as a concrete example. With this choice each J α=1,2,3 has to have endpoints on z ∈ {±ρ, ±s} since this is where h(z) decays to −∞, while K α=1,2,3 can only have endpoints on z = ∞ as it is the only place where h(z) → +∞. This alone fixes the shape of the paths of steepest descent and ascent uniquely up to contour deformations. We illustrate them in Fig. 3.
FIG. 3: Morse-Smale complex associated to the Morse function h(z) = Re(log u(z)) with eq. (63) and ρ>s>0, γ>0. The set of filled dots corresponds to P ω = {±ρ, ±s, ∞} removed from X. Empty dots at z * (α) represent critical points of the Morse function, with paths of steepest descent J α (solid lines) and ascent K α (dashed lines) extending from them. They give a triangulation of X = CP 1 −P ω . The arrows indicate the direction of the flow towards lower values of h(z).
The critical points together with paths of steepest of descent and ascent triangulate the manifold X into what is known as a Morse-Smale complex. Denoting the number of q-dimensional elements of this complex by b q (called the Betti number) we have (72) 8 Likewise, the paths of steepest ascent of h(z), Kα are integration cycles along which the dual integral Kα u(z) −1 ϕn converges the most rapidly and {Kα} n α=1 can be used a basis of H −ω n .
For more background on Morse theory, see, e.g., [49,50] and in the context of twisted geometries [18,20,21,33]. Now, with the help of eq. (8) and using (16) for the computation the individual multi-variate (here 2 and 3forms) intersection numbers, we determine the coefficients c i in (77).