The SAGEX Review on Scattering Amplitudes, Chapter 3: Mathematical structures in Feynman integrals

Dimensionally-regulated Feynman integrals are a cornerstone of all perturbative computations in quantum field theory. They are known to exhibit a rich mathematical structure, which has led to the development of powerful new techniques for their computation. We review some of the most recent advances in our understanding of the analytic structure of multiloop Feynman integrals in dimensional regularisation. In particular, we give an overview of modern approaches to computing Feynman integrals using differential equations, and we discuss some of the properties of the functions that appear in the solutions. We then review how dimensional regularisation has a natural mathematical interpretation in terms of the theory of twisted cohomology groups, and how many of the well-known ideas about Feynman integrals arise naturally in this context. This is Chapter 3 of a series of review articles on scattering amplitudes, of which Chapter 0 [arXiv:2203.13011] presents an overview and Chapter 4 [arXiv:2203.13015] contains closely related topics.

Here γ E = −Γ ′ (1) is the Euler-Mascheroni constant, ν = (ν 1 , . . . , ν p ) ∈ Z p , and we work in dimensional regularisation with D = D 0 − 2ǫ dimensions and D 0 ∈ N. Depending on the context, we will use D or ǫ interchangeably to denote quantities that depend on the dimensional regulator. The factors 1/(m 2 j − q 2 j − iε) are called propagators, and we use the usual Feynman iε-prescription to deform the integration contour away from the propagator poles. Hence note the distinction between the dimensional regularisation parameter ǫ and the infinitesimal ε for the Feynman prescription for the propagator. The integral is a function of the propagator masses m 2 1 , . . . , m 2 p and of the external momenta p 1 , . . . , p E , constrained by momentum conservation E j=1 p j = 0 (we assume without loss of generality that all external momenta are incoming). The squared propagator masses are assumed to be positive, and the external momenta are real Minkowski momenta which we also assume to be D-dimensional (we assume that, aside from momentum conservation, they are all linearly independent). The numerator N is a polynomial in the dot products involving at least one of the loop momenta k i . The momenta q i flowing through the propagators have the form We define |ν| as the sum of the exponents of the denominators, The Feynman integral in eq.
As a consequence, the integral only depends on the external scales We will often write I(x; ν; D) instead of I(p 1 , . . . , p E ; m 2 1 , . . . , m 2 p ; ν; D), in order to make the Lorentz invariance manifest. We note that because of momentum conservation not all {p i · p j } 1≤i,j≤E are independent, so it should be understood that x contains only an independent set of such dot products. Remark 1. It is possible to relax the condition ν i ∈ Z. This leads to the notion of generic Feynman integrals considered in ref. [3,4]. In particular, for certain classes of Feynman integrals, an L-loop Feynman integral can be written as an (L−1)-loop generic Feynman integral.

Some properties of Feynman integrals
Feynman graphs. Feynman integrals as presented in eq. (1) are associated to Feynman graphs with L loops, where for every propagator there is an internal edge e j labelled by (q j , m 2 j , ν j ), and momentum must be conserved at every vertex. The external edges are labelled by the inflowing external momenta.
Mass dimension of a Feynman integral. Consider the action on Feynman integrals of a rescaling of the external momenta and the propagator masses, (p j , m j ) → (λ p j , λ m j ), λ ∈ R * = R \ {0}. It is straightforward to check that all external scales in eq. (5) rescale like x → λ 2 x. Assuming that the numerator is homogeneous (as is always the case in physics applications), N (λ 2 {k j · k l , k j · p l }; D) = λ α N N ({k j · k l , k j · p l }; D) , (6) where α N = [N ] is the mass dimension of the numerator, we have I(λ 2 x; ν; D) = λ α I I(x; ν; D) , where the mass dimension of the integral is Note that [I(x; ν; D)] = m − 2Lǫ for some integer m ∈ Z, which implies that the mass dimension of an integral is never zero in dimensional regularisation.
Proof. Since there is no external scale, the integral is a constant. On the other hand, we must have I(λ 2 x; ν; D) = λ α I I(x; ν; D) , for all λ ∈ R * .
Proof. If ν j ≤ 0 for all 1 ≤ j ≤ p, then the integrand is a polynomial. We can then write the integral as a linear combination of scaleless integrals.

Parametric representations
The representation of Feynman integrals in eq. (1) is the one that is the most straightforwardly associated with Feynman graphs and Feynman rules. However, it is not the most convenient to evaluate the integrals and it obscures some of their properties. In this subsection we present alternative representations that can be helpful in addressing these questions, mostly focusing on the Feynman integrals with unit numerator, N = 1. The parametric representations of the integrals that will be discussed here can also be used as the definition of a Feynman integral in dimensional regularisation, where D is simply another parameter the integrals depend on, instead of having the physical interpretation of being the non-integer dimensional space in which the momenta live. We will not discuss all existing parametric representations for Feynman integrals, but instead restrict ourselves to the ones that will feature in the rest of this review. We refer the reader to other references for a more detailed discussion, e.g., refs. [5][6][7].
Schwinger parametrisation. The Schwinger-parameter representation is obtained by using: 1 Here, the α j are called Schwinger parameters, and U(α) and F (α; x) are polynomials that are associated with the corresponding Feynman graph (i.e., they can be determined directly from the Feynman graph associated with the Feynman integral). Their precise form is not important for us here, but we highlight some of their properties that will be used later: • U(α) is the determinant of an L × L matrix whose entries only depend on the Schwinger parameters α j . It is an homogeneous polynomial of degree L in the α j .
• F (α; x) depends on both the Schwinger parameters α j and the external scales x.
It is an homogeneous polynomial of degree L + 1 in the α j of mass dimension 2. It has the form: where F(α; x) is independent of all internal masses m 2 j . The polynomials U and F are also known as the first and second Symanzik polynomials.
Feynman parametrisation. This is perhaps the most well known parametric representation, and can be derived either directly from eq. (1), or from the Schwingerparameter representation. Starting from eq. (14), we insert and after some manipulation we obtain where the α j are now called Feynman parameters. The Feynman-parameter representation in eq. (17) is not as general as it could be. Indeed, it can be shown that it is in fact a projective integral over a simplex in real projective space of dimension p − 1, and eq. (17) is only a particular realisation of such an integral. Other realisations can, for instance, be obtained by restricting the sum in the delta function to be over only a subset of the Feynman parameters (usually referred to as the Cheng-Wu theorem in physics [8]), or even any other linear combination thereof [6].
The Feynman-parameter representation has been widely used to compute Feynman integrals by direct integration, both numerically and analytically. It also allows us to make an observation on Feynman integrals that is completely obscured in the momentum-space representation of eq. (1) and not as clear in the Schwinger parameter representation of eq. (14): we see that the powers of the denominators ν i and the dimension D appear in a very similar way in the integrand of eq. (17), namely they are both in the exponents of the U(α) and F (α; x) polynomials. We will see later that we can find relations between integrals with different values of the ν i , but also with different values of D.
Cutkosky-Baikov parametrisation. The final parametric representation we discuss is the Cutkosky-Baikov representation. We start by noting that the integrand of the Feynman integral in eq. (1) depends on K = L + E − 1 momenta: the L loop momenta k 1 , . . . , k L and the E − 1 independent external momenta, say p 1 , . . . , p E−1 . Let us denote these momenta as t 1 , . . . , t K , in the order specified above. We then consider all dot products τ ij = t i · t j for 1 ≤ i ≤ L and 1 ≤ j ≤ K that involve at least one loop momentum. It is straightforward to check that there are N = L(L + 1)/2 + L(E − 1) such dot products. We can change variables from the components of the loop momenta to the τ ij to find: where G(a 1 , . . . , a n ) is the Gram determinant of the momenta a 1 through a n . The integration domain ∆ is bounded by the surface G(t 1 , . . . , t K ) = 0. We next note that the inverse propagators z a = (m 2 a − q 2 a ) are linear in the τ ij , see eq. (2). Provided there is an invertible transformation from the z a to the τ ij , we can then trivially change variables from the τ ij to the z a . It is very common that there are fewer propagators than the number N of τ ij . This situation can be dealt with simply by considering a larger class of integrals with extra propagators, so that the transformation from the z a to the τ ij is invertible (this defines a complete family of Feynman integrals, see section 2.2 below), and then set the power ν i associated with the extra propagators to zero. Assuming there is an invertible transformation between the z a and the τ ij , we can write: where B(z) is the Baikov polynomial associated with the Feynman integral I(x; ν; D) (strictly speaking, with the complete family of Feynman integrals defined by I(x; ν; D)): Remark 2. (Cut Feynman Integrals) The Cutkosky-Baikov representation in eq. (19) is often not the most practical for computing I(x; ν; D). However, this representation is well suited for studying certain properties of Feynman integrals. For instance, it is a natural starting point for defining and studying so-called cut Feynman integrals, where a subset of the inverse propagators are set to zero. (The representation in eq. (19) was in fact first introduced in ref. [9] precisely to study the Landau singularities of Feynman integrals [10], which are closely related to their cuts; see also refs. [11][12][13] for discussions about computing cut integrals from the Cutkosky-Baikov representation.) In the representation of eq. (19), this cut condition can be imposed by taking residues where the corresponding z a are zero. We note that there is a close connection between evaluating residues and choosing a contour that encircles the associated poles, so that cut Feynman integrals correspond to Feynman integrals evaluated on modified integration contours (see e.g. ref. [14] for a discussion in the context of one-loop integrals). Of particular interest are contours that encircle the poles associated with all active propagators (those whose ν i are positive), corresponding to a maximal cut of the integral. Beyond one loop, maximal cuts might not be unique, and there may be inequivalent choices of how to integrate over the variables that are not localised by the maximal-cut conditions.
Given the relation between cut Feynman integrals and residues, it follows that: If ν i ≤ 0 for any of the cut propagators of a cut Feynman integral, then the cut integral vanishes.

Linear relations among Feynman integrals
Instead of seeing the Feynman integral as a function of the external scales x for fixed values of the propagator exponents ν and the space-time dimension D, we can also interpret it as a function ofν = (ν 0 , ν 1 , . . . , ν p ) ∈ Z p+1 , with ν 0 = D 0 /2, for fixed x. In other words, if we fix the external scales x, we can associate to every pointν on the lattice Z p+1 the Feynman integral I(x; ν; D). In this way we obtain an infinite family of Feynman integrals labelled by the lattice points. Note that Corollary 1 implies that the integral vanishes unless at least one of the ν 1 , . . . , ν p is strictly positive. The point of this section is to show that these integrals are not independent, but there are linear relations among integrals for different lattice points, and one can always identify a finite basis of integrals that generates this infinite family. We start by discussing the linear relations relating integrals with different values of ν ∈ Z p but for the same space-time dimension ν 0 , and we comment on relations between integrals in different dimensions at the end of this section.

Total derivatives in dimensional regularisation
We start by presenting a theorem first used in refs. [15,16] to study linear relations among integrals. In order to state the theorem, it is useful to introduce the notation with Proposition 3. In dimensional regularisation, we have for every D-dimensional vector v µ .
Proof. We follow the argument given in ref. [17] and use the simplified notation We expect the integral to be invariant under general linear changes of variables, k µ i → λ k µ i + v µ , with λ a non-zero real number and v µ a D-dimensional vector independent of k µ : • Under an infinitesimal translation, k µ i → k µ i + εv µ , we have Hence, we must have: • Under an infinitesimal rescaling, k µ i → e ε k µ i , we have Hence, we must have Remark 3. We have stated Proposition 3 for a specific loop momentum k i , but it clearly also holds for any of the other loop momenta. Note that the regularisation scheme plays an important role, and Proposition 3 may be false in another scheme. For example, if we use a cut-off as a regulator, and Proposition 3 does not necessarily hold as it might need to be corrected by boundary terms.

