Evaluating single-scale and/or non-planar diagrams by differential equations

We apply a recently suggested new strategy to solve differential equations for Feynman integrals. We develop this method further by analyzing asymptotic expansions of the integrals. We argue that this allows the systematic application of the differential equations to single-scale Feynman integrals. Moreover, the information about singular limits significantly simplifies finding boundary constants for the differential equations. To illustrate these points we consider two families of three-loop integrals. The first are form-factor integrals with two external legs on the light cone. We introduce one more scale by taking one more leg off-shell, $p_2^2\neq 0$. We analytically solve the differential equations for the master integrals in a Laurent expansion in dimensional regularization with $\epsilon=(4-D)/2$. Then we show how to obtain analytic results for the corresponding one-scale integrals in an algebraic way. An essential ingredient of our method is to match solutions of the differential equations in the limit of small $p_2^2$ to our results at $p_2^2\neq 0$ and to identify various terms in these solutions according to expansion by regions. The second family consists of four-point non-planar integrals with all four legs on the light cone. We evaluate, by differential equations, all the master integrals for the so-called $K_4$ graph consisting of four external vertices which are connected with each other by six lines. We show how the boundary constants can be fixed with the help of the knowledge of the singular limits. We present results in terms of harmonic polylogarithms for the corresponding seven master integrals with six propagators in a Laurent expansion in $\epsilon$ up to weight six.


