Lectures on differential equations for Feynman integrals

Over the last year significant progress was made in the understanding of the computation of Feynman integrals using differential equations. These lectures give a review of these developments, while not assuming any prior knowledge of the subject. After an introduction to differential equations for Feynman integrals, we point out how they can be simplified using algorithms available in the mathematical literature. We discuss how this is related to a recent conjecture for a canonical form of the equations. We also discuss a complementary approach that allows based on properties of the space-time loop integrands, and explain how the ideas of leading singularities and d-log representations can be used to find an optimal basis for the differential equations. Finally, as an application of the differential equations method we show how single-scale integrals can be bootstrapped using the Drinfeld associator of a differential equation.


Introduction
Feynman integrals are ubiquitous in quantum field theory. They occur when quantities such as correlation functions of local operators, or scattering amplitudes are computed in perturbation theory. Beyond the lowest order in perturbation theory, integrals over Ddimensional space time over (rational) propagator factors have to be evaluated.
Our ability to evaluate Feynman integrals in perturbation theory is critical for connecting quantum field theory to experiment. For example, in order to make precise theoretical predictions for collisions at the Large Hadron Collider (LHC), higher-order Feynman integrals for complicated kinematical processes are a crucial ingredient, see e.g. [1]. This not only applies to virtual amplitudes, but, via the optical theorem, also to phase space integrals and to cross sections [2,3]. Other applications include statistical field theory, where critical exponents can be computed; higher order calculations are also important for instance for a better theoretical description of the anomalous magnetic moments.
From the mathematical point of view it is an interesting question to ask what class of functions Feynman integrals give rise to. These are typically functions of a set of parameters, such as the particle momenta and masses in the case of scattering amplitudes, or a set of points in the case of correlation functions. Experience shows that typical functions occurring in Feynman integrals are certain classes of iterated integrals, elliptic functions, and possibly generalizations thereof. It is an interesting open problem in general to predict, for a given Feynman graph, what class of functions it is described by. Understanding the above questions also implies, for fixed values of the external parameters, an understanding of what periods Feynman integrals give rise to. Periods are a class of numbers that is situated between the rational and irrational numbers [4].
Some Feynman integrals have divergences that can come from different regions of integration, either ultraviolet (short distances) or infrared (long distances). In that case, integrals can be defined in D = 4 − 2 dimensions, and one is usually interested in the solution as a Laurent series in . Physically relevant quantities can be defined via a renormalization procedure, such that divergences ultimately cancel, but it is useful to be able to compute the expansion, in principle to any order in . From the mathematical point of view, the question about periods and functions can be posed in a well-defined way for each term in the expansion.
In these lectures we present the method of differential equations (DE) for computing Feynman integrals. It has a long and successful history [5][6][7]. Earlier reviews include [8]. Despite the success of the method, so far it was mostly used on a case by case basis, with considerable amount of work needed by hand to solve the equations. Over the last year a more systematic picture has emerged, and it is the goal of these lecture notes to present this in a coherent fashion.
One of the key new ideas of ref. [9] is to choose an optimal basis of integrals that leads to a system of differential equations in a canonical form. We present two different, complimentary approaches to this end, one of which is more algebraic, the other one more geometric in nature.
Functions coming from Feynman integrals have regular singularities as kinematic vari-ables approach a given singular point. At the level of the differential equations, the canonical form should make the behavior of the Feynman integrals near singular points manifest. This can be done using known algorithms available in the mathematical literature. Having put the kinematic dependence into a canonical form, the next step is to simplify the dependence on the dimensional regularization parameter . For example, in cases where the answer is given by iterated integrals, the idea is to make this iterative structure manifest at the level of the differential equations. We review the features of this algebraic approach, and explain how elliptic and more complicated functions arise from this viewpoint. The approach mentioned above is an algebraic one, which to a large extent could be applied to any system having regular singularities only. There is a complimentary approach that allows to gain insights on an appropriate basis choice already at the level of the Feynman integrand. The idea is to analyze the structure of (generalized) unitarity cuts of the latter [10][11][12][13], where, roughly speaking, a number of propagators is replaced by delta functions. It had been observed earlier that the analysis of generalized cuts (and in particular, leading singularities) is a natural way to choose a basis for loop integrands that make expressions for scattering amplitudes simple, and e.g. exhibit their properties in soft limits [12,14]. In ref. [9] it was proposed to use these ideas in the context of differential equations in order to simplify the latter. In that context, considering cuts of integrals is very natural, since the cut integrals satisfy the same differential equations [3], albeit with different boundary conditions. Moreover, the cuts can be viewed as projection operators onto sectors of the DE, which is a useful feature in practice. This explains why cuts are an efficient guide to finding a canonical form of the DE.
A closely related tool are certain d-log representations, either at the level of parameter integrals, or directly at the level of loop integrals. In the case where the answer is given by iterated integrals, the former can be used to prove transcendental weight properties of the answer, and guarantee that a canonical form of the DE can be found. In the latter case, the transcendental weight property is expected conjecturally [13].
Over the last year, these tools have been applied successfully to the calculation of many classes of Feynman integrals. In practice, often the most efficient strategy is to choose the basis integrals according to their cut properties, which is most natural for Feynman integrals. This usually leads to a form of the differential equations that is very close to the desired canonical form, and the remaining algebraic transformations are then easily found.
These lecture notes are intended for students and researchers alike. We hope that they are useful to phenomenologists desiring to get acquainted with this method, in order to apply it to loop integrals relevant for their calculations, as well as to more formal theorists seeking a general understanding of multi-loop scattering amplitudes, and last but not least to mathematicians interested in new ideas that do not rely on the usual Feynman parametrization approach.
The style we have tried to adopt is one where the main ideas are explained with the help of simple examples. In this way, we avoid the burden of a general but possibly unintelligible notation, and are able to keep the lectures at a digestible length 1 . We frequently summarize general conclusions.
These lecture notes are organized as follows. In section 2, we recall basic definitions of Feynman integrals and mention important properties that follow from them. We continue in section 3 to introduce notions that allow to understand algebraic relations between different Feynman integrals and to derive differential equations for them. In section 4 we show how general properties of Feynman integrals allow to transform the differential equations into a canonical form. In section 5 we discuss solutions in the case where the answer is given by iterated integrals. In section 6, we explain the analysis of generalized cuts / leading singularities of Feynman integrals, and d-log representations to choose an optimal basis. As an example, we fully explain the basis choice made in ref. [9] using these concepts. The final section 7, which can be read independently, combines the ideas of the previous sections for a sample application, namely the computation of single-scale integrals via differential equations.

Definition and basic properties of Feynman integrals
Let us begin by recalling the main definitions and introduce the notation that we will use in the following. Our aim will be to be brief, as more details can be found in standard textbooks, e.g. [15,16]. The main conclusions of this section are summarized at the end for the benefit of the reader already familiar with this material.

