From Infinity to Four Dimensions: Higher Residue Pairings and Feynman Integrals

We study a surprising phenomenon in which Feynman integrals in $D=4-2\varepsilon$ space-time dimensions as $\varepsilon \to 0$ can be fully characterized by their behavior in the opposite limit, $\varepsilon \to \infty$. More concretely, we consider vector bundles of Feynman integrals over kinematic spaces, whose connections have a polynomial dependence on $\varepsilon$ and are known to be governed by intersection numbers of twisted forms. They give rise to differential equations that can be obtained exactly as a truncating expansion in either $\varepsilon$ or $1/\varepsilon$. We use the latter for explicit computations, which are performed by expanding intersection numbers in terms of Saito's higher residue pairings (previously used in the context of topological Landau-Ginzburg models and mirror symmetry). These pairings localize on critical points of a certain Morse function, which correspond to regions in the loop-momentum space that were previously thought to govern only the large-$D$ physics. The results of this work leverage recent understanding of an analogous situation for moduli spaces of curves, where the $\alpha' \to 0$ and $\alpha' \to \infty$ limits of intersection numbers coincide for scattering amplitudes of massless quantum field theories.


Introduction
One of the more surprising developments in the recent study of scattering amplitudes has been the introduction of scattering equations [1], which allow for writing tree-level amplitudes and loop-level integrands-such as those of Yang-Mills or gravity theories-in terms of certain localization integrals on moduli spaces of punctured Riemann spheres [2,3]. Scattering equations can be understood as critical-point conditions for a certain "potential" function W , determined by the vanishing of its first derivative, dW = 0.
It was later understood that such localization formulae are not at all specific to moduli spaces and can be broadly extended to more general cases [4]. To be specific, let us consider a complex manifold M written as a complement of a finite number of hypersurfaces in CP m with inhomogeneous coordinates (z 1 , z 2 , . . . , z m ), as well as a potential function W (z 1 , z 2 , . . . , z m ) with logarithmic singularities on those hypersurfaces. To two top holomorphic forms, ϕ − = ϕ − d m z and ϕ + = ϕ + d m z, we associate a pairing, which following [4,5] we state as a Grothendieck residue around the critical points, where d m z is the measure form and ∂ i = ∂/∂z i . More geometrically, it should be understood as a self-duality pairing of the cohomology of the Koszul complex (Ω • M , dW ∧). When applied to the moduli space of Riemann spheres with n punctures, M =M 0,n , the pairing (1.1) coincides with the Cachazo-He-Yuan formula [2].
The first correction to (1.1) is given by the higher residue pairing [6] Res dW =0 1 2 Notice that it has m+1 powers of W in the denominator, compared to just m in (1.1). This motivates an introduction of a book-keeping parameter τ and sending W → τ W (in quantum field theory τ is proportional to the inverse of Planck's constant −1 , in string theory it is the inverse string tension α , while for Feynman multi-loop integrals it becomes the dimensionregularization parameter ε). All higher residue pairings may be compactly written as a τ −1 expansion of a single object, ϕ − |ϕ + dW = (1.1) + τ −1 (1.2) + · · · , which in fact gives a compact expression that generates all-order corrections. Geometrically (1.3) is the intersection number of cohomology classes associated to the twisted de Rham complexes (Ω • M , d±τ dW ∧), see, e.g., [26,27], which will be reformulated in terms of aČech-de Rham double complex later in the text.
The physical meaning of intersection numbers on M =M 0,n (with τ = α ) is that they compute tree-level scattering amplitudes of quantum field theories with a finite spectrum of masses, m 2 ∈ Z/α [4,5], which are rational functions of kinematic invariants with simple poles of the form 1 p 2 +Z/α . As a matter of fact, they were used to resolve a long-standing puzzle regarding scattering equations, which-despite computing low-energy physics [2]-determine worldsheets dominating in the high-energy limit of string theory [28,29]. 1 On the one hand, in the α → 0 limit intersection numbers coincide with the low-energy limit of string-theory scattering amplitudes. On the other hand, in the α → ∞ limit they reproduce the localization on scattering equations. One may ask when the two limits agree. This clearly happens when the intersection number is independent (or homogeneous) of α in the first place [4,5], as then it does not matter if we send α → 0 or α → ∞! Physically, this property corresponds to propagators of the form 1 p 2 , i.e., when intersection numbers compute amplitudes of massless quantum field theories.
The main goal of this paper is to leverage this new understanding to other problems in scattering amplitudes.
In particular, we focus on multi-loop Feynman integrals in D = 4−2ε space-time dimensions, as those have a known interpretation in the same geometric language [30]. The reason for employing dimensional regularization is that such integrals most often do not converge in strictly four dimensions. They can be written as where Γ is some integration cycle and the potential W is determined in terms of so-called Symanzik polynomials that specify the topology of a given graph G. In this case the role of τ is played by ε and M =M G is the moduli space of Riemannian metrics on G with coordinates given by Schwinger parameters. Alternatively, Feynman integrals can be understood as sections of vector bundles over the kinematic space, defined by the solutions of the system of differential equations (D − Ω∧) I = 0, (1.4) where I is a vector of integrals of the type (1.3), D is the differential on the kinematic space, and Ω is a matrix-valued one-form, subject to integrability constraints, that needs to be determined. Together with boundary conditions, which we assume are known, (1.4) fully characterizes the behavior of Feynman integrals in a given family around ε → 0. Thus the problem amounts to finding the matrix Ω. It was recently realised that fibers of the vector bundle can be described by the cohomology of (Ω • M , d+τ dW ∧) and hence the entries of Ω can be computed by the same intersection numbers (1.3) described above [30,31].
As a matter of fact, on physical grounds Ω must be a polynomial in ε, as it can be shown that any pole in ε must be spurious (see, e.g., [32]), (1.5) Because of this we can expand the matrix Ω around either ε → 0 or ε → ∞ and still obtain the exact result with a finite number of terms! We use the latter option, which combined with the expansion (1.3) allows us to compute Ω in terms of higher residue pairings (1.1), (1.2), and their further corrections. Notice that critical points contributing to these computations correspond to places on the moduli space M G that normally dominate the ε → ∞ physics. This is yet another example of what seems to be a more general moduli space localization phenomenon [5], in which physical quantities in one limit can be extracted from the exact opposite one. We illustrate this new idea by performing explicit computations for two families of integrals. We start with arguably the simplest case of a single-box massless integral and follow with a two-loop sunrise diagram with masses running in the loops.
Outline. In Section 2 we review the geometric setup based on twisted de Rham andČechde Rham cohomologies, which leads to explicit expressions for higher residue pairings. In Section 3 we formulate Feynman integrals as twisted periods and describe how to obtain their differential equations from higher residue pairings. Explicit examples are given in Section 3.3 for the one-loop box diagram and in Section 3.4 for the two-loop massive sunrise diagram. We conclude in Section 4 with a discussion of future directions. This paper comes with Appendix A, where we clarify the computation of α → ∞ asymptotics of genus-zero string amplitudes.

Geometric Setup
In this section we briefly review the geometric setup underlying the remainder of the paper. The understanding of Sections 2.1-2.3 is not needed to compute higher residue pairings in practice, but rather is meant to give an intuition about where they come from. Explicit expressions for higher residue pairings are given in Section 2.4, and the way of relating them to integrals over middle-dimensional cycles is explained in Section 2.5.

From Koszul to Twisted de Rham Complex
The formula (1.1) can be understood geometrically in the following way [5]. Let us consider Integrals defined on such spaces are ubiquitous in physics, e.g., in Feynman multi-loop integrals or string perturbation theory.
We introduce a holomorphic function W on the covering space M of M with logarithmic singularities along each H i . For instance, if H i are defined by equations {f i = 0} then for generic constants α i , with k i=1 α i = 0, is a valid choice of W . 2 We will often call W a potential to use the same nomenclature as in the literature on mirror symmetry [33]. Let us consider a single-valued holomorphic one-form dW . Defining Ω k M to be the space of smooth k-forms on M (with O M := Ω 0 M being the space of functions), we introduce the following sequence: Here each map is given simply by wedging the element ϕ k ∈ Ω k M from the left with dW ∧, that is for k = 0, 1, . . . , m−1. Since dW ∧ dW = 0, the sequence (2.2) is exact, meaning that image of each map is equal to the kernel of the following one (applying two consecutive maps gives zero). This allows us to construct cohomology groups H k (M, dW ∧) associated to (2.2), which are given by kernel of each dW ∧ modulo the image of the preceding dW ∧, or in other words One can show that only the case k = m is non-trivial, provided the constants α i are generic [34]. In addition, since dW ∧ defines a rank-1 flat connection we have which allows us to compute the dimension of the above cohomology group purely topologically in terms of the Euler characteristic χ(M ) of M . From now on we assume that (W ) is a Morse function [35] with isolated and nondegenerate critical points. It is easily seen that such critical points are given by dW = 0, i.e., coincide with the critical points of the potential function W . For later convenience let us introduce notation for the critical locus of W : (2.6) which by the above assumptions is a finite set. Since W is holomorphic, all critical points have the same Morse index, i.e., the same number of independent upwards and downwards directions extending from it. This tells us that [36,37] which combined with (2.5) allows one to compute the dimension of H m (M, dW ∧) by counting critical points. One can define a self-duality pairing of H m (M, dW ∧), which is given by [4] (ϕ − |ϕ + ) dW,0 := Res dW =0 Here the hat denotes stripping an overall differential from a form, ϕ d m z := ϕ and d m z = ∧ m i=1 dz i . Alternatively we can think of the hatted function as being defined in the ring of functions modulo the ideal generated by ∂ i W = 0, (2.9) The symbol Res dW =0 denotes a sum over Grothendieck residues around each critical point [38], which is simply given by Here the contour is oriented by d(arg ∂ 1 W ) ∧ · · · ∧ d(arg ∂ m W ) > 0 and has support only on small tubular neighbourhoods are each critical point. It will be evaluated directly in many situations later in the text.
It is important to note that the pairing (2.8) is not unique. As a matter of fact, we can embed it into the following formalism. Consider the exact sequence called a twisted de Rham complex, where we introduced a differential which defines an integrable connection, or equivalently a flat line bundle, since ∇ 2 dW = 0. As was the case before, we can define cohomology groups based on this complex 3 which are spaces of ∇ dW -closed modulo ∇ dW -exact forms. We will call H k dW := H k (M, ∇ dW ) for short from now on. As in the case of (2.4), only k = m gives a non-trivial cohomology [40] and hence we have dim H m dW = (−1) m χ(M ). (2.14) Let us introduce a dual twisted cohomology H m −dW defined in the same way as H m dW but with W → −W . We will often refer to cohomology classes [ϕ ± ] ∈ H m ±dW as twisted cocycles and specific representatives ϕ ± as twisted forms. Duality of the two cohomologies is induced by the intersection pairing and called an intersection number. The overall normalization is chosen for later convenience.
Since M is non-compact, for this definition to make sense one needs to use a compactlysupported form ϕ c + (i.e., one which vanishes in infinitesimal neighbourhoods of the boundary divisor ∂M ) in the same cohomology class as ϕ + for direct computations. This not only makes the result depend on the potential W , but also means the above integral localizes on ∂M , as in the bulk of M we have ϕ − ∧ ϕ c + = ϕ − ∧ ϕ + = 0 for two top holomorphic forms, see, e.g., [5]. Intersection numbers are rational functions of α i 's and τ . Different ways of evaluating them in practice were given in [4,5,27,[41][42][43][44][45][46][47][48][49]. For spaces M admitting a fiber bundle decomposition (or, more precisely, such that the connection decomposes generically on the fibers), the most efficient computation method is currently given by recursion relations [5]. When M = M 0,n is the moduli space of genus-zero curves with n marked points, intersection numbers have an intrinsic interpretation as computing tree-level scattering amplitudes of quantum field theories [4,5]. Worldsheet models that might underlie these computations were studied in [50][51][52][53]. In those cases the corresponding potential W can be given an interpretation as an electrostatic potential for a system of particles [54].
In the context of Landau-Ginzburg models, intersection numbers between basis elements compute entries of a metric g ij = ϕ i |ϕ j dW on the space of physical operators, see, e.g., [12]. A particular problem, related to the theory of Frobenius manifolds, is finding bases that make this metric flat, i.e., g ij = δ ij .
Since (2.4) looks like a limit τ → ∞ of (2.13) it is natural to expect that the residue pairing (ϕ − |ϕ + ) dW,0 should be related to the limit of the intersection number ϕ − |ϕ + dW . Before deriving this result in full generality, let us consider the one-dimensional case dim C M = 1 to gain some intuition about this relationship.

Example: One-Dimensional Case
Let us consider a one-dimensional case, where each hypersurface H i is a single point, say {z = p i }, removed from CP 1 , and the corresponding potential is W = k i=1 α i log(z−p i ) with α i 's adding up to zero. The Euler characteristic is simply χ(M ) = 2 − k, which is the same as the number of critical points, as the solutions of dW = 0 are roots of a degree-(k−2) polynomial in z.
The intersection number of two twisted forms ϕ − ∈ H 1 −dW and ϕ + ∈ H 1 dW is defined by Let us construct ϕ c + explicitly as which is manifestly cohomologous to ϕ + . Here Θ(x) is a step function equal to one for x>0 and zero otherwise. In this way, each term in the sum has support on an infinitesimal disk around p i with radius ε. The inverse differential ∇ −1 dW is defined such that ∇ dW ∇ −1 dW η = η. In particular, ∇ −1 dW ϕ + is a zero-form. Let us check that the resulting form has compact support by expanding the above expression: (2.20) The first term vanishes inside small disks around each p i and the second term has only support on the circles with radii ε imposed by Dirac delta functions δ(x). Therefore ϕ c + has compact support. One could have performed the same computation in a smooth way with bump functions instead of step functions, leading to the same final result, see, e.g., [41], but will not do it here for the sake of clarity.
Plugging (2.20) back into (2.18), the first term wedges to zero and only the second contribution survives, which straightforwardly expresses the intersection number as a sum of k residues around each point removed from M : Note that because of the residue, ∇ −1 dW ϕ + needs to be computed only locally as a holomorphic expansion around the each p i to some finite order in z−p i (depending on the order of the pole of ϕ − ). Simple power counting reveals that a given boundary at {z = p i } gives non-zero contribution only when the orders of poles of ϕ − and ϕ + at this point add up to at least two. Of course, we could have imposed compact support on ϕ − instead, which after an analogous computation gives a different representation It will later turn out to be convenient to symmetrize between the two type of expressions, but for the time being let us stick with (2.21).
In order to make connections to the residue pairings, let us expand the inverse of the twisted differential ∇ −1 dW in powers of τ −1 , where ∂ z := ∂/∂z. This expression can be confirmed by imposing ∇ dW ∇ −1 dW ϕ + = ϕ + orderby-order in τ −1 . Substituting this expansion into (2.21), we can notice two facts. The first one is that k−2 new poles, at the positions of each critical point, have been introduced to the argument of the residue. Secondly, argument of each residue is now the same one-form. This allows us to deform the original contour from enclosing ∂M = ∪ k i=1 {z = p i } to enclosing the set of critical points Crit(W ) by the residue theorem. Therefore we obtain where ∇ −1 dW is understood as an expansion in (2.23). We can now start collecting terms proportional to different powers of τ −1 , where we assumed for simplicity that ϕ ± themselves are independent of τ . For example, the leading term is given by which coincides with (2.8) for m = 1. As a matter of fact, we have an infinite number of corrections (ϕ − |ϕ + ) dW,k given in (2.25), which can be straightforwardly obtained by expanding (2.24) to higher orders. These are the simplest examples of higher residue pairings [6]. The above example motivates looking for generalizations to higher-degree forms. In principle one should be able to carry out a similar derivation, starting with the integral expression (2.16) and showing that it localizes on Crit(W ) using a global residue theorem, though this path requires an involved computation. Fortunately, we can circumvent it by a change of perspective, by considering an extension of the twisted de Rham complex.

TwistedČech-de Rham Complex
In this section we give an alternative, though equivalent, definition of intersection numbers that evaluates directly to the localization formula (2.8) and all its τ −1 corrections. We follow the work of Saito [6][7][8][9] (for reviews see, e.g., [55,56]) in a language adapted to the present context.
We start by introducing a (locally-finite) open cover It allows us to define a double complex (C • (U, Ω • M ), δ, ∇ dW ) called the twistedČech-de Rham complex, which is an extension of (2.11), Here each C p (U, Ω q M ) denotes the space of p-cochains of the cover U with coefficients in q-forms on M , such that their elements [57] for a textbook reference. For instance, in the two extreme cases p = 0 and p = m−1, which will be of our main interest, we have Each vertical line in (2.28) then becomes copies of a twisted de Rham complex with a differential ∇ dW . In the horizontal direction we have aČech coboundary operator δ satisfying δ 2 = 0, which acts as where the hat denotes an omitted index. One can check that δ ∇ dW − ∇ dW δ = 0. Let us group terms along the anti-diagonal of (2.28) by defining followed by an introduction of the differential operator D : K p+q → K p+q+1 given by One can check that it satisfies which allows us to define a cohomology of the complex (K • , D) often called the hypercohomology of the cover U with coefficients in the double complex (2.28) and denoted by H p (U, (Ω • M , ∇ dW )). As before, replacing W → −W at all steps allows us to define a dual hypercohomology.
The intersection number (2.16) can be re-stated in this formulation as [6] ϕ The equivalence to (2.16) follows from the fact that both definitions satisfy Saito's uniqueness theorem [6], up to an overall constant. This constant is fixed by matching the leading τ → ∞ asymptotics of the intersection number to (2.8), which was done independently in [4], and determines the prefactor τ m on the right-hand side of (2.34) in our conventions.
We can compute ψ + as follows. Let us first use the embedding defines an element in the top-left corner of the twistedČech-de Rham complex (2.28). Explicitly, (2.36) Then ψ + , in the bottom-right corner, can be obtained by solving DΨ + = ϕ + and extracting ψ + as the p = m−1 component of Ψ + . Concretely this can be done by tracing a zig-zag path through the diagram: This gives us an expression for ψ + , which involves applying the inverse operator ∇ −1 dW m times and δ m−1 times in alternating order: Since we are interested in explicit formulae, let us show how to compute ψ + step-by-step. It will be convenient to introduce the following notation for each component of ∇ dW , from which the general patter should be clear (the offset by 1 is simply a consequence of our conventions for the covering U = {U i+1 } m−1 i=0 whose index traditionally starts from 0 instead of 1). After m steps we find where the inverses ∇ −1 i are understood in terms of their expansion around the τ → ∞ limit. Note that ∇ −1 i 's commute and hence we do not need to specify the order in which they are applied. We will evaluate each order in τ −1 in the following subsection.
As a matter of fact, we can derive a dual formula for intersection numbers given by whose evaluation is entirely analogous. Defining In general, one can interpolate between the two definitions (2.34) and (2.42) by considering a paring with p+r = m−1 and q+s = m, however we will not do so here.

Higher Residue Pairings
To summarize, we found that intersection numbers of the cohomology classes [ϕ ± ] ∈ H m ±dW can be expressed in terms of Grothendieck residues around the critical points Crit(W ) as follows: , and the residue is defined as in (2.10). We also have a dual formula for the same object, obtained by switching the roles of ϕ − and ϕ + : We are interested in expanding such intersection numbers in powers of τ −1 , as follows: where the leading order starts at τ 0 because of the overall normalization of intersection numbers. The coefficients (ϕ − |ϕ + ) dW,k are called higher residue pairings [6]. They have the symmetry property For instance, when ϕ + = ϕ − all the odd higher residue pairings vanish identically.
In the following subsections we extract higher residue pairings directly from the expressions (2.45) and (2.46) by expanding each inverse derivative operator according to which can be shown by requiring that ∇ i ∇ −1 i η = η order-by-order. Likewise we have which is obtained simply by replacing τ → −τ in (2.49). In general the two types of expansions will involve similarly-looking terms that could cancel out upon averaging between (2.45) and (2.46). We exploit this fact in deriving explicit expressions for k = 0, 1, 2 below.

Leading Order
At the leading order we see straightforwardly that the two expression evaluate to which was in fact shown previously in [4] using complex Morse theory. In order to evaluate (2.51) explicitly, let us make use of the m × m Hessian matrix Φ with entries which, by the assumption on non-degeneracy of the critical points, has maximal rank and hence is invertible. To compute the residue we first change the form variables from (z 1 , z 2 , . . . , z m ) to (∂ 1 W, ∂ 2 W, . . . , ∂ m W ), at a cost of dividing by the Jacobian det Φ, so that we obtain Here we also used the fact that ϕ ± do not have poles on the critical locus Crit(W ). It is known that if both twisted forms ϕ ± are logarithmic, their intersection number is homogeneous in τ and (ϕ − |ϕ + ) dW,0 is the only non-vanishing residue pairing, see, e.g., [5,42].

Subleading Order
Expanding (2.45) to order τ −1 we find at subleading order (2.54) Using the expression (2.46) leads to a similar expression which is related to the above one by exchanging ϕ − ↔ ϕ + and an overall minus sign. Thus after symmetrizing the result the final two terms cancel out and we are left with which matches the subleading higher residue pairing [6]. Repeating the steps from the previous subsection we obtain: where the first line is a result of performing a residue around the double pole in ∂ i W , while in the second line we changed the variables back to (z 1 , z 2 , . . . , z m ) in order to evaluate the derivative explicitly.

Subsubleading Order
Expanding (2.45) to order τ −2 , we find that the subsubleading correction to the intersection number is given by the higher residue pairing Due to the complicated combinatorics we will not explicitly expand the derivatives above. However, let us comment on some features of the above expression. The subsubleading term has poles in ∂ i W of maximal order 5, which come from the first line of equation (2.57). After expanding the derivatives one finds multiple sums, the largest being a 4-fold sum with O(m 4 ) terms. Equation (2.46) also leads to a similar expression for the subsubleading pairing, which is related to (2.57) by exchanging ϕ − ↔ ϕ + , and hence there would be no cancellations after symmetrizing.
Given that the dimension of H m dW is |χ(M )|, an arbitrary twisted form ϕ can be expanded into a basis {ϕ a } |χ(M )| a=1 of this cohomology group. This can be done explicitly by is a basis of the dual cohomology group H m −dW , which is orthonormal in the sense that ϕ ∨ a |ϕ b dW = δ ab . Naturally, the above equality implies a relation between integrals (2.58). As a matter of fact, similar decomposition can be achieved in the homology basis, leading to a |χ(M )| × |χ(M )| period basis of integrals, but we do not use it as in our applications Γ is always kept constant.
There is one caveat, however, in that an orthonormal set of bases might not be easily found, as is in fact generically the case. This can be alleviated by introducing an auxiliary basis which follows from a simple linear algebra exercise. In this way we can use (2.60) to perform expansion of an arbitrary integral into a basis. This simple property, when used together with higher residue pairings, turns out to be quite powerful.

From Infinity to Four Dimensions
In this section we apply the formalism reviewed above to extract the information about analytic properties of Feynman integrals. After reviewing their representation as twisted periods in Section 3.1, we discuss how to construct vector bundles of such integrals over the kinematic space in Section 3.2, followed by a determination of their connections in terms of intersection numbers and higher residue pairings. We finish with explicit examples in Sections 3.3-3.4 for one-and two-loop diagrams.

Feynman Integrals as Twisted Periods
Let us start by reviewing how to translate a given Feynman integral from its momentum-space form into a representation using Schwinger parameters. An L-loop integral with P propagators {D a } P a=1 is given by where the integration contour is ( 1 , 2 , . . . , L ) ∈ (R 1,D−1 ) L in Lorentzian signature and the overall constant is for later convenience. We use mostly-plus conventions for the metric. Each integral is labelled by integers (ν 1 , ν 2 , . . . , ν P ) ∈ Z P that specify the powers to which the corresponding denominators are taken. There is no substantial difficulty in repeating our analysis for multi-loop scattering amplitudes (sums over multiple Feynman integrals of the above type), by allowing the set {D a } P a=1 to be large enough, however we focus on individual integrals in order to make universal statements that do not depend on a specific quantum field theory.
We employ a "Schwinger trick" in which each denominator D a is expressed as an integral over a variable x a representing Schwinger time associated to the corresponding edge of the Feynman diagram: For the time being let us not worry about a possible divergence of the gamma function and treat ν a as formal parameters. Applying (3.2) to the above Feynman integral P times will involve the following combination in the exponent: P a=1 x a D a =: which defines the L × L matrix Q, the length-L vector L, and the scalar c in terms of kinematic invariants and x a 's. Since our goal is to integrate out the loop momenta, we first complete the square in this combination: The Gaussian integral over i 's gives (iπ D/2 ) L /(det Q) D/2 thus cancelling the prefactor in (3.1) (the factors of i come from a Wick rotation to Euclidean time) and the resulting expression becomes The integration contour is now independent of the space-time dimension D and from now on we employ dimensional regularization by setting D = 4 − 2ε.
Let us briefly comment on the meaning of the space M with coordinates (x 1 , x 2 , . . . , x P ). Given a graph G whose internal edges are prescribed by the set of propagators {D a } P a=1 , each Schwinger parameter x a parametrizes proper length the edge associated to D a . For this reason we will refer to M as the moduli space M G of Riemannian metrics on G.
The above integral is already in a form similar to (2.58), however for our purposes we would like to massage it into an integral where the potential W is proportional to ε, i.e., identify it with the expansion parameter τ in (2.58). To this end we follow a standard procedure by inserting 1 = R + δ(ρ − M a=1 x a )dρ into (3.5), followed by rescaling x a → ρx a . Collecting all the Jacobians this leaves us with where |ν| := P a=1 ν a . For brevity of notation we also expressed the result in terms of the so-called Symanzik polynomials: One can recognize that the ρ integral is of the form (3.2) and hence evaluates to On the other hand we can rewrite the right-hand side of (3.8) using Feynman parametrization as (3.9) Followed by rescaling x a → x a /ρ and undoing theρ integration with R + δ(ρ − M a=1 x a )dρ = 1 we obtain the representation [73]: This family of integrals is almost of the form (2.58), if it were not for the following fact. Taking the potential function suggested by the above representation, τ W = ε log(F +U), leaves us with forms ϕ of the type M a=1 x νa−1 a dx a /(F+U) 2 . These might have poles on {x a = 0} and/or {x a = ∞} depending on the values of ν a 's, which violate the assumption that the pole divisor of ϕ is contained within the divisor of the integration space M (similar issue appears in the definition of the integration cycle R P + ). This is actually a physical effect, since such singularities of ϕ correspond to propagators pinching, and thus cannot be removed. It signals that one should have instead considered a twisted homology relative to such singularities, see [72] for a formulation of Feynman integrals in this setup.
Nevertheless, the philosophy of the present paper is that one should study properties of integrals on M globally and in particular without worrying about stratification of its boundaries and related issues. Thus we follow a different path and redefine (3.10) by infinitesimally deforming the integer parameters ν a to ν a → ν a + εδ a , (3.11) where each δ a is a generic infinitesimal variable. We send each δ a to zero at the end of the computation and assume that the resulting regulated integrals I ν 1 ,ν 2 ,...,ν P are smooth in this limit. With this deformation we have W = log(F+U) + P a=1 δ a log(x a ), (3.12) and since M is defined as CP m minus the pole divisor of dW , boundaries of M now include {x a = 0} and {x a = ∞} for all a. We define twisted cohomologies H P ±dW with this potential and identify τ = ε. This leads to a family of twisted forms: where |δ| := P a=1 δ a . We will often set all δ a 's to be equal, δ a = δ. Note that the deformation also regularized the gamma functions. In this language (regulated) Feynman integrals become twisted periods: where the middle-dimensional integration cycle is Γ := R P + and the hat denotes the fact that we expect (3.14) to agree with (3.1) only after taking the limit δ a → 0. Closely related ways of rewriting Feynman integrals as twisted periods were introduced in [30], see also [31,68,70,72,[74][75][76][77].
To simplify ε-power counting we will normalize basis forms by appropriate powers of ε such that they start at ε 0 . In addition, rescaling the integration variables x a → βz a by a constant β typically allows one to remove one mass-scale outside of the integral, which is what we will do in the explicit examples below.
The dimension of twisted cohomology groups H P ±dW is the absolute value of the Euler characteristic |χ(M )| of M , which is given by where C × := C−{0}, in agreement with [31,78]. Physically it counts the number of linearlyindependent Feynman integrals that involve the set of propagators {D a } P a=1 over Q(K, ε, δ a ), where K in the set of kinematic variables appearing in Q, L, c. 4 It is the most convenient to compute |χ(M )| by invoking Morse-theory arguments, which for sufficiently generic W imply that it is equal to the number of critical points Crit(W ) determined by the condition dW = 0. From this point of view the regulators δ a are needed to ensure that the Morse flow is transverse to the divisors {x a = 0} and {x a = ∞}. Explicitly, dW = 0 gives a system of equations ∂ a (F+U) F+U + δ a x a = 0 (3.16) for a = 1, 2, . . . , m. Recall that in our work we always assume the critical points are isolated and non-degenerate. Since we will be interested only the limit δ a → 0, it is sufficient to solve the above constraints as an expansion of z a in δ a , which greatly simplifies finding the critical points. (One should not expand the equations (3.16) themselves, as the number of solutions is generically discontinuous in such a procedure.) According to [73], the number of critical points that solve (3.16) with strictly δ a = 0 computes the number of top-level Feynman integrals in a given family.
There is no agreed-upon way of finding bases of Feynman integrals. There is, however, a criterion for what constitutes a "good" basis based on its behaviour near the ε → 0 limit. This leads us to the following discussion.

Vector Bundles Over the Kinematic Space
Feynman integrals are functions on the kinematic space K with coordinates given by the kinematic variables K, such as Mandelstam variables or particles' masses. For our purposes, however, it is more convenient to think of a basis of Feynman integrals as a section of a vector bundle V over the kinematic space. To make this concrete, let us split the differential operator on the total space as where the matrix-valued one-form Ω is given by intersection numbers [30], as a special case of (2.60): Since the integration domain Γ is always kept constant, each fiber of V is isomorphic to a twisted cohomology H P dW , where W is determined by a point K ∈ K on the base space. Connections obtained in the above way are always integrable (flat), meaning that D − Ω∧ squares to zero, i.e., DΩ − Ω ∧ Ω = 0. (3.20) In physics parlance we are considering a non-abelian gauge theory, with zero curvature and gauge group GL(|χ(M )|, C), on the kinematic space K. Fixing a gauge corresponds to choosing a basis of Feynman integrals. The equation (3.19), together with a specification of boundary conditions, can be understood as an alternative definition of a family of Feynman integrals [79][80][81]. As a matter of fact, solving such differential equations provides one of the most efficient ways of evaluating Feynman integrals in practice, see, e.g., [32,82,83] for reviews. Our goal will therefore be to derive the connection in (3.19).
In physically relevant situations the matrix Ω can be expanded as a polynomial in ε [32], To be more precise, if Ω had any pole in ε, it would have to be spurious [32] and here we assume that such poles do not appear. This means we can evaluate intersection numbers that comprise the entries of Ω either as an expansion around ε → 0 or ε → ∞ and both of them truncate. The latter expansion can be consistently carried out using higher residue pairings. As demonstrated in [30,31,68] the intersection numbers in Ω can be also computed exactly in ε, though the techniques used there rely on the knowledge of either stratification of M or its fibration properties. Using higher residue pairings allows us to forget about these technicalities. In many cases one can bring (3.21) into a so-called ε-form, where only a single term Ω (1) is non-vanishing [84]. In those cases iterating the system of differential equations becomes particularly simple (an additional simplification would be if the matrix was triangular). A basis leading to an ε-form of Ω is called canonical. Literature on these aspects of differential equations includes [84][85][86][87][88][89][90][91][92][93][94][95][96].
There are two points we should discuss before diving into explicit computations. Firstly, the twisted form (D + εDW ∧)ϕ a in (3.18) is manifestly non-homogeneous in ε, and therefore we should consider the two terms Dϕ a and εDW ∧ ϕ a separately to be consistent with counting powers of ε. Secondly, as remarked before the dual orthonormal basis {ϕ ∨ b } is typically not accessible a priori. Instead, we can use a more complicated expression where C bc := ϑ b |ϕ c dW as in (2.61) and D ca := ϑ c |Dϕ a dW + ε ϑ c |DW ∧ϕ a dW . of H P −dW , which we take to be independent of ε.
In order to make ε-power counting simple, we normalize each ϕ a such that it contains only positive powers of ε up to the global maximum ε n ϕ a = n k=0 ε k ϕ (k) a (3.24) In terms of higher residue pairings we have simply: Taking into account non-homogeneity of twisted forms in D we have (3.27) To be concrete let us finish by giving an expression for the connection matrices in terms of the above quantities: ) All computation of Euler characteristics χ(M ) and differential equations below have been double-checked with computational software Macaulay2 [97] and FIRE6 [98] respectively.