Introduction
A new strategy of solving differential equations (DE) for Feynman integrals [1][2][3][4][5][6] was recently suggested [7]. It is based on choosing a convenient basis of master integrals that are Q-linear combinations of iterated integrals [8][9][10] of uniform weight, i.e. pure functions of uniform transcendental degree. The strategy was then successfully applied to the evaluation of all the three-loop four-point massless planar diagrams with all four legs on the light cone [11] and to two-loop planar diagrams relevant to Bhabha scattering [12]. In the present paper, we develop this strategy further and obtain new results at the three-loop level.
We pointed out in [11] that, as a by-product of the evaluation of four-point massless planar diagrams, we also obtained analytic results for planar single-scale three-point form-factor integrals, although, formally, DE written for one-scale integrals are trivial and express only the homogeneity of the integrals. Nevertheless, the solution of a more complicated problem, with one more scale, provided, in a purely algebraic way, the solution of the one-scale problem, in agreement with the results of [13][14][15][16][17][18].
In this work, the finiteness of planar integrals in the u-channel as u = −s − t → 0 played a decisive role because these boundary conditions turned out to be very restrictive. However, in the non-planar case, there are no such simple boundary conditions. One of the goals of the present paper is to argue that DE, within the strategy of [7], can systematically be applied to single-scale Feynman integrals also in this case.
As we will see, one of the key reasons why this is possible has to do with the fact that the differential equations contain valuable information about singular limits of Feynman integrals. To illustrate this, let us take the case of a set f (x, ǫ) = {f 1 (x, ǫ), . . .} of master integrals that depend on a dimensionless variable x and where D = 4−2ǫ. When applicable, the method of [7] produces a system of differential equations of the Fuchsian type, with a set of constants x i and constant matrices A i . The perturbative solution in ǫ is given by iterated integrals built from the alphabet of differential forms {d log(x − x i )}. Note that at order k in the ǫ expansion one has Q-linear combinations of iterated integrals of uniform weight k.
When one approaches one of the singular points x i , which are often of particular physical interest, f has logarithmic singularities. A typical problem is that the limits x → x i and ǫ → 0 in general do not commute. Here we point out that the knowledge of eq. (1.1) allows to resolve this order of limits ambiguity. It is easy to see from eq. (1.1) that the leading behavior of f as x → x i is where g(ǫ) = {g 1 (ǫ), . . .} represents the boundary constants. Therefore, the singular behavior is governed by the eigenvalues and eigenvectors of the matrix A i . The crucial point is that eq. (1.2) allows us to control the non-commutativity of the limits, since both limits can be generated from the same boundary information g(ǫ). In practice, it is often the case that some limit is particularly simple, or can be related to a previously solved problem.
In that case, one can determine g(ǫ) in that limit, and then use it in the other limit. 1 This simple observation leads to numerous applications. It can be used to determine the asymptotic behavior of Feynman integrals from fixed-order calculations. This applies to physically important singular limits of Feynman integrals and amplitudes such as threshold expansions, soft limits or Regge limits, to name a few examples. 2 Usually such limits are analyzed using the strategy of expansion by regions [20][21][22][23], where the different scalings in eq. (1.2), corresponding to different eigenvalues of A i , are related to various contributions in asymptotic expansions of Feynman integrals. The applications we pursue in this paper use the information about limits that is provided by the DE to compute single-scale and non-planar integrals.
To illustrate our strategy we use the example of a family of three-loop form-factor master integrals with two external legs on the light cone -see Fig. 1.
We introduce one more scale 3 by turning to the corresponding family of integrals with one more leg off-shell, p 2 2 = 0, i.e. depending on two non-zero external momenta squared. After solving DE for the master integrals we obtain analytic results in a Laurent expansion in ǫ = (4 − D)/2. Then we show how to obtain analytic results for the corresponding one-scale integrals in an algebraic way. As mentioned above, this is made possible by eq. (1.2), which allows us to match solutions of differential equations in the limit of small p 2 2 to our results at p 2 2 = 0. Another application is to three-loop four-point massless integrals with all four external legs on the light cone. We previously computed all planar integrals of this type [11], and there are various motivations for extending this to the non-planar case. It would allow to compute complete three-loop scattering amplitudes in supersymmetric Yang-Mills and supergravity theories which currently are only known in un-integrated form [25]. In particular, this would shed light on the infrared properties of gauge and gravity theories. Another motivation is to find out whether one can obtain an equation of the form of eq. (1.1), or whether there is some obstruction due to the non-planar nature of the diagrams.
The non-planar case is, however, much more complicated, for various reasons. Our second goal in the present paper is to evaluate, within the strategy of [7], a particularly interesting subfamily of this class corresponding to the complete four-vertex graph K 4 consisting of four external vertices which are connected with each other by six lines -see Fig. 2(b). It can be considered as a part of the family C in [25]. These integrals have fifteen indices: we associate the first ten of them to the edges of the graph C shown in Fig. 2(a) and the last five to numerators. Explicitly, we have Here s = (p 1 + p 2 ) 2 and t = (p 1 + p 3 ) 2 denote the Mandelstam invariants and the causal prescription −k 2 −→ −k 2 − i0 is implied. For the subfamily associated with the graph K 4 , we have not only a 11 , . . . , a 15 ≤ 0 but also a 1 , a 2 , a 5 , a 6 ≤ 0. When we are interested in K 4 Feynman integrals only without negative indices, we will use the notation K a 1 ,a 2 ,...,a 6 = F C 0,0,a 1 ,a 2 ,0,0,a 3 ,a 4 ,a 5 ,a 6 ,0,...,0 . (1.4) It will be convenient to choose elements of the uniformly transcendental basis which can have the following non-zero indices: a 2 , a 3 , a 4 , a 7 , a 8 , a 9 , a 10 , with a 2 ≤ 0. In this case, we will use the notationK a 1 ,a 2 ,...,a 6 ,a ′ = F C 0,a ′ ,a 1 ,a 2 ,0,0,a 3 ,a 4 ,a 5 ,a 6 ,0,...,0 , where a ′ is always non-positive. The massless graph K 4 was recently discussed [26] in the context of the strategy of evaluating Feynman integrals by iterative integrations over Feynman parameters [27], using multiple polylogarithms. 4 The graph K 4 was discovered not to be linearly reducible [26], i.e. it is impossible to find an order of integration over Feynman parameters such that the dependence on a current integration parameter would be linear so that every iterative integration could be performed in terms of multiple polylogarithms. It was also claimed [26] that the presence of K 4 as a subgraph is crucial for the linear irreducibility at higher loops. The kinematics of the corresponding Feynman integral was considered to be the simplest one at which the linear irreducibility holds, i.e. all the legs were assumed to be on the light cone.
In the present paper, we show that the K 4 Feynman integrals can be evaluated in terms of harmonic polylogarithms, in spite of the fact that the graph is linearly irreducible. To do this, we apply the strategy of [7] to evaluate all the seven master integrals for K 4 with six propagators in a Laurent expansion in ǫ up to weight six, which is the typical order for three-loop calculations.
The paper is organized as follows. In section 2, we derive the system of DE satisfied by the class of non-planar form-factor integrals and discuss the behavior of the integrals in singular limits. We then use this information to analytically determine all integration constants. In section 3 and 4, we discuss the family of non-planar four-point K 4 integrals. We present two ways of computing them. In the first method, presented in section 3, we directly derive a system of differential equations in x = t/s and determine the boundary conditions from symmetry properties and certain asymptotic limits computed via expansion by regions. In the second method (section 4), we introduce one more scale by making one of the external legs off shell. This allows us to fix the integration constants without additional computations, and to solve this three-scale problem analytically. We conclude in section 5.
Together with the paper, we present ancillary files which contain our results (not only presented in the text) with explanations.
2. Evaluating single-scale diagrams by differential equations 2.1 Analyzing asymptotic behavior with differential equations Let us first show how to use differential equations in order to determine the asymptotic behavior of Feynman integrals.
As an example, we consider a family of massless form-factor integrals with two legs off-shell, see Fig. 1. We have where a 1 , . . . a 9 can take any integer values, while a 10 , a 11 , a 12 correspond to potential numerators and therefore can only take non-positive integer values. Moreover, we use the notation k 123 = k 1 + k 2 + k 3 , etc. This is a two-scale problem, with kinematic invariants (p 1 + p 2 ) 2 = s and p 2 2 , while p 2 1 = 0. We denote the dimensionless ratio by x = p 2 2 /s.
Let us write down the α (or Feynman) representation 5 where a = L i=1 a i , h = 3 is the number of loops and U and W are basic polynomials which are given by well-known graph-theoretical formulae, see e.g. [31]. The main point is that where W s and W p 2 2 are positive polynomials in the α i . From this and eq. (2.2) it follows that all integrals are real when s < 0, p 2 2 < 0, i.e. for x > 0. The same fact allows us to absorb the i0 from the Feynman prescription of the propagators in the kinematical variables, −s → −s − i0 and −p 2 2 → −p 2 2 − i0. Setting s = −1 without loss of generality, this means that x acquires a small negative imaginary part, x → x − i0. We will leave this implicit in the formulas presented below.
Using IBP relations with the help of the c++ version of FIRE [32,33], we find that this family can be spanned by a basis of 39 master integrals f = {f 1 , f 2 , . . . , f 39 }. For example, one of the most difficult nine-propagator integrals we take as basis elements is We have normalized all integrals such that they are dimensionless functions, and so that their ǫ expansion starts at ǫ 0 . In particular, they cannot have branch cuts starting at x = 1, and this information will be useful when determining boundary constants for the integrals. We find the following system of differential equations where A and B are constant 39 × 39 matrices. The singular points x = 0, 1, ∞ of eq. (2.5) correspond to the on-shell limit p 2 2 = 0, the two-point function limit p 1 = 0, and the on-shell limit s = 0. It turns out that x → 1 is an excellent limit for determining boundary conditions, because many integrals either vanish or are known simple functions at that point. For a few integrals we have also used the limit x → ∞. These limits, together with simple analytic expressions for propagator-type integrals that are known in terms of gamma functions, completely determine the boundary constants, and therefore allow us to obtain the full solution. 5 Eq. (2.2) is for a10 = a11 = a12 = 0. Let us note that integrals with numerators (negative indices) can be considered by the same formula (2.2) where auxiliary α-parameters are introduced for the negative indices. A differentiation of order −ai for such indices is implied and then they are set to zero, so that a resulting α-parametric integral has the same structure as with nine positive indices but the powers of the two basic functions become shifted and an extra polynomial in the integrand appears. This extra polynomial comes from the differentiation and therefore is only in the numerator and does not alter the analytic structure.
The solution at any order ǫ k is given by a linear combination (with rational coefficients) of harmonic polylogarithms H a 1 ,a 2 ,...,an (x) [34] of weight k. The latter are iterated integrals built from the alphabet of differential forms d log x, d log(1−x), d log(1+x). More precisely, and at least one of the indices a i is non-zero. For all a i = 0, one has From eq. (2.5) we see that in fact only the first two letters of this alphabet are required, or, in other words, only the indices 0 and +1. We explicitly expanded the solution to weight eight. In the remainder of this section, for the sake of readability, we will often truncate formulas at some lower order in ǫ.
To evaluate massless form-factor integrals, let us consider the limit x → 0. This limit does not in general commute with the ǫ expansion, and therefore naively one cannot use the result for fixed order in ǫ to compute the massless form-factor integrals. However, the missing information can be obtained by the differential equation (2.5). It is easy to solve it for small x and for finite ǫ, by neglecting the B/(x − 1) term. The solution in that regime takes the form where g(ǫ) are boundary constants, to be determined. Note that x ǫA is a matrix exponential, which can easily be computed for a given constant matrix A. In a typical situation, where the matrix A is non-diagonalizable, we find x −α j ǫ log k (x), where α j are eigenvalues of A. So, the solution in that regime looks like with the α j are eigenvalues of A, and the c ijk (ǫ) are determined by g(ǫ). (One could make this relationship more concrete by referring to the eigenvectors and in general power vectors of the matrix A.) Expanding this formula for small ǫ, we can determine the matching coefficients c ijk by comparing to our results in a Laurent expansion in ǫ. Then, we can return to formula (2.11), keep only the terms with α j = 0 and, therefore arrive at the form-factor integrals with the external momentum p 2 on-shell, i.e. x = 0. Indeed, the integrals considered on-shell, or at a threshold are, by definition, obtained from integrals at general values of a given external momentum by setting it on-shell or a threshold under the integral sign either in integrals over loop momenta or in parametric integrals. On the other hand when we consider the limit x → 0 we can apply the strategy of expansion by regions [20][21][22][23] (see also Chapter 9 of [31]). According to this strategy, the expansion in a given limit is given by a finite number of series corresponding to so-called regions which are scalings of certain components of the loop momenta in terms of a given parameter of expansion, e.g. x. One of the regions corresponds to all the components of the loop momenta (or all the parameters in alpha representation) to be of order x 0 . It is usually called hard. Its contribution is given by Taylor expanding the integrands in x. The leading term in this contribution is obtained just by setting x = 0 under integral sign. Obviously, this contribution corresponds to α j = 0 in eq. (2.11). Other regions typical for Sudakov and Regge limits are called collinear (with α j = 1 per loop) and ultrasoft (with α j = 2 per loop).