Definitions and Feynman parametrization
We will discuss Feynman integrals in D-dimensional quantum field theory. In the momentum space language, we have integrals over D-dimensional space d D k, with the integrand consisting of propagator factors like 1/[−(k + p) 2 + m 2 − i0]. The Feynman i0 prescription allows one to perform a Wick rotation from Minkowski space with metric + − ...− to Euclidean space, see e.g. [15,16]. In most of the following we drop the i0 from our formulas for simplicity of notation.
As an example, let us start with a momentum-space box integral at one loop. This example will recur frequently in these notes.
Here p 1 , p 2 , p 3 , p 4 are D-dimensional momenta satisfying momentum conservation 4 i=1 p i = 0 and the on-shell conditions p 2 i = 0. Poincaré invariance implies that the integral depends on the Mandelstam invariants s = (p 1 + p 2 ) 2 and t = (p 2 + p 3 ) 2 only. Moreover, the integral is covariant under dilatations, so that after normalizing it with some appropriate power of s or t it becomes a function of the dimensionless variable x = t/s only.
A very convenient notation for planar integrals as the above are dual (or region) coordinates. This constitutes a change of variables. We can rewrite the above integral as where now s = y 2 13 and t = y 2 24 , with y ij = y i − y j . This notation is very convenient, as it allows us to give a very simple expression for the Feynman parametrization of such integrals. Let us consider a general one-loop integral in the dual notation, where we have allowed the propagators to depend on masses m i , and to be raised to arbitrary powers a i . Introducing Feynman parameters in the usual way, performing a Wick rotation to Euclidean space, and carrying out the D-dimensional integration using a Gaussian integral, one arrives at This formula is derived for integer dimension D and exponents a i , in a region where the Feynman integral converges. We will use eq. (2.4) as the definition for integrals with other values of these parameters, via analytic continuation. At L loops, there is a generalization of eq. (2.4), with U and V being homogeneous polynomials in the α parameters that have a graph theoretical definition, see e.g. [15]. Example: One-loop on-shell box integral in D = 6 dimensions. In this simple case, we can directly carry out the Feynman parameter integrations. Applying the general one-loop formula (2.4) with D = 6, we obtain .
Here we assume s < 0, t < 0, which allows us to drop the Feynman i0 prescription. Changing variables according to with Jacobian w 3 z(1 − z), this is easily evaluated, with the result (2.8)

Infrared and ultraviolet divergences
Feynman integrals occurring in quantum field theory can have two types of divergences having a specific physical origin. The first are ultraviolet (UV) divergences, related to the region of large loop momentum k. The degree of UV divergence of an integral can be easily determined via power counting. For example, in the large momentum region the one-loop integral I n can be approximated by for R = |k| 1, so that a divergence occurs if D < 2a. Therefore, bubble integrals with a 1 = a 2 = 1 in four dimensions have a logarithmic divergence, while triangle and box integrals are UV finite, for example. This explains why the box integral I D=6 box considered above is UV finite.
Feynman integrals in on-shell kinematics can have another type of singularity, infrared divergences, that originate from different loop integration regions. By infrared divergences we mean both soft and collinear divergences. The former originate from integration regions where k 1, while the latter come from regions where the loop momentum becomes collinear to an on-shell external momentum, k ∼ p.
Both UV and IR divergences can be regulated using dimensional regularization.

Behavior near singular points
We will be interested in computing Feynman integrals as a function of kinematic invariants (e.g. s, t, or masses) and of the space-time dimension D = 4 − 2 . In order to know what class of functions to expect it is important to think about the asymptotic behavior of the integrals in (potentially) singular limits. A typical example is a Regge limit, where s t, where one expects in general divergences of the form s −p log q s/t. Cf. eq. (2.8) for an example.
The fact that the functions we are computing have an integral representation such as eq. (2.6) restricts the type of singular behavior that we can obtain. E.g., in the Regge limit, the leading behavior can be predicted by finding the dominant region in the space of Feynman parameter integration [10]. It is rather obvious from such an analysis that an integral will be bounded in a limit by a power with a certain exponent. This property means that we are dealing with integrals with regular singularities only (i.e., no essential singularities). Moreover, the special dependence of the Feynman parameter integral formula (2.6) suggests that the scaling exponents that can appear are linear in the space-time dimension. In fact, the analysis of regions can be used to determine the different scaling exponents [17].
In summary, we have the following important properties of Feynman integrals, • Feynman integrals only have regular singularities in the kinematic variables.
• The scaling exponents near a singularity are linear in the space-time dimension D.
In the section 4, we will analyze differential equations satisfied by Feynman integrals. We will see that these properties imply that a certain canonical form of the differential equations can be reached in an algorithmic way.

Integral families and differential equations: an invitation
In this section we explain the following concepts: Given a Feynman graph (and corresponding Feynman integral), we define a family of Feynman integrals associated to it. This family consists of, roughly speaking, all Feynman graphs with the same propagator structure, but arbitrary powers of the propagators. This includes cases with fewer propagators, i.e. subgraphs. We derive linear identities between elements of such family, which imply the notion of a basis. Finally, we explain how to derive differential equations in the external invariants for the basis integrals.

Integral families and basis
To illustrate the ideas we will proceed with the example of the one-loop box integral considered in eq. (2.2). The first step consists in generalizing it to arbitrary (integer) powers of the propagators, see Fig. 1. We recall that due to the on-shell conditions, we have y 2 12 = y 2 23 = y 2 34 = y 2 41 = 0, and the integral depends on s = y 2 13 and t = y 2 24 (and on the dimension D). For positive a i , this is a box integral, with propagators raised to general powers. If one of the a i is zero, we have a triangle integral, etc. Negative values of the a i correspond to numerator factors. We call the set of G for arbitrary integer powers of the a i an integral family (associated to the box diagram). We will see presently that this notion is useful to understand the structure of the differential equations satisfied by this integral.
The integrals in a given family are in general not independent. There are linear relations that have a very simple origin, namely integration by parts (IBP) relations [18]. These identities follow from the fact that total derivatives vanish in dimensional regularization. 2 Therefore we have e.g.
where ξ is some vector, e.g. y − y 1 . Acting with the differential operator on the rest of the integrand, and performing a little amount of algebra, one observes that the r.h.s. of eq. (3.2) can be expressed in terms of members of the integral family, albeit with different values of {a 1 , . . . a 4 }. The coefficients in this relation are rational functions of s, t, and D. Introducing the operators Y ± i , with Y ± 1 G a 1 ,a 2 ,a 3 ,a 4 = G a 1 ±1,a 2 ,a 3 ,a 4 , etc., we have the IBP relations as well as similar relations obtained by cyclic symmetry. We note that the relations are linear in the integrals G, as promised.
The IBP identities relate integrals having different indices a i . For integer values of the a i , one can devise a simple strategy to reduce any integral to a known basis [5]. This is done by noticing that the relations (3.3) can be used to reduce the value of a = 4 i=1 a i . This is done until one of the indices is zero, and then repeated for the remaining indices. In this way one sees that one can relate any integral in the family to a basis. In the present case, this basis consists of three elements, which can be chosen to be the s-and t-channel bubble integrals, and the box integral, with unit powers: G 0,1,0,1 , G 1,0,1,0 , G 1,1,1,1 .
Let us give two examples of the integral reduction.
The first relation demonstrates how higher powers in the denominators can be removed, as explained above. The second example shows that the on-shell triangle is trivially related to the bubble integral, which explains the absence of triangle integrals in the above basis. In general, the IBP relations imply the existence of a finite basis of integrals for the family under consideration. In the literature the basis integrals are often referred to as master integrals. We prefer the word basis, as it reminds one that in linear algebra, there is a freedom in the choice of basis. This will be very important in the following sections.
Here we considered a simple example, but the properties observed turn out to be completely general. Let us summarize the main points. When studying a Feynman integral, it is useful to consider the family of Feynman integrals associated to it. Integrals in this family satisfy relations following from integration-by-parts identities. The latter are linear in the integrals, and rational in the kinematic variables and the dimension D. They imply that the family has a finite basis. At one loop it is possible to solve the IBP relations directly and to write down a relation between arbitrary integrals and the chosen set of basis integrals. At higher loops, such a general solution is not known, so that for each case under consideration one usually constructs a sufficient number of IBP relations. Although straightforward in principle, this is rather tedious to do by hand, so that one usually employs some appropriate computer algebra program. Various implementations exist, see [19][20][21][22].

