Intersection Numbers from Higher-order Partial Differential Equations

We propose a new method for the evaluation of intersection numbers for twisted meromorphic $n$-forms, through Stokes' theorem in $n$ dimensions. It is based on the solution of an $n$-th order partial differential equation and on the evaluation of multivariate residues. We also present an algebraic expression for the contribution from each multivariate residue. We illustrate our approach with a number of simple examples from mathematics and physics.

A wide class of integral functions, such as Aomoto-Gel'fand integrals, Euler-Mellin integrals, Gel'fand-Kapranov-Zelevinsky integrals, to name a few, which embed and generalise Feynman integrals, can be considered as the pairing of regulated integration domains, known as twisted cycles, and of regulated forms, known as twisted cocycles [18].Within this definition, the integrand appears as the product of a multivalued function, called twist, and of a differential form.The twist carries information on the integral regularisation: for the case of dimensionally regulated Feynman integrals, the space-time dimensionality appears in the exponent of the twist.In this fashion, the algebraic properties of the integrals can be thought as coming more fundamentally from the algebraic properties of the corresponding cycles and cocycles.In particular, evaluation of intersection numbers for twisted differential forms becomes a crucial operation to derive linear and quadratic relations for integrals mentioned above, and to systematically derive differential and difference equations the latter obey [14][15][16][17][19][20][21][22][23][24].See [25][26][27][28][29][30][31], for recent reviews, and [24,32] for applications to multi-loop calculus.Applications of intersection theory and co-homology to diagrammatic coaction have been presented in [33][34][35][36], and to other interesting physical contexts in [37][38][39][40].
The evaluation of the intersection numbers for twisted forms is based on the twisted version of Stokes' theorem [41].In particular, for the case of logarithmic n-forms, intersection numbers can be computed by applying the algorithm proposed in [2] or by means of the global residue theorem [8].For generic meromorphic n-forms, the evaluation procedure can become computationally more demanding, and it can be performed by means of an iterative approach, as proposed in [16,17,42], elaborating on [10].The iterative approach has been further refined in [43], by exploiting the invariance of the intersection numbers for forms belonging to the same cohomology classes.In [22,23], this algorithm has been extended to account also for the relative cohomology cases [44], required when dealing with singularities of the integrand which are not regulated by the twist.
As an alternative to the evaluation procedure based on the Stokes' theorem, intersection numbers can also be computed by solving the secondary equation built from the Pfaffian systems [9,45,46].Within this algorithm, the determination of Pfaffian systems obeyed by the generators of the cohomology group is required, and efficient methods for their derivation have recently started to be proposed by means of Macaulay matrix in [46].
In this article, we propose a new algorithm for the computation of the intersection number of twisted n-forms, based on a novel way of applying Stokes' theorem that requires the solution of a higherorder partial differential equation and application of the multivariate residues.The computational algorithm hereby proposed can be considered as a natural extension of the univariate case [41], and, just as the latter, its application requires the solution of a (partial) differential equation around each intersection point, that belongs to the set of zeroes of the twist.In this work, we show that the solution of the differential equation can be found analytically by multiple Laurent series expansions, and that each residue admits a closed expression in terms of the Laurent coefficients of the two forms entering the pairing and of the twist.
The structure of the paper is as follows: in Section 2 we discuss aspects of twisted cohomology theory and intersection numbers; we introduce a new method for computation of multivariate intersection numbers as multivariate residues using a higher-order partial differential equation, and discus its solution locally around each intersection point.Section 3 contains application of our new approach to integrals and functions of interests for physics and mathematics.In Section 4 we give a closed, algebraic expression for each residue, contributing to the multivariate intersection number.Section 5 contains our concluding remarks.The paper includes four appendices: Appendix A contains the link of our new approach to Matsumoto's algorithm, explicitly shown in the simple case of 2-forms; Appendix B contains further details of the examples discussed in Section 3; Appendix C contains the derivation of the algebraic expression given in Section 4.
2 Intersection numbers for twisted n-forms