Example
As an example of this, let us discuss the solution for one of the most non-trivial integrals, f 38 . Solving the system of differential equations with the appropriate boundary conditions as discussed above, we find is the result for integral f 38 in the small ǫ limit. When taking in addition x → 0, divergences appear from the logarithms in that formula. As explained above, we can understand these divergences from the general solution of eq. (2.5), namely eq. (2.11). In the present case, we find 6 (2.13) For ǫ → 0, the x −αǫ terms lead to the logarithms we observed above. We can also see that there is an order of limits issue when considering ǫ → 0 and x → 0. On the other hand, we arrive at this general structure from the analysis of the system of differential equations. Moreover, we determine the expansion coefficients a i (ǫ) by comparing against the solution for small ǫ. Indeed, let us compare the small x limit of eq. (2.12) to the small ǫ limit of eq. (2.13). This leads to a constraint on the c i (ǫ). Of course, we have an analog of eq. (2.13) for all integrals, which implies that we can completely determine the coefficients a i (ǫ) by this procedure. For notational brevity, we present the results only up to order ǫ 4 , 14) Having obtained these 'matching coefficients', we can return to eq. (2.13) and consider the opposite order of limits, i.e. x → 0 for finite (negative) ǫ. In this case, we ignore the terms x −αǫ with positive α and we are left with This is nothing but the value for the massless form-factor integral. Comparing to eq. (22) of ref. [18], we find perfect agreement to order ǫ 8 . We have calculated, in a similar way, all the master integrals for the family of Feynman integrals (2.1) at p 2 1 = p 2 2 = 0 and found agreement with the results [18] up to weight eight and earlier results up to weight six [15,16]. We stress that the calculation performed here was done entirely from first principles, using only algebraic steps.
Let us comment on the general structure of the asymptotic expansion. The x −αǫ terms correspond to contributions of certain regions [20]. Although we did not carry out a detailed analysis, one would expect that the term x −ǫ corresponds to regions where one of the loop momenta is collinear and two loop momenta are hard, the term x −2ǫ corresponds to regions where two loop momenta are collinear and one loop momentum is is hard, and the term x −3ǫ corresponds to regions where all the loop momenta are collinear. As discussed above, the term x 0·ǫ corresponds to the region of the three hard momenta. Although intuitive, in general it is hard to find these regions in momentum space systematically. In contrast, revealing regions in the space of Feynman/alpha parameters [22,31] can be made automatic.
To do this, one can apply an open code asy.m [35,36]. In this particular example, this code reports about seven contributions corresponding to certain regions. The power dependence on x in these contributions exactly corresponds to the exponents present in (2.13). More specifically, there is one region with x 0·ǫ (this always takes place), one region with x −ǫ , two regions with x −2ǫ , and three regions with x −3ǫ .