Example I: One-Loop Massless Box
Let us study arguably the simplest example that illustrates the idea behind this paper in a straightforward manner. We consider a one-loop scalar massless box graph G box , where P =4 and the set of propagators is given by where all external momenta are massless, i.e., p 2 i = 0. The Symanzik polynomials (3.7) are given by They depend on the two Mandelstam invariants s = (p 1 +p 2 ) 2 and t = (p 2 +p 3 ) 2 . We can factor out one of these mass-scales by rescaling x a = z a /(−s) and defining y := t/s to be the only kinematic variable. This leaves us with the family of integrals where we factored out an overall kinematics-independent normalization for later convenience. Setting δ a = δ, the potential is given by Twisted forms are defined through . (3.34) The moduli space of metrics on the box graph, M G box = (C × ) 4 − {G = 0}, has the Euler characteristic |χ(M G box )| = 3, as can be checked with Macaulay2 [97], and therefore we need to choose three basis forms {ϕ a } 3 a=1 to span H 4 dW . Following [32], we take: The powers of s, y, and ε are chosen such that the resulting basis elements depend only on the ratio y and in particular be independent of s and ε to leading orders in δ. Explicitly, we have For simplicity we also choose the same basis for the dual cohomology {ϑ a } 3 a=1 ∈ H 4 −dW , ϑ a = ϕ a , which guarantees that the intersection matrix C, as in (3.25), starts at order ε 0 in the expansion around ε → ∞. In this case, equations (3.25) and (3.26) simplify to As the first step in computing intersection numbers we find the critical points given by dW = 0, whose positions to the order O(δ 2 ) are given by Let us remark that the fact that only the last critical point remains at a non-singular position as δ → 0 signals that there is only one top-level diagram (e.g. given by ϕ 3 ) in agreement with the program Mint provided in [73]. Computing matrix C from (3.25) using higher residue pairings gives to leading orders: Vanishing of the diagonal entries of C (−1) is actually an exact-δ statement since we used the same bases for the two cohomologies and the subleading higher residue pairing is antisymmetric. In order to compute the matrix D from (3.26) we first need to evaluate how the kinematic space differential D = dy∂ y acts on forms. We have as well as Plugging into the definition (3.26) we find to order O(δ 2 ): (3.42) Finally, we put everything together using (3.28) and (3.29), which can now be truncated to the finite order in δ and read Note that even though some entries of the inverse matrix C −1 (0) contain terms of order O(δ −1 ) coming from O(δ 4 ) of C (0) not given above, they drop out in the final contraction C −1 (0) D (1) . Assuming no subsubleading terms contribute, which would violate ε-polynomiality, we can write the differential equations matrix as in agreement with FIRE6 and [32]. The connection has simple poles at y = 0, −1, ∞, which physically correspond to singularities at s=0, u=0, and t=0 respectively. This gives a definition of the vector bundle (3.19) for the family of box integrals.
Therefore we must choose a basis of seven forms {ϕ a } 7 a=1 to span H 3 dW , which we take to be where we have chosen normalization constants that will turn out to bring the connection into the form Ω = Ω (0) +εΩ (1) . Explicitly, to leading orders in δ we have (3.53) For simplicity, we also take ϑ a=1,2,3 := ϕ a ε , ϑ 4 := ϕ 4 1−2ε , ϑ a=5,6,7 := ϕ a 2−3ε (3.55) as our basis for the dual cohomology {ϑ a } 7 a=1 ∈ H 3 −dW , where the normalization is taken to remove the dependence on ε in the leading δ-order compared to (3.53) and (3.54), which simplifies power counting. In particular, equations (3.25) and (3.26) become  for k ≥ 0 since it is not needed for this example.
At this stage let us comment on counting the size of the basis. In the δ → 0 limit, only the first 4 critical points from (3.58) remain at non-singular positions implying that there are 4 top-level diagrams (e.g. given by {ϕ a } 7 a=4 ) in agreement with Mint [73]. (In the general-mass case the behavior of critical points around δ remains unchanged and the same conclusion holds.) We also chose our bases of twisted forms in (3.51) and (3.52) such that {ϕ 1 , ϕ 2 , ϕ 3 } as well as {ϕ 5 , ϕ 6 , ϕ 7 } are related by relabelling symmetry given by permuting the masses (m 1 , m 2 , m 3 ). Including this non-linear relation would mean there are only 3 independent integrals to compute, and only 2 in the top-level. In addition, in the equal-mass case the integrals over {ϕ 1 , ϕ 2 , ϕ 3 } as well as {ϕ 5 , ϕ 6 , ϕ 7 } are equal even without invoking symmetry relations. Let us stress that none of these facts contradicts the result that dim H 3 dW = 7. This is because even though the integrals might evaluate to the same function, the corresponding twisted forms are not cohomologous (the reason why such a possibility exists is that the integration contour in (3.48) is permutation-symmetric and kept constant). The fact that {ϕ a } 7 a=1 provide a basis of twisted cohomology associated to W with (3.49) can be checked by confirming that the matrix C has full rank.
Following the procedure outlined in Sections 3.2 and 3.3, we expand C and D to O(ε 0 ). We found it sufficient to keep their expressions up to order O(δ), but given that they are large and not illuminating, we will not spell out the details here. Their contractions give rise to the expansion of Ω as in (3.28) and (3.29). Assuming polynomiality, the result reads in the equal-mass case: where in agreement with FIRE6 [98]. The connection has simple poles at y = 0, −1, −9, ∞ which correspond to singularities at s = 0, −m 2 1 , −(3m 1 ) 2 and m 1 = 0 respectively. In the generic-mass case, where (m 1 , m 2 , m 3 ) are all distinct, the resulting connection Ω becomes more complicated and would not fit within the margins on this paper. Nevertheless, since the kinematic space is 3-dimensional, we can perform a non-trivial check on integrability of the connection. To be precise, Ω takes the form ( Ω (0,i) + ε Ω (1,i) )dy i (3.63) and hence the constraint (3.20) gives more explicitly, for i, j = 1, 2, 3, , Ω (0,j) ] = 0, (3.64) (1,i) , Ω (0,j) ] = 0, (3.65) [ Ω (1,i) , Ω (1,j) ] = 0, (3.66) at orders ε 0 , ε 1 , and ε 2 respectively. We checked that the result of computing Ω with higher residue pairings satisfies these integrability conditions and agrees with FIRE6.