Integration-by-parts relations
We now show that the vanishing of the surface terms in Proposition 3 implies that Feynman integrals satisfy linear recursion relations in the propagator exponents ν. We first note that, when the differential operator ∂ ∂k µ i v µ acts on the propagator (m 2 j − q 2 j ) −ν j , it shifts the value of the exponent. For example, if v = k i and q j = k i + p 1 , we have: We see that we recover the same propagator, with the power raised by one unit. However, we have also introduced a numerator, which may not be present in our original Feynman integral. If v is chosen to be a loop or external momentum (or any linear combination thereof), then the numerator is a polynomial in the dot products involving loop and/or external momenta. We may then express these dot products as inverse propagators. For example, if we choose v = p 1 in the above example, we have We see that we have indeed obtained a linear combination of propagators raised to different powers (we assume without loss of generality that there is a propagator with momentum q 1 = k i and mass m 1 ). In general, the coefficients of the linear combination are polynomial functions in the external scales (and the dimensional regulator ǫ).
Numerator factors involving loop momenta appear as propagators raised to negative integer powers. However, we had to assume that we can write all scalar products in terms of inverse propagators. This may not always be the case, and we say that a family of Feynman integrals is complete if it contains enough propagators to express all dot products involving at least one loop momentum in terms of inverse propagators (we also use the word topology to refer to a complete family of Feynman integrals, while elsewhere in the literature it is sometimes used for families that may not be complete). Note that this is the same assumption we made in section 1.2 when discussing the Cutkosky-Baikov representation and, as we argued there, every family of Feynman integrals can be completed by adding enough propagators. We therefore always assume from now on that our families of Feynman integrals are complete. This discussion can then be summarised as follows: Proposition 4. Every complete family of Feynman integrals (that is, every topology) satisfies linear recursion relations in the propagator exponents ν ∈ Z p , called integrationby-parts (IBP) relations [15,16]. The coefficients of the linear combinations are rational functions in the external scales x and the dimensional regulator ǫ.
IBP relations allow one to relate different Feynman integrals from the same family to each other, and therefore to reduce considerably the number of Feynman integrals to be computed. These relations involve only rational functions of the scales and ǫ. It is possible to solve the recursion relations and to express every member of a given (complete) family in terms of a basis of integrals. Elements of such a basis are conventionally called master integrals.
Proposition 5. The number of master integrals is always finite.
The proof of this result can be found in refs. [18,19]. It is also possible to predict the dimension of the basis (up to some caveats) via the number of critical points [20] or a certain Euler characteristic [19]. We stress nevertheless that there is considerable freedom in how one can choose a basis of integrals for the complete family (the dimension of the basis, however, is fixed). In fact, different choices may lead to improved methods for their evaluation, as will be discussed in section 3.
To close this subsection, let us introduce some concepts that will be useful later on. We define a map where θ(x) denotes the Heaviside step function: We say that two Feynman integrals I(x; ν; D) and I(x; ν ′ ; D) belong to the same sector if ϑ(ν) = ϑ(ν ′ ). Roughly speaking, a sector is the collection of all Feynman integrals of a family that share the same set of active propagators (where we qualify a propagator as active if it is raised to a strictly positive power ν i , i.e., θ(ν i ) = 1). Integrals from the same sector may, however, differ by the choice of the numerator factors. There is a natural partial order on sectors, and we say that If we derive IBP relations starting from Proposition 3 for a given integral I(x; ν; D), then these relations will only involve integrals from the sector ϑ(ν) or from lower sectors. It can also happen that a Feynman integral can be expressed as a linear combination of integrals from lower sectors only. We call such an integral reducible. If all integrals from a given sector are reducible, we call the sector reducible. ‡ Remark 4. In applications it is often not necessary to solve the IBP relations in all sectors. Indeed, one often encounters the situation that certain propagators only enter with negative exponents (e.g., because they were added to obtain a complete family of integrals). In such a scenario it is then only necessary to solve the IBP relations in the subsectors that describe the active propagators needed for the application one has in mind (with a caveat discussed in section 2.6).
Remark 5. IBP relations for cut Feynman integrals (see Remark 2) follow from those for uncut integrals. In an IBP relation, the 'cut' only acts on Feynman integrals and not on the rational functions. According to Proposition 2, some of the Feynman integrals in the IBP relation might be set to 0. In particular, this implies that the maximal cut of a reducible integral always vanishes.
Remark 6 (Lorentz-invariance identities). At the heart of the existence of linear relations among Feynman integrals is the fact that first-order differential operators act via shifting the propagator exponents. This is of course not restricted to derivatives with respect to loop momenta, as considered in Proposition 3, but it is easy to see that the same conclusion holds for derivatives with respect to external momenta (or even propagator masses). In particular, this means that the fact that the generator of a Lorentz transformation annihilates a Feynman integral (cf. eq. (9)) translates also into a linear relation among Feynman integrals with shifted exponents. These so-called Lorentz-invariance identities are not independent from the IBP relations generated by Proposition 3 [17].