Twisted cohomology
Let B i , with i = 1, . . ., m , be complex homogeneous polynomials in the homogeneous coordinates Z = (Z 1 , . . ., Z n+1 ) of the complex projective space CP n .We introduce a manifold X = CP n − m i=1 S i , where the hypersurfaces S i are identified by the equations: (2.1) In the following we work in the chart Z 1 ̸ = 0 with the local coordinates (see Appendix B.1 for details) We introduce the Aomoto-Gel'fand integrals, defined as twisted period integrals, where: u is a multivalued function called the twist, which regulates the integral; Γ (n) is a regularised cycle called twisted or loaded cycle, i.e. an n-chain with empty boundary on X (usually Γ (n) is denoted as Γ (n) ≡ Γ (n) ⊗ u to separate the integration domain Γ (n) and a specific choice of the branch of multivalued u along it); φ (n) is a meromorphic differential n-form defined on X, called the twisted cocycle.In general u is a multivalued function that "vanishes" on the integration boundary: u(∂Γ (n) ) = 0.The latter property ensures that for any generic (n − 1)-form φ (n−1) the integral of the total differential is zero: where we introduced the covariant derivative: with d = n i=i d zi , where d zi = ∂ zi dz i , and ωi = u −1 ∂ zi u , using the shorthand notation ∂ zi ≡ ∂/∂z i .When dealing with individual integration variables, it might be convenient to consider the decomposition of the full covariant derivative: with the partial covariant derivatives defined as: ∇ ωi := ∇ωi dz i , and ∇ωi = ∂ zi + ωi . (2.7) Aomoto-Gel'fand integrals represent a wide class of special functions, such as Gauß hypergeometric functions, Lauricella functions, and their generalizations, Euler-type integrals, and Feynman integrals [14].Integrals of this class are invariant under shifting of the differential n-form by a covariant derivative: 1) , explicitly: (2.8) Similar results are obtained also for the so called dual integrals, obtained from the integrals (2.3) by replacing u → u −1 and ω → −ω in definition (2.5).
In the case of Feynman integrals in the Baikov representation [47], the twist u admits the following factorized form: where the exponents γ i are non-integer, and the factors B i are polynomials build out of the kinematical and Baikov variables (corresponding to propagators).For this set of functions, analyticity, unitarity, and algebraic structure are related to the geometry captured by the Morse function h := Re(log(u)), see [48,49].
The multivalued twist u carries information on the regularization: for dimensionally regulated Feynman integrals, it depends on the integration variables as well as on the external scales, such as the Mandelstam invariants and masses (all appearing in the polynomials B i ), and on the space-time dimensionality d (appearing in the γ i ).The topological information of integrals and dual integrals is contained in ω that is a differential form with zeroes and poles 1 , collected in the respective sets: Z ω = zeros of ω , and P ω = poles of ω . (2.10) The invariance of integrals and dual integrals under the transformation (2.8) can be used to expose the algebraic structure of Aomoto-Gel'fand integrals.Let us introduce two vector spaces of twisted cocycles: the twisted n-th cohomology group, which is the quotient space of closed n-forms 1) }; and the dual twisted cohomology group: (H n ω ) ∨ := H n −ω .These spaces are isomorphic, and their dimension: can be determined by counting the number of critical points of B i , namely ν = dim(Z ω ) [49], or equivalently from the Euler characteristics χ(P ω ) of the projective variety generated by the poles of ω, as ν = (−1) n (n + 1 − χ(P ω )) [16], see also [50,51], or by the Shape Lemma [17].
We denote the elements of the twisted cohomology (2.11) , and use them to define the following natural twisted Poincaré pairings: • Integrals: (2.13) • Dual integrals: • Intersection numbers for twisted cocycles: where φ L,c is a compactly supported cocycle equivalent to φ L [2] (see also Section 2.3 below and Appendix A).

Integrals and Relations
Let us now briefly review some practical applications of the twisted cohomology theory to the study of Feynman integrals, focusing on the IBP decomposition method.Consider the bases generating the cohomology groups introduced in eq.(2.11): and {|h i ⟩} i=1,...,ν ∈ H n −ω .These two bases can be used to express the identity operator in the cohomology space [14,15] as follows: where we defined the metric matrix : whose elements are intersection numbers of the twisted basis forms.Linear relations for Aomoto-Gel'fand-Feynman integrals, the differential equations, and the finite difference equation they obey are consequences of the purely algebraic application of the identity operator (2.16), see also [14].
In the context of Feynman integral calculus, the decomposition of scattering amplitudes in terms of master integrals (MIs), as well as the equations obeyed by the latter, are derived by means of IBPs [52,53] and of the Laporta method [54].These relations emerge naturally from the algebraic properties of twisted cocycles.
Indeed, generic twisted cocycles can be projected onto the bases in the corresponding vector spaces as: The latter formula is dubbed the master decomposition formula for twisted cocycles [14,15].It implies that the decomposition of any Aomoto-Gel'fand-Feynman integral as a linear combination of MIs is an algebraic operation (similarly to the decomposition/projection of any vector within a vector space), which can be carried out by computing intersection numbers of the twisted de Rham differential forms.
Using the master decomposition formula, a Feynman integral I can be decomposed in terms the MIs J i as: with the decomposition coefficients c i given by eq.(2.18).Let us remark, that the metric matrix (2.17), in general, differs from the identity matrix.The Gram-Schmidt algorithm can be employed to build orthonormal bases from generic sets of independent elements, using the intersection numbers as scalar products.But more generally the coefficients appearing in the formulas (2.18, 2.19) are independent of the respective dual elements.Therefore, exploiting this freedom in choosing the corresponding dual bases may yield striking simplifications [15,22,23].The decomposition formulas hold also in the case of the relative twisted de Rham cohomology, which allows for relaxation of the non-integer condition for the exponents γ i that appear in eq.(2.9), see [22,23,44].