Evaluating K 4 integrals
Here we evaluate the K 4 integrals defined in the introduction, see eq. (1.3) and Fig. 2(b). As mentioned in the introduction, these integrals are a non-planar version of the three-loop integrals solved for in [11]. In that reference, we used a simple boundary condition of the absence of singularities in the u-channel in order to determine the boundary constants. For non-planar integrals, we do not have a similar condition. However, the setup discussed in section 2 gives us control even in the case where the limits are singular, and this will help us in fixing the integration constants.
There is another complication related to the non-planar nature of the integrals that can be seen by looking at the α representation. The polynomials U and W in eq. (2.2) are given by We see from eq. (3.2) that W does not have a definite sign (for some region of s, t), and as a consequence we cannot treat the i0 prescription in eq. (2.2) simply as a complex deformation of s and t. This can also be seen simply by looking at bubble integrals in the s, t and u-channel, which give rise to logarithms where we used u = −s − t. We therefore need to be careful about those different i0 prescriptions. As in section 2, we define dimensionless functions of one variable, x = t/s. In the calculation below we will assume x > 0, unless otherwise stated. Solving IBP relations [37] with the help of the c++ version of FIRE [32,33], for the family of integrals (1.3), we find that, at a 1 , a 2 , a 5 , a 6 , a 11 , . . . , a 15 ≤ 0, there are three trivial master integrals with four propagators and seven master integrals 7 with six propagators. In different situations, it is reasonable to choose different bases of the master integrals.
An eight-fold MB representation for K 4 -integrals without numerator is presented in Appendix A. However, the straightforward procedure of evaluating these integrals by the MB representation works well only up to weight three in the epsilon-expansion, while we have an implicit obligation to obtain results up to weight six. Therefore, we will now use differential equations, and the results obtained with the MB representation will be used for checks.