Solving IBP relations
In all but the simplest cases, it is not known how to find a solution to the IBP relations in closed form. That is, we do not know how to find an expression for an integral I(x; ν; D) for generic ν in terms of a basis of master integrals. In practical applications, however, we are only interested in solving the relations for a finite set of ν (e.g., the ones appearing in the calculation of an amplitude). As noted in ref. [21], instead of solving a recursion problem, this can be formulated as a linear algebra problem: one explicitly writes down all IBP relations satisfying some criterion (typically one puts a bound on the sum of the negative and positive exponents ν i ; we note that care must be taken when choosing this bound, as if it is too constraining one might miss some relations and find a number of master integrals that is larger than it should be), which includes all the integrals one wishes to rewrite, and then solves for the 'complicated' integrals in terms of the 'simple' integrals. The notion of 'complicated' and 'simple' integrals requires introducing an ordering criterion, and it is natural to favour integrals with fewer propagators (see footnote 2). The approach to solving the IBP relations proposed in ref. [21], commonly called the Laporta algorithm, has since been implemented in a number of public programs such as AIR [22], REDUZE2 [23], LiteRed [24], FIRE [25], and Kira [26] (we only refer to the latest releases of each code).
For integrals with several loops and many scales, solving the IBP relations is still a very challenging problem and often poses a bottleneck in their evaluation. The main difficulty comes from the fact that the number of IBP relations to solve grows very fast, and the complexity of the rational functions involved can become unmanageable when the number of scales increases. Many approaches have been proposed to overcome these challenges, starting with the implementation of sophisticated algorithms to solve such systems of equations in the programs mentioned above. Aside from these algorithms, there have also been new developments which are more specifically targeted to IBP relations. First, one can try to construct simpler IBP relations by choosing vectors v µ in eq. (23) with certain properties. For instance, one might choose v µ so that the IBP relations do not involve integrals with propagators raised to higher powers (this would be natural for reducing amplitudes to a basis of master integrals, since in that case one is in general trying to reduce numerator factors) [27]. It is not trivial to find such vectors v µ , but it is a question that can be answered with tools from computational algebraic geometry, namely by solving syzygy equations [28][29][30][31], which have had a strong impact in modern approaches to the calculation of scattering amplitudes [32][33][34][35]. Second, one might try to choose a basis for which the coefficients in the IBP relations will be simpler. It has been observed that working with a so-called canonical basis (see Proposition 10 below) can be very beneficial as the singularity structure of the coefficients greatly simplifies [36][37][38], and in particular the dependence on D and on x completely factorises. Combined with multivariate partial-fractioning techniques [37,[39][40][41][42], one obtains much more compact solutions to the IBP relations. Finally, and closely connected to the point above about the simplicity of the coefficients, one can solve the IBP relations for numerical values of x and D, and use this numerical data to determine the analytic expressions [43][44][45][46][47]. This approach is particularly powerful when combined with the second point mentioned above: using numerical evaluations bypasses the intermediate analytic complexity generated in the solution of the IBP system, and one directly reconstructs the much simpler rational expressions appearing in the solution of the IBP relations. The numerical evaluations are most commonly done in finite-field arithmetic [43,44], for which there exist very efficient linear algebra algorithms and where all numerical calculations are exact, thus removing any questions about numerical precision and stability.
Let us conclude by mentioning that an alternative approach is being developed to achieve the reduction to master integrals without the need to solve the IBP relations. We will review this new approach in section 5. Currently, however, this new approach is still in its infancy, and solving the IBP relations using the techniques described in this section remains the method of choice for all applications.

Example: The one-loop bubble integral
As an example of the application of IBP relations, we will consider the one-loop bubble integral: Note that any polynomial N ({k 2 , p 1 · k}; D) can always be written as a polynomial in the propagators and in the variables in x = (p 2 , m 2 1 , m 2 2 ), so without loss of generality we consider the integral with a trivial numerator.
The IBP relations can be constructed from with v µ = k µ and v µ = p µ . For v µ = k µ , we obtain and for v µ = p µ we obtain where we simplified the notation to only keep the dependence of the integrals on the ν i , i.e., we set I(ν 1 , ν 2 ) = I(p 2 ; m 2 1 , m 2 p ; ν 1 , ν 2 ; D). To simplify the discussion, let us first consider the case where m 2 1 = m 2 2 = 0. Before exploring the IBP relations, we note that in this limit I(ν 1 , ν 2 ) = 0 if either ν 1 ≤ 0 or ν 2 ≤ 0, as under those conditions the integrals become scaleless. The IBP relations above can be rewritten as: This means that any integral I(ν 1 , ν 2 ) is either 0 or can be related to I(1, 1), which is trivial to compute (e.g., using Feynman parameters): We then find that this family of Feynman integrals contains a single master integral. We also note that for this particular case it is not hard to compute the Feynman integral for arbitrary powers of the propagators: From this representation, we recognize the IBP relations in eq. (38) as consequences of the well-known recurrence relations between Γ-functions, Γ(1 + z) = z Γ(z). If instead m 2 1 = 0 and m 2 2 = 0, then I(ν 1 , ν 2 ) = 0 if ν 1 ≤ 0. The IBP relations for this case are obtained from eqs. (36) and (37) evaluated at m 2 2 = 0, and we find that there are two master integrals, which we can choose to be I(1, 0) and I(1, 1), to which any I(ν 1 , ν 2 ) can be related. The tadpole integral I(1, 0) is while the bubble integral I(1, 1) is where 2 F 1 is Gauss' hypergeometric function: Just as for the bubble with massless propagators, the IBP relations for this case follow from the recurrence relations (known in this case as contiguous relations) satisfied by the 2 F 1 hypergeometric function [48].

Dimension-shift relations
So far we have only considered linear relations connecting Feynman integrals with different propagator exponents ν, but with the same values for the space-time dimension D and external scales x. The Feynman-parameter representation of Feynman integrals in eq. (17), however, makes it manifest that there is no substantial difference between the dimension D and the exponents ν i : they both appear as exponents of the polynomials in the integrand. More precisely, the space-time dimension enters the exponents of the Feynman parameter integral only through the combination ν 0 = D 0 /2. It is therefore natural to expect that there are linear relations relating Feynman integrals with different values of ν 0 . This was worked out in detail in refs. [49,50].
Proposition 6. For Feynman integrals depending on generic and non-zero propagator masses, we have: Before we give a derivation of this relation, let us make some comments. The differential operator in the right-hand side involves the U-polynomial that appears in the Feynman and Schwinger parametrisations, but with the Feynman/Schwinger parameters α i replaced by the differential operators ∂ ∂m 2 i . For this reason we need to consider integrals with generic propagator masses. The action of this differential operator produces on the right-hand side a linear combination of Feynman integrals in D dimensions with shifted propagator exponents. Proposition 6 asserts that this linear combination equals the Feynman integral in (D − 2) dimensions (up to an overall factor). Moreover, once all derivatives have been carried out, we can set the masses to non-generic (possibly zero) values, and we obtain a relation between integrals in different dimensions also for non-generic masses.
Proof. We start from the Schwinger parametrisation in eq. (14). The dependence of the integrand in eq. (14) on the space-time dimension and the masses is particularly simple: the space-time dimension only enters through the exponent of the U-polynomial, and the masses only appear in the exponent in the integrand, see in particular eq. (15). We then have We see that the application of the differential operator amounts to multiplication by U(α) in the integrand, changing the exponent from −ν 0 to −(ν 0 − 1).

Proposition 7.
An integral in D + 2 dimensions can be written as a linear combination of integrals in D dimensions as: where (x) L is the Pochhammer symbol, B is the Baikov polynomial defined in eq. (20), and the b i are operators that lower the value of the exponent ν i , that is Proof. We start from the Cutkosky-Baikov parametrisation in eq. (19), and note that the space-time dimension appears in a very simple way. In particular, in the integrand it only appears in the exponent of the Baikov polynomial. The dimension-shifting relation of eq. (48) then follows simply by noting that a polynomial on the inverse-propagator variables in the numerator gives a linear combination of integrals with shifted powers of the propagators.
where we used the same compact notation as in eqs. (36) and (37), but this time keeping the dependence on D explicit. The factor of λ(p 2 , m 2 1 , m 2 2 ) is in this case a singularity related to the integral on the left-hand side; however, IBP relations can generically introduce spurious poles as well. We can also relate the (D + 2)-dimensional bubble to the master integrals in D dimensions with eq. (48). Using the same basis as above, we find: The dimension-shift relations for the cases where m 2 1 = 0 and/or m 2 2 = 0 are obtained from the above by simply setting the masses and the scaleless integrals to zero. As for the IBP relations, the dimension-shift relations correspond to recurrence relations of the special functions the integrals evaluate to, either gamma functions or (generalised) hypergeometric functions (see eq. (40) for an explicit example). These recurrence relations are now with respect to the parameter ν 0 = D/2, but, as already pointed out, there is no substantial difference between the exponents ν and ν 0 .

Other relations among Feynman integrals
Let us conclude this section with some comments about to what extent the IBP and dimensional-shift relations capture all relations between Feynman integrals. Since IBP and dimensional-shift relations are linear, we need to discuss linear and non-linear relations separately.
Let us start by discussing linear relations among Feynman integrals. Currently, there is no indication of homogeneous linear relations among Feynman integrals in the same family, in dimensional regularisation, that do not follow from the IBP and dimensional-shift relations, and conjecturally we have Conjecture 1. All linear relations among Feynman integrals from a given complete family follow from IBP and dimension-shift relations.
As an example, we already pointed out in Remark 6 that the Lorentz-invariance identities follow from the IBP relations. In applications, however, there may be various caveats: • We have already mentioned that it is useful to separate Feynman integrals into sectors, and in concrete applications one only needs to consider those sectors where a certain subset of propagators is active. Consequently, one only needs to solve the IBP identities in those sectors. It can happen, however, that certain relations among integrals from lower sectors are only found if IBP relations involving integrals in higher sectors are considered. Hence, if IBP relations from higher sectors are neglected, it may appear that there are relations among Feynman integrals in lower sectors that seem not to follow from IBP relations (while in fact they do follow from IBP relations in higher sectors).
• It can also happen that certain new relations arise in limits where the external scales take degenerate values. Those additional relations arise from IBP relations once the degeneracy among the scales is resolved, cf., e.g., ref. [51,52]. Similar types of relations were derived in ref. [53] from functional equations, although the examples shown there relate integrals of different families and so fall outside the scope of the conjecture above.
• IBP relations detect linear relations of Feynman integrals in generic space-time dimensions D. In applications one is usually interested in external momenta that lie in 4 space-time dimensions. Since at most 4 vectors can be linearly independent in 4 dimensions, it may happen that certain combinations of scales vanish when the external momenta are chosen four-dimensional. This may lead to new relations for four-dimensional external momenta, which were not present for D-dimensional external momenta. We note that, due to momentum conservation, this situation arises for the first time for 6 external momenta in 4 space-time dimensions.
Finally, let us comment on non-linear relations among Feynman integrals. Much less is known about such relations, though over the last couple of years several instances of nontrivial quadratic relations between Feynman integrals (or their maximal cuts) have been discovered [54][55][56][57][58][59]. By nontrivial, we mean both that it is not possible to obtain these relations as a consequence of linear relations, and that some of the loop integrations do not trivially factorise. As we will mention in section 5, the existence of nontrivial quadratic relations seems to be very general, and it would be interesting to explore them more generally in the future. Currently, there is no example of nontrivial relations of higher degree (cubic or higher).

The method of differential equations
The IBP relations reviewed in the previous section allow one to reduce the problem of computing a given family of Feynman integrals to the computation of a (finite) set of basis integrals, usually called master integrals. The master integrals must be evaluated by other means, and there are various well-established methods to compute them (cf., e.g., ref. [5] and references therein for a review). Some of these methods rely on direct integration of parametric representations such as the ones in eqs. (14) and (17), using various methods to perform the integrals (see, e.g., refs. [60][61][62][63][64][65][66][67][68]), but over the last decades the method of differential equations [69][70][71][72][73][74] has established itself as one of the most powerful. There are several good reviews and lecture series on how to use differential equations to compute Feynman integrals, see, e.g., refs. [75,76]. Here, we attempt to give an overview of the general strategy in a broader mathematical context.

Differential equations satisfied by Feynman integrals
In the following we denote by T a vector of basis integrals depending on the scales x = (x 1 , . . . , x s ), and we assume that the entries of I(x, ǫ) are ordered in a way that is compatible with the natural order on the sectors. Let us consider the derivative of I(x, ǫ) with respect to an external scale x i . We can exchange the derivative ∂ x i = ∂ ∂x i and the loop integration, and we act with the derivative on the loop integrand. If x i is a propagator mass, x i = m 2 j , the action on the integrand is easy to compute. If x i is a scalar product between external momenta, x i = p j · p k , then we can use the chain rule to express ∂ x i in terms of the differential operators p µ j ∂ ∂p µ k , whose action on the loop integrand is straightforward and very similar to the calculation done in eq. (31), see also Remark 6. The expression one obtains after acting with the differential operators can then be rewritten in terms of master integrals using the IBP identities described in section 2. In summary, the derivative of master integrals with respect to an external scale x i can be expressed as a linear combination of master integrals. That is, we can write: where the A x i (x, ǫ) are N × N matrices. Since IBP relations involve only rational coefficients, the entries of A x i (x, ǫ) are rational functions in x and ǫ. Moreover, if the entries of I(x, ǫ) are ordered such that they respect the natural partial order on the sectors (as we are assuming), then the matrices are block upper-triangular.
Example 2 (The differential equations for the one-loop bubble). Let us return once more to the example of the one-loop bubble integral with two massive propagators. We already established in section 2.4 that there are three master integrals associated with this topology. We set where we use the propagator powers in eq. (34) to distinguish the three different master integrals. From now on we will drop the dependence on all other quantities. The derivatives with respect to p 2 , m 2 1 and m 2 2 are given by where λ(a, b, c) was defined below eq. (45). In the equations above, we first present the action of the differential operator on the integrand of eq. (34), and then what one obtains after using IBP relations to rewrite the expressions in terms of master integrals. Note in particular that ∂ p 2 I(0, 1) is only explicitly zero after accounting for IBP relations.
It is convenient to package the different partial derivatives into a total differential with respect to all external scales, where is a matrix whose entries are rational one-forms. Remark 7. The matrices A x i are in fact not independent. Indeed, the total differential must satisfy d 2 = 0, and we have It follows that A(x, ǫ) must satisfy the integrability condition where '∧' denotes the wedge product between differential forms. Equation (58) gives a set of differential relations between the matrices A x i (x, ǫ) that can serve as a useful check of the correctness of the differential equations. It is straightforward to verify, for instance, that the matrices in eqs. (53) to (55) satisfy these relations.
Remark 8. The differential operators ∂ x i are not all independent, even if the scales x i are. Indeed, for every integral I(x, ν, ǫ) (basis integral or not) we have the relation: where α I = I(x, ν, ǫ) was defined in eq. (7). To see why eq. (59) holds, we note that the differential operator on the left-hand side is the infinitesimal generator of the dilatations where D was defined in eq. (11). The eigenvalues of the dilatation operator are the mass dimensions (the additional factor of 1/2 comes from the fact that α I is the mass dimension of the integral, and the scales x i themselves have mass dimension equal to 2). Equation (59) implies that the nontrivial functional dependence of I(x, ν, ǫ) can only be in the ratios so that eq. (59) implies It is therefore sufficient to consider the derivatives with respect to the ratios y i , 1 ≤ i ≤ s − 1, rather than in the individual scales x i (or, equivalently, we may put x s = 1 during the computation). We note that this implies that the differential equation satisfied by a one-scale integral (s = 1) is trivial, and the method of differential equations is not suitable for computing such integrals. These integrals must then either be computed by other means, e.g., direct integration or the dimensional recurrence and analyticity method [50,77]. If one still wishes to use differential equations, it is often possible to introduce an additional scale x 2 . The method of differential equations may then be applied to the ratio y 1 = x 1 /x 2 , and we recover the original integral in the limit y 1 → ∞. We refer to ref. [78] for details. Equation (59) leads to the following differential relation between the master integrals where I(x, ǫ) = diag I 1 (x, ǫ) , . . . , I N (x, ǫ) is the diagonal matrix whose entries are the mass dimensions of the basis elements. Equation (63) provides another nontrivial check of the correctness of the differential equations. For instance, using the matrices in eqs. (53) to (55), we verify that in the case where I denotes the master integrals for the two-mass bubble discussed in section 2.4.
3.1.1. Change of basis The differential equation (56) can be very hard to solve in closed form, including the full dependence on the dimensional regulator ǫ. In applications we are only interested in the Laurent expansion of I(x, ǫ) up to some finite order in ǫ, and this order is typically relatively low. The basis I(x, ǫ) is not unique, and we may use this freedom to change basis to a new basis J (x, ǫ) which brings the differential equation into a form that can be more easily solved. We are therefore interested in determining how the differential equation (56) behaves under a change of basis.
Since I(x, ǫ) and J (x, ǫ) are bases for the same vector space, there must be an invertible matrix M(x, ǫ) such that The new basis J (x, ǫ) satisfies the differential equation where the matrix A ′ (x, ǫ) is related to the matrix A(x, ǫ) from eq. (56) by At this point we have to make an important comment: we have already mentioned that, due to the rational nature of the IBP relations, the original matrix A(x, ǫ) in eq. (56) is a matrix of rational one-forms. However, we have not specified what the functional dependence of the matrix M(x, ǫ) is, and no-one forces us to choose M(x, ǫ) to have rational entries. On the other hand, if the entries of M(x, ǫ) are not rational, the new basis J (x, ǫ) will not consist of Feynman integrals of the form I(x, ν, ǫ) multiplied by rational functions of x and ǫ. In the following, we will refer to a basis related to one consisting of integrals of the form I(x, ν, ǫ) via a rational transformation matrix M(x, ǫ) as an IBP-basis. We will only consider changes of bases M(x, ǫ) that depend rationally on ǫ, and we will call the transformation in eq.
We will later see that it is particularly convenient to consider the change of basis with λ(a, b, c) as defined below eq. (45). This is an example of an algebraic transformation.
For the case where one of the propagators is massless, say m 2 2 = 0, the basis is two-dimensional. In section 2.4 we chose the basis For this example, the change of basis equivalent to that of eq. (69) is This is an example of a rational transformation. In both cases above, we can compute the matrices A ′ (x, D) associated with the bases J (x; D) using eq. (67), and we find that That is, the D-dependence completely factorises from A ′ (x, D), and for D = 2 − 2ǫ the matrix is proportional to ǫ. As we will later see, such bases have particularly interesting properties.

Solving the differential equations
In this section we discuss a strategy for solving the differential equations satisfied by the master integrals (for an alternative strategy in the case of one-variables problems, see, e.g., ref. [79]). We focus here on analytic approaches, though we point out that it is also possible to solve the differential equations numerically, cf., e.g., refs. [80][81][82][83][84][85][86][87][88].
We are interested in the first few terms in the Laurent series around ǫ = 0 of I(x, ǫ). Typically, some basis integrals will have poles at ǫ = 0, but it is generally easy to determine the first non-zero order in the Laurent expansion. More generally, we have: Proposition 8. For every IBP-basis, there is a rational transformation to an IBP-basis whose elements are finite and non-zero at ǫ = 0.
A basis with the property that every basis element is finite and non-zero at ǫ = 0 is called ǫ-regular [89] (see also ref. [90,91] for the closely related notions of ǫ-finite and quasi-finite basis). The proof of this statement is constructive and provides an algorithm to determine the ǫ-regular basis [89]. Note that, since the transformation obtained from this algorithm is rational, the ǫ-regular basis will also be an IBP-basis. In the following we therefore assume without loss of generality that the starting IBP-basis I(x, ǫ) is ǫregular. This implies that the matrix A(x, ǫ) in the differential equation (56) is finite for ǫ = 0, and we define Let us now return to the discussion of how to solve he differential equation. We recall that we assumed that the ordering among the basis elements in I(x, ǫ) is compatible with the natural order on the sectors. We can then write The θ i ∈ {0, 1} p label the N sec irreducible sectors, and the irreducible basis integrals in the sector θ i are collected in . The irreducible basis integrals in I θ i (x, ǫ) have all the propagators of the sector (ϑ(ν i j ) = θ i ), and no linear combination of them can be reduced to a lower sector via IBP relations. Since the matrix A(x, ǫ) is block upper-triangular, we can split the homogeneous system of differential equations in eq. (56) into N sec inhomogeneous systems for the basis in a given sector: where N θ i (x, ǫ) collects contributions from lower sectors. Written in this form, the equations can be solved sector by sector. The procedure involves two steps, which we describe separately.
The associated homogeneous equation. We start by considering the homogenous equation for ǫ = 0 associated to eq. (75): Since this is a linear system of N θ i homogeneous first-order differential equations, the general solution takes the form: where the c k are unknown complex constants and the I (k) θ i ,h (x) form a basis for the solution space of eq. (76). We can then form a matrix where the columns are these vectors, called the Wronskian matrix of eq. (76) (sometimes also called the fundamental solution matrix ): The general solution in eq. (77) can then be cast in the form where c = (c 1 , . . . , c N θ i ) T is determined by the initial condition, and the Wronskian matrix satisfies the equation We see that the entire information about the general solution to the homogeneous equation (76) is encoded into the Wronskian matrix, which can be interpreted as a matrix solution to the homogeneous equation at ǫ = 0. Since the columns of W θ i (x) form a basis for the solution space, the determinant of W θ i (x) must be non-zero (at least for generic values of x). The determinant satisfies the differential equation Although the determinant is often an algebraic function, the individual entries of W θ i (x), are in general not algebraic, but rather transcendental functions. Determining the entries of the Wronskian matrix is often a very complicated task and no closed form for the solutions is known, especially if N θ i > 2. In that case it can be useful to turn the system of N θ i first-order equations for the vector I θ i ,h (x) into a system of differential equations of higher order for an individual entry of I θ i ,h (x). It may then be possible to obtain local power-series solutions from the Frobenius method using standard techniques, which may be continued to multi-valued solutions for all values of x. For recent applications in this direction in the context of Feynman integrals, see refs. [59,[92][93][94]. In the following we will assume that we know the expression for W θ i (x), which is typically the case for N θ i ≤ 2.
Remark 9. The entries of the Wronskian matrix W θ i (x) have an interpretation in terms of maximal cuts. Indeed, the cuts of a Feynman integral satisfy the same differential equations as the full, uncut, integral, where we have to put to zero all cut integrals where we would need to cut propagators raised to non-positive powers (cf., e.g., refs. [95,96]). It then follows that if we cut all propagators in the sector θ i , the resulting maximal cuts of I θ i (x, ǫ), denoted I max-cut θ i (x, ǫ), satisfy the homogeneous equation associated to eq. (75) [11-13, 97, 98]: Note that eq. (82) holds in D = D 0 − 2ǫ dimensions, not just for ǫ = 0.
The number of independent maximal cuts in the sector θ i always agrees with the number N θ i of irreducible basis elements in that sector. We can then interpret the Wronskian matrix W θ i (x) as the matrix of all independent maximal cuts at ǫ = 0 in the sector θ i . This implies that we can also obtain the entries of the W θ i (x) by evaluating maximal cuts [11-13, 97, 98]. It is in fact possible to say a bit more about the homogeneous equation in eq. (82), at least at a conjectural level. In ref. [58] it was conjectured that it is always possible to find an algebraic transformation M(x, ǫ) such that A θ i (x, ǫ) has the form for every irreducible sector θ i .

Iterative solution of the inhomogeneous equation.
We have established that we can associate the Wronskian matrix W θ i (x) to equation (75). Let us now consider the change of basis is invertible, so this is a valid change of basis). We obtain: with Note that, since the entries of W θ i (x) will generally be transcendental functions, this change of basis will in general be a transcendental transformation, and the new basis will no longer be an IBP-basis. Since our basis is ǫ-regular, all quantities admit a Taylor series expansion close to ǫ = 0: Equation (84) can now easily be solved order by order in ǫ. At O(ǫ 0 ) we find: This equation is solved by quadrature: where J (0) θ i ,x 0 is a constant vector of initial conditions and the integration runs over a (fixed) path from an initial point x 0 to a generic point x. Note that the integrability condition in eq. (58) ensures that the solution only depends on the choice of the end points x 0 and x, but it does not depend on the details of the path. We will discuss how to fix the initial conditions in section 3.3.1. At O(ǫ 1 ) we find the equation: The solution is We can of course continue in this way and write down the solution at every order in ǫ.
We see that at every order we have one extra integration compared to the previous order, which naturally leads us to a class of functions called iterated integrals. They can be defined as follows: consider a path γ : [t 0 , t] → X in some space X, and ω 1 , . . . ω n differential one-forms on X. We define their iterated integral along γ by [99] γ ω 1 · · · ω n := where f i : C → C are complex functions which are the pull-backs of the ω i along γ, γ * ω i = dt i f i (t i ). We will study iterated integrals and their properties in detail in section 4. Our arguments lead to the following conclusion: Proposition 9. At every order in ǫ, the basis integrals I(x, ǫ) involve rational functions in x, the maximal cuts for ǫ = 0, and iterated integrals of these functions.

Differential equations in canonical form
The strategy spelled out in the previous section is very general, but it is often not the most efficient approach to find an explicit expression for the Feynman integrals. It has been suggested in ref. [73] that there is a class of distinguished bases in which the differential equations take a particularly simple form. The existence of this form for the differential equation is still conjectural, though the conjecture is supported by all existing computations of multi-loop Feynman integrals. The idea is to find a transformation matrix M(x, ǫ) such that the matrix A ′ (x, ǫ) in eq. (67) is as simple as possible. The main conjecture is the following, first formulated in ref. [73] (in the slightly restricted setting discussed below): The initial condition J θi,x0 can be thought of as the value of J θi (x) is singular, in which case the situation requires special attention, because in that case the value at x = x 0 , as well as the integral in eq. (88), are ill-defined. We will discuss this in more detail in section 4.3. For now, it suffices to assume that J Conjecture 2. For every IBP-basis I(x, ǫ) satisfying the differential equation there is a (possibly transcendental) transformation to a new basis where A(x) is a matrix of one-forms with at most logarithmic singularities.
If the new basis satisfies eq. (93) we say that the differential equations are in canonical form, and we call the basis a canonical basis. Here, a differential one-form with logarithmic singularities should be thought of as a differential one-form ω i such that, in some appropriate choice of local coordinates ξ = (ξ 1 , . . . , ξ s ), all singularities take the form where p is an algebraic function of ξ and the dots denote power-suppressed terms that are regular at p(ξ) = 0.
The simplest examples of canonical differential equations are then those with where the A i are constant matrices and p i (x) are algebraic functions of the external scales. We refer to this case as a canonical dlog-form. This is the setup originally discussed in ref. [73]. It is known that there are Feynman integrals for which one cannot obtain a system of differential equations in canonical dlog-form as in eq. (95). To see this, we start by arguing that, if there is a canonical dlog-form, it can necessarily be reached by an algebraic transformation. Indeed, there must be a transformation M(x, ǫ) such that We know that A(x, ǫ) is a matrix of rational one-forms. The matrix A(x) in eq. (95) is a matrix of algebraic one-forms, i.e., its entries are of the form s k=1 dx k a k (x), with the a k (x) algebraic functions (simply take a k (x) = ∂ x p(x)/p(x)). This shows that M(x, ǫ) must itself be algebraic (unless there are some really unnatural cancellations). Since we know that the solutions to the homogeneous differential equations for the maximal cuts (cf. eq. (82)) generically involve transcendental functions already for ǫ = 0, it follows that the differential equations cannot always be cast in the form (95). A necessary condition for a canonical dlog-form to exist is thus that the maximal cuts for ǫ = 0 (which are the solutions to the homogeneous equations (82)) are all algebraic functions.
Remark 10. It is an interesting question if the previous condition is also sufficient and not only necessary. For all known examples for which the homogeneous equations can be solved in terms of algebraic functions, a canonical dlog-form also exists. From a purely mathematical standpoint, however, this is not automatic, because the matrix A(x) may involve one-forms without any singularities. It would be interesting to understand if, for Feynman integrals, the fact that all maximal cuts for ǫ = 0 are algebraic is equivalent to the existence of a canonical dlog-form. 3.3.1. Solving canonical differential equations One of the main advantages when working with differential equations in canonical form is that they are particularly easy to solve as a Laurent expansion in ǫ. To see this, consider the differential equation and assume that we know the value of J 0 (ǫ) = J (x 0 , ǫ) at some point x = x 0 . For simplicity, we assume that J (x, ǫ) is regular at x = x 0 . We will return later to the situation where J (x, ǫ) develops singularities at this point. In order to get the value at some other point x, consider a path γ from x 0 to x. The value of J (x, ǫ) is then obtained by parallel transporting the solution from x 0 to x along γ: where W(γ, ǫ) is given by the path-ordered exponential: Note that W(γ, ǫ) is actually a Wronskian matrix for the differential equation (that is, a matrix built out of a basis of general solutions to the differential equation like W θ i (x) in eq. (78) ; unlike W θ i (x), however, W(γ, ǫ) solves the full differential equation and not just the homogeneous part). In particular, it satisfies eq. (97): It is straightforward to obtain the ǫ-expansion of the path-ordered exponential. As an example, consider the case where A(x) = A 1 ω 1 + A 2 ω 2 , with ω i one-forms and A i constant matrices. We have: We see that at every order in ǫ the Laurent coefficients are iterated integrals, in agreement with the observation made in Proposition 9.
We can actually say a bit more if we work with a canonical basis: since the matrix A(x) has at most logarithmic singularities, the iterated integrals arising from the pathordered exponential W(γ, ǫ) will diverge at most logarithmically, and in particular they have no poles or power-like singularities. More precisely, if the differential one-forms ω i diverge at most like in eq. (94), then their iterated integrals will behave in the limit p(ξ) → 0 like γ ω 1 · · · ω n ∼ n k=0 c k log k p(ξ) + · · · , where the c k are constants, and the dots indicate contributions that are powersuppressed in the limit. Iterated integrals of this type are often called pure functions in the physics literature [105]. We thus conclude Proposition 10. Every order in the ǫ-expansion of a canonical basis involves iterated integrals that are pure functions.
Let us point out a major difference between Propositions 9 and 10. The fact that only pure functions arise in the ǫ-expansion is tightly linked to a canonical basis, whose existence is still conjectural in the general case. Proposition 9 is not conjectural, but it follows from direct considerations of how to solve the system of differential equations. By itself, however, Proposition 9 does not allow us to conclude that the iterated integrals are pure.
Remark 11. If a system of differential equations is in canonical form, then the integrability condition in eq. (58) takes the form: Since this equation must hold for all values of ǫ, the integrability condition takes a simpler form for systems in canonical from: In particular, this shows that for a system of differential equations in canonical form, the matrix A(x) must be a matrix of closed one-forms. This is manifestly the case for systems in canonical dlog-form (because all dlog-forms are closed). It was also observed to hold for non-dlog cases, see refs. [64,106,107].
Remark 12. If the matrix A(x) satisfies the condition (104), the value of the pathordered exponential does not depend on the choice of the path γ (note that this is in fact true for any basis: the integrability condition of eq. (58) is sufficient to guarantee path-independence of the solutions). More precisely, if γ and γ ′ are two paths that can be continuously deformed into one another without crossing any singularities, then W(γ, ǫ) = W(γ ′ , ǫ). We may then choose a path which makes the iterated integrals easy to evaluate in terms of known classes of special functions. In particular, we may write γ as a composition of segments where all but one variable are constant. We refer to such a path as piecewise-constant. On each segment the iterated integrals may be evaluated in terms of iterated integrals in one variable, which are easier to evaluate. The value of W(γ, ǫ) can be recovered as follows: If γ 1 and γ 2 are two paths such that the end-point of γ 2 coincides with the initial point of γ 1 , we denote by γ 1 γ 2 the path obtained by first traversing γ 2 and then γ 1 . We then have the relation: Let us conclude this discussion by briefly mentioning how we can fix the initial conditions. Clearly, one way is to know the value of the integral in one point. It is often useful to consider the value of the integral at some singular point, for instance where some scales vanish and the integral is expected to simplify. If the differential equation is in canonical form, then all singularities are logarithmic. For example, close to the singularity described by eq. (94), the differential equation takes the form where A p is a constant matrix and the dots indicate terms that are power-suppressed in the limit p(ξ) → 0. The solution to this equation is simply where in the last step we have inserted the Jordan decomposition for the matrix A p = RJR −1 and J p (ǫ) is related to the value of J (ξ, ǫ) in the limit p(ξ) → 0. The matrix J is block diagonal: where the λ i are the eigenvalues of the matrix A p and r i is the size of the i th Jordan block. The matrix exponential is then easy to carry out for each Jordan block J i : We see that in a canonical basis we can easily predict the leading-power logarithmic behaviour close to the singularities. The asymptotic behaviour of Feynman integrals can also be analysed with independent techniques, like Mellin-Barnes integrals or the Strategy of Regions [108][109][110][111]. By matching the two perspectives, we can fix the initial condition at the singular point. We will not say more about this approach here, and refer instead to the literature (cf., e.g., refs. [78,112]), because in many cases the initial condition can be fixed purely from physics input, without the need to explicitly evaluate the integrals in a limit. The idea is that the differential equation exposes all the singularities that its solutions may have. The Feynman integrals, however, are very special solutions, and the structure of their singularities is very constrained from physical principles such as unitarity. In particular, a Feynman integral can only be singular if we go to a kinematic configuration where intermediate particles go on-shell. This implies that several singularities, present in the differential equation and therefore also in the general solution, must be absent in the particular solution we are looking for. It has been observed in many cases that by imposing this condition we are in fact able to determine all the initial conditions (up to some simple one-scale integrals, that are simply a choice of normalisation of the solution).
Example 5 (One-loop bubble with a single massive propagator). As an example of how to solve a differential equation in canonical form, we return to the example of the oneloop bubble integral. For simplicity, we fix m 2 2 = 0 and we will consider the canonical dlog basis J (p 2 , m 2 1 ; ǫ) obtained by setting D = 2 − 2ǫ in eq. (71). We start by changing variables from (p 2 , m 2 1 ) to (u = p 2 /m 2 1 , m 2 1 ) and defining: The differential equation satisfied by J (u; ǫ) is found to be where we have made it explicit that it is in canonical dlog form. Since the singlescale tadpole integral only depends on m 2 1 , this differential equation does not constrain J 2 (u; ǫ). This integral is considered trivial compared to the bubble, and was already quoted in eq. (41). Adjusting the normalisation according to the basis in eq. (71), it is: Let us now focus our attention on solving the differential equation (111) order by order in ǫ, and let J i,u 0 (ǫ) be the initial condition for J i (u; ǫ) at u = u 0 . Equation (112) then implies J 2,u 0 (ǫ) = e γ E ǫ Γ(1 + ǫ), for all values of u 0 . We define The first two orders of J 1 (u; ǫ) are: Let us now note that eq. (111) is singular at u = 0, corresponding to p 2 = 0, hence any general solution to eq. (111) should have a logarithmic divergence at that point. This is clearly the case in eq. (114). The one-loop bubble, however, should only have singularities at the thresholds m 2 1 = 0 and p 2 = m 2 1 (see, e.g., ref. [113] for a textbook discussion of the singularities of Feynman integrals, and refs. [114,115] for a more modern perspective that applies to Feynman integrals with massive propagators). It is clear that this implies that J 2,u 0 , that is the leading term of the bubble integral should be equal to the leading term of the tadpole integral, which we have already computed above. It is in fact particularly convenient to choose u 0 = 0, and we will now use the strategy outlined above to determine the boundary condition to all orders in ǫ (since this is a singular point of the differential equation, we refer the reader to section 4.3, in particular eq. (138), for a discussion on the regularisation of the integrals). In the neighbourhood of u = 0 we can expand the differential equation as in eq. (106): The solution is of the form of eq. (107), that is where The requirement that J (u; ǫ) should be regular at u = 0 order by order in the ǫ-expansion implies The solution for the bubble integral is then: After accounting for the different normalisation, this expression is found to agree with the expansion of the all-order solution in eq. (42).

Finding a canonical form
The discussion of the previous sections makes it clear that it is of great advantage to bring a system of differential equations into canonical form. The existence of a canonical basis is however conjectural, and there is no general algorithm to find a transformation that brings a given system into canonical form. In fact, it is only known how to find a canonical basis for an arbitrary number of external legs and for arbitrary propagator masses in the case of one-loop integrals [100][101][102][103][104].
Over the last couple of years substantial progress has been made in understanding how to find a canonical dlog-form, if it exists. We have already argued that a necessary condition for a canonical dlog-form to exist is that all maximal cuts at ǫ = 0 evaluate to algebraic functions. This is usually equivalent to requiring that the loop integrand can be completely localised via a sequence of residues. The algebraic functions obtained in this way are closely related to the maximal cuts of the integral at ǫ = 0, and they are often referred to as leading singularities [116]. It was conjectured in ref. [105] that an integral evaluates to a pure function if an only if it can be normalised in such a way that all its non-vanishing leading singularities are ±1. This observation was at the heart of the original conjecture of ref. [73], which states that a basis is canonical if and only if it has unit leading singularities. A very important strategy for finding a canonical dlog-form therefore consists in finding a basis of master integrals with unit leading singularity. Closely related to integrals with unit leading singularity are Feynman integrals whose integrals can be written in terms of dlog-forms with algebraic arguments, cf. ref. [117]. Since leading singularities are also closely related to maximal cuts, which themselves solve the associated homogeneous system of differential equations (cf. section 3.2), many strategies to finding a canonical basis start by solving these homogeneous equations (cf., e.g., [118]). There are several public computer codes that can be used to find a canonical basis, for example Fuchsia [119], Canonica [120], Libra [121] and Initial [122], all of which are tailored to slightly different situations. A method to find a canonical basis based on the computation of intersection numbers (see section 5) has also been proposed [104,123]. For large systems of integrals depending on many scales, we find that a combination of all approaches is often necessary to find a canonical basis, see e.g. refs. [124][125][126][127][128].
In cases where the maximal cuts cannot be evaluated in terms of algebraic functions for ǫ = 0, no canonical dlog-form can exist. While several examples exist where it was possible to find a canonical form nonetheless [106,107,[129][130][131], there is no systematic understanding of how to do so. In the future, it would be interesting to clarify the connection between generalisations of the concept of leading singularity to such cases [132] and the existence of the canonical basis.

Iterated integrals
It follows from Propositions 9 and 10 that the natural class of functions that arise from the Laurent-expansion of dimensionally-regulated Feynman integrals are iterated integrals. In this section we present this class of functions in detail.

General definitions
In the following we consider a geometric space X, and we always fix a set of local coordinates ξ = (ξ 1 , . . . , ξ s ). If γ is a path in X, and ω 1 , . . . , ω n are one-forms on X, we have already defined the iterated integral of ω 1 · · · ω n along γ in eq. (91). In the following it will be useful to refer to the one-forms ω i as letters and to ω 1 · · · ω n as a word of length n. The set of all (independent) letters is called the alphabet. It will be useful to consider linear combinations of words (unless stated otherwise, we always consider linear combinations with rational numbers as coefficients), and integration is linear: By convention, it is useful to define the integral of the empty word to be γ () = 1. Iterated integrals satisfy several general properties [99]: (i) Shuffle product: where in the right-hand side we have introduced the shuffle product on words, defined recursively by and (ω 1 · · · ω n ) ¡ () = () ¡ (ω 1 · · · ω n ) = (ω 1 · · · ω n ).

Homotopy invariance
We have seen in section 3 that whenever the matrix A satisfies the integrability condition in eq. (58), then the solution does not depend on the details of the path γ, but only on its endpoints (as long as we do not cross any singularity). As a consequence, the iterated integrals that arise in the ǫ-expansion should have the same property. The iterated integrals defined in section 4.1, however, will in general depend on the details of the path γ, and not just on its endpoints. For example, consider the case X = C 2 with coordinates ξ = (ξ 1 , ξ 2 ) and ω i = dlog ξ i . We consider the line segments with 0 ≤ t ≤ 1. Then γ 12 = γ 1 γ 2 and γ 34 = γ 3 γ 4 are two paths from the point (1, 1) to the point (ξ 10 , ξ 20 ). These two paths clearly have the same endpoints, but the iterated integrals of the word ω 1 ω 2 depend on the details of the paths. Indeed, using the path composition formula (123), one finds: ω 1 ω 2 = log ξ 10 log ξ 20 and If instead we consider the linear combination of words w = ω 1 ω 2 + ω 2 ω 1 = ω 1 ¡ω 2 , then we have This illustrates the fact that the ǫ-expansion of the path-ordered exponential can only furnish special linear combinations of words such that the iterated integrals only depend on the endpoints of the path, provided that the differential-equation matrix satisfies the integrability condition in eq. (58). In the remainder of this section we describe a necessary and sufficient condition that these special linear combinations need to satisfy. Two paths γ 1 and γ 2 that can be continuously deformed into each other while keeping the endpoints fixed are called homotopic. Homotopy defines an equivalence relation on paths, and the equivalence classes are called homotopy classes. An (iterated) integral that does not depend on the details of the path, but only on its endpoints, is called homotopy invariant. The iterated integrals that arise from the ǫ-expansion of the path-ordered exponential are homotopy invariant, provided that the differentialequation matrix satisfies the integrability condition (58).
In order to understand the criteria for homotopy invariance, let us start by understanding the case of length n = 1. Consider a one-form ω 1 and two homotopic paths γ 1 and γ 2 with the same endpoints. It is easy to see that the integrals of ω 1 along the two paths agree if and only if It is clear that γ 1 γ −1 2 is a closed curve. Consider now a surface D whose boundary is ∂D = γ 1 γ −1 2 . Stokes' Theorem implies If we want this identity to hold for all paths, we see that the integral of ω 1 along a path is homotopy invariant if and only if ω 1 is closed, dω 1 = 0. Hence, for n = 1, the criterion for homotopy invariance reduces to requiring that the form is closed, i.e., it is annihilated by the total differential. This criterion can in fact be generalised to higher length [99]: if w is a linear combination of words, then the iterated integrals of w are homotopy invariant, and thus independent of the details of the path, if and only if Dw = 0, where D is the differential that acts on words of one-forms ω 1 · · · ω n via A linear combination of words w that is annihilated by this differential is called integrable. We note that this is just another incarnation of the integrability condition of eq. (58). This discussion leads us to the following conclusion [99]:

Proposition 11. An iterated integral is homotopy invariant if and only if the linear combination of words is integrable.
Remark 13. In the case where all one-forms in an integrable word w are dlog-forms, ω i = dlog a i (ξ) for some algebraic functions a i (ξ), then this integrable word can be identified with the symbol of the iterated integral γ w, which has prominently appeared in the context of Feynman integrals and scattering amplitudes, cf. refs. [99,[133][134][135][136].
Remark 14. If all the one-forms ω i are closed, dω i = 0, then eq. (130) takes the simpler form: This is in particular the case when all one-forms are dlog-form, in which case the condition Dw = 0 reduces to the well-known integrability condition for the symbol w, see for example ref. [114]. More generally, this is the case for iterated integrals that arise from differential equations in canonical from, cf. eq. (104). A special case is when X is one-dimensional and all one-forms ω i are holomorphic. Then not only are all one-forms on X automatically closed, but we also have ω i ∧ ω i+1 = 0. Hence, if X is one-dimensional, all words are integrable, and so all iterated integrals are homotopy invariant.

Regularisation
We have argued that in applications it might be useful to choose a singular point to fix the initial condition of a system of (canonical) differential equations for Feynman integrals. However, if an endpoint of the path γ is a singular point, then the iterated integrals arising from the expansion of the path-ordered exponential will typically be divergent. In this scenario, we need to replace the integrals that arise from expanding the path-ordered exponentials by a suitably regularised version. The regularisation should satisfy certain conditions. For example, it should agree with the naive definition of iterated integrals in section 4.1 whenever the integral converges, and all algebraic properties (like the shuffle product, path composition formula, etc.) should apply also to the regularised version. In the remainder of this section we present such a regularisation, sometimes called shuffle-regularisation or tangential base-point regularisation (cf., e.g., ref. [137]). We restrict the discussion to the case where all singularities are logarithmic (which is automatically the case for systems of differential equations in canonical form) and where the space X is one-dimensional. This case is typically sufficient for applications, as we may always choose a path that is piecewise constant, and on each segment we can reduce the problem to a one-dimensional case. We therefore assume from now on, without loss of generality, that the path γ goes from the origin ξ = 0 to the point ξ = x, and some of the one-forms ω i have a logarithmic singularity at the origin, We also assume that there is no other singularity on the integration contour. The regularised value of x 0 ω 1 · · · ω n = γ ω 1 · · · ω n is then defined as follows: (i) We introduce a small cut-off ε, i.e., we replace the integral x 0 ω 1 · · · ω n by x ε ω 1 · · · ω n . This integral is convergent for all ε = 0.
(ii) Since all singularities are logarithmic, the integral behaves in the limit ε → 0 as (iii) The regularised value of x 0 ω 1 · · · ω n is then defined to be the constant term I 0 (x). It is easy to check that this regularisation satisfies all our requirements: • If the original integral is convergent, then it agrees with the regularised version: x 0 ω 1 · · · ω n = lim ε→0 x ε ω 1 · · · ω n = I 0 (x) .
• The regularisation is consistent with the shuffle product, etc. We can perform all algebraic manipulations in a naive way, without having to worry about the regularisation. Indeed, x ε ω 1 · · · ω n is convergent for ε = 0, and so all naive manipulations apply. Moreover, the projection onto the constant term I 0 (x) is consistent with multiplication (the constant term of a product is the product of the constant terms), and so the regularisation and multiplication operations commute.
Example 6. Let us illustrate this on the example of the integral x 0 dξ ξ . Clearly, this integral is divergent. We introduce a cut-off, and consider the integral There is one important issue we need to address: The regularised value obtained in this way depends on a 'choice of regularisation scheme'. Indeed, we are free to rescale the cut-off by some non zero constant, ε → v ε, and the regularised value will depend on the choice of v. This rescaling factor can be interpreted as the angle and speed of approach to the singularity. It is therefore natural to identify v with a tangent vector of γ at the point x 0 = 0 (called a tangential base-point), and we will denote it by v x 0 = v 0 . The regularised value x 0 ω 1 · · · ω n with the tangential base-point v 0 at the origin is then denoted by x v 0 ω 1 · · · ω n . Note that, if x 0 ω 1 · · · ω n converges, its value is independent of the choice of tangent vector.
Example 7. Let us return to Example 6. If we choose as a cut-off v ε, we have Hence, the shuffle-regularised value with respect to the tangent vector v 0 is Remark 15. The previous discussion implies that whenever we want to fix the initial condition of a system of differential equations at a singular point, then we have to interpret the iterated integrals that arise from the ǫ-expansion of the path-ordered exponential as their regularised versions. This regularisation then depends on a 'schemechoice', namely the choice of the tangent vector v x 0 at the initial point x 0 , and the Laurent coefficients of the path-ordered exponential will depend explicitly on log v. The final result for the integral can of course not depend on the scheme choice: the initial condition J 0 (ǫ) will also depend on the choice of v in such a way that the dependence cancels. It is customary to choose v = 1, so that no explicit logarithms of v appear.
The regularisation we just described seems to be very hard to implement in practice, because it requires one to introduce a cut-off, expand in the regulator ε and then project onto the constant term. It is actually possible to obtain a closed formula for the shuffleregularised version. We only present here the result, and refer to ref. [138], section 4, for details. It is useful to define ω ∞ i = a i dlog ξ, so that ω i = ω ∞ i + · · · (cf. eq. (132)). We then have [138] x v 0 where the map R is defined by The map R replaces the word ω 1 · · · ω n by a linear combination of words so that their iterated integral is convergent over the whole range [0, x].
Remark 16. Equation (136) takes a very simple form if there is a single one-form that diverges at ξ = 0. In that case eq. (136) is equivalent to unshuffling all the occurences of this one-form. For example, for ω 1 = dlog ξ and ω 2 = dlog(1 − ξ), eq. (136) is equivalent to the well known shuffle-regularisation formula: x v 0 where Li 2 (x) is the dilogarithm function, defined for |x| < 1 by Remark 17. The discussion in this section only applies to logarithmic singularities, which is sufficient to cover systems of differential equations in canonical form. If also one-forms with higher-order poles are present, the regularisation becomes more involved. For a discussion of this more general case, see ref. [139].

Linear independence of iterated integrals
In applications one is often interested in knowing that the special functions introduced to express the answer are (linearly) independent. Indeed, working with an independent set of objects often leads to shorter analytic expressions that are free of hidden cancellations. The goal of this section is to state a linear independence result for iterated integrals (in particular those from systems of differential equations in canonical form). The mathematical background can be found in refs. [99,140]. Let us consider a system with A(x, ǫ) = ǫ A(x) = ǫ i A i ω i (we follow the notation of section 3, but we do not explicitly require the ω i to have only logarithmic singularities for the purpose of this section). At every order of the ǫ-expansion, the path-ordered exponential will involve iterated integrals of words in the one-forms ω i . Clearly, if we work with arbitrary one-forms ω i , these iterated integrals will not be independent. For example, every linear relation between the ω i will induce linear relations among the iterated integrals. Moreover, if a one-form is a total derivative, ω i = df for some function f , then we can integrate out ω i , again leading to relations among the iterated integrals.
Here we have to make an important comment: of course, we can always locally find a primitive f for every one-form ω, and this f will in general be a transcendental function. In applications, however, we are typically interested only in primitives taken from some restricted subalgebra C of functions such that in local coordinates ω i = j c ji (ξ) dξ j with c ij (ξ) ∈ C. For example, in the case of dlog-forms ω i = dlog a i (ξ) with algebraic a i (ξ), we would take C to be the field of algebraic functions in ξ (in the case of non-dlogforms non-algebraic functions may also appear, such as modular forms). More generally, we will then say that the ω i are linearly dependent with respect to C if there exists a function f ∈ C and constants α i not all zero such that Linear independence with respect to C is defined in the obvious way, such that every identity of the form (140) implies α i = 0 for all i. It turns out that this linear independence of the one-forms is sufficient to guarantee that iterated integrals are independent as functions.
Proposition 12. The iterated integrals arising from the path-ordered exponential are linearly independent over C as functions if and only if the one-forms ω i are linearly independent with respect to C.
The proof of this statement can be found in ref. [140]. ¶ Remark 18. It is important to understand that we consider the linear independence of the iterated integrals as functions. If we consider evaluations at special values, there may be additional relations. For example, if ω 1 = dlog ξ and ω 2 = dlog(1 − ξ), then the iterated integrals x 1 0 ω 1 ω 2 and x 1 0 ω 2 ω 1 are linearly independent as functions of x. If we evaluate these integrals at x = 1, we find

Iterated integrals and Feynman integrals
We have already seen that Propositions 9 and 10 imply that iterated integrals naturally arise from differential equations satisfied by dimensionally-regulated Feynman integrals. In applications it is often desirable to relate the results of a computation to known definitions of special functions which have been studied independently in the literature, e.g., there may be computer codes and algorithms for their manipulation and evaluation.
There are several classes of iterated integrals that are well studied in mathematics and physics, which we briefly review in this section. Before doing so, however, we note that in some applications it is beneficial to consider classes of iterated integrals that are specific to a given kinematic configuration, such as pentagon functions for five-point Feynman integrals [141][142][143] or hexagon functions for specific six-point integrals [144,145], even if these very specialised iterated integrals can be related to more general classes of iterated integrals.
Note that if a n = 0, this integral is divergent, and we need to regularise it. This is typically done by introducing a tangential base-point 1 0 for the lower integration boundary. This is equivalent to the special definition (see eq. (136)) MPLs contain the ordinary logarithm and the classical polylogarithms as special cases, e.g., if a = 0, G(a; x) = log 1 − x a and G(0, . . . , 0 It is easy to check that MPLs are in fact the prime examples of pure functions (see section 3.3). The properties of MPLs, in particular their algebraic structure and the relations they satisfy, are well understood. An important tool in the study of MPLs is their so-called symbol and their coaction/coproduct, cf. refs. [99,[133][134][135][136][157][158][159][160][161]. Moreover, there are several publicly available computer tools to manipulate and evaluate MPLs [63,[162][163][164][165][166][167][168][169][170][171][172]. It has been known from the early days of quantum field theory that not all Feynman integrals can be expressed in terms of MPLs. Indeed, A. Sabry discovered new functions of elliptic type in the calculation of the two-loop corrections to the electron propagator in QED with massive electrons already in 1962 [173]. Since then, many other Feynman integrals have been identified that cannot be expressed in terms of MPLs [80,81,[129][130][131]. There is a natural class of iterated integrals that generalise MPLs to functions of elliptic type, called elliptic multiple polylogarithms (eMPLs) [199,200]. eMPLs have first appeared in physics in the context of one-loop scattering amplitudes in string theory [201][202][203][204], but they have also been used to express several multi-loop Feynman integrals that cannot be expressed in terms of eMPLs [68,[205][206][207][208][209] (see refs. [189-192, 210, 211] for alternative definitions of elliptic generalisations of polylogarithmic functions that are closely related to eMPLs). Closely related to eMPLs are iterated integrals of modular forms [138,212], including meromorphic modular forms [139], which have also appeared in the context of Feynman integrals [107,[213][214][215][216][217]. Iterated integrals of holomorphic modular forms and eMPLs are examples of pure functions [205], and the corresponding Feynman integrals satisfy differential equations in canonical form [106,107,131] (though it is not possible to bring them into canonical dlog-form). The properties of these functions are not yet as well understood as in the case of MPLs, although substantial progress was made in recent years to understand their algebraic properties and numerical evaluation, see, e.g., refs. [218][219][220][221].
MPLs, eMPLs and iterated integrals of modular forms are not the only classes of iterated integrals that have been observed to arise from Feynman integral computations. Other classes of special functions have been discovered [89,129,130,167,196,206,210,[222][223][224][225][226][227][228][229][230][231][232][233], even though for some of them it is known by now that they are related to the classes of functions discussed above. A completely new class of special functions has recently been discovered that is related to higher-dimensional generalisations of elliptic curves known as Calabi-Yau manifolds [93,94,194,[234][235][236][237], with new classes of iterated integrals which generalise the iterated integrals of modular forms to Calabi-Yau varieties of higher dimension [59].
MPLs and systems in canonical dlog-forms. Let us conclude this discussion of iterated integrals that arise from Feynman integrals by commenting on iterated integrals of dlog-forms, i.e., iterated integrals of words of letters of the form dlog a i (ξ), with a i (ξ) an algebraic function. So far in all cases relevant to Feynman integrals these algebraic functions only involve square roots. + It is clear that iterated integrals of dlog-forms are closely related to MPLs. Indeed, if the a i (ξ) are linear functions, then we immediately reproduce the definition in eq. (142). In many cases it is still possible to evaluate the resulting iterated integrals in terms of MPLs, even if the a i (ξ) are more general rational functions or even involve square roots.
• If all the a i (ξ) are rational functions of ξ, we can evaluate the iterated integrals in terms of MPLs in an algorithmic fashion. Indeed, we can assume without loss of generality that the a i (ξ) are polynomials. We can choose a piecewise constant path to evaluate the path-ordered exponential (cf. section 3.3.1). On the segment γ i 0 where all elements of ξ but ξ i 0 are constant, we can factor each polynomial a i (ξ) into linear factors in ξ i 0 , a i (ξ) = c i k (ξ i 0 − ξ i 0 ,k ), and we obtain only dlog-forms depending on linear arguments: • If the a i (ξ) involve square roots, the situation is different. We have to distinguish two cases. First, if we can find a change of variables ξ = ϕ(χ) that rationalises all square roots, then the functions a i (ϕ(χ)) are rational in χ, and we reduce the problem to the previous case. In some cases it is possible to find such a change of variables in an algorithmic manner, or to show that is does not exist cf., e.g., refs. [238][239][240][241][242][243][244]. If not all square roots can be rationalised, then one cannot conclude a priori if a representation in terms of MPLs (with algebraic arguments) exists. Indeed, in ref. [245] an explicit example of a double-iterated integral of dlog-forms with a square root was constructed that cannot be expressed in terms of MPLs. So far, however, all instances of systems in canonical dlog-form with + It is interesting to ask if algebraic functions not involving square roots can also arise in Feynman integral computations.

Intersection theory for Feynman integrals
The previous sections have discussed the evaluation of Feynman integrals by fairly direct approaches. Recently there has been some interest in applications of the mathematics going by the name of intersection theory, which deals with certain multivalued integrals. One area of application is an alternative approach to computing IBP reduction coefficients, which might help to overcome some of the shortcomings of current approaches. Others involve an algebraic operation called the coaction, which can expose the behaviour of functions under differentiation or operations computing discontinuities.
In this section, we review the formalism of intersection theory for a class of integrals that includes Feynman integrals in dimensional regularisation, as well as the Euler-type hypergeometric integrals that arise upon their evaluation. A principal classic reference for the mathematics is ref. [249]. For recent treatments including applications to physics, see refs. [7,[250][251][252].
An integral is a pairing between two objects, the integrand and the contour of integration. Denoting these objects respectively by ω and γ, the integral is represented as We will consider cases where the differential form ω is a holomorphic differential form of degree n, and the real dimension of the integration contour is n, within the complex projective space P n (C). The integral remains unchanged if we add a surface term to the integrand, and ω is properly understood as a cohomology class, or cocycle. Likewise, the integration contour γ is a homology class, or cycle. In the types of integrals with which we are concerned, the integrand involves multivalued functions. This can be seen in the parametric representations of section 1.2. For instance, let us consider the Baikov representation of a Feynman integral in eq. (19). The Baikov polynomial is raised to the power D−K−1 2 = n + µ 2 − ǫ, where n is an integer and µ ∈ {0, 1}. Thus the integrands of Feynman integrals in dimensional regularisation are always multivalued functions. Because of this multivaluedness, we need to work in the frameworks of twisted cohomology and homology, which we review below.
An example that we will keep in mind and follow through this section is the integral representation of the hypergeometric function 2 F 1 defined in eq. (43). While this function is not a Feynman integral by itself, some simple Feynman integrals such as the bubble studied in the previous sections are naturally expressed in terms of 2 F 1 functions (see eqs. (42) and (44)), and the additional prefactors do not affect the main points of intersection theory. Feynman integrals, as presented in eq. (17) and eq. (19), can be viewed in terms of generalisations of this integral, i.e., polynomials in the kinematic invariants raised to complex powers. For the purpose of the present discussion, we will simplify this integral even further by dividing through by the gamma functions. Furthermore, we relabel the exponents to make their independence manifest. Thus we write our basic reference integral as

Twisted cohomology and twisted homology
Twisted cohomology. We begin by considering the integrand ω, which we take to be a holomorphic n-form ω = du where u = (u 1 , u 2 , . . . , u n ) and du = du 1 ∧ . . . ∧ du n , where the P I are polynomials in the integration variables u i and additional kinematic variables x j , and with complex exponents, α I ∈ C. In the case of Feynman integrals in dimensional regularisation, we would assume that the exponents take the form α I = n I + a I ǫ, with n I ∈ Z, a I ǫ ∈ C * , I a I = 0, and where ǫ can be taken to be infinitesimally small. We restrict ourselves to the case where the coefficients in the Laurent expansion of the integral γ ω can be given in terms of multiple polylogarithms, which requires the integer values of n I and some further restrictions on the form of the polynomials P I (u).
We decompose the integrand as ω = Φϕ, where Φ is a multi-valued function and ϕ is a single-valued differential form, Φ = I P I (u) a I ǫ and ϕ = du We define the twist dlog Φ and consider the covariant differential We then have d(Φξ) = Φ ∇ Φ ξ, where ξ can be any smooth differential form. Stokes' Theorem implies that for an arbitrary smooth (n − 1)-form ξ we have Thus the integrand is only defined up to adding a total covariant derivative, and we are therefore considering elements of the (twisted) cohomology groups One might consider the twisted cohomology group with the differential of eq. (150) acting on differential forms of different degrees, but a theorem of Aomoto [253] states that n is the only dimension with nonvanishing twisted cohomology.
Equation (151) provides an alternative way to generate the IBP identities discussed in section 2. The expansion of these covariant derivatives in terms of the original polynomials P I (u) has the effect of shifting the integer parts n I of the exponents in eq. (148), while leaving the complex part a I ǫ unchanged. The space generated by the integrands with different values of n I gives the full space of twisted cohomology.
Example 8. Let us return to our example of the integral in eq. (147). Comparing eq. (147) with eq. (148), we see that n = 1, |I| = 3, and the set of polynomials can be taken to be {P 0 = u, P 1 = 1 − u, P 1/x = 1 − xu}. We have set u 1 = u for simplicity, not to be confused with the multi-index u in general formulas such as eq. (148).
The space of integrands modulo the IBP relations is known to be two-dimensional. We will discuss this fact along with a choice of basis shortly.
Twisted homology. The integration contour γ is a n-dimensional (relative) cycle in Strictly speaking, {P I (u) = P I (u 1 , . . . , u n ) = 0} is an affine variety in C n . We use the same notation for the affine variety and its lift to projective space. In other words, γ is a domain whose boundary is contained in the union of the varieties defined by P I (u) = 0 (or equivalently Φ = 0). Note that this is in particular the case for the Baikov representation in eq. (19), where the integration cycle ∆ is bounded by the variety where the Baikov polynomial vanishes. Since a I = 0, then Φ vanishes on the boundary of γ, at least for some ranges of values of ǫ, and thus for all values by analytic continuation. Therefore, if a I = 0, there are no boundary contributions when performing integration by parts.
The concept of a twisted cycle, as an element of twisted homology, is essentially a version of the equivalence class of γ that is taken together with a specific choice of branch of Φ to deal with the multi-valuedness. The ordinary Stokes' Theorem holds on the branch and corresponds to the version with the covariant derivative given above.
The effect of the twist is to replace the boundary of the cycle by a contour that encircles the singularity. The simplest twisted cycle is a loop around a singular point. Precise definitions and explanations may be found in section 2.3 of ref. [249]. When used as an integration contour, the result is the same as the physicists' notion of the so-called 'plus-prescription' that regularises the integral in the ǫ expansion (see e.g. [254]).
In the case of the integral in eq. (147), we are given the contour γ = [0, 1], which is consistent with the description in eq. (154). We will not write the twist explicitly, but it can be constructed with the knowledge of Φ, and its use is implicit in regularising divergences if we perform a series expansion in ǫ before evaluating the integral. We note that there are various other cycles that can be constructed with boundaries contained in the union of varieties defined by Φ = 0, namely cycles whose endpoints belong to the set {0, 1, 1/x, ∞}. It turns out that at most two of them can be linearly independent (up to boundaries), as we discuss below in the context of a basis choice.
Bases of twisted cohomology and twisted homology. It will be useful to identify explicit bases for the twisted cohomology and homology groups associated to a particular integral. The problem of constructing bases algorithmically appears to be difficult to solve. In practice, it is helpful to first identify the dimension of these groups, then select a set of elements of that cardinality, and finally test for linear independence. It is not obvious that the dimension of twisted cohomology should equal the dimension of twisted homology, but in the case where a universal coefficient theorem holds, there is a natural isomorphism between these groups, leading to the bilinear pairing expressed by integration (see Lemma 2.3 of ref. [249]). For all cases of interest in this review, the dimensions of the two groups coincide. Even in such cases, however, determining the dimension of the groups is not straightforward in general. There is an upper bound given by the number of critical points of the function Φ, i.e., the number of independent solutions to the equation and in many cases this bound is saturated, for example under conditions outlined in ref. [249] (see also refs. [20,255,256]).
Example 9. In the example of the integral in eq. (147), we see from eq. (153) that the equation dlog Φ = 0 has two solutions generically, since we can factor out a quadratic polynomial in u after multiplying the denominators through. Accordingly, the dimension of co/homology in this case is 2. To choose a basis of cohomology, we can try two different sets of integers {n 0 , n 1 , n 1/x } in ϕ = u n 0 (1 − u) n 1 (1 − xu) n 1/x du and check their linear independence as cohomology classes. To choose a basis of homology, we can choose two sets of endpoints among {0, 1, 1/x, ∞}. In order to follow this example through the following sections, let us make a concrete choice and take and Here we have chosen differential forms that have logarithmic singularities exactly at the endpoints of the cycles. This type of differential form is called a canonical form associated to the cycle [257]. * In the following subsection, we use pairings of bases of twisted co/homology to construct the period matrix and matrices of intersection numbers in twisted * The notion of canonical form used here is a differential form and should not be confused with the differential equations in canonical form introduced in section 3.
co/homology. The linear independence of a putative basis can be confirmed by checking the rank of these matrices.
Remark 19. One of the assumptions made in the mathematical treatment of intersection theory is that none of the exponents α I is an integer, and moreover that their total sum is also not an integer. This assumption is not valid for typical Feynman integrals and their associated hypergeometric functions, but it is found in practice that physical results can be obtained by taking suitable limits of the generic case.

Pairings
Now that we have introduced the twisted cohomology and homology groups associated to an integral γ ω, we consider functions that pair their elements through complexvalued bilinear maps. The first pairing is the familiar integral map, with a possible physical interpretation of a Feynman integral or a cut Feynman integral. Other pairings can be constructed between elements of the same group, giving intersection numbers of the twisted cohomology or twisted homology.
Integral pairings. The integral pairing brings us back to our starting point, namely the integral γ ω obtained from the pairing of γ and ω. Now that we have identified these objects as elements of twisted homology and cohomology groups, respectively, we can consider the full space obtained from this pairing, giving integrals related to the original one.
Let γ denote a basis of the (twisted) homology group, and let ϕ denote a basis of the (twisted) cohomology group. The integral pairings of the cycles γ l ∈ γ with the forms ϕ k ∈ ϕ give entries of the so-called period matrix, where each row is associated to a differential form, and each column is associated to a cycle. Here we have also introduced the bra-ket notation of ref. [258], which will also be used for the other types of pairings. This notation should properly be accompanied by a subscript Φ to keep track of the twist, but since the twist remains constant throughout our discussion, we will suppress it.
In the context of Feynman integrals, each column of the period matrix comes from a basis of contours associated to generalised cuts, while each row of the period matrix comes from a basis of integrands associated to master integrals. The period matrix can in fact be interpreted as (the transpose of) a fundamental solution matrix for the system of differential equations satisfied by the master integrals, and is thus closely related to the Wronskian that was extensively used in section 3.
The period matrix P is a square matrix whose dimension is given by the dimension of the (co)homology group. It follows that any integral pairing of elements of the twisted homology and cohomology groups (with the same twist) can be written as a linear combination of the elements of the period matrix, The algebraic properties of any integral of this type can then be studied from the entries of the period matrix.
Example 10. For the bases given in eqs. (156) and (157), we find the following entries of the period matrix: Cohomology intersection numbers. A less obvious pairing can be constructed between two differential forms. Consider two bases ϕ and ψ, not necessarily distinct, of the same twisted cohomology group. We can then compute cohomology intersection numbers ϕ i |ψ j between these forms. To be more precise, we must first construct a dual twisted cohomology group, which is also generated by ψ but for which the covariant differential is ∇ Φ −1 . For dimensionally regularised integrals, this dual operation corresponds to taking ǫ → −ǫ in Φ. We can then pair generators ϕ i | of the cohomology with elements |ψ j of the dual cohomology [249,258], where X(C) is as defined in eq. (154), and ι Φ is the map that associates to a form ψ j a form ι Φ (ψ j ) in the same cohomology class but with compact support, so that the integral is well defined [258,259]. That is, the form ι Φ (ψ j ) vanishes in a neighborhood of the space Φ = 0 in P n (C). Unlike integral pairings, the cohomology intersection numbers are single-valued. The intersection numbers of basis elements can be arranged in the matrix C kl ( ϕ, ψ) = ϕ k |ψ l which has the same dimensions as the period matrix P , and they similarly form a basis of all intersection numbers in the twisted cohomology. The definition of the intersection numbers in eq. (160) is not always the most convenient for practical calculations, even though it has recently been used in refs. [103,260], so alternative formulations have been found. One version pertains to the case where n = 1, and the ϕ i and ψ j are dlog-forms, i.e. wedge products of 1forms dlog α k , where the singularities of these forms are contained within the boundary variety Φ = 0. In this case, a more explicit formula for the intersection numbers is given by [250,258,261] where P(Φ) is the set of poles of dlog Φ. This formula has been generalised to the case where n > 1 in ref. [262]. When the ϕ i and ψ j are not necessarily dlog-forms, other alternative formulas were proposed in ref. [258]. For n = 1, and setting u 1 = u, where the sum is over the critical points, i.e., the points u * satisfying dlog Φ(u * ) = 0, and ϕ i = ϕ i du and similarly for ψ j . In the case n = 2, with (u 1 , u 2 ) = (u, v), where the sum extends over the critical points (u * , v * ) satisfying and these formulas generalise naturally to higher values of n. Reference [263] presents a different algorithm using a Gröbner basis for generating a basis of twisted cohomology, with the feature that algebraic extensions such as square roots are avoided. Cohomology intersection numbers are algebraic functions of the kinematic variables x j and exponents α I . For a dlog basis, the leading order in ǫ of the period matrix agrees with the matrix of intersection numbers [261], which gives yet another method of computing these numbers.
Example 11. In the example of the integral in eq. (147), and the basis of eq. (157), it is straightforward to see from eq. (162) that the matrix of cohomology intersection numbers is Homology intersection numbers. For completeness, we comment briefly on homology intersection numbers, although these have not yet been applied as widely in the context of physics. In the definition of cohomology intersection numbers above, we saw that one of the differential forms needed to be given a compact support. The dual construction for twisted cycles is to use a locally finite chain group. Then the intersection numbers [γ k |γ l ] can be computed directly as geometric intersections of the cycles γ k and γ l . We refer to section 2.3.3 of ref. [249] for a fuller explanation and an illustrative example. Given a pair of bases γ, δ of twisted homology, the intersection numbers can be collected in matrix form, Here the bra elements [γ l | are the twisted cycles described above, while the ket elements |δ l ] are generators of a locally finite version of twisted homology.
Twisted period relations. The pairings expressed in the period matrix co/homology intersection numbers give rise to the following quadratic relations [264]: The subscripts on the period matrix refer to two possible pairings: P −kl ( ϕ, γ) = ϕ k |γ l ] is the usual integral γ l Φϕ k , while P +kl ( ϕ, γ) = [ϕ k |γ l is a pairing between the corresponding compactified cocycles and locally finite cycles, which can simply be evaluated as the integral γ l Φ −1 ϕ k . For dimensionally regularised Feynman integrals, this duality is simply the transformation ǫ → −ǫ.
As applied to the Gauss hypergeometric function, the twisted period relation can be expressed as the following quadratic relation: Since twisted period relations are a generic feature of twisted cohomology and homology theories, they imply that there are generically quadratic relations among dimensionallyregulated Feynman integrals. It would be interesting to explore these relations more systematically, and to relate them to the quadratic relations among Feynman integrals discussed in section 2.6.

Applications of intersection theory
This section presents a brief review of a few different applications of the intersection theory introduced above to the computation of Feynman integrals. We would also like to mention that intersection theory has found physical applications in the context of amplitudes in string theory and related theories, as initiated in ref. [258,265], but we will not discuss those ideas here. 5.3.1. Reduction and differential equations Intersection theory was applied to compute the coefficients in a reduction to master integrals in refs. [103,256,260,261]. Suppose that the integral of interest can be represented as The reduction to master integrals can be represented as where ϕ i is a basis of twisted cohomology, and the reduction coefficients are denoted by c i . This equation is fundamentally an operation on the cocycles, which we can write as By projecting this equation onto the basis of cocycles, one finds that the reduction coefficients can be computed by cohomology intersection numbers: Thus the computation of cohomology intersection numbers can be substituted for reduction techniques. In particular, this approach avoids having to solve large systems of equations, addressing one of the bottlenecks identified in section 2.3. However, and despite much progress in the last few years, the difficulty associated with computing intersection numbers still makes the Laporta algorithm the method of choice in practical calculations. When applied to the example of the 2 F 1 integral, we find a version of the Gauss contiguous relations, which express the property that the space generated by integer shifts of the parameters α, β, γ in 2 F 1 (α, β, γ; x) is two-dimensional. In this case, it is equivalent to the statement that the twisted cohomology group is two-dimensional.
Intersection theory can also provide the coefficients of the differential equations of eq. (56) [261]. If the differential operator is applied to the integral pairing, and if the integration contour does not depend on the kinematic variables, then we obtain with the covariant derivative ∇ Φ,x i = ∂ x i + ∂ x i log Φ∧, where here we take Φ = B(z) . The right-hand side can be expanded in master integrals as in eq. (169), and the same reasoning followed to obtain the coefficients in terms of intersection numbers. Thus we can identify the matrix of coefficients as In the Baikov representation, uncut (or less than maximally cut) Feynman integrals have singularities that are not regulated by the parameter ǫ. These singularities can be treated with analytic regularisation [3], which then allows the application of intersection theory [262,266]. A complementary approach is given in refs. [103,260], which focus on the cohomology that is dual to the cohomology of Feynman integrands, for the purpose of constructing intersection numbers as above. These papers find that the correct mathematical framework is that of relative twisted co/homology as introduced recently in [267], which extends the scope of intersection theory to allow relative boundaries.
Remark 20. The decomposition into a basis of master integrals with the help of intersection numbers has a dual application for integration cycles. Assume that we have fixed a basis of cut integrals, i.e., we have fixed a basis of integration cycles γ. Then every other integration cycle |γ] can be decomposed into this basis using the approach we have just described, but using the intersection numbers among cycles instead of intersections numbers among differential forms: In this way we can decompose a given cut integral γ|ω] into the basis of cut integrals γ i |ω].

Coaction
MPLs are well known to be equipped with a coaction [150,157,158,160], a mathematical operation that is naturally compatible with the actions of differential operators and taking discontinuities across branch cuts. As such, it can be employed as a computational tool. There is also a natural mathematical coaction on Feynman integrals [268][269][270]. Motivated by these observations, it was conjectured that there exists a diagrammatic coaction on Feynman integrals which, for the class of integrals admitting a Laurent expansion in ǫ in terms of MPLs, corresponds precisely to the coaction on MPLs. This conjecture was expressed in refs. [100,101,271], which proposed a general form for a coaction on integrals:♯ where the c ij are rational or algebraic coefficients, the {ω i } are a basis of the (twisted) cohomology group associated with the integral on the left-hand side, and the {γ j } are are a basis of the corresponding homology group. † † It is straightforward to check that eq. (175) satisfies the algebraic properties of a coaction, such as coassociativity. For oneloop integrals, the diagrammatic coaction was made fully explicit by providing bases of forms and contours for arbitrary one-loop integrals [100,101]. It was further conjectured in refs. [276,277] that eq. (175) can be applied to generalized hypergeometric functions in their integral representations, of the type ♯ This form of a coaction is also believed to apply to the disk integrals in string tree-level amplitudes [272,273], as well as the more general genus-zero integrals in ref. [274]. † † See also ref. [275] for a separate proposal of a diagrammatic coaction. considered in this section, of which the prototypical example is the function 2 F 1 that we have been following. The coaction operation proposed for hypergeometric functions is claimed to be compatible with the diagrammatic coaction, as well as to commute with taking the Laurent expansion in ǫ, at least in the case where there is an expansion in terms of MPLs. For this class of functions, a rigorous (motivic) mathematical treatment has been developed in ref. [278] and applied specifically to Lauricella F D functions. It is hoped that compatible coactions could be formulated for more general Feynman integrals, but this is more speculative, as it requires formulating coactions on the new classes of functions arising in the ǫ expansion. This has been done for elliptic eMPLs in ref. [214], and an example of a corresponding diagrammatic coaction is given in ref. [279].
If one can choose generators of the (co)homology such that C(Ω( γ), ϕ) = δ ij , then the coaction takes a particularly simple form: Example 12. Throughout this section, we have followed the example of the hypergeometric function 2 F 1 . Any choice of bases will give an expression for the coaction using the general formula eq. (176). But let us see how to arrive at a particularly elegant expression.
For the basis of contours, let us take γ 1 = [0, 1] from the usual Euler integral, and γ 2 = [1/x, ∞). One benefit of using γ 2 is that its pairing with the general integrand can be recognised as a 2 F 1 function of x after a simple change of integration variable from u to v = 1/(xu).
It may be desirable to choose a second basis of differential forms ϕ i , such that the matrix of intersection numbers C( ϕ, ψ) has a minimum number of nonvanishing offdiagonal elements. Indeed, C will be diagonal if each ϕ i is taken to be a dlog-form whose singularities overlap with the boundary components of γ j if and only if i = j. In this case, our chosen cycles γ 1 and γ 2 have no common boundaries. The choice ϕ = (ϕ 1 , ϕ 2 ) with ϕ 1 = a 0 a 1 a 0 + a 1 ǫ du u(1 − u) , ϕ 2 = (a 0 + a 1 + a 1/x )a 1/x a 0 + a 1 ǫ x du 1 − xu , uses dlog forms normalised so that the matrix of cohomology intersection numbers is the identity.