Partial Differential Equation
By elaborating on the method proposed in [2], we hereby propose to evaluate the intersection number for n-forms, using the multivariate Stokes' theorem, yielding (see also Appendix A): where: • ψ is a function (0-form), that obeys the following n-th order partial differential equation (nPDE): • p = (p 1 , p 2 , . . ., p n ) ∈ P ω is a pole of ω, i.e. an intersection point of singular hypersurface S i defined in eq.(2.1), at finite location or at infinity.
• The residue symbol stands for: where the integral goes over a product of small circles ⟲ i , each encircling the corresponding pole z i = p i in the z i -plane, see [55].
Representation (2.20) can be derived by rewriting the intersection number integral as a flux of a certain local form η: (2.23) Working term-by-term in the sum on the RHS, let us temporarily denote by (z 1 , . . ., z n ) the local coordinates centered at the intersection point p.As the integration domain we take the polydisc where hi := 1 − h i and h i is the Heaviside step-function: By requiring that the auxilary 0-form ψ is the solution of the following nPDE: we obtain: where the compactly supported n-form φ L,c is defined as: The middle expression here is equivalent to the φ L,c introduced by Matsumoto in [2] and, therefore, the integration of eq.(2.23) can be carried out via iterated residues.Indeed, since φ R is a holomorphic n-form, in eq.(2.26)only the last term gives a non vanishing contribution: where the product of small circles ⟲ 1 ∧ . . .∧ ⟲ n (i.e. an n-dimensional torus) is the distinguished boundary of the polydisc D p .The last equation above2 reproduces the result shown in eq.(2.20).
For more details we refer the interested reader to the discussion in Appendix A.
Finally, let us once again highlight the crucial eq.(2.27) and write it as: This nPDE, equivalent to eq. (2.21), is the natural extension of the equation L presented in [2] for the single variable case.Equation (2.31) constitutes the first main result of this communication, as it offers a new algorithm for the direct determination of the scalar function ψ, hence a simpler strategy for the evaluation of the intersection numbers between twisted n-forms.

Solution
The solution of eqs.(2.21, 2.31) can be formally written as3 : (2.32) A crucial observation is in order.The global solution ψ is, in general, a transcendental function.The evaluation of the intersection number in (2.20), though, because of the calculation of residues, requires the knowledge of ψ only locally around each of the contributing pole, say ψ p = ψ| z→p .To determine the local expression of ψ around the point p, we propose two equivalent strategies: 1. solution by ansatz4 : whose coefficients ψ p,a can be determined by requiring the fulfillment of (2.31) around the pole p.The values of the a min and the a max depend on the Laurent expansion of φ around p. The determination of the Laurent coefficients ψ p,a can be carried out by solving the triangular systems of linear equations that arises after inserting the ansatz in the multivariate differential equation.

solution by series expansion:
which can be directly obtained by a series expansion under the integral sign, and a subsequent multifold integration.Taking p = 0 without loss of generality, we observe that as z → 0, the twist (2.9) admits a factorized expansion u| z→0 = z γ ′ • i≥0 u i z i with noninteger exponents γ ′ .A similar expansion holds for the cocycles: φL Using these expansions integration in eq.(2.34) can be done term by term via z dx x a = z a+1 /(a + 1).
Let us remark that the second strategy allows for a straightforward determination of the intersection numbers bypassing any linear system solving procedure, and it constitutes the second main result of this study.

Simplified formulas from Rescaling
In the previous sections we saw how the intersection number was given by where formally It is not hard to see that the above expressions are invariant under the following simultaneous rescalings: where q may be any5 function of z.Since each individual residue of eq.(2.35) possesses this invariance, the rescaling may be done locally at each individual p.Such rescalings are of interest since certain specific choices for q introduce simplifications.Two special choices of q are of particular interest: Right rescaling.This is defined as the choice q = φR .It corresponds to the rescaling: In this case, the nPDE eq.(2.31) becomes where ω R := d log(u R ), and the argument of the residues appearing in eq.(2.35) is directly given by f .Therefore, upon the right rescaling, it is sufficient to know the solution of the differential equation in (2.38) up to the simple pole term, because the higher order terms, from O(1) on, would not contribute to the residue.This observation turns in a computational advantage, and, for this reason, eq.(2.38), which is a special form of eq.(2.31), will become in Section 4 the starting point for deriving an algebraic expression of Res(f ).
Left rescaling This form of rescaling is defined by the choice q = 1/ φL .It corresponds to Upon this rescaling eq.(2.31) becomes and we also define ω L := d log(u L ), and where the argument of the residues appearing in eq.(2.35) becomes g φL φ R .This equation is conceptually simpler to solve than the case with a generic form on the RHS.Finally let us remark, that another source of simplifications is the shift invariance of intersection numbers: which is valid for generic (n − 1)-forms ξ L and ξ R , but discussion of this falls out of our scope.Now we move on to applications of the formalism that was introduced here.