Differential equations
One result of the previous section was that for a given family of Feynman integrals there exists a finite-dimensional basis f . In practice, the latter is found by writing down a sufficient number of IBP relations. Given such a basis, any integral in the family can be written in terms of a linear combination of basis integrals, with rational prefactors in the kinematic variables and D. Therefore it is sufficient to compute the basis integrals. We will use differential equations in the kinematic variables to that end.
In the example of the one-loop box integral family, the kinematical variables are s and t. We implement differential operators ∂ s and ∂ t acting on the integral representation (3.1), which depends on the vectors y i (or, equivalently, on the p i ), via the chain rule. When doing so one has to make sure that the differential operators commute with the on-shell and momentum conservation constraints.
For an analogy, think of a two-dimensional space with constraint x 2 +y 2 = 1, i.e. points lying on the unit circle. In that case, only the operator y∂ x − x∂ y is allowed. Of course, one could also introduce radial coordinates x = r cos θ, y = r sin θ, so that the constraint becomes r 2 = 1, which means that one can freely vary θ. For on-shell massless scattering amplitudes, likewise, one can solve the momentum conservation and on-shell constraints using momentum twistor variables [23]. Here we will use the momentum-space variables and construct differential operators commuting with the constraints.
Example: differential operators Let us use the momentum-space notation for the family of one-loop box integrals, and eliminate p 4 using momentum conservation, Let us construct a differential operator for ∂ s that can act on the r.h.s. of this equation. This is easily achieved by making the ansatz Imposing that this operator should commute with the on-shell conditions p 2 1 = 0 and (p 1 + p 2 + p 3 ) 2 = 0, and imposing the normalization condition ∂ s (p 1 + p 2 ) 2 = 1 fixes the parameters in eq. (3.7) to be . (3.8) When acting with such differential operators on the Feynman integral representation, we have to perform algebraic manipulations similar to those when deriving the IBP relations.
It is clear that one obtains integrals within the same family of integrals. The fact that there is a basis means that we can rewrite the result of the differentiation as a linear combination of basis integrals. In other words, we have where A s and A t are N by N matrices, with N being the number of basis integrals f . By construction, they contain only rational functions of s, t, as entries.
In other words, Feynman integrals satisfy first-order systems of (partial) differential equations. The matrices A i can be computed algorithmically, as outlined in this section.
Example: Differential equations for the family of one-loop 2 → 2 integrals. We already saw that in this example there are three basis integrals. Integral reduction suggests the following basis choice, With this choice, we find the following matrices in eq. (3.9), We can make the following observations.
• Computing sA s +tA t = diag(− , − , −2− ), the scaling dimensions of the integrals are correctly reproduced. We can set them to zero by choosing appropriate dimensional normalization factors, so that we only have one non-trivial variables x = t/s.
• The equations for the bubble integrals f 1 and f 2 are trivial, and indeed being singlescale integrals, their functional dependence follows from dimensional analysis.
• The equations have the singular points s = 0, t = 0, s = ∞, t = ∞, and s = −t (i.e. u = 0). The latter singularity may be surprising for planar integrals, and as we will see occurs only after analytic continuation.
As a preview of the general method to be discussed in the following sections, let us make the following educated basis choice (to be justified later), with c = e γ E being a normalization factor, and with γ E being Euler's constant. The g i are chosen to be dimensionless, such that they depend on x and only. Implementing the derivative ∂ s as explained above, and using the chain rule, we find The system (3.14) can be solved easily in an expansion in . One sets 16) and plugging this into eq. (3.14) it becomes clear that at each order in , the r.h.s. of that equation is known and can be integrated. Let us discuss the boundary conditions for the equations. As already discussed, the bubble integrals are trivially known: a short calculation using the formulas of section 2 shows that they are given by with a = a 1 + a 2 . In application to our case, we have Finally, we need a boundary condition for g 3 . We can use the fact that planar integrals should not have u-channel singularities, which implies that g 3 should stay finite as x → −1, despite the presence of the matrix b in eq. (3.14). This fixes the solution to all orders in the expansion. The first few orders are given by where Li n is a polylogarithm, defined by and Li n (0) = 0. In section 5 we will discuss a more general class of functions that is useful for writing the solutions to such differential equations.
We note that one can associate a notion of weight to the above functions, corresponding to the number of integrations. For example Li n has weight n. Moreover, if one also associates weight −1 to , then one observes that the g i have uniform weight 0 for all terms in the expansion.
The main conclusion of this section is that the basis integrals satisfy a linear system of differential equations with rational coefficients. This follows simply from the existence of a basis, and the way the derivative operators act on integrals of a given family.

Singular points of the differential equations
In section 2 we observed general properties of Feynman integrals near their singular points. Let us study the one-loop box example from this point of view. We see that eq. (3.14) has the singular points x = −1, 0, ∞. (The analysis of the point at infinity can be treated as the other ones after performing an inversion.) What is their physical meaning? It is easy to see that they correspond to the limits u = −s − t → 0, t → 0, and s → 0, respectively.
The differential equations tell us how the solutions behave near those points. Let us consider the limit x → 0. Keeping only the leading term on the r.h.s., we find the solution where f 0 ( ) is a boundary vector. The matrix exponential evaluates to This illustrates the statements at the end of section 2. The solutions are linear combinations of different terms x α , where α are linear in the dimension. We note that the α are the eigenvalues of the matrix multiplying the singular point, in this case a. This example illustrates that in general there are two sources of logarithms in the expansion around a singular point. The first is the matrix exponential itself, c.f. eq. (3.22), and the second is the expansion for small of exponentials such as x − .
In general, the structure of singularities will not always be as manifest as in this example. To understand this better, consider the scalar differential equation which is more singular at x = 0 compared to the previous case. Indeed, the solution f (x) = e −a/x f 0 has an essential singularity at x = 0. Such a function cannot appear in individual Feynman integrals. Can one conclude therefore that the differential equations contain simple poles only? Unfortunately, this conclusion would be premature, since the matrix nature of the equations allows for 'spurious' terms to occur. This fact can be seen from the following example, (3.24) Here the 1/x 2 term is spurious. It can be removed by a simple change of basis, which leads to The reason the singularity structure was manifest in the example of the one-loop box integral crucially had to do with the basis choice made above, namely eq. (3.13) vs. eq. (3.11). Obviously, the two basis choices are related by a transformation T . The next sections are devoted to a better understanding of this point. In particular, section 4 discusses this from the point of view of the DE, while section 6 relates it to properties of the loop integrands.