Discussion
In this work we studied a surprising phenomenon in which the information about Feynman integrals in dimensional regularization around D = 4 can be fully extracted from a finite expansion around saddle points on the moduli space of graphs M G . This behavior is in contrast with a more conventional 1/D expansion of Feynman integrals, previously studied in the context of gravity [105][106][107], which in principle requires an infinite number of corrections to reach D = 4. We nonetheless hope that a deeper connection between the two approaches can be made in the future. One of the more intriguing questions is whether there exists an intrinsic property of a basis of twisted forms that could determine if the connection matrix Ω in (3.18) is homogeneous in ε without direct computations. For instance, it is known that intersection numbers of logarithmic forms are always homogeneous in ε and only the leading higher residue pairing is non-vanishing, see, e.g., [5,42]. In the present context we are looking a property of a basis of twisted forms, rather than an individual one. It would be fascinating to understand a similar geometric condition that leads to an ε-form differential equations, or decide whether such a basis could even exist. Although we used a representation in terms of Symanzik polynomials, as in (3.12), there is no substantial difficulty in repeating our analysis in other ways, e.g., using the original loop-momentum variables [72] or Baikov representation [30,31,68], where the answer to this question might prove easier.
One of the byproducts of our investigation, which has not been given the attention it deserves, is a new connection between intersection numbers in scattering amplitudes and Landau-Ginzburg models. Indeed, for a given potential W one can define a Hilbert space of such a theory with states given by twisted cohomology classes on M . Two-point functions in this model are computed by the intersection numbers ϕ − |ϕ + dW and can be expanded in 1/τ using higher residue pairings, see, e.g., [12]. It would be very interesting to construct concrete realizations of such models in the case of M =M 0,n or M =M G .
Given that a motivation for the present paper partially came from the scattering equations formalism [2], let us comment on why we have not studied higher residue pairings on M =M 0,n in more depth. It is known that intersection numbers compute tree-level amplitudes with poles of the type 1 p 2 +Z/α . For concreteness, let us give an example in a simple case of massive cubic scalar theory with m 2 = 1/α , whose 4-pt amplitude can be written as in the notation of [5]. Clearly, to leading order in the massless limit α → ∞, the intersection number computes the 4-pt amplitude of massless scalars, as expected since the leading higher residue pairing coincides with the Cachazo-He-Yuan formula [2]. Alas, the expansion in 1/α does not truncate in the presence of massive propagators and higher residue pairings do not seem terribly useful in this context. Combined with the fact that scattering equations do not have algebraic solutions for n > 5, suggests that in order to compute scattering amplitudes on M 0,n one should instead employ much more efficient recursion relations [5] that are exact in α . Finally, the theory of primitive forms, which gave rise to higher residue pairings, has been developed in order to generalize the classic theory of elliptic integrals to more general spaces [6,8]. We expect it to play a crucial role in recent developments connecting Feynman integrals to Calabi-Yau geometries [108][109][110][111][112][113][114][115].
but distinct phases. As we shall see, this allows for an infinite number of saddles-point contributions to be resummed into a simple oscillatory term.
In order to analyze the asymptotic behavior of (A.1) we first analytically continue its integrand to a complex variable z, The key observation is that the integrand of (A.2) is no longer defined on the moduli space M 0,4 = {z ∈ CP 1 | z = 0, 1, ∞}, but rather on its universal cover M 0,4 . 5 This is the case because the part of the integrand z α s (1−z) α t is multi-valued and hence defines a infinitelysheeted surface M 0,4 . For example, going around the branching point z=0 in an small anti-clockwise circle p times the integrand changes to e 2πiα ps times its original value, and likewise going around z=1 q times multiplies it by e 2πiα qt . As a result we find infinite number of saddle points characterized by the number of windings (p, q) around z=0 and z=1 respectively (note that winding r times around z=∞ is equivalent to (p, q) = (−r, −r) and thus not independent). There are several ways of consistently computing contributions from all the saddles. For example, we can consider a change of variables In order to impose the last constrain we introduce a Lagrange multiplier w, so that (A.2) becomes α yields an infinite number of solutions given by for every (p, q) ∈ Z 2 which are the winding numbers introduced above. Indeed, it is easily seen that undoing the change of variables (A.5) each (u * , v * , w * ) is mapped to the same z * = s/(s+t), which is why it "looked" like a single saddle point to begin with. In order to discern which saddles contribute to the α → ∞ limit, one should first find all 5 The space described by the number of windings (p, q) ∈ Z 2 is in fact the maximal Abelian cover of M0,4, which is smaller than the universal cover M0,4 (covering group of the former is the 1-st homology group H1(M0,4, Z) as opposed to the fundamental group π1(M0,4), cf. [118]), but is sufficient for our purposes.
the cycles of steepest descent (Lefschetz thimbles) J p,q and ascent K p,q of the Morse function (W ) = (su + tv + 2πiw(e u + e v − 1)) (A.7) emanating from each saddle-point labelled by (p, q). Clearly, if the integral (A.4) was over such J p,q , it would have been dominated by the (p, q)-th saddle in the α → ∞ limit. Thus, it remains to translate the original integration over Γ into those over the steepest descent paths, which can be done by using the relation Here Γ|K p,q ∈ Z is the homology intersection number of the corresponding cycles, see [119,120] for standard references. Using saddle-point approximation and separating the (p, q)-dependent terms, in the high-energy limit we find Crucially, the shape of each J p,q and K p,q might change drastically depending on values of the parameters s and t. For instance, when s, t > 0, the original integration domain Γ is already homologous to J 0,0 and the sum f (s, t) in (A.8) has only one term equal to 1. As a consequence, only a single saddle contributes to the high-energy limit. Note that it is exactly the one that in the original variables lies on the integration contour, i.e., z * ∈ (0, 1).
Nevertheless, in a generic kinematic region, none of the steepest descent cycles equals to Γ and the sum in (A.8) generically involves an infinite number of terms and consequently an infinite number of saddles contribute to the high-energy limit. As we will see, in those cases f (s, t) can be resumed into a concise expression. 6 There is however a different approach to this problem, using homologies with local coefficients (or twisted homologies, as they were called in the main text), which computes f (s, t) in one go. In a nutshell, it allows to "collapse" the information about all the branches of M 0,4 by endowing each integration cycle with a coefficient of the form e 2πi(ps+qt) for a given 6 Finding Kp,q and their respective intersection numbers Γ|Kp,q is a generalization of a similar problem considered in [121,122] for the gamma function Γ(s). As a matter of fact, we could have used the result already on the Veneziano amplitude on the right-hand side of (A.1), though this approach would not give us any intuition about generalizations to higher-point or higher-genus amplitudes that cannot be expressed in terms of gamma functions.
p, q. This allows for computations directly on M 0,4 , which are much easier than those on the covering space. In this language f (s, t) becomes a single twisted intersection number, which can be easily computed [5]. In fact, following a computation from Appendix A of [5] we find in the physical region s<0, t, u>0 f (s, t) = e −2πiα t − e 2πiα s 1 − e 2πiα s = e −2πiα t which in the second equality we rewrote as a sum over the lattice (p, q) ∈ Z 2 , as in (A.9), from which one can read-off integer coefficients of each e 2πiα (ps+qt) . Performing similar computations in other kinematic regions it is easily confirmed that in all of them, except for s, t>0, an infinite number of saddles contribute. Note that f (s, t) is in general not real, but works out so that the whole right-hand side of (A.9) remains real (i.e., compensates for the fact that log(−x) = log(x) + iπ for x>0 in the exponential). One can check that in the region t<0, s, u>0 the factor f (s, t) is obtained from (A.11) by exchanging s ↔ t, which is a consequence of crossing-symmetry of the Veneziano amplitude.
As remarked before, there are many other indirect ways of obtaining the same result. For instance, one can notice that the no matter which ray in the kinematic space we are considering, the critical point x * is always on the integration contour of some partial amplitude. For instance, in the physical region we have 1 < x * < ∞. We can then use two independent monodromy relations [123] to solve for the integral over (0, 1) in terms of that over (1, ∞). In doing so one recovers the result (A.12), where the oscillatory phases come from solving the monodromy relations. See [124] for related discussion. This approach, however, does not seem to scale well to higher-point amplitudes, as Lefschetz thimbles generically do not coincide with open-string integration cycles. The discussion of closed-string scattering is almost identical, except for an additional step at the beginning. One starts by homologically splitting the corresponding complex integral into two copies of integrals over Lefschetz thimbles. In doing so one encounters a homological intersection number of a similar oscillatory type as that in (A.11). Since Lefschetz thimbles themselves depend on the direction in which the α → ∞ is taken, so does the asymptotic behavior of amplitudes, which is generically dominated by an infinite number of saddles. We refer the reader to Appendix A of [5] for details of this computation.