Applications
In this section we apply the nPDE method for computing intersection numbers of twisted 2-forms (n = 2 variables).Explicit results presented below agree with the iterative method of [16,17,42].

Two-loop massless sunrise diagram
The first example is the massless sunrise diagram on the maximal cut.Using the Baikov parametrization, we denote the two irreducible scalar products as z 1 and z 2 .The twist u is built from the Baikov polynomial on the maximal cut.In the three coordinate charts U z , U x , and U y of the projective plane CP 2 (see Appendix B.1 for a brief review) it looks like this: Chart Point Coordinate transformation(s) Table 1: Intersection points P ω (middle) from the three coordinate charts of CP 2 (left) that contribute to the intersection numbers between forms (3.3) with the massless sunrise twist (3.1).On the right we show the coordinate transformations used to compute the residues (2.20).We sum over all the contributions of each displayed transformation.
where γ 0 := γ 1 + γ 2 + γ 3 , and the polynomial factors (2.9 The corresponding singular hypersurfaces S i defined in eq. ( 2.1) intersect at 3 + 2 + 1 = 6 points in CP 2 (see Figure 1a): In this example we consider intersection numbers between generic φ L and φ R cocyles: where • ∈ {L, R} , and n i ∈ Z .Our task is to compute the multivariate residues shown in the main formula (2.20) at each of the points (3.2).One way to do this is to find coordinate transformations that localize on these points P ω , and allow for direct series expansion and solution of the nPDE (2.21) either using the algebraic formula (4.8) or the integral formula (2.32) (see also an alternative algebraic formula in Appendix B.4).We collect such transformations in Table 1 and refer the interested reader to Appendix B.2 for further details.

Example 1 Consider the intersection ⟨φ
Out of the six points collected in eq.(3.2) only two give non-trivial residues, yielding: Example 2 We can also compute ⟨φ L | φ R ⟩ of meromorphic forms containing the B 3 factor: (3.6) Out of the six points shown in eq.(3.2) the three points from the U z chart give non-zero residues: where we separated contributions of the transformations shown in Table 1 into individual terms.By adding them up, the intersection number becomes:

Two-loop planar box diagram
Now we turn to the massless 2-loop planar box on the maximal cut in the Baikov representation with the twist: where γ 0 := γ 1 + γ 2 + γ 3 , and the factors (2.9) are: There are 3 + 1 + 1 = 5 intersection points in CP 2 of the singular hypersurfaces S i (see Figure 1b), their location is as follows: Here we only consider intersection numbers between monomial φ L and φ R cocycles: with the same notation as in eq.(3.3).Some of the intersection points (3.13) turn out to be degenerate: they have three singular hypersurfaces S i passing through them.An example of this is the point x = (0, 0) located on the line at infinity, as can be seen from eq. (3.11) and Figure 1b.One way to amend this issue is to employ the resolution of singularities procedure (closely related to the sector decomposition algorithm, see [56][57][58]).In Appendix B.3 we give further details and collect the full list of coordinate transformations used for computation of intersection numbers in this setup.
Example 3 Let us consider the intersection number ⟨φ L | φ R ⟩ between the two logarithmic forms: Only the origins of the three charts (3.13) contribute to this intersection number producing: . (3.16) Example 4 Similar to Example 1, we compute the intersection between: (3.17) We find the three points from P ω contribute to this intersection number, resulting in: (3.18)
Example 5 In continuation of Example 3 we calculate the intersection of the logarithmic forms: Once again, only the origins of the charts contribute to this intersection number giving the total: This concludes our collection of examples, for intersection numbers in n = 2 variables, involving Feynman integrals and other mathematical functions.Let us remark, that the presented method is also applicable to meromorphic n-forms in higher dimensions.

Algebraic solutions
In this section, we will introduce an algebraic solution for the contribution to the intersection number coming from each individual residue.The main formula of this section, eq.(4.8), is derived by solving the nPDE using an ansatz.We will start the section by introducing an essential combinatorial object, that of vector compositions.

Vector Compositions
Vector compositions are a convenient tool for manipulating algebraic expressions involving multivariate monomials, in particular when dealing with product of multivariate series.

Definition 1 Vector compositions
are a function of n non-negative integers τ 1 to τ n .Its output S is a set of ordered lists σ, each containing a number of n-vectors s, with the property s∈σ s = τ (s vectors contain non-negative integers; the zero-vector is excluded, corresponding to6 all s > 0).Each s can be represented as a step in Z n , therefore σ, being a collection of steps, can be mapped to a path in Z n , going from 0 to τ .Within this representation, S can be viewed as the set of all paths Z n , connecting 0 to τ , in one or in multiple Example 7 (n = 2 case) The two dimensional case, where the vector τ has two components, is the first non-trivial case for applying vector composition.For illustration, let us consider τ = (1, 2), yielding where the three rows correspond to entries with |σ| being 1, 2, and 3, respectively.Each element of VC(τ ) admits a representation in terms of lattice paths in Z 2 , as shown in Figure 2.
More details on vector compositions can be found in e.g.[59].In the current context, they turn out to be a useful tool to express the solution of the nPDE eq.(2.38).

Algebraic expression for residues
In Section 2 we saw that multivariate intersection numbers can be expressed as a sum over contributions from a number of points, with each contribution expressed as a residue.In the "right rescaling" framework of Section 2.5 the relation is given as where f is defined as the solution of eq.(2.38), namely and where ∇ω R,i := (∂ zi + (∂ zi w)) , with w := log(u/ φR ) .
In Section 2.4 it was discussed how to solve systems of the form of eq.(4.5) by making a series ansatz for f and solving it locally around each point, in what is essentially an analytic use of the "Frobenius method".In this section we are going to pursue that direction a step further, to obtain an algebraic solution for each term in eq.(4.4).The first step is to write w and φ as series around each point, where the appropriate variable changes have been made in order to ensure that the point is located at z = 0 and in addition that the series expansions are well defined so no further variable changes (or blowups) are needed.With this we may write where In this notation µ ϕ should be interpreted as the powers corresponding to the leading term in the expansion of φL φR .The w-expansion begins at powers 0 (lower ones cannot appear when u is of the form i B γi i , as in our case).Using an ansatz for f as a Laurent series expansion around each intersection point, the nPDE constraints the coefficients, so that the residue of f can be obtained in closed form as, where and U is a function dependent only on the terms in the expansion of w, which is conveniently defined in terms of an auxiliary function, R, as whose general definition is reported in Appendix C.1.2.Here, for illustration purposes, we showcase its expressions for a few small values of n: • n = 1 : R(α, β) = (β+v)δ α,0 + αw α (4.11) (α 1 −j 1 )j 2 w α−j w j (4.12) • n = 3 : R(α, β) = (β 1 +v 1 )(β 2 +v 2 )(β 3 +v 3 )δ α,0 For a general n-variate expression of R, see Appendix C.1.2.The residue formula in eq.(4.8), using vector compositions, constitutes the third main result of this work.Further details on its derivation can be found in Appendix C.
The arXiv version of this paper is supplemented with a Mathematica implementation of eq.(4.8), consisting of a package algebraic residue.m and a notebook file load algebraic residue.nb that contains examples of its use.

Example
Here we show an application of the algebraic formula eq.(4.8).

Example 8 Using the setup of Section 3.1, let us compute ⟨φ
In the "right rescaling" scheme, the first step is to compute Φ = φ L φR at each of the points in Table 1 and identify the leading term corresponding to the µ ϕ -vector.The result is given in Table 3.  1, and the additional labels (a,b,c) refer to the adopted coordinate transformation.
We see from the general expression (4.8) that a given point will only contribute if −µ ϕ − 2 ≥ 0. So we realize that there only will be contributions from three of the intersection points: 1, 4, and 6.In each case eq.(4.8) will contain one term only, of the form We observe that ϕ −2,−2 = 1 at each of the three points, while the values of the v i depend on the considered point, and are computed locally.The non-vanishing expressions of Res(f ), are: therefore, the intersection number is obtained by adding them up, as: which is the correct expected results.

Conclusion
In this work we proposed a new computational method for the evaluation of intersection numbers for twisted meromorphic n-forms through Stokes' theorem, which is based on the solution of a n-th order partial differential equation (nPDE).Our finding can be summarised simply as: (5.1) The evaluation of the last integral is performed by multivariate Laurent series expansions and multivariate residues.The analytic properties of the intersection numbers and of the nPDE yielded the algebraic determination of the contributing residues, which we were able to cast in closed forms for an efficient evaluation.
The presented method requires the knowledge of the intersection points as input.In the case of a twist corresponding to a hyperplane arrangement and normal crossing at the intersection points, this information can be easily extracted from the integrand.Instead, for generic configurations, when more hypersurfaces pass through a given intersection point, the desingularization procedure constitutes a challenging, hence interesting, mathematical problem.
The formulation of the intersection number for twisted cohomology presented here applies to the case of regulated singularities.We are confident that it can be extended to the relative cohomology case, where the singularities of the integral are not regulated within the twist.
The new method presented here can be applied to derive linear relations, differential equations, difference equations, and quadratic relations for Feynman integrals as well as for a wider class of functions, such as Aomoto-Gel'fand integrals, Euler-Mellin integrals, and GKZ hypergeometric functions.Therefore the method can be useful for computational (quantum) field theory, and computational differential and algebraic topology.ϵ 2 respectively, with ϵ 1 < ϵ 2 , such that V i ⊂ U i .Following [2], we then write an explicit formula for the compactly supported 2-form: where h j is a smoothened version of the Heaviside step function: h j = 1 in V j , h j = 0 outside U j , and it smoothly interpolates between the two boundary values 0 ≤ h j ≤ 1 in U j \ V j .These properties of h j imply that the two forms coincide φ L,c = φ L outside of the singular neighborhood ∪ j U j ⊃ ∪ j S j .Furthermore, it is not difficult to show that eq.(A.2) has compact support, i.e. the 2-form φ L,c vanishes on the singular neighborhood ∪ j V j ⊂ ∪ j U j , if we impose certain constraints on the auxiliary 0-and 1-forms appearing on the RHS of eq.(A.2).For our purposes here it is enough to investigate these constraints locally7 around a given intersection point p ij for some fixed i and j, where we introduce local coordinates (z 1 , z 2 ), such that the hypersurfaces S i and S j are parametrized by z 1 = 0 and z 2 = 0 respectively.Then the auxiliary 1-forms ψ i and ψ j , and the 0-form ψ (we will drop the subscripts in the following) from eq. (A.2) must satisfy the following differential equations [2]: with the z-dependence shown explicitly: Then on the polydisc U i ∩ U j eq.(A.2) reduces to: which vanishes on the inner polydisc V i ∩ V j .Similarly one can show that eq.(A.2) vanishes on the rest of the singular tube V j \ ∪ i̸ =j U i for each j, which proves that φ L,c indeed has compact support.In the integral (A.1) the compactly supported 2-form φ L,c is wedged against the holomorphic form φ R , hence out of eq.(A.2) only the anti-holomorphic part will contribute to the intersection number.The only source of this is the last term in eq.(A.5), so that the intersection number integral localizes onto neighborhoods of intersection points p ij as: where the integration domain is a difference of two polydiscs: Table 4: The three coordinate charts (left) that cover CP 2 with their local coordinates (middle).
On the right we show the coordinate transformations from the U y and U x charts back to U z .In this section we focus on the U z chart and its line at infinity, i.e. the projective line CP 1 given by Z 1 = 0.
The last factor here is problematic: its series expansion around the new origin z 1 = z 2 = 0 is ill-defined.For example, expanding 1/(z 1 + z 2 ) in the following two ways we get different results: To fix this we further change (z 1 , z 2 ) → (z 1 , z 1 z 2 ), so that the twist (B.2) becomes regular: meaning that all the polynomial factors now have constant terms.Series expansion of eq.(B.4) no longer depends on the order of operations and we may proceed with evaluation of the multivariate residue in eq.(2.20) as discussed in Section 2 and Section 3.1.Putting together the transformations that led us to the twist (B.4) gives: where we recognize one of the elements listed on the second row of Table 1.
In general, coordinate transformations such as (B.5) are also known as resolutions of singularities of algebraic varieties.They can be computed algorithmically with the help of, for example, pySecDec [56,57].

B.3 Two-loop planar box diagram
Here we provide further details for Section 3.2, in particular in Table 5 we collect the coordinate transformations needed for evaluation of intersection numbers between forms (3.14) with the twist (3.11).
Consider the point z = (−t, 0) appearing on the second row of Table 5.This point is located at the intersection of the two singular hypersurfaces S 2 ∩ S 3 , i.e. satisfies B 2 = B 3 = 0 (see eq. (3.12)).We look for a pair of local coordinates that roughly looks like (z we may keep z 2 as one of our new coordinates.To derive the transformation rule for z 1 solve: and expand this solution near z 2 = 0 producing: Table 5: Intersection points (middle) from the three coordinate charts of CP 2 (left) that contribute to the intersection numbers between forms (3.3) with the massless 2-loop planar box twist (3.11).
On the right we show the coordinate transformations used to compute the residues (2.20).Some transformations are parametrized by N ∈ Z ≥0 .We sum over all the contributions of each displayed transformation.
Truncating this expansion up to the N th order and an additional change (z 1 , z 2 ) → (z 1 , z 1 z 2 ) reproduces (up to the z 1 (z 1 z 2 ) N term) the coordinate transformation shown on the second row of Table 5.
In these new coordinates, the twist becomes: where each polynomial factor now has a constant term, so that the series expansion around z 1 = z 2 = 0 is well-defined.A generic monomial 2-form (3.14) transforms as follows: φ The pole order of the form shifts by N in both variables, which implies that for large enough N the non-vanishing condition (B.19) will no longer be satisfied and the corresponding transformation rule will not contribute to the intersection number.Hence we conclude, that even though the Table 5 collects an infinite number of coordinate transformations, only a finite number of them contributes to a given intersection number.

B.4 Series expansion of the integral formula
Here we derive the expression of a given residue appearing in the main formula (2.20) that utilises the series expansion method presented in Section 2.4.Throughout this section we assume that the pole of the residue is located at the origin z = 0. Using the multi-index notation introduced in Section 2.4, we write the z-expanded cocycles as: Examples of the domains of summation appearing here are shown in Figure 3a.In the following we will also refer to such domains of z-expansions as supports 8 .Expansion of the twist u takes the following form:  where u 0 is a z-independent constant, the vector γ is a linear combination of exponents appearing in (2.9), ũ collects terms of the expansion that have at least one power of z inside, and we denote: We show the supports of u and ũ in Figure 4a and Figure 4b respectively.In the following we implicitly assume an additional condition i ≥ j whenever we write (i − j) • 1 ≥ m for some vector shift j.
According to the integral formula (2.32), our goal is to compute the following residue: The key idea is to expand the inverse of the twist u(z) −1 in terms of ũ(z) using the geometric series9 : where in ρ m we gather the residues with m insertions of ũ(z).The summation bounds in the second term follows from the fact that the z-expansion of ũ(z) m has terms with at least m powers of z inside, as reflected in Figure 4b and Figure 4c for m = 1 and m = 2 respectively.The residue ρ 0 in eq.(B.14) is straightforward to evaluate: where the fraction stands for 1/(i + γ) := 1/ j (i j + γ j ) , and the n-dimensional summation domain is shown in Figure 3b.
To compute ρ m in eq.(B.14) we separately expand its three factors, multiply them, and take the residue.The first factor is given in eq.(B.10).The second factor is a power of a series, so we may write its i th non-zero coefficient ũ(z) m = i (ũ m ) i z i using the Kronecker symbol δ i,j as: One can come up with various evaluation strategies for this multiple sum using, for example, the nested parametrization as in eq.(C.32), or the (precomputed) multivariate Bell polynomials [60].For the i th non-zero coefficient of the last factor of ρ m we get  Finally the sum representation of the ρ m residue reads: where examples of the outer sum's domain are shown in figures 3c and 3d (see also Figure 4d for an illustration of the ≤ notation in the inner sum).Alongside the non-vanishing condition: we obtain the algebraic expression for a given residue: This corresponds to where we have used the obvious notation ∂ i := ∂ zi .Inserting the expansions we get Changing the summation variables we may rewrite this as where 12∆ s := which reduces to the expression given in eq.(4.12) in the case where µ w = 0.In general this function R may be interpreted as R(α, β) being the coefficient in ∆ of ψ β z α+β−1 .
Having ∆ = 0 in general requires that each term in the z-expansion of eq.(C.9) is zero independently.This means that our task now is to solve ∆ s = 0 for each value of s.
This can be done.If we define μψ := µ L − nµ w + 1 (C.12) we find the recursive solution where and where the exact definitions of ≥ and other binary operators applied to the index vectors are given is Section C.4.The quantity μψ may be interpreted as the "correct" value for µ ψ , or at least as a value for which µ ψ ≰ μψ would make it impossible for the above derivation to go through.The recursive solution given by eq.(C.13) can be useful on its own to find ψ-solutions to use as arguments of the residue function using for instance eq.(C.5), but let us continue the derivation in order to find a closed expression.
Evaluation of the residue.The next step is to re-express eq.(C.13) in a form in which only the coefficients of the individual φL,i are recursive.This can be done as where W (η, λ) := Q(η+λ, η) or correspondingly Eq. (C.17) is a recursive expression that only involves the coefficients of the expansion of w.The recursion can be solved with the result where VC refers to the vector compositions defined and discussed in Section 4.1.By combining eqs.(C.15) and (C.5), and changing the summation variables we get and finally we may insert the expression for Ψ from eq. (C.19) and we get the most general form of the algebraic expression: where We notice that τ is not really needed as an argument to Y since τ = j σ j .
which after a similar derivation yields which reduces to eq. (4.13) in the case of µ w = 0. Going through the derivation for various values of n allows us to realize the pattern and find the general expression for R which is given in the following.

General expression for R
For R in the n-variate case, we may write where R n,k contains the terms with k w-factors and correspondingly k−1 sums.We find it convenient to use distinct expressions for R n,k for k = 0, k = 1, and13 k ≥ 2. Defining we have We should point out that it is possible to write eq.(C.32) nicer and more symmetrically at the cost of introducing an extra sum and a delta function, that is14 where evaluating the i 1 -sum would bring back eq.(C.32).

C.2 Algebraic expression in the rescaling frameworks
In Section 2.5 it was discussed how to simplify the computation of the intersection number by rescaling the quantities on which it depends, locally at each intersection point.Two different cases were presented, right rescaling and left rescaling, which we will discuss below.

Right rescaling
Right rescaling is defined by the rescalings and renamings given in eq.(2.37) Eq. (C.21) may then be used for these rescaled variables, with φL,i , φR,i , w i , and v j interpreted as the coefficients in the expansions of the rescaled variables φ, 1, and log(u R ), and µ L , µ R , µ w as the corresponding powers of the leading coefficients.But the expansion of 1 is of course trivial, corresponding to φR,i = δ i,0 and µ R = 0.This allows us to do the h-sum of eq.(C.

Left rescaling
One may go through the same steps for the left rescaling framework, which is defined in terms of the renamings and rescalings defined in eq.This fact is useful for implementation purposes, since it tells which coefficients to extract before applying the algebraic expression.Applying one of the rescaled versions of the algebraic expression does not change these limits if expressed in terms of the original φ L and φ R .

The univariate case
In the univariate case there are two simplifications taking place which deserve discussion.The first is that the vector compositions VC reduce to standard compositions as discussed in Section 4.1.The second is that only one ω is present which makes it desirable to express the result in terms of expansions of that object, replacing the w-expansion of eq.(C.4) with its separately treated log-terms.That is Combining this with eq.(2.20) we get an expression for the intersection number similar to an expression from refs.[2,8].Additionally it is worth noticing that that the intersection number is linear in both φ L and φ R which may be seen from the integral definition in eq.(2.15).That property is easy to see in the algebraic expression as it is given by eq.(C.21).But in the rescaling framework linearity becomes hidden.For instance looking at the solution in the right-rescaling framework as it is given by eq.(4.8), linearity in φ L is obvious.Linearity in φ R , however, becomes obscured as the φ R,i -terms mix with the w i -factors, which enter the expressions in a very non-linear fashion.

C.4 Comparing index vectors
In this appendix and in Section 4 we used the notation for the n-variate formula eq.(4.8)where the sums and products are parameterized by index vectors written in bold, of length n and containing integers.Such vectors can be compared using the so-called product order, let us briefly review how it is done.If we let a and b be two such index vectors of length n, the relations we use in this paper are  6, and some are illustrated in Figure 5.We see that ≰ and > are not equivalent, for instance it is true that 1  2 ≰ 2 0 while it is false that 1 2 > 2 0 .We notice that the index vectors form a partially ordered set under the ≤ operator.

2 Figure 1 :
Figure 1: Singular hypersurfaces S i of the massless sunrise twist (3.1) on the left, massless 2-loop planar box twist (3.11) in the middle, and the 3 F 2 twist (3.19) on the right.The red sphere (with identified antipodal points) depicts the real slice of the projective plane CP 2 , whose equator is the line at infinity.The blue dots represent the intersection points collected in eqs.(3.2, 3.13, 3.20).

Figure 4 :
Figure 4: Supports of the expansions (B.11) in (a, b, c), where we implicitely assume an extra i ≥ 0 condition.In (d) we illustrate the notation of eq.(B.18) (an additional i ≤ 0 here is implied as well).

17
) and i ≥ µ L + 1. Integrals here are done via the same technique as discussed in Section 2.4.

47 )Figure 5 :
Figure5: Illustration of some of the binary operators defined in Table6.The shaded areas cover the integer values for which the illustrated relation is true.

Table 6 :
a = b ∀i: a i = b i a ≤ b ∀i: a i ≤ b i a ≥ b ∀i: a i ≥ b i a < b a ≤ b and a ̸ = b a ̸ = b not(a = b) a ≰ b not(a ≤ b) a ≱ b not(a ≥ b) a > b a ≥ b and a ̸ = bThe definitions needed to compare index vectors.defined as in Table .25) so that the differential dh i is localized on the circle |z i | = ϵ .The action of the partial derivatives in eq.(2.23) gives: d z1 . . .d zn η = h1 . . .hn u ∇ ω1 . . .∇ ωn ψ + . . .+ (−1) n

Table 2 :
Intersection points P ω (middle) from the three coordinate charts of CP 2 (left) that contribute to the intersection numbers between forms (3.21) with the 3 F 2 twist (3.19).On the right we show the coordinate transformations used to compute the residues (2.20).We sum over all the contributions of each displayed transformation.

Table 3 :
Values for Φ and µ ϕ at the twelve sub-intersection points for Example 8.The points are defined in Table