Asymptotic behavior and boundary conditions
The first three functions g 1 , g 2 , g 3 , are simple functions that can be given analytically in terms of Γ functions. Therefore we only need to specify boundary values for the remaining seven functions.
As discussed above, the eigenvalues of ǫA, ǫB, and ǫ(A + B) characterize the three singular limits of f . We have the eigenvalues {−3ǫ, 0, ǫ} for t → 0 and {−4ǫ, −3ǫ, 0} for s → 0. (Notice the rescaling by (−s) −3ǫ in eq. (3.4).) Finally, for u → 0, we have {−3ǫ, 0, ǫ}. We see that some of the eigenvalues are positive, e.g. ǫA has eigenvalues +ǫ. This is slightly surprising, for the following reason. All integrals in our basis f are UV finite, and IR divergent. Therefore, they can be defined for ǫ < 0, and in particular they should stay finite if we take a limit such as x → 0 (with ǫ finite). However, a term like x ǫ corresponding to the eigenvalue +ǫ would diverge when x → 0. Therefore we expect that the coefficients of such terms must vanish. This requirement fixes some of the integration constants.
In order to have further analytic boundary conditions, we computed asymptotic expansions using the computer code asy.m [35,36] which is now included in FIESTA3 [42]. This code uses the information about the propagators of a given Feynman integral as an input and produces the corresponding set of regions relevant to a given limit, in the language of Feynman parameters. The search of regions reduces to finding faces of maximal dimension of the Newton polytope associated with the two basic polynomials in the alpha representation. This code provides various contributions to a given asymptotic expansion as parametrical integrals and performs as many explicit integrations in these integrals as possible. In the case of the limit t → 0, these are the contributions (called hard and collinear) characterized by exponents x 0 and x −3ǫ of the expansion parameter, in agreement with the discussion above. Starting from these parametric integrals we derived a one-fold MB representation for the collinear-type contributions (with the exponent x −3ǫ ). The evaluation of the corresponding MB integrals in the ǫ-expansion is straightforward. It reduces to summing up one-fold series, which can be done with the help of public computer codes [43,44]. The results obtained can be expressed in terms of multiple zeta values.
Let us illustrate this procedure using the master integral K 2,2,1,1,1,1 . The code asy.m reveals one hard contribution and two collinear contributions. To deal with the two collinear contributions individually, one has to introduce an auxiliary analytic regularization. This can be done by shifting the first index 2 by an analytic parameter λ, i.e. we will consider a 1 = 2 + λ. Then the code asy.m produces the following expression for the sum of the two collinear contributions: As is well known, one can choose any subset of the parameters in the argument of the delta functions involved.
The second of the two integrals can be evaluated as follows. We choose the argument of the delta function as α 4 −1 so that we set α 4 = 1 and obtain an integral from 0 to ∞ over the remaining three parameters. We turn to the new variables by α 1 = ηξ, α 3 = η(1 − ξ) and integrate explicitly over η. Then we separate the two terms in (−(1 − ξ) + ξα 2 − i0) by introducing a MB integration and take explicitly integrations over α 2 and ξ in terms of gamma functions. One can proceed similarly with the first integral. Then we can take the limit λ → 0 in the sum of the two integrals to obtain the following expression for the sum of the two collinear contributions in terms of a one-fold MB integral: One can evaluate this integral in a Laurent expansion in ǫ by the standard procedures [45][46][47][48], first, resolving the singularities of the integral in ǫ and then converting the integrals obtained into a series and summing it up [43,44]. We obtain the following result for the sum of the leading order collinear contributions to K 2,2,1,1,1,1 in the limit x → 0: For the hard contribution (i.e. terms with the exponent x 0 ), it was possible to derive a four-fold MB representation. Moreover, we could simplify the evaluation taking into account the fact that the corresponding contributions are given by Feynman integrals depending on two, rather than three external momenta because the kinematics t = 0 implies p 3 = −p 1 . Therefore, we could apply an IBP reduction to such integrals. The number of the corresponding master integrals drastically reduces: in the sector with six positive indices, it is equal to two instead of seven. We used the four-fold MB representation to evaluate these master integrals in lower orders of the ǫ expansion. However, as we will discuss in the next subsection, we obtained sufficient boundary information in other ways, so that we used such results only for checks.