Algebraic simplifications of differential equations
Here we discuss how the expected singularity structure can be made manifest, leading to a canonical form of the differential equations. In the physics literature, some of these ideas were presented e.g. in [9,[24][25][26][27]. We also discuss how iterated integrals, and more complicated functions such as elliptic functions occur.
We saw in the previous section that integrals for families of Feynman integrals satisfy first-order coupled systems of differential equations. Here we will focus on the case where the integrals depend on one kinematical variable, x, for simplicity. This is the case of the example discussed above, since one can always normalize the integrals to remove one of the scales. Then, for some choice of basis f we have where A is an N × N matrix. From the structure of the IBP relations it follows that A depends on x and in a rational way 3 .

Simplifying the dependence on x
The singularities of the differential equations (4.1) have to correspond to singularities of the original Feynman integrals. If we denote the singular points by x k , we expect the leading behavior of the integrals to be given by terms that grow like ∼ (x − x k ) α , for some values of α. This is to say that (4.1) should only have regular singularities, i.e. it should be a Fuchsian system of differential equations.
To specify what this means in terms of the system of DE (4.1), we first have to recall that there is a gauge degree of freedom in the choice of the basis. Indeed, starting from eq. (4.1), we can switch basis according to for some invertible matrix T , which transforms (4.1) into an equivalent system, Let us now investigate the DE near a singular point. We take x = 0 without loss of generality. The matrix A in eq. (4.1) has the expansion 4 for some value of p. The system of DE is regular singular in x = 0 if there exists some gauge transformation T , for which the matrix A appearing in the DE has p ≤ 1 (if p < 1 the solution has no singularity at x = 0). In other words, the fact that (4.1) is a Fuchsian system implies that for each singular point, one can find a gauge transformation T such that the gauge equivalent system has a matrix which has the leading behavior i.e. with p ≥ 1 in eq. (4.5). As a consequence, near each singular point the solution behaves as x A 0 ( ) . The degree of singularity of a system of differential equations was studied by Moser [28,29]. It was shown that under certain conditions on the two leading matrices in the expansion (4.5) near a singular point, the order of the singular term can be reduced. It is important to note that the necessary transformation is rational in x, see e.g. [30,31]. This implies that removing spurious singularities at one singular point does not influence the behavior at other points (except possibly at infinity). Therefore one can algorithmically construct a rational matrix T such that the DE system reads where p(x, ) is polynomial in x. If p(x, ) = 0, this means there is still an undesired spurious singularity at infinity. The question whether p(x, ) can be removed without introducing further singular points is related to the Riemann-Hilbert problem, see e.g. [32]. It is an interesting question how to decide whether this is possible and to construct such a transformation matrix in the positive case. A more pragmatic solution consists in introducing another singular point, not yet present in the list {x k }, to "balance" the transformation at infinity. In this way, we can use the above algorithm to obtain the form (with different matrices compared to eq. (4 .7)) This is a system with manifestly only regular singularities.

Simplifying the dependence on
We have just explained how for Feynman integrals one can put the differential equations in the form of eq. (4.8), using the algorithms of refs. [28][29][30]. (We restricted ourselves to the single variable case for the sake of the presentation, but similar results can be obtained for the multi-variables case.) This equation makes the behavior of f near the singular points in x manifest. Indeed, the solution takes the form (we again take x = 0 without loss of generality), where f 0 ( ) is a boundary vector, independent of x, at x → 0, and is a matrix polynomial in x, whose expansion coefficients P m ( ) can be determined recursively from the information in eq. (4.8) [29].
In other words, this form already contains all the information about the scaling behavior of the integrals near the singular points. In particular, as discussed in section 2, the eigenvalues of a k ( ) are related to different regions in integration space, and one expects them to be linear in .
This raises the question whether the dependence on in eq. (4.8) can be simplified. By construction, we know that only appears in a rational way. Moreover, it is clear that poles in in a k must be spurious, and these can be removed similarly to the removal of spurious divergences in x, see e.g. [29,33].
For a polynomial dependence on , we can distinguish between two cases: if the r.h.s. of eq. (4.7) is O( ), the solution at each order in can be obtained in terms of iterated integrals. If the r.h.s. starts at order 0 , the solution may be more complicated.
So we see that a crucial question is whether we can construct a transformation that removes the 0 part of the matrix on the r.h.s. of eq. (4. 8), and what the nature of this transformation is. This is best explained via a few examples.
• Case 1: integrating out the 0 term amounts to choosing a rational normalization factor.
As an example, imagine choosing a normalization factor s 2 instead of s t for g 3 in eq. (3.13). In that case, one obtains a 0 in the DE. The latter can be removed (by construction) by a simple rational transformation T .
• Case 2: integrating out the 0 term can be done using algebraic functions; sometimes a change of variables leads to a rational dependence. This is something that occurs typically for integrals involving masses. As an example let us choose the 2 × 2 system for a massive bubble and tadpole integral [34,35]. The bubble integral depends on an external invariant s, as well as on an internal mass m.
We set s = x and m 2 = 1 without loss of generality. Before choosing appropriate normalization factors, the differential equation in x reads Here, integrating out the constant term in amounts to choosing a transformation matrix T = diag(1, 1/ 1 − 1/x). Note that contrary to all transformations discussed so far, this transformation is not rational (in the chosen variables). Under f −→ T f , the system of DE becomes (4.12) The r.h.s. is ∝ , as promised. Finally, we note that in this case, one can recover a rational form of the equations by employing the change of variables x = −(1 − y) 2 /y, with the resulting system having regular singularities in y = ±1, 0, ∞.
• Case 3: integrating out the 0 term leads to elliptic or more complicated functions.
A simple example of this kind is the system of DE satisfied by complete elliptic integrals. The difference w.r.t. the previous example lies in the fact that this is a coupled 2 × 2 system of differential equations, even at = 0. The fact that the appearance of elliptic functions can be seen in this way was also pointed out in the algorithm of ref. [25].
For Feynman integrals, the simplest example where this can occur is perhaps the twodimensional two-loop sunrise integral with equal masses, see e.g. [36][37][38][39][40]. There are singular points at x = 0, −1, −1/9, ∞, where x = m 2 /(−p 2 ). Making an appropriate basis choice, one obtains a system of (4.13) One may verify that all eigenvalues are linear in , as expected. Analyzing this system at = 0, and writing it as a second-order eq. for one of the integrals, one recovers the Picard-Fuchs equation discussed e.g. in refs. [38,39]. In other words, integrating out the 0 piece of the equation leads to elliptic functions. This could obviously be generalized to cases where more than two equations are coupled at = 0, leading to higher-order Picard-Fuchs equations for the individual integrals.
Note added: After these lectures were given, ref. [41] appeared, which has some overlap with the material presented in this section. This reference proposes a refined version of balancing the transformations of the Moser algorithm, and, for the case of a Fuchsian system with eigenvalues already normalized to O( ), it gives a formula for constructing a transformation matrix to simplify the dependence of the differential equations. The methods presented here and in ref. [41] allow for a systematic, albeit not fully algorithmic, simplification of differential equations. Challenges for an algorithmic implementation include reaching a Fuchsian form without spurious singularities, i.e. the step from eq. (4.7) to eq. (4.8), and dealing with non-rational dependence on variables. As we have seen above, the latter can appear when integrating out 0 terms, or, equivalently, when normalizing eigenvalues [41]. Another open problem is the extension to multi-variable cases. In section 6, we present an alternative method for constructing an integral basis. This method, based on insights from the singularity structure of Feynman integrands, is very simple to use and has already been successfully applied to complicated multi-variable cases with algebraic dependence on the variables. Also, it typically gives much simpler expressions for the integral basis. Of course, the ideas of this section and section 6 are complementary, and having both at one's disposal allows one to solve very complicated problems.

Iterated integrals from differential equations
In the previous section we discussed how differential equations can be simplified systematically both in x and in in an algebraic way. Before discussing a completely orthogonal approach to this, that makes direct use of properties of the Feynman loop integrand, we wish to discuss properties of the solutions of the differential equations.

Canonical form of differential equations
Here we will discuss the case where the differential equations can be put into the form suggested in ref. [9]. Let us start with the case of one non-trivial variables x. The form of the equations proposed in ref. [9] reads In the case where this representation can be reached using only rational transformations (for some choice of variables, see the examples at the end of the previous section), the α k are rational functions, i.e.
where the x k are the locations of the singularities. In the general case they can depend algebraically on x.
The generalization to the multi-variable case is straightforward, we can simply replace x by x in eqs. (5.1) and (5.2).
The canonical form of eqs. (5.1) and (5.2) has the virtue that the information about the functions f is encoded in a minimal way: the alphabet α specifies which class of generalized functions will be required for writing down the answer. As we will see presently, the coefficient matrices for each letter determine which linear combination of those functions is required. Indeed, the general solution to eqs. (5.1) and (5.2) in the multi-variable case can be written as Chen iterated integrals [42] f where P stands for path ordering along the integration contour γ, and f 0 ( ) is a boundary value. This formula is to be understood in an expansion in , where the k-th term in the expansion is a k-fold iterated integral (along γ). Let us be more specific about the notation, following closely the recent lecture notes [25,43,44] on iterated integrals. We denote by M the space of kinematical variables, and let ω i be some differential one-forms (corresponding to entries of dÃ). Moreover, define the pull-back of the differential forms to the unit interval [0, 1] via Then, an ordinary line integral is given by The iterated integral of ω 1 , . . . ω n along γ is defined by Iterated integrals have many nice properties, see [43,44]. Moreover, the iterated integrals appearing in eq. (5.4) are homotopy invariant (on M with the set of singularities removed). This property makes them very flexible, and we will see later that how to rewrite them in terms of more familiar functions, if desired. The freedom in choosing the integration contour γ is also useful e.g. for analytic continuation, or for writing integral representations that have certain desired properties. It is important to emphasize the simplicity of the solution (5.4). Recall the notion of weight introduced earlier, corresponding to the number of iterated integrations. It is obvious from eq. (5.4) that each term in the -expansion is a Q-linear combination of iterated integrals of the same weight. This is a remarkable property that is not true in a generic integral basis, where results would look far more complicated, with terms of different weights being mixed, and where prefactors are in general algebraic functions of the kinematic variables. Finally, note that if is assigned weight −1 (and defining the weight of a product as the sum of individual weights), then one can say that (5.4) has uniform weight zero.
These remarkable simple properties will in fact be a guiding principle for finding an appropriate integral basis in section 6. As we will see there, it is possible to anticipate these properties by inspection of the Feynman integrand, i.e. even before carrying out any integrations.

Fixing the boundary conditions
The system of first-order differential equations (5.1), (5.2) determines the answer up to an integration constant. In principle, the latter can be fixed by an independent (and simpler) calculation at a preferred kinematic point. However, in practice the differential equations themselves, together with insight coming from the original Feynman integral representation, allow one to determine boundary conditions without calculation.
One idea is that Feynman integrals are multi-valued functions, and not all of the singularities present in the differential equations can appear on the first sheet of the functions. For example, for planar integrals, there cannot be a u-channel singularity or branch cut, and this information was used in ref. [45] to determine all non-trivial boundary constants (the only integrals that had to be computed were integrals that trivially evaluated to Gamma functions). We have seen this in the example considered in section 2. In [46] it was shown that the same idea can be used for non-planar integrals, at the cost of introducing temporarily one additional scale.
The idea of introducing an additional variable, together with the idea of "bootstrapping" the boundary information was further used in [46] to compute single-scale integrals. Here the crucial observation was that the differential equations provide the exact form of the answer, e.g. x a 0 f 0 ( ) in the x → 0 limit, valid for any . This is important since the matrix exponential contains in general terms for which the x → 0 and the → 0 limit do not commute. The knowledge of this exact term allows one to translate between the two orders of the limits. Moreover, terms with different scaling behavior can be cleanly separated. We give a pedagogical example of this method in section 7. We note that very similar ideas were also recently applied in [47].
In the next subsections we show various examples of alphabets appearing in practice, and comment on different ways of representing the answer.

Examples of function alphabets
Let us give some examples of alphabets that appear in practice in the computation of Feynman integrals: • For on-shell massless four-point integrals to three loops [7,45,46] only the letters {x, 1 + x} are needed.
• For one-variable functions, the alphabet {x, 1 − x, 1 + x} often makes an appearance. This is the case e.g. for vector boson fusion via a top quark loop to two loops [48] (at three loops, additional letters appear), Wilson line integrals [49,50] forming a cusp to three loops, and certain contributions to Higgs cross sections at N 3 LO [27,47].
Cases with several mass scales typically have a larger alphabet of the order of 10 to 20 letters, see e.g. [25,26,[56][57][58]. In general, there can be also polynomials of higher degree, e.g. 1 + x 3 [59], or even algebraic dependence on the kinematic variables, as already mentioned in one example. It should be stressed that the statements about the alphabets occurring in certain types of integrals are only firmly established up to a given loop order, and in general the alphabet is not stable when increasing the loop order, i.e. additional letters appear, or one may even leave the space of iterated integrals.
Many cases are related to a sphere with n marked points [43]. It is in general a difficult problem to know whether two alphabets are related because of the freedom to redefine variables.
It is interesting to note that some alphabets have been observed to be related to cluster algebras [60][61][62].

Iterated integrals and hyperlogarithms
The solution to eq. (5.1) is given by Chen iterated integrals [42]. In the case where the alphabet can be written in terms of rational functions (in at least one variable), one can write the answer in terms of Goncharov polylogarithms (also called hyperlogarithms). These functions go back to Lappo-Danilevsky [63], who introduced them precisely as solutions to differential equations of the kind that we are discussing. The Goncharov polylogarithms can be defined iteratively as follows, G(a 1 , . . . a n ; z) = z 0 dt t − a 1 G(a 2 , . . . , a n ; t) , For a 1 = 0, we have G( 0 n ; z) = 1/n! log n (z). See refs. [43,64] and the lecture notes [44,65] for a discussion of their properties. For a given class of Feynman integrals, one will need only a subset of allowed indices a i . One case that appears particular often in practice is that corresponding to the alphabet {x, 1 − x, 1 + x}, i.e. indices drawn from 0, ±1. This case is referred to as harmonic polylogarithms (HPL) in the physics literature [66,67]. They are denoted by H(a 1 , . . . a n ; x). 5 The example considered in section 2 can be solved, to any order in , in terms of these functions.

Different representations of the answer
There are various ways in which the answer to the differential equations can be written. As already discussed, it is always possible to write the answer in terms of Chen iterated integrals. If the function alphabet is rational in at least one variable, the latter can be represented in terms of Goncharov polylogarithms.
In practice, we are often not interested in the whole expansion, but can truncate the expansion at a given order. For example, for calculations at next-to-next-to leading order in perturbation theory, we are typically interested in the answer only up to weight four. In such a case, one can try to rewrite the answer in terms of a minimal function basis. According to a conjecture by Goncharov, all weight four functions can be written in a basis spanned by the following functions (for certain arguments, to be determined), {log x log y log z log w, log x log yLi 2 (z), Li 2 (x)Li 2 (y), log xLi 3 (y), Li 4 (x), Li 2,2 (x, y)} , (5.10) where The rewriting in terms of a minimal function basis can be done with the help of projections onto different terms in (5.10). Moreover, the necessary arguments in that equation can be found in a systematic way. For these purposes, the notion of the "symbol" of iterated integrals is useful [51,68]. We note that in the differential approach, the symbol of the answer is entirely manifest, as it is encoded in the matrixÃ of eq. (5.1). We will illustrate the different representations using the example of a one-loop on-shell box integral [25], called g 6 there, which depends on the variables u = 4m 2 /s and v = 4m 2 /t. This integral is a weight 2 function. As already discussed, it is always possible to write the answer in terms of Chen iterated integrals, in this case, The boundary condition is that g 6 vanishes as u, v → ∞. The monodromy invariance allows many choices of γ, but a particularly simple parametrization is (u(t), v(t)) = (u/t, v/t), with t ∈ [0, 1].
In oder find a representation in terms of Goncharov polylogarithms, we first need to find a change of variables that rationalizes the function alphabet. This is achieved by setting (5.14) which allows us to write (This particular case only requires Goncharov polylogarithms with indices 0, ±1, so they could be replaced by HPL.) Finally, we can also write the answer in terms of the weight two minimal function basis, namely {log x log y, Li 2 (x)}, for certain arguments. In the present case, we have [69] Let us discuss advantages and disadvantages of the various representations.
• The Chen iterated integral formulation is usually the most compact way of writing the answer. Its monodromy invariance makes it very flexible, which is of great use e.g. when computing limits, or for analytic continuation. Examples for numerical integration can be found in [25]. Finally, we stress that the "symbol" of the answer is completely manifest.
• The Goncharov polylogarithm version of the answer amounts to fixing a specific integration contour. This usually leads to a proliferation in the number of terms, and obviously the monodromy invariance is no longer manifest, which makes e.g. analytic continuation more difficult. It should also be noted that this representation is by no means unique, since e.g. different choices of integration contour lead in general to different expressions in terms of Goncharov polylogarithm. On the positive side, there exist dedicated numerical integration routines for these functions [70,71].
• The representation in terms of a minimal function basis also usually leads to a larger number of terms (compared to the Goncharov polylogarithms, one trades simple arguments of the functions for simple indices); the representation is not unique due to identities of functions involving different arguments, so that it is difficult to find an "optimal" representation; However, since the functions involved are well studied, this representation may be advantageous for numerical evaluation. Finally, in some cases, e.g. when expressing a final result that is expected to have simple properties, this might be made manifest by an astute choice of function arguments. 6 Finding an optimal basis using d-log forms and generalized unitarity cuts In the previous section we discussed an algebraic approach to simplifying the system of differential equations. Here we discuss another (in our opinion, more natural) approach that is based on exploiting the original integral representation, and d-log representations.
The canonical form of the differential equations in section 5 lead to answers that have two important properties: • They are given by iterated integrals of uniform weight, i.e. where the k-th coefficient in the -expansion has weight k. Many examples of individual Feynman integrals with these properties (at least, up to some order in the expansion) were previously known in N = 4 SYM. However, even understanding the transcendental weight of the leading term in the expansion of a given integral is not immediately obvious. For example, at the one-loop order it is known that four-dimensional Feynman integrals give at most weight two functions. However, in momentum space the number of initial integrations is higher (four), and likewise in the Feynman representation, the number of integrations is proportional to the number of propagators minus one, so that in general neither of these representations makes the weight properties of the answer manifest.

d-log representations
Sometimes one can bring the integral representation to a form where the weight properties of the answer are manifest. This is very desirable, since such an understanding usually implies also an algorithmic way of computing the answer.
Let us discuss as an example the Wilson line integral shown in Fig. 2, following [49]. (This is also relevant to scattering integrals, since the latter are in some cases dual to Wilson line integrals.) The integral over the line parameters s and t can be written as where on the r.h.s. we have dropped differentials involving dx because they do not contribute to the integral, and where the integration region Λ is s ∈ [a, S] and t ∈ [b, T ]. In the Wilson line language, the weight properties are easy to see: at one loop, one has two line integrations to carry out, and the result is a weight two function, as expected . This is a first example of a d-log representation. It has the virtue of making the fact that the answer is a Q-linear combination of weight k functions manifest, where k is the number of integrations. This can be seen algorithmically, and the same algorithm in fact computes the function [49]. We comment that additional factors, such as X , that can occur in dimensional regularization, do not change this conclusion.
Sometimes one can find similar d-log representations by manipulating the Feynman representation. Let us take a one-loop s-channel bubble integral as an example, with integral mass m. For general powers a 1 , a 2 of the propagators, it is given by We see that for a 1 = 1, a 2 = 2 and D = 4 − 2 , we have a d-log representation (up to some normalization factor), with some additional, but inconsequential factor (. . .) . Inspecting further the Gamma functions, one can easily prove that upon including a normalization factor, e.g. e γ E , this yields a uniform weight one function. This method of identifying uniform weight functions works in practice in many cases, especially for cases with bubble and/or triangle (sub)integrals.

Generalized cuts
In the previous section we discussed ideas to prove that a given integral has the desired uniform weight properties discussed above. In practice, sometimes a necessary condition can also be very valuable.
The idea we want to pursue is to consider discontinuities across branch cuts of the functions under consideration. It is rather obvious that a uniform weight function will yield another uniform weight function under this operation. On the other hand, discontinuities of Feynman integrals are sometimes easier to compute than the integrals themselves. Consider the one-loop on-shell box integral for example. If we apply both an s-and t-channel discontinuity, we can replace its four propagators by four delta functions. This localizes the four-dimensional integration. Computing the Jacobian, (and ignoring the -dimensional part of the integration), we obtain 1/s/t as the answer for the cut integral. This explains the normalization factor st in eq. (3.13) of the uniform weight basis discussed there.
We now generalize the above idea from unitarity cuts to generalized cuts, where any number of propagators can be cut (replaced by delta functions). Such integrals satisfy the same IBP relations and hence the same differential equations as the original functions, albeit with different boundary conditions. Roughly speaking, thinking about a system of DE, each way of choosing cuts projects onto a smaller, but nontrivial subsystem. If we are looking for a canonical form of the full system, the same canonical form has to be present also at the level of the generalized cuts.
Therefore, we can impose as a criterion for a putative basis of uniform weight functions that any generalized cut (i.e. replacing any number of propagators by delta functions) should give a pure uniform weight function. This is a very powerful test that can be easily performed. In particular, as we will see, at higher loops one can "recycle" the knowledge obtained previously from lower loops. Let us give two examples.
Example: One-loop pentagon integral. Consider the one-loop pentagon integral with some numerator N (k), There are five possibilities for computing a maximal cut (cutting four propagators). The answer on a given cut will be of the form where k * is a solution to the cut conditions, J is the Jacobian, and P is the uncut propagator factor, e.g. P (k) = (k − p 5 ) 2 , if the first four propagators are cut. One may verify that for a constant (i.e. k-independent) N , there is no choice of N such that one obtains a Q number on each cut. Therefore the scalar pentagon integral in four dimensions is not a pure uniform weight integral. We prefer not to use it in our basis. However, it is now clear how to proceed. We can allow a more general numerator N (k), and impose the condition that all cuts should yield numbers in Q only. We refer to reference [13] where this philosophy was introduced and many integrals with these properties were constructed. This reference also gives an introduction to momentum twistors that allow one to visualize more easily the geometry of the cuts, and in terms of which the cut solutions are simple.
Example: Two-loop planar and non-planar box integrals. For a multi-loop example, we may consider a massless double box integral. Here we wish to illustrate how the knowledge of lower-loop information can be recycled. Consider an ansatz of the form where p 12 = p 1 + p 2 . Consider cutting all propagators of the k 2 integral, as well as two adjacent propagators. Localizing the k 2 integration we obtain a Jacobian factor of 1/s/(k 1 − p 4 ) 2 , so that we are left with (ignoring the -dimensional part of the integration) (6.6) We recognize this as a one-loop box integral, with a numerator N (k 1 ). Our knowledge of one-loop integrals now tells us that choosing N (k 1 ) = s 2 t makes this a uniform weight box integral. However, there is a second possibility: choosing N (k 1 ) = s 2 (k 1 − p 4 ) 2 , we obtain a triangle integral that is also a uniform weight function (since triangle and box integrand simply differ by the fact that for the triangle one of the points is at infinity). This simple analysis explains the choice of integral basis for the double box integrals in ref. [9]. For a non-planar example, consider the two-loop non-planar double box integral. It contains both a pentagon and a box subintegral. Let us consider the generalized cut that localizes the box subintegral. The Jacobian 1/(P 1 P 2 ) gives two propagator factors that depend on the remaining loop momentum. This means that the cut integral contains a scalar pentagon integral. We have seen above that this integral does not have the desired uniform weight properties that we are looking for. However, we already know how to remedy this problem, by using as the starting point a non-planar double box with appropriate numerator factors. In complete analogy to the analysis above, suitable candidate numerator factors are P 1 , P 2 , or P 1 P 2 , multiplied by an appropriate overall normalization.
Example: Two-loop 2 → 2 integrals of ref. [9]. As an application to the discussion of this and the last subsection we are now in a position to fully understand the basis choice made e.g. in ref. [9], see. Fig. (3). The propagator-type integrals can all be directly evaluated in terms of (iterations of the one-loop) formula (3.17). It is easy to see that for D = 4 − 2 dimensions, this formula gives results of homogeneous transcendental weight for the choice of doubled lines as shown in Fig. 3.
Let us now discuss the remaining integrals shown in that figure. Integrals having bubble sub-integrals are also immediate: integrating out the bubble integral yields a Feynman integral already known from the analysis of the previous loop order [45]. E.g., the first integral in the second line is related to a box integral. The latter is known to be of uniform weight, with known normalization factor (the latter can be determined e.g. through a maximal cut). Similarly, the third integral in the first line is related to a one-loop triangle. Triangle integrals in four dimensions have the same properties as box integrals. This can be seen by employing momentum twistor variables [13,23], or in the embedding space formalism [72]. Roughly speaking, the triangle is a box with one point at infinity.
The second integral in the second line can be understood in two different ways. One possibility is to introduce one Feynman parameter to combine two propagators adjacent to an incoming on-shell momentum. In this way, one obtains a one-parameter integral (with a d-logarithmic integration kernel) over an integral it is easily seen to be of uniform weight, in complete analogy to what was discussed above. Another possibility is to use the embedding formalism, in which case the cut analysis as explained above is straightforward.
Finally, for the most complicated, seven-propagator sector, the basis choice was explained above using cuts. One numerator is just the overall normalization N = s 2 t, while the other one contained a loop momentum dependence, s 2 (k 1 − p 4 ) 2 , which is indicated a dashed line in the figure.
This explains how to choose a basis of uniform weight integrals for this problem. The critical reader may wonder about the absolute weight for the individual basis elements. The latter are indeed different: for example, the double box integrals have weight 4, whereas e.g. the double bubble integrals as defined here have weight 2. This overall difference in weight is compensated by powers of included in the definition of the basis in ref. [45]; this final adjustment leads to the canonical form quoted there.
We leave it to the interested reader as an exercise to reproduce the basis choice for higher-loop integrals [45], or for partially off-shell kinematics [46,57], or cases with masses [25]. Further explicit examples can be found in section 7.
We note that this method typically gives a more compact form of the basis compared to the algebraic method of section 4, and it is more likely to make physically important properties, such as e.g. infrared behavior, manifest.
We close this section with a number of remarks.
• While the cut analysis in this section only provides consistency checks for an integral to be of uniform weight, we wish to emphasize that this hypothesis can be immediately tested (and proved!) by deriving the system of differential equations for it.
• In practice, one often finds that after this analysis, the system obtained is close to the desired canonical form. It may happen that small transformations, in the spirit of section 4, are required to fix e.g. terms that vanish on the cuts considered. This combined analysis is what we have found most useful in practice.
• The geometry of the cut equations for massless on-shell Feynman integrals is particularly transparent and simple when using momentum twistor variables [13,23]. As was already mentioned, the same variables, or also the embedding formalism [72], also have the virtue of making the point at infinity manifest, so that triangle integrals and box integrals are in fact treated on the same footing.
• We also wish to mention that the above criteria could be tested algorithmically, for a given graph. In practice, we have not found this necessary since the application by hand is very simple.
• As we have seen in the examples, the cuts relate different loop orders, and one can therefore "recycle" information about an optimal basis of lower loops when going to the next loop order; also, once the behavior of certain subgraphs is understood, they can be used at higher loops as "lego blocks", i.e. without having to revisit the justification for using them.
• We also mentioned that cuts very naturally project the DE onto subsectors; this allows one to gradually construct the full system of DE, dealing with smaller matrices only at a given time. This can be of practical importance, since it means that one has to set up the entire N × N system only when the optimal basis is essentially known, and since the DE are much simpler in the optimal basis, this leads to a significant speed-up in the algebraic manipulations needed.
• Finally, it is important to point out that the initial motivation for considering integrals of this type is to make infrared properties of scattering amplitudes especially simple by choosing basis integrals that are either infrared finite or "almost" finite. This is related to generalized cuts since the latter can probe regions of loop integration that produce infrared divergences. It was observed that choosing integrals defined in this way made the answers particularly simple, even before carrying out the loop integration [12][13][14].

Momentum-space d-log representations
We wish to mention that there exists another type of d-log representation that is closely related to the cut properties discussed in the last section. Namely, it was observed that certain loop integrals/integrands can be algebraically put into a momentum-space d-log form.
For example, for the integrand of the one-loop box integral this representation is [73] where k * is one of the quadruple cuts solutions of the box. More examples of d-log representations of this kind, for example for the planar and non-planar double box integrals discussed in the previous subsection are given in ref. [73]. One virtue of this representation is that it is straightforward to take generalized cuts, since no Jacobians have to be computed. See also the discussion in [60] regarding the expected transcendental weight properties of the integrated result.
For progress in the direction of a direct integration in momentum space starting from such expressions see refs. [74,75].  Figure 4. Non-planar massless form factor integral. q 2 1 = 0, q 2 3 = 0, and we will consider both q 2 2 = 0 and q 2 2 = 0.
7 Bootstrapping single-scale integrals using the Drinfeld associator In this section we give a simple pedagogical example of the application of differential equations to single-scale Feynman integrals. We explain how to find a basis of integrals of uniform weight in practice, and discuss the structure of the differential equations, placing particular emphasis on the singular points and asymptotic limits. This is based on refs. [9,46]. Naively, one cannot use DE to compute single-scale integrals, since their scale dependence is trivial. The main idea is to introduce a parameter x in a natural way, so that for x = 0 the original integrals are recovered. The system of DE w.r.t. x turns out to have another singular point at x = 1, i.e.
The point x = 1 turns out to yield a simple boundary condition, without computation, similar to the example considered in section 2. One can then transport this boundary information back to x = 0, via the DE. This is precisely what is computed by the Drinfeld associator. The latter depends only on the matrices a and b, and we can compute its expansion in , and hence the original single-scale integrals, to any desired order in . 6 Moreover, it follows from eq. (7.1) that the associator contains only multiple zeta values. Consider the family of massless non-planar form factor integrals at two loops, The recent ref. [76] also uses the idea of introducing an additional parameter that is expected to have simple boundary values.
where q µ 1 +q µ 2 +q µ 3 = 0 and a 5 ≤ 0. We are interested in this for the on-shell case q 2 2 = q 2 3 = 0. There is one non-trivial crossed ladder integral, see Fig. 4, that we would like to evaluate using differential equations. In order to do this, we introduce another scale by letting q 2 2 = 0. In that case, there are seven basis integrals. They are shown in Figs. 4 and 5. The integrals now depend on two scales, q 2 1 and q 2 2 , and we can derive differential equations in x = q 2 2 /q 2 1 . On general grounds it is known that we will obtain a first-oder system of differential equations for the basis integrals. In the following, we will explain how to simplify that system by making an appropriate choice of basis integrals.

Choosing basis integrals that have uniform weight
We would like to choose basis integrals that have uniform weight, i.e. that are built from linear combinations of functions of a given weight. We will call such functions UT (loosely for uniform transcendental degree). There are several ideas for how to find such functions. Let us apply a few of them in the present case.
and we indeed see that each term in the expansion has uniform weight. This is already obvious from the first line of eq. (7.6), if one takes into account that The situation is more interesting for the triangle integrals shown in the second line of Fig. 5. Consider f 4 for example. We already know that we can always integrate out a massless bubble insertion. The resulting integral will be a triangle integral, where one line has a shifted power a 3 + a 4 − D/2. We will argue that the choice a 1 = a 2 = a 3 = 1, a 4 = 2 is a good one. This can be see e.g. from the Feynman representation, which reads (7.8) We will see that this integral is UT, without having to evaluate it. The main point is that its integrand is a rational form with logarithmic singularities. Solving the δ function e.g. by α 1 = y, α 2 = (1 − y)z, α 3 = (1 − y)(1 − z), we see that the y integration is elementary, and the remaining integral becomes . From this is is clear that there is a unique normalization factor, q 2 2 − q 2 1 . Moreover, upon integration, the function at = 0 will have weight 1. It is easy to convince oneself that expanding to k will increase the weight of the function by k.
Hence we conclude that is a UT basis integral (the other factors were chosen for later convenience). In the same way, one arrives at the choices f 5 = 3 (−q 2 1 ) 2 (−q 2 1 + q 2 2 )G 0,2,1,1,0,1,0,0,0 , (7.11) f 6 = 4 (−q 2 1 ) 2 (−q 2 1 + q 2 2 )G 1,1,1,1,0,1,0,0,0 . (7.12) Finally, for the last basis integral, it is more convenient to use a different approach based on unitarity cuts. The idea is that, if an integral is UT, this property is preserved by unitarity cuts. In practice, it is convenient to perform maximum cuts that localize the loop momentum (in four dimensions). We have already seen in the above analysis that often "O( )" terms do not influence the UT property, and therefore we will ignore them in the first approximation. (If they do matter, one can always perform a more careful analysis.) In the present situation, we can perform a maximal cut. We proceed in two steps: cutting all visible propagators allows us to localize, say the k 2 integration. The resulting Jacobian produces further propagator factors, which we cut as well. In this way, we obtain a leading singularity of 1/(−q 2 1 + q 2 2 ) 2 . We are thus led to define the candidate integral A comment is in order here. Unlike the analysis of the Feynman parametrizations, the cut analysis did not lead to a proof that the integral in question is UT. However, we will be able to prove this presently from the differential equations. The leading singularity method is invaluable for finding candidate integrals that can then be rigorously proven to be UT.

Differential equations, boundary conditions, and solution
The candidate integrals f = {f 1 , . . . , f 7 } chosen above are all dimensionless functions of the single variable x = q 2 2 /q 2 1 , and . We can derive differential equations in x by differentiating the Feynman integrals w.r.t. q 1 or q 2 , and using integration-by-parts identities to re-express the result of the differentiation in terms of f .
In this way, we find the following system of differential equations (DE) , We see that the DE have three singular points, x = 0, 1, ∞. They correspond to different physical limits, namely q 2 2 = 0, q 3 → 0, q 2 1 = 0. At this point, it is useful to recall the analytic structure of the functions under consideration. From the Feynman parametrization it is obvious that they are real in the unphysical region q 2 1 < 0, q 2 2 < 0, i.e. for any x > 0. As a consequence, they cannot have a singularity or branch cut at x = 1. We can use this information as a boundary condition for the integrals. In fact, thanks to their normalization factors, the integrals f 4 , f 5 , f 6 , f 7 have to vanish at x = 1. Since f 1 , f 2 , f 3 are elementary, this provides the complete boundary information, without any calculation.
We are now in the position to prove that the functions f are indeed UT, as expected. In fact, this is obvious: expanding them in , according to 16) we see that the DE decouples at each order in . It follows that the solution at order k will + 6 2567 9 ζ 2 3 + 59 1512 π 6 + O( 7 ) . (7.25) This is in perfect agreement with results found in the literature. It is straightforward computer algebra to extend this to higher orders in the expansion.

Discussion and conclusion
We have given a review of the differential equations technique for the computation of Feynman integrals. We pointed out general properties of the equations that allow one to make the singularity structure manifest. Moreover, we discussed simplifications in the dependence on the dimensional regularization parameter . In the case where the answer is given by iterated integrals, we discussed a particularly simple canonical form of the equations and various aspects of their solution. We discussed two systematic strategies for obtaining this form, the first one being algebraic in nature, using algorithmic ideas available in the mathematical literature, and the second, more geometric one, using properties of the original loop integrand.
Although some of the algorithmic ideas of the first approach were already used in some form or another in the physics literature, we hope that this presentation clarifies in particular which aspects of the proposal of [9] are always true, and which are conjectural. In particular, we explained how elliptic and more complicated functions appear from this point of view.
The second approach is based on choosing loop integrals with simple generalized cuts (and in particular leading singularities), as proposed in [13], and applied to the DE method in [9]. Although unitarity cuts are frequently used in other contexts, e.g. for finding coefficients of loop integrals in unitarity-based calculations, see e.g. [77][78][79][80] and references therein, to our surprise [13] does not seem to be widely known in the QCD community, and we hope that these lecture notes help to disseminate these useful ideas.
In practice, we have found that a combination of the two methods is most promising: after choosing a basis according to the unitarity cut analysis, the DE are usually either already in the canonical form, or very close to it, so that the remaining transformations are simple to find.
There are several interesting avenues for exploration, and related work we wish to mention.
We have mostly focused around expansion near D = 4−2 dimensions. The same ideas outlined here can also be applied for expansions around other values of the dimension. In particular, it is easy to see from the Feynman representation that Feynman integrals in D and D ± 2 dimensions are related [81]. Therefore, the existence of a uniform weight basis in D = 4 − 2 dimensions implies the existence of a similar basis in dimensions related by multiples of 2. This fact should imply non-trivial constraints on the matrices appearing in the canonical form (4.7) of the differential equations. Similarly, one can apply this method to study Feynman integrals in D = 3 − 2 dimensions, e.g. with applications in ABJM theory.
The dimensional shifts mentioned above can also be helpful in finding a uniform weight basis. For example, the four-dimensional pentagon integral with a numerator discussed in section 6 is equivalent to scalar pentagon integral in six dimensions. Likewise, the fourdimensional bubble integrals with a doubled propagator can be equivalently understood as two-dimensional bubble integrals with standard propagators.
It is important to emphasize that despite the simple dependence on , the differential equations are valid for any dimension, not just for small . Therefore, they can also be used as a starting point for writing down a solution for general dimension. In this case one obtains hypergeometric functions and their generalizations. It would be interesting to connect this to the dimensional recurrence method of refs. [81][82][83].
We wish to mention that in the case of integrals having neither UV nor IR divergences, additional simplifications occur when evaluating the differential equations directly in D = 4 dimensions. In fact, it is possible to set up the IBP technique directly for the case of finite integrals. For an example and more ample discussion, see ref. [25].
We have briefly discussed how elliptic and more complicated functions appear in this setup. The question how to treat such functions systematically is an important conceptual (see e.g. ref. [84]) and practical question, e.g. for scattering amplitudes involving top quarks.
The question whether a given class of Feynman integrals can be evaluated in terms of iterated integrals is crucial to an approach based on direct integration in Feynman parameter space [85,86], where it is related to the question of linear reducibility. The differential equations method also applies to cases that are outside of this class.
Given a Feynman integral/graph, it should be possible in principle to determine what class of functions arises (and in the case of iterated integrals, what the alphabet is) from its propagator structure. See e.g. [87] for work in this direction. The transcendental weight properties of certain Feynman graphs have also been studied in [88].
As we have discussed, it is straightforward to obtain series representations from the differential equations. It would be interesting to explore the relationship to representations in terms of nested sums, see e.g. [89].
The criteria for finding an optimal basis as presented here had the main purpose of transforming the DE into a canonical form, which in turn makes e.g. their singularity structure manifest and makes it obvious which class of functions is needed for their solution. The main tool for this is the analysis of generalized cuts and leading singularities, as proposed in ref. [9]. We already mentioned that leading singularities were also used to construct loop integrands having desirable physical properties, such as good behavior in the infrared [12][13][14]. Such properties are usually lost when employing IBP, since the IBP relations mix for example IR and UV poles. We wish to emphasize that the requirement of a certain IR or UV behavior at the level of the loop integrand, i.e. disregarding IBP relations, is in general a stronger one, and we expect this to be of use e.g. when constructing a natural basis for higher loop integrands. This is closely connected to the generalized unitary approach for computing coefficients of loop integrals in a certain basis, see e.g. [78][79][80] and references therein.