Crossing symmetry
The K 4 integrals have a huge amount of symmetry under exchange of external momenta. For example, studying the exchanges p 1 ↔ p 2 and p 1 ↔ p 3 , we find the following relations, and 17) and similarly for f 1 , f 2 , f 3 . This can be seen by inspecting the α representation, and in particular the polynomial W of eq. (3.2) for those integrals. In this way one sees that the relations that interchange x and −1 − x are valid for arbitrary real x, while the relations that interchange x ↔ 1/x above are only valid for x > 0. For x < 0, the i0 prescription leads to the following modification, where * stands for complex conjugation and ǫ is supposed to be real. We used these symmetries in order to fix the boundary constants remaining from the discussion in subsection 3.1. The remaining relations served as a check of our calculation. K 1,1,1,1,1,1 Using the results for the boundary conditions of subsection (3.1) we solved the differential equations (3.10) for f to order ǫ 6 , i.e. weight six, which is the typical weight required for three-loop computations. The explicit results are given in attached text files.

Result for
Here, as an example and an application of these results, we will define the integral which is of special interest, as discussed in the introduction. It is related to the basis above via As a consequence, it has the same uniform weight properties as the f i . The first terms of its expansion in ǫ are given by The terms of order ǫ 5 and ǫ 6 are presented in the appendix, for completeness. Note that in eq. (3.21), we have chosen to represent the answer in such a way that the branch cuts of all functions involved lie on the negative real axis. 9 Finally, let us mention that integral K 0 of eq. (3.19) is completely crossing symmetric (see subsection 3.2), i.e.
As was discussed in subsection 3.2, for x < 0 eq. (3.23) is to be replaced by One may explicitly verify eq. (3.22) for −1 < x < 1 (cf. footnote 9) and (3.23) for x > 0 and this is a non-trivial test of our result (3.21). We have done so using the convenient Mathematica implementation of harmonic polylogarithms [49]. We stress that, as a consequence of the form of the differential eqs. (3.10), all integrals are pure functions of uniform weight 10 .

Further analytic and numerical checks
The terms up to weight three are in agreement with analytical results which we obtained with the MB representation presented in Appendix B.
We have checked our results for the master integrals numerically by FIESTA [42,50,51]. To do this we evaluated with FIESTA the canonical master integrals because evaluating integrals with an index equal to 2 is preferable to evaluating integrals with an index equal to −1. Then the numerical results for the elements of our uniformly transcendental basis could be obtained because we have an IBP reduction at hand. The agreement between our analytical and numerical results was achieved at least at the level of three digits in the ǫ expansion up to weight six.
The alphabet of differential forms we obtain is the same as the one occurring at the previous loop order. (We have verified that the same form of the equations (4.14) and (4.15) holds at two loops, by choosing a slightly different basis compared to [5,6].) Leaving the issue of the boundary conditions aside for the moment (this will be dealt with in the next subsection), eq. (4.14) allows us to solve for f to any desired order in the ǫ expansion. It is clear that each term in the answer will be given by iterated integrals over the differential one-forms shown above. This class of functions forms a subset of multiple polylogarithms and was studies in [5,6].
More generally, we can write the solution to eq. (4.14) in the beautiful language of Chen iterated integrals. Let M be a (in general complex) manifold describing the kinematical data, in this case s, t, u. Each element of the matrix dÃ is a one-form on this manifold. The integration contour is then a path on this manifold. We can parametrize it by defining a map γ :  Using these iterated integrals, we can write down the general solution to eq. (4.14). It is given by where h(ǫ) represents the boundary condition. Expanding the exponential in eq. (4.18) perturbatively in ǫ, one obtains at order ǫ k (linear combinations of) k-fold iterated intervals. The latter are homotopy invariant line integrals, with γ connecting the base point (s 0 , t 0 , u 0 ) to the argument of the function, (s, t, u). Upon choosing a specific contour of integration, one can recover expressions in terms of multiple polylogarithms. In the next section, we will provide the information for determining the boundary constants.

Determining the boundary conditions
There are several boundary conditions that we can use, as we discuss presently: • Elementary integrals: integrals f 1 , f 3 , f 5 are trivial bubble-type integrals that can be expressed in terms of Γ functions.
• Branch cut structure: for massless integrals, we expect branch cuts to start only at positions p 2 i = 0, (p i + p j ) 2 = 0, etc. Inspecting the terms on the r.h.s. of eq. (4.14), we see that only the logarithms on the first line have arguments of this form. Imposing the absence of branch cuts coming from functions like log(−s − t) then imposes constraints on the answer. We have verified this expectation using the computer code asy.m mentioned earlier.
• Asymptotic limit and UV behavior: just as in the previous section, we can determine some of the coefficients in the asymptotic expansion by requiring the absence of UV divergences in the basis f . For example, in this way it can be seen that this implies that f 8 ∝ K 2,1,1,1,1,1 → 0 as t → 0. • Simple limits: In general, the limits at the singular points of the DE do not commute with the ǫ expansion. However, for some integrals the situation is simpler. For example, for integral K (0) one expects soft limits such as p 1 → 0 to commute with the ǫ expansion, since they do not change the divergence structure of the integral.
(The integral is IR finite, and stays IR finite in the limit. The UV divergences are unchanged by the limit.) At p 1 = 0, however, K 1,1,1,1,1,1 becomes a known planar form-factor integral, and we can use its value as boundary condition. 11 We have found the above requirements to be sufficient to determine all boundary constants for all 16 integrals, order by order in ǫ. In fact, one can see that the first three elements on the list above are sufficient to fix all integration constants. We stress that these conditions do not require any integrations and can be implemented in an algebraic way.

Analytic solution and on-shell limit
Here we present the analytic solution for the first orders in the ǫ expansion of the integrals. As discussed above, it can be written as (4.18), with the boundary conditions following from the considerations in the previous paragraph. In (4.18), one has the freedom of choosing a base point for the iterated integral. One reasonable choice would be s = t = u = −1, since this stays away from all potentially singular surfaces. However, this leads to rather awkward integrations, such as log 3 at weight one. In the literature, results for the alphabet (4.15) are usually represented in terms of so-called two-dimensional harmonic polylogarithms, a subset of Goncharov polylogarithms (GPL). GPL are defined as follows.
G(a 1 , . . . , a n ; z) = z 0 dt t − a 1 G(a 2 , . . . , a n ; t) , (4.19) 11 Note that the form-factor integral can itself be determined by bootstrap arguments [11]. For a i = 0, we have G( 0 n ; x) = 1/n! log n (x). The total differential of a general Goncharov polylogarithm is dG(a 1 , . . . , a n ; z) = G(â 1 , a 2 , . . . a n ; z) d log z − a 1 a 1 − a 2 + G(a 1 ,â 2 , a 3 , . . . , a n ; z) d log a 1 − a 2 a 2 − a 3 + . . . + + G(a 1 , . . . , a n−1 ,â n ; z) d log a n−1 − a n a n , (4.21) whereâ means that this element is omitted. The subset of two-dimensional harmonic polylogarithms is obtained by specifying labels to be from the set {0, −1, −1 − y, −y} and argument z = x. One can easily convince oneself that this set of functions, together with HPL of argument y, is sufficient in order to represent the solution (4.18). It essentially corresponds to choosing s = −1, t = 0, u = 1 as base point (after separating the logarithmic divergences as t → 0.) Here we followed a slightly different approach, by first integrating the differential equation in x, and then in y. The procedure is almost the same as in [12].
The boundary constants are determined from the conditions discussed in section 4.2. When approaching the various limits discussed there, one sometimes encounters spurious divergences of the Goncharov polylogarithms. This is a slight disadvantage of having gone from the language of Chen iterated integrals to the latter. Such divergences can be extracted before taking limits by using the shuffle product formula of Goncharov polylogarithms, G(a 1 , . . . , a n ; z)G(b 1 , . . . , b m ; z) =G( a; z)G( b; z) = c∈ a ⊎ b G( c; z) , (4.22) where a ⊎ b is the shuffle product of two ordered sets, i.e. all combined sets where the relative order of the elements of a and b is preserved. In this way, the complete solution can be obtained algorithmically. We give an example for illustration. For integral f 8 up to order ǫ 2 , we have and similarly for the other integrals. Of course, the results up to weight two can easily be rewritten in terms of logarithms and dilogarithms, and probably in a more compact way. We prefer to use the language of GPL and HPL because it is valid at any order in ǫ. Finally, we wish to outline how to recover the on-shell case discussed in section 3. In this way, one sees that the boundary constants for those integrals also follow from the general considerations made here. The limit p 2 4 → 0 is governed by the term ǫA 4 log(−s − t − u) in the differential equation. Analyzing ǫA 4 , we find that it has two possible eigenvalues −2ǫ and 0. One can then proceed as explained in section 2 in order to resolve the order of limits ambiguities, and obtain results for K 4 on-shell.
This completes our discussion of the K 4 integrals. In summary, we have seen that the latter are completely determined by the differential equations discussed here, and that the boundary constants follow from simple physical considerations. In our setup, it is clear that the results are pure functions of uniform weight in the ǫ expansion. We have also outlined how to recover results in the on-shell case from this setup. Note that in this approach, in contrast to section 3.1, no integrals have to be calculated in order to determine the boundary conditions.

Discussion and outlook
In this paper we observed that the differential equations for master integrals can be used to infer the structure of asymptotic expansions of the master integrals. We applied this information to the computation of single-scale and non-planar integrals. Although we mainly had in mind to give examples showing the scope of this method, many of the results derived in this paper are new, and the integrals computed are highly non-trivial. In this paper, we mainly focused on using the information on singular limits in order to provide simple boundary conditions for the DE. Of course, one can also go the other way. We expect that this will have many applications for the analysis of physically interesting limits.
In order to use the method of DE for single-scale integrals, we first generalized the problem by introducing an extra scale. On the one hand, this leads to an increase of the number of master integrals needed, from 14 to 39. On the other hand, it allows to use the powerful DE technique. In fact, once the basis of master integrals is chosen appropriately [7], the number of integrals does not play an important role in the structure of the equations. Moreover, the additional scale gives access to new limits where the boundary constants can be determined easily. In this way, we solved the more general twoscale problem, using only algebraic means. Finally, the knowledge of the precise scaling behavior of the Feynman integrals, also inferred from the differential equations, allowed us to relate the two-scale problem to the one-scale problem we started with.
The procedure leading from the Feynman integrals to the set of differential equations for master integrals is entirely algebraic. The differential equations, especially when written in the simple form of [7], make it clear which class of special functions is needed to describe the Feynman integrals, to all orders in ǫ. In particular, this also determines what types of transcendental constants can appear at special values of these functions. These periods of Feynman integrals are heavily studied in the mathematical literature. The approach proposed here, which consists of solving a more general problem via differential equations and then to obtain the periods as a corollary, has also been used in the mathematical literature [52]. Feynman integrals depending on a parameter, called graphical functions there, were computed there with the help of the second-order differential equations of [53,54]. The desired periods where then obtained at special values of those functions, and this was used to prove a conjecture made in ref. [55].
The second class of Feynman integrals computed in this paper is a family of non-planar on-shell four-point functions. We consider integrals corresponding to the graph K 4 . They were found not to be linearly reducible in the framework of ref. [26]. Here we studied them in the context of DE and found that the DE have the same form as in the previously studied planar case [11], and in particular, lead to the same class of multiple polylogarithms and transcendental constants.
We did two calculations for these integrals, the first one, in section 3 being a direct one. In order to determine the integration constants, we used information from the asymptotic expansion of the integrals. This was possible thanks to the control the DE give over such expansions. In section 4, we performed the calculation for the same integrals with one external leg off-shell. As in the case of the form-factor integrals, the number of master integrals increased, here from 10 to 16. On the other hand, this allowed us to fix the integration constants in a clear and simple way. We outlined how the results for the onshell integrals can be recovered. This was done mainly as a proof of principle that no integrals have to be performed in order to find the integration constants.
The method and results presented here for the K 4 integrals strongly suggest to us that the planar results of [11] can be carried over to the non-planar case. The results presented here and in [11] already allow the computation of non-planar scattering amplitudes in φ 4 models. Completing the calculation for all non-planar master integrals will allow the computation of non-planar scattering amplitudes in super Yang-Mills and supergravity theories that are currently only known at the integrand level [25]. This will give valuable insights into the generic structure of infrared divergences in gauge and gravity theories. The methods developed in the present paper should be extremely helpful in determining the boundary constants for the required non-planar integrals.
However, at concrete integer indices, there is usually the possibility to take two integrations by means of the first Barnes lemma. In particular, we obtain To evaluate this and other above mentioned master integrals one can apply public computer codes [45][46][47][48]. This straightforward procedure works well only up to weight three in the ǫ-expansion.

B. Explicit results for K (0) up to weight six
Here we present analytical results for the integral (3.19) which is expresed in terms of the master integrals discussed in the main text via eq. (3.20) and therefore has the same uniform weight properties as the f i . We have K (0) (x, ǫ) = 2ζ 3 ǫ 3 + ǫ 4 3iπζ 3 + 3π 4 20 + 2iπH −3 + 1 2