Iterative structure of finite loop integrals

In this paper we develop further and refine the method of differential equations for computing Feynman integrals. In particular, we show that an additional iterative structure emerges for finite loop integrals. As a concrete non-trivial example we study planar master integrals of light-by-light scattering to three loops, and derive analytic results for all values of the Mandelstam variables $s$ and $t$ and the mass $m$. We start with a recent proposal for defining a basis of loop integrals having uniform transcendental weight properties and use this approach to compute all planar two-loop master integrals in dimensional regularization. We then show how this approach can be further simplified when computing finite loop integrals. This allows us to discuss precisely the subset of integrals that are relevant to the problem. We find that this leads to a block triangular structure of the differential equations, where the blocks correspond to integrals of different weight. We explain how this block triangular form is found in an algorithmic way. Another advantage of working in four dimensions is that integrals of different loop orders are interconnected and can be seamlessly discussed within the same formalism. We use this method to compute all finite master integrals needed up to three loops. Finally, we remark that all integrals have simple Mandelstam representations.


Introduction
In perturbative quantum field theory, beyond the leading order interesting quantities such as scattering amplitudes or correlation functions involve Feynman integrals. They are multivalued functions of the kinematical invariants and masses. One finds that one-loop integrals in four dimensions can be expressed in terms of logarithms and dilogarithms. At higher orders, generalizations to other special functions are typically required. Understanding the structure of loop integrals, in particular precisely which special functions are needed, and their analytic computation is an important problem.
In the literature on Feynman integrals one finds various generalizations of the logarithms and dilogarithms encountered at one loop order. Most of them fall into the rather general class of multiple polylogarithms [1,2], also referred to as Goncharov polylogarithms or hyperlogarithms. 1 As we will see, our method will naturally suggest to work with a slightly more general class of functions, Chen iterated integrals. They enjoy beautiful mathematical properties that we are going to benefit from.
A variety of different approaches for computing Feynman integrals can be found in the literature, and often there are interrelations between them. Here we wish to focus on the method of differential equations [3][4][5][6], which has provided many results in the past, in particular for virtual corrections to scattering processes. Recently, a proposal has appeared to render this technique more systematic. The suggestion is that a standard form for the differential equations can be obtained, by choosing a particularly convenient basis [7]. The choice of the latter is guided by the properties of the functions expected to appear in the answer, and in particular their transcendental weight properties.
Following the approach of [7] we write a system of first-order differential equations for the required integrals, and bring the latter into a simple form by choosing a convenient basis. Denoting the set of basis integrals by f , we seek to express their differential as where the differential acts on the kinematical variables s, t, m 2 . The special feature of eq. (1.1) is that its dependence on is simple, namely the matrixÃ is independent of . Moreover, it contains only logarithms of combinations of the kinematical variables, in agreement with the expected singularities of Feynman integrals. More precisely, where a k are certain constant matrices, and the α k are algebraic functions of s, t, m 2 . They are called letters, and the set set of letters appearing inÃ specifies the alphabet required for the integrals f . Once the form (1.1) is reached, one can immediately write down a solution to eq. (1.1) in terms of Chen iterated integrals [8], to all orders in the expansion. Depending on the properties of the alphabet {α k }, one may choose a representation in terms of more familiar functions. The method was successfully applied to many integrals, including [9][10][11][12]. Although the conjecture in [7] is still unproven, that the special form (1.1) can always be achieved, this was verified in all cases.
A convenient feature of (1.1) is that the equations give all orders in the expansion. This feature will not, however, be so desirable for us here, for several reasons. First of all, since all the integrals we will consider are already convergent in D = 4, the expansion represents unneeded information. Possibly the differential equation would become even simpler if it did not need to carry this information. Second, for particular applications, one may have a symmetry enhancement at D = 4. In particular, the integrals we will be interested in enjoy (dual) conformal symmetry, in the sense of ref. [13], which is broken away from D = 4. A similar situation occurs for example for correlation functions in conformal field theories. Clearly, it is desirable to set up the computation in a way which preserves the symmetries. Third, in integer dimensions integrals with different loop orders can be related to each other [14][15][16], and this allows us to discuss them in a seamless way.
We will thus develop a simplified version of the approach of [7], applicable directly in D = 4 dimensions. In a first step, we will assume that, by using suitable integral identities that are valid among (convergent) integrals defined in strictly D = 4, a basis of "master integrals" has been identified such that the derivatives of all integrals with respect to external parameters can be expressed within the basis. We will then seek a change a basis where the differential takes the special form d g(s, t, m 2 ) = (dA) g(s, t, m 2 ) , (1.3) where the matrix A will be required to be block triangular. By this we mean that there should exist a grading ("weight") such that the derivative of a weight n integral is expressed in terms of only weight (n − 1) integrals. In particular, all diagonal elements of A vanish. The weight will turn out to correspond the so-called transcendental weight, although at this stage we are merely discussing the matrix structure of A.
Just as in the D-dimensional case, one can immediately read off from eq. (1.3), and more precisely from the entries of the matrix A, what class of functions the answer will be expressed in. It defines the integration kernels that can appear in the iterated integrals that the solution is expressed in. Moreover, the block triangular structure of A implies that the integrals can be organized according to their (transcendental) weight, and can be computed in an iterative way, as will be explained in detail in section 3.
In this paper we will compute families of one, two, and three-loop three-scale integrals defined in D = 4 and giving the complete light-by-light scattering in planar N = 4 super Yang-Mills [17]. Anticipating other applications, we will also compute the one and two loop families at general value of . This set of integrals contains a number of two-scale integrals that appear in massive form factor calculations. They were computed previously (in a different basis) to some order in in ref. [18]. We confirm these results. Moreover, the formulation we give can be trivially expanded to any order in , where the result is given by a homogeneous expression in terms of harmonic polylogarithms. This paper is organized as follows. In section 2, we give definitions of the loop integrals that appear in the paper. Then, in section 3, we briefly review the differential equations method for loop integrals, and discuss simplifications in the four-dimensional limit. We use the one-loop integrals as a pedagogical example. In section 4, we explain how to systematically set up the differential equations directly in four dimensions, and present an algorithm for putting the latter into a canonical block-triangular form. We give the differential equations at two and three loops and discuss the iterative structure of the analytic solution. In section 5, we discuss the analytic properties of the functions to three loops and show that they satisfy a Mandelstam representation. We discuss checks of our results in section 6. We conclude in section 7. There are three appendices. In appendix A we apply the differential equation method to the set of two-loop master integrals in D = 4 − 2 dimensions and compute them using the method of ref. [7]. Additional material on our method for writing down identities and differential equations for four-dimensional loop integrals can be found in appendix B. Appendix C contains expressions for the oneand two-loop box integrals in terms of multiple polylogarithms. Appendix D contains the differential equations up to three loops. We supply several ancillary electronic files together with the arXiv submission of this paper.

Definitions and kinematical preliminaries
Let us introduce the integral families that we are going to discuss in this paper. As was already mentioned, the motivation for these particular integrals is that they can be used to study light-by-light scattering in planar N = 4 super Yang-Mills [17].
We define the inverse propagators Note that we changed the convention for the metric w.r.t. [13] to the mostly minus metric (+ − . . . −). The integrals we discuss include the one-loop family where the a i can take any integer values. At two loops we discuss the family where a 4 ≤ 0, b 2 ≤ 0 represent possible numerator factors. See Fig. 1(a) and Fig. 1(b) for the one-and two-loop cases, respectively. Finally, defining the two three-loop families we discuss are the following. For the first family, a 4 ≤ 0, b 2 ≤ 0, b 4 ≤ 0, c 2 ≤ 0, f ≤ 0 represent possible numerator factors, and the second family is defined in the same way, but instead with a 1 ≤ 0, a 4 ≤ 0, b 3 ≤ 0, c 1 ≤ 0, c 2 ≤ 0. They are represented in Fig. 2(a) and Fig. 2(b), respectively. Each G is a function of the dimension D, the Mandelstam variables s = (p 1 + p 2 ) 2 and t = (p 2 + p 3 ) 2 , and the internal mass m. The on-shell conditions for the external legs read A one-loop integral in the above family will be said to be dual conformal invariant, or DCI for short, if The name is motivated by a natural generalization to a multiple mass case, where one can define dual conformal transformations under which the integral is covariant, see [13] and the review [19]. We use a similar definition at higher loops, equating to 4 the sum of indices associated with each loop variable. For example, a two-loop integral is DCI if (a 1 + a 2 + a 3 + a 4 + c) = 4 and (b 1 In the following, we will define all integrals in the Euclidean region s/m 2 < 0, t/m 2 < 0, where all functions are real-valued. One then defines the functions elsewhere by analytic continuation, using the Feynman prescription. This implies giving the kinematical variables a small imaginary part, according to Let us have a first look at these integrals, and take I 1 as an example. It is given analytically by (the form below is due to [20]), Here we introduced dimensionless variables 2 and the following abbreviations, The functions appearing in eq. (2.6) are examples of polylogarithms. For these and more general classes of integral functions that we will discuss one can define a "symbol" [1,2,21]. Roughly speaking, the symbol contains information about the integration kernels leading to those functions, while forgetting about boundary constants at each integration step. We note that the symbol of the above formula is very simple, and visibly more compact compared to eq. (2.6), This foreshadows a simple structure of the integral under the action of differential operators.
In the next two sections, we will see how this structure arises in a systematic way. It will however not be necessary to restrict the analysis to the level of the symbol, rather the observations will apply to the functions directly. Indeed, we will be able to write compact formulas in terms of iterated integrals that make the simplicity manifest, while at the same time keeping track of the integration constants.
In particular, the complete information specifying the multi-loop integrals we will discuss will be contained in simple formulas similar to eq. (2.9). In the next two sections, we will first see how to reproduce this formula from differential equations, and then proceed to compute the required integrals at two and three loops.

Differential equation at one loop: 4 versus D dimensions
Computing loop integrals via differential equations is by now a fairly standard procedure [3][4][5][6], so we will only briefly outline the main steps. For a given class of integrals under consideration, defined by the propagator structure, one defines integral families for any possible powers of the propagators. E.g. at one loop, the family we consider is defined by eq. (2.2), see also Fig. 1(a), for any integer values of the parameters a i . The integrals in a given family are in general not independent. They satisfy integration-by-parts (IBP) identities [22]. The latter are linear in the integrals. In practice, this means that one can relate the calculation of any given integral to a finite number of master integrals. The algebraic manipulations needed are straightforward but cumbersome, so that one usually performs them using computer algebra codes. We have used FIRE [23,24], as well as an own implementation, to be described below. Other public computer codes include [25,26]. The independent integrals obtained in this way are usually referred to as master integrals. The basis choice is not unique. We have found [27] helpful to verify the number of master integrals needed in a given sector.
For a given set of master integrals, one sets up differential equations in the kinematical invariants. The derivative w.r.t. m 2 is straightforward to implement. The derivatives w.r.t. s and t can be written down via the chain rule, in terms of derivatives w.r.t. the external momenta, e.g. 3 When acting with the operator (3.1) on integrals such as (2.2), it is clear that the resulting terms will still be within the same family of integrals. Therefore, any derivative of a master integral can be expressed in terms of a linear combination of master integrals. This implies that the master integrals satisfy a linear system of first-order equations. Let us discuss this system in the one-loop case. It is well known that the set of master integrals at one-loop, in the presence of internal masses, consists of the box, triangles, bubbles and tadpole topologies. Thanks to symmetries, for the present configuration this amounts to six master integrals. As we will see, the corresponding 6 × 6 system of differential equations can be put in the form of eq. (1.1) by choosing the following basis of master integrals, where c = m 2 /Γ(1 + ) is a convenient overall factor. We start by setting up the system of first-order differential equations for the basis functions f i . which takes the form, in our basis choice, 4 The crucial simplification of our basis choice lies in the fact that the partial derivative matricesÃ s ,Ã t andÃ m 2 are independent of . Moreover, as we will see presently, they can be expressed in terms of derivatives of logarithms only. The f i are normalized to be dimensionless functions and therefore satisfy which is a useful consistency check. Other consistency checks are obtained from requiring that partial derivatives commute, which implies where i, j can be s, t, or m 2 . The structure of the differential equations, and hence that of the f i , is greatly clarified by combining equations (3.3) by writing them in differential form This is the form quoted in the introduction. HereÃ was obtained by integrating the partial differential equations This integration is very similar to the one that is needed to integrate f (at the first order in the expansion). Performing it at this stage allows us to completely elucidate the properties of the functions f . In fact, we find thatÃ can be written in terms of logarithms only, We have 'boxed' the elements which survive in the → 0 limit, to be discussed presently.
To fully define the answer we need to specify a boundary condition. For the integrals under consideration a particularly natural boundary point is m = ∞, where Before commenting on the solution to the differential equations, let us first discuss a simplification for finite integrals. We will be particularly interested in the → 0 limit. Since all the above integrals are both ultraviolet and infrared convergent (thanks to the internal masses, and thanks to the choice of representatives we made), it is natural to remove the powers of in eq. (3.2) in the limit. We thus let (g 1 , g 2 , g 3 ,g 4 ,g 5 , g 6 ) = lim The powers of reflect the transcendental weights of the functions g i . We have placed 'tilde' ong 4 andg 5 to distinguish them from two-loop integrals g 4 and g 5 which will be introduced in the next section. The integrals g 1 , g 2 , g 3 and g 6 will carry the same meaning throughout this paper. The differential equation for the g i takes the similar canonical form (1.3), with the A-matrix obtained fromÃ in eq. (3.8) by retaining only the boxed elements.
Two observations will be important for us.
• The integrals g i have uniform transcendental weights (0, 1, 1, 2, 2, 2), respectively. This follows immediately from the block-triangular structure of the A-matrix: the derivatives of g 6 are expressed in terms of g 2 , g 3 only, while the derivatives of the latter are expressed in terms of g 1 only. Another way of seeing this is to note that the functions f i have uniform weight zero, and the rescaling by powers of , which can be assigned weight −1, leads to the weight of the g i given above.
• The box g 6 and trianglesg 4 ,g 5 are fully decoupled in D = 4.
The block triangular structure of the matrix and its implication for the weight of the functions is visualized in a different way in Fig. 3. There one can also see the decoupling of the box and triangle integrals. The decoupling of the triangles implies that, if we are interested only in evaluating the box integral in D = 4, we could consistently work with the truncated basis g 1 , g 2 , g 3 and g 6 . Specifically, the system of differential equations in this case reduces to, cf. eq. (3.8), 11) in agreement with the general form of eq. (1.3). Such truncations will be exploited extensively in the next section. Let us now discuss the general solution in D = 4 − 2 dimensions, and then come back to the simplifications as → 0. With the differential equations in the form (3.6) we can immediately write down the analytic answer in terms of Chen iterated integrals [8]. We have f (s, t, m 2 ; ) = Pe γ dÃ h( ) . (3.12) Here the integration contour γ is a path in the space of kinematical variables, which begins at a base point, in our case m → ∞, where we have the simple boundary condition (3.9).
Let us be more specific about the notation, following closely the recent lecture notes [28,29] on iterated integrals. We denote by M the space of kinematical variables, here (u, v) ∈ R 2 , and let ω i be some differential one-forms (corresponding to entries of dÃ or d A). Moreover, define the pull-back of the differential forms to the unit interval [0, 1] via Then, an ordinary line integral is given by (3.14) Then, the iterated integral of ω 1 , . . . ω n along γ is defined by Iterated integrals have many nice properties, see [28,29]. Moreover, the iterated integrals appearing in eq. (3.12) 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. We will also discuss expansions in limits and their numerical evaluation. Expanding eq. (3.12) to second order in and noting that the box integral I 1 in eq. (2.6) is simply 1 2 2 βuv f 6 , it is trivial to reproduce its symbol given in eq. (2.9) starting from thẽ A-matrix. With a little bit more work, it is also possible to verify that I 1 indeed obeys the claimed differential equation and boundary condition.
More generally, one may expand the solution to any desired order in . With the normalization of the master integrals of eq. (3.2), the term multiplying k will be a Qlinear combination of k-fold iterated integrals, i.e. integrals of weight k. Such functions are sometimes called pure functions.
Let us now discuss the simplifications that occur for the four-dimensional basis (3.11). The solution can be again be written as in eq. (3.12), with the difference that is removed, and that only a finite number of terms have to be kept in the exponential.
The block diagonal form of d A implies the particularly simple iterative structure of the integrals that we already alluded to earlier when discussing the symbol. The iterative solution starts with the weight zero function 1, represented by the tadpole integral g 1 . To obtain the weight one integrals, we simply have to follow the arrows shown in Fig. 3 back up (starting from the tadpole). To each arrow corresponds an element of the d A-matrix, and the contribution is given by the integral over the latter. For example, As was explained before, the integration is along a contour γ in (u, v) space, starting at the boundary point u, v → ∞ corresponding to m 2 → ∞. Of course the above integral just evaluates to log βu−1 βu+1 , but this notion will be important when generalizing to higher weight functions. Next, the formulas for the weight two functions are obtained by summing over all paths from the tadpole to the integral under consideration. The contribution of each path is an iterated integral over the differential forms corresponding to each arrow. In this way, one arrives at Note that only the sum of the two iterated integrals in eq. (3.17) is homotopy invariant, not the individual terms. We can compare this to eq. (2.9), keeping in mind that g 6 = lim →0 βuv 2 I 1 . As was already mentioned, the difference is that, unlike the symbol, eq. (3.17) completely determines the answer.
We see that the solution is specified by the set of differential forms ω i allowed to appear in the iterated integrals. This alphabet can be read off directly from the dÃ or d A matrices. Of course, for multivalued functions depending on many kinematical variables, one can expect this alphabet, which is related to possible singularities of the functions, to be relatively complicated. As we will see, one can often choose specific contours γ that express the Chen iterated integrals in terms of more common multiple polylogarithms [1,2], at the cost of giving up the homotopy invariance and the compactness of the expressions. An example is given in appendix C. The general formulas are also very useful for discussing analytic continuation and for studying simplifying limits, where again one often finds expressions in terms of multiple polylogarithms. Many examples can be found in [17]. Discontinuities of the functions are also very transparent, see section 5.
In this section we have explained how to obtain a simplified version of the differential equations in the → 0 limit, starting from the D = 4 − 2 dimensional case. We saw that this implied that we could work with a smaller set of integrals which satisfy block-triangular systems of differential equations. In the next section, we will explain how to arrive at such results working directly in four dimensions.

Differential equation in D = 4: two and three loops
In D = 4 a natural subset of integrals to consider are the dual conformal ones, cf. eq (2.5). As we will see shortly, the derivatives of dual conformal integrals with respect to kinematic variables are themselves dual conformal. Therefore, one expects to be able to find a closed set of differential equations among dual conformal integrals. This has several advantages over the D-dimensional approach. 5 By restricting to a subset of integrals closed under differential equations, the number of integrals we need to consider in total decreases. This, and the possibility of setting D = 4 from the start results in much faster computer implementations of the algebra needed to derive the differential equations. This allowed us to extend the analysis to the three-loop level.
In this section we describe the method. The first step, to generate tables of integrationby-part identities among dual conformal integrals, is discussed in subsection 4.1. An interesting feature of working in D = 4 dimensions is that the identities sometimes contain contact terms, allowing different loop orders to be seamlessly merged. Once some basis of dual conformal master integrals has been found, our next step was to find a change of basis that puts the differential equation into the canonical block triangular form. It turns out that in D = 4 this problem is much simpler compared to the = 0 case, and can be solved systematically. We detail the procedure in subsection 4.2. We conclude this section with the main result for the two-and three-loop integrals.

Integration by part identities among dual conformal integrals
It is natural to expect that among all integration-by-parts identities, there exists a nonempty subset of identities, which hold exactly in D = 4, and which relate only dual conformal integrals to one another. Identifying this subset starting from a list of D-dimensional reduction identities among non-conformal integrals might seem like a daunting task. It turns out that this admits a rather elegant solution in terms of the embedding formalism.
The embedding formalism is ideally suited to this problem, because writing a nonconformal expression in this formalism requires a deliberate effort. The main idea is to write propagators as linear functions of a where the Y i 's obey certain equivalence relations. The formalism is described in appendix B. As an example, the one-loop integral of eq. (2.2) is rewritten in the form The numerator is necessary to ensure that the integrand has uniform homogeneity degree D under rescaling of Y 1 . The important feature here is that dual conformal integrals are precisely those for which the "point at infinity" I is absent. By simply never introducing this point, one automatically stays within the class of conformal integrals.
Since the (D+2)-vectors Y i are a trivial change of variable away from the original loop momenta k i , it is straightforward to carry through all the usual operations (generation of integration by parts identities, derivatives with respect to external variables, etc.) in terms of these variables. The upshot is that we automatically generate identities which involve only dual conformal integrals.
In order not to distract from the main course of this paper, these technical matters are discussed extensively in appendix B. Here we record only the most salient features in which D = 4 differs from D = 4: • Integration-by-part identities in D = 4 sometimes contain contact terms, which appear through the identity (B.13) This identity effectively relates integrals with different loop orders, since the δ-function removes one integration variable.
• When working without a regulator it is important to restrict attention to convergent integrals, which converge in all integration regions. To this aim we imposed a simple set of sufficient conditions on the set of integrals we have considered, also listed in appendix B.
• Derivatives of dual conformal integrals with respect to kinematic variables are themselves dual conformal, see e.g. eq. (B.18).
By generating identities among dual-conformal integrals in D = 4 and following standard integral reduction algorithms, we have reduced all integrals which appearing in derivatives to a minimal basis. Before we describe this basis, it is convenient to perform a certain basis change.

Algorithm for obtaining the differential in canonical form
Perhaps the most striking feature of working in D = 4 is that the A matrix, in a suitable basis, is expected to admit a block-triangular form. This reflects the existence of master integrals which have uniform transcendental weights. Assuming that such a basis exists, a constructive algorithm can find it! The method is actually rather straightforward: first we put the system into triangular form, by finding the proper normalization of (suitable combinations of) the master integrals. This already introduces a first approximation to the notion of weight. Then we recursively shift integrals by lower weight ones, until all differentials can be written in "d-log" form; the matrix will then automatically be block-triangular. 6 We explain the three steps of our method, illustrating them with simple examples.
Step 1. Obtaining a triangular form.
We decompose the integrals into two sets: correctly "normalised" and "remaining" ones. We assume that the constant integral (g 1 ) is part of the basis; by definition it is correctly normalised. We then successively normalise the remaining integrals as follows.
First we find the smallest subset of remaining integrals, whose derivatives involve only themselves and normalised ones. As a trivial example, in the two-loop case after we had normalised the integrals g 1 , g 2 , g 3 , g 4 (see below for their definitions), the smallest subset consisted of just one element: We integrate this equation, working modulo already normalised integrals: Since the right-hand side involves only other, already normalised integrals, we declare the parenthesis to be correctly normalised and we add it to the list. We then repeat the same procedure for the remaining integrals.
Sometimes the smallest closed set of remaining integrals contains more than one element. This happened to us at three loops; for example at some point we had the "entangled" pair of integrals h =(G 0,1,2,0,2,0,0,0,2,0,1,0,1,1,0 , G 0,2,1,0,2,0,0,0,2,0,1,0,1,1,0 ) with the differential equation The right-hand side can be removed by considering the ansatz I = h.(c 1 , c 2 ) and solving the differential equation d ds (c 1 , c 2 ) = −M (c 1 , c 2 ). 7 The solution gives that These two combinations can now be added to the list of normalised integrals, and the procedure repeated for the next subset. 8 It is noteworthy that even at three loops we never had to integrate anything more complicated than a 2 × 2 matrix. 9 This procedure successfully puts the differential equation into strict triangular form, which is already a tremendous simplification. We stress the critical role of D = 4 for this to be possible (see for example appendix A for a discussion of the D = 4 case). The integrals do not, yet, however, have uniform transcendental weight.
Step 2. Obtaining a block-triangular form.
We start from the triangular form of the differential equation and again split the integrals into two sets, "UT" (for uniform transcendental degree, or weight) and "remaining" ones. The UT integrals are assigned definite weights. At the beginning only the constant integral g 1 is UT, and its weight is defined to be 0.
We pick one of the remaining integrals whose derivatives involve only integrals in the UT set. Such integrals always exist, since the system is already strictly triangular. We must figure out the weight of that integral. To do this we integrate the coefficients of the maximal weight UT integrals on the right-hand-side of the differential equation and ask whether they contain transcendental parts or not. If they do, the weight of the integral will be one more than the weight on the right-hand side. If they do not, we remove the highest weight on the right-hand-side by a judicious redefinition, and repeat.
Let us illustrate this with a simple example. After g 1 , g 2 , g 3 , g 4 were put in UT form, we had the equation The presence of g 4 on the right-hand side, which has weight 2, would suggest that the parenthesis has weight 3. However, when we integrate the coefficient of g 4 , we find that it contains no logarithm: where c 1 is an (s-independent) integration constant. This shows that g 4 can be removed by a simple redefinition of the integral. After we do this, we are led to integrate the coefficient of g 2 , which is now the highest weight integral on the right-hand side (it has weight 1). This time we find that its primitive contains a logarithm: this cannot be removed by a simple redefinition. Including an integration constant for g 2 , and setting c 1 = 0 since there is no use for it anymore (the integral under consideration turned out to have the same weight as Finally we set c 2 = −2 to remove the last term, which would have the wrong weight. The parenthesis is now a UT integral of weight 2, one more than the weight of g 2 . We add it to the list and repeat the procedure for the next integral. The general procedure is as follows. Once we have identified the weight of the integral, we have two tasks: i) we must remove non-logarithmic terms in the highest-weight part of the derivatives ii) we must remove all terms involving lower-weight integrals (such as g 1 in the example). Both of these tasks are carried out order by order in decreasing weight. Non-logarithmic terms at high weights are removed by shifting the definition of the integral (like we did for g 4 above), while at lower weights it is also necessary to adjust integration constants from the previous weight. We believe that this procedure will always successfully find a UT basis if it exists.
Step 3. Optimizing the solution.
The UT master integrals generated by the above procedure can be quite lengthy and far from optimal. This is especially so for the highest weight UT integrals, which can end up expressed as combinations of all the other integrals in the original basis.
The good news is that once some UT basis has been found, it is very easy to find other, more compact, ones.
By construction, each of these may be reduced to linear combinations of UT master integrals. One can then search, among these integrals, for combinations which reduce directly to the parenthesis in eq. (4.7). That way we found the rather compact combination recorded as g 5 in eq. (4.8) below.
Contrary to the previous steps, there is of course no systematic procedure for finding "the nicest possible" representative. In most cases, we could find satisfactory combinations by trying a rather obvious set of candidates.

Result at two and three loops: UT basis and differential equation
The complete basis of four-dimensional DCI integrals for the family (2.3) is given as follows. First we have the four one-loop integrals   Figure 4. Hierarchy of master integrals up to two loops. The integrals are classified according to their (transcendental) weight, shown in the leftmost column. Each arrow corresponds to one non-zero element of the derivative matrix A, cf. eq. (4.11). The fact that arrows only link integrals in adjacent rows is the statement that the matrix is block triangular. The result for an integral can immediately be written down by summing over all paths leading up from the tadpole integral g 1 = 1. Each path gives a contribution to an iterated integral, with the integration kernels being specified by the 'letters' written next to the corresponding arrows. Solid and dashed lines denote massive and massless propagators, respectively. Note that the pictures are intended to give an idea of how the integrals look like, but omit details such as e.g. numerator factors.

They obey the differential equation in the canonical form (1.3), with the A-matrix
We have abbreviated each nonzero entry of the matrix by omitting the logarithm; for example A 2,1 ≡ log βu−1 βu+1 . The differential equation (1.3), with the A-matrix (4.11) and boundary condition (4.10), completely determines the integral. The block-triangular nature of the matrix, highlighted by its presentation, reflects the uniform transcendental weights of the integrals. The weights can be read off from the matrix structure: g 1 ; g 2,3 ; g 4,5,6,7 ; g 8,9 and g 10 have uniform weights 0,1,2,3 and 4.
As in the preceding section, the differential equation admits a natural graphical representation, see Fig. 4. The solution to the differential equation can be read off directly from this figure, as was explained in section 3.

Numerical evaluation
We would like to stress that, for practical purposes, Fig. 4 contains all relevant information about the double box integral.
Consider, for example, its numerical evaluation. As discussed in the end of section 3, each path between g 10 and g 1 translates to an iterated integral: The integration variables are ordered: 1>t 4 >t 3 >t 2 >t 1 >0. The shown term corresponds to the rightmost path in Fig. 4. In the end g 10 is independent of the choice of path going to (u, v) starting from infinity, but one simple choice is ( There are in total 10 paths in the figure, and the sum of the 10 associated terms represents g 10 exactly. Alternatively, the 10 terms can be generated by solving formally the differential equation dg i = (dA) ij g j , with the A-matrix given in eq. (4.11) and boundary condition lim u,v→∞ g i = δ i1 . Although it can be evaluated numerically, the four-fold integral (4.12) is obviously not optimal. Many integrals can be done analytically. For example, the t 1 and t 4 integrations obviously only produce logarithms, so at most two integrations really need to be done numerically. In fact, since analytic expressions for the weight 2 functions g 4,5,6,7 are easy to obtain, the t 2 integration can be also done analytically. In this way Fig. 4 immediately produces a one-fold integral representation for the double box: where everything in the square brackets is evaluated at (u(t), v(t)). In addition to g 6 = βuv 2 I 1 which was given previously in eq. (2.6), the functions are (x ≡ βu−1 βu+1 ): (log x) 2 + ζ 2 , (4.14) and g 7 = g 4 (v). We found the representation (4.13) to be very well suited for numerical evaluation, with the number of digits of accuracy increasing linearly with computing time (when implemented in Mathematica).
It is also possible to convert Fig. 4 automatically into expressions involving Goncharov polylogarithms. This exploits a judicious parametrization of the kinematic variables, as discussed in appendix C. These can then be evaluated numerically, using, for example [33,34]. In practice, we found the representation (4.13) to be both more direct to implement, easier to analytically continue and generally faster.
Before closing this section, we wish to emphasize that Fig. 4 is also a useful starting point for expanding the integrals in various kinematical limits. Often, the alphabet simplifies in these cases, and a representation in terms of simpler standard functions can be given. For example, let us consider the Regge limit, s → ∞, t fixed. Choosing the parametrization −t/m 2 = x/(1 − x) 2 for convenience, after a little algebra one sees that the alphabet required for all functions including the three-loop ones simplifies to This implies that the Regge expansion of all integrals, to any order in 1/s, can be given in terms of powers of s and log s, multiplied by harmonic polylogarithms [35] of argument x. Such expansions, as well as their physical interpretation, will be discussed in greater detail in ref. [17].

Result at three loops
A remarkable feature of the preceding computation is the surprisingly small amount of CPU power and memory usage that it required, as well as the fact that all manipulations are essentially algorithmic and with little need for trial-and-error. Furthermore the final form of the differential equation is pleasingly compact. These findings encouraged us to go to three loops. Indeed, the algorithms described above carried seamlessly through the three-loop integrals of the considered family, if only that more CPU power was required, as detailed in appendix B. For the two families in Fig. 2.4 combined (and including all one-and two-loop subtopologies), our final result for the UT basis contains 48 master integrals, 38 of them being relevant to the computation of the three-loop light-by-light amplitude [17] (i.e. the remaining 10 integrals decouple from the problem). In particular, the triple ladder integral, and the 'tennis court' integral with a single numerator are part of our basis, 1,1,0,1,0,1,0,1,0,1,1,1,1,0 , (4.16)  The system of differential equation, which takes a hierarchical form as found previously at one-and two-loops, is given for the 38 first integrals in appendix D. The full system, including the definitions of all the integrals in the basis, is attached as an ancillary file with the arXiv submission of this article. As before, this results contains all the information needed for efficient numerical evaluation of the integrals, the study of their limits, analytic continuation, etc. Here we simply note that the alphabet at three loops is not the same as at two loops. In addition to the entries in eq. (4.11) and their u ↔ v symmetric values, the A-matrix, and hence the alphabet, contains the new entries

Analyticity properties and Mandelstam representation
The differential equation can be most readily integrated in the Euclidean region u, v > 0, but this is not the only region of physical interest. To carry out their analytic continuation to other regions, it is important to know where the singularities lie.
A useful way to encode the analytic properties of an amplitude is through dispersivetype representations. The one-loop box integral (2.6) is known to admit the following remarkable integral representation [20,36] The integration region is the triangle ∆(ξ, η) = {ξ, η : ξ ≥ 0, η ≥ 0, ξ + η ≤ 1}. This representation makes the analyticity properties of the integral in the complex u and v planes completely manifest, because all the singularities originate from the denominator. The integral is therefore analytic for (u, v) ∈ (C\[−1, 0]) 2 . Furthermore the discontinuity with respect to u enjoys analytic properties with respect to v.
Representations of this type are known as Mandelstam representations, following an early proposal by Mandelstam [37]. 10 To our knowledge, the existence of the Mandelstam representation has never been conclusively proved nor disproved, not even to all orders in perturbation theory. There is a known perturbative obstruction at one-loop, namely the presence of so-called anomalous thresholds. These are singular branch points that are unrelated to the physical thresholds s, t = 4m 2 and which occur when the internal and external masses in a problem violate certain inequalities. In the equal mass case this obstruction is absent, in agreement indeed with the existence of the formula (5.1). Mandelstam proposed that if such a representation exists at one-loop for a given process, it should continue to exist at all loop orders for this process. We can test this hypothesis using the obtained differential equation.
We checked the analyticity properties of the integrals as follows. Starting from the point u, v = ∞, we integrate the differential equation along u with v = ∞. We find a single branch cut along u ∈ [−1, 0], as expected (in this interval s > 4m 2 ) and we evaluate the discontinuity across it. (In fact the u dependence at v = ∞ can be expressed analytically in terms of harmonic polylogarithms, similarly to the Regge limit as discussed around eq. (4.15), so at this stage everything is analytic.) We then integrate the discontinuity along the v direction, keeping u fixed. Although the differential equation has in principle many loci of singularity, we find (numerically) for all integrals in our three-loop basis only two branch points, located at v = 0 and v = −1 − u, indicating that the other singularities are spurious. This implies that all integrals, in our two-and three-loop families, possess the analyticity properties implied by the Mandelstam representation.
We should note that, for these statements to be true, it is important to consider integrals with rational coefficients, e.g. we have to divide the g i 's in eq. (4.8) by the square roots which appear in their normalization. Computing the discontinuity of the discontinuity, for v ∈ [−1 − u, 0], constructively produces the integrand of the Mandelstam representation.
Up to two loops we obtained closed-form analytic formulae: All of these integrals are manifestly analytic in the cut complex plane (u, v) ∈ (C\[−1, 0]) 2 . The square roots β u , β uv factored out in front of the integrals are precisely those which appear in the corresponding normalization factors in eq. (4.8), as just mentioned. The remaining integrals in our two-loop basis depend on only one variable at a time and admit the more familiar single-variable dispersion relations: .
The analyticity properties which we have verified for the three-loop integrals imply that they also admit Mandelstam representations, but we did not study these explicitly.
These dispersion relations can be written in a more familiar form using the change of variable u = 4m 2 /(−s), ξ = 4m 2 /s : Here we have included Feynman's i0 prescription, implicit in all preceding equations, but required to make sense of the integral when s > 4m 2 or t > 4m 2 . We recognize the righthand-side as a once-subtracted dispersion relation. The subtraction constants happen to vanish for all integrals in eqs. (5.3), because the left-hand-sides vanish at s = 0. Similar comments apply to the Mandelstam representations (5.2), where the one subtraction for g 9 is due to its nonvanishing value at t = 0.
We mention that at an early stage of this work, we actually successfully guessed an analytic expression for the two-loop integral g 10 , by simply assuming that it admitted a Mandelstam representation! The reason is that its integrand (i.e. the function h 10 above) has only transcendental weight 2, so by assuming that all entries of its symbol are taken from the finite alphabet (u, v, β u , β v , β u + 1, β v + 1, β uv + 1, β u + β uv , β v + β uv ) as suggested by the one-loop computation, one arrives at a relatively compact ansatz. We could then fix all coefficients using reality considerations, symmetries, Regge limits, and numerical evaluation at one point. However, we must warn the reader that such a procedure is far from systematic. At three loops our ansatz would most likely have missed the new entries of the symbol shown in eq. (4.18), and would thus have led us to a wrong answer, if any. The differential equation approach in contrast is completely systematic and reliable.

Checks
We have performed a number of checks on our calculations. First of all, there are many internal consistency checks for the differential equations, for example the integrability conditions (3.5) for the matrices appearing in different partial derivatives.
Moreover, as discussed in appendix A, the integrals we computed contain simpler twoscale form-factor type integrals that were previously computed [18]. We have performed tests against the results given in that reference. Our results also include (a subset of) two-loop integrals occurring in top quark pair production cross sections [38]. They can be obtained from our formulas by taking an sor t-channel discontinuity and going to certain forward limits. In particular, in this way we checked that the singularity structure of the integrals shown in Fig. 4(b)-(d) of the above reference can be understood by taking forward limits of the alphabet we found here in the general kinematical case. Moreover, we have done a detailed check for the cut double box integral shown in Fig. 4(b) of that reference, which corresponds to Disc v [g 10 ] (u, −u), up to a normalization factor. 11 Perhaps the easiest check can be done using the Mandelstam representation discussed in section 5, cf. eq. (5.2), which gives In this way, we numerically verified the result of [38] for that integral. It is also possible to give an analytic answer for this (and other) integrals, starting from the Chen iterated integral representation. In eq. (C.1) we give a convenient parametrization of u and v for rewriting those integrals in terms of multiple polylogarithms. Here it will suffice to remark that by specializing to the v = −u case gives u = 4z(1 − z 2 )/(1 + z 2 ) 2 , in which case the alphabet relevant for Disc v [g 10 ] (u, −u) reduces to z, 1 ± z, 1 + z 2 , i.e. a particularly simple case. Our results can also be used to derive similar formulas for other integrals in the same limit v → −u, including the three-loop ones, albeit using a larger alphabet.
Very strong tests were also provided by expansions of the integrals in various physical limits that will be discussed in detail in ref. [17]. For example, we have successfully reproduced the known Regge limit of the three-loop ladder and tennis court integrals g 37 and g 38 . Because of the iterative structure of the solution, this also tests to some extent the whole system of differential equations.
Finally we have performed numerical tests. The iterative structure of the differential equation allows for a simple numerical evaluation, as detailed in section 4.4. In the threeloop case, by working out analytic expressions for all weight 3 g i functions in terms of Li 3 functions, and for the t 5 , t 6 integrals in terms of dilogarithms, we actually succeeded in obtaining a one-fold integral representation similar to eq. (4.13), which also converge rather fast. We then performed numerical comparisons against values obtained from FIESTA [39,40]. While it is difficult to estimate the expected error of multidimensional integration routines, based on gradually increasing the number of sampling points and observing stable results, we arrived at the following values and error estimates We have also performed successful checks at other kinematical points, and for other integrals.

Conclusions
We have shown in this paper how to implement the differential equations method [3][4][5][6] for a closed subset of finite (and dual conformal) integrals. The restriction to D = 4 turns out to be a significant simplification of the method of ref. [7], which was originally developed for arbitrary D dimensions. We expect that our approach can also be applied for other classes of convergent fourdimensional integrals, such as for example position space correlation functions. For instance, in N = 4 super Yang-Mills, expressions for the integrands of four-point correlation functions are known to relatively high loop order [41], and our method could be used to shed light on their functional dependence and to evaluate them analytically. Currently, the planar integrals up to three loops are known [42]. Another set of finite and dual conformal functions that this method could be applied to are ratio and remainder functions in N = 4 super Yang-Mills [43]. 12 . Using a suitable four-dimensional regulator, such as for example differential renormalization [46,47], it seems reasonable to anticipate applicability of the method to correlation functions in an arbitrary quantum field theory.
It is important to stress that this approach to calculating integrals is systematic. We find it extremely satisfying that the algorithm developed in section 4 for the two-loop problem worked without any modification at three loops. However we do not have a proof that the algorithm will always succeed, and it will be interesting to try it on other problems, and also maybe generalize it to the D-dimensional case.
Our final result for the double box in D = 4 reveals a simple iterative structure, summarized in Fig. 4. In fact the figure contains all needed information about the integral! For example, it can be evaluated numerically almost directly, as explained in section 4.4. It also provides an efficient starting point for expanding around various kinematic limits, as will be pursued in a companion paper [17]. Finally, it can be automatically converted to a sum of standard special functions, for example Goncharov polylogarithms, as shown in appendix C. It is tempting to call Fig. 4 the answer for the double-box; it is the most concise and useful presentation which we can imagine on this day. The similar figure at three-loops is given in matrix form in appendix D.
It is tempting to expect such an iterative structure of the differential equation to be a general feature of Feynman integrals in integer dimensions. We expect this regardless of whether or not the integrals are expressible in terms of polylogarithm-type functions. For example, a similar iterative structure is also seen for integrals involving elliptic-type functions (see, for example, section 5 of [48]).
We comment that the method automatically finds the 'alphabet' relevant to a problem. This is important since the alphabet may change from one loop order to the next. The three-loop octagon Wilson loops of ref. [49] are a case in point; this paper provides another example, where the letters (4.18) would likely have been very difficult to guess. The alphabet specifies the class of two-variable functions describing the problem. Although the latter are rather complicated, we have demonstrated in this paper that they can be easily dealt with. In particular, we have shown how to evaluate them numerically, and discussed their analytic structure.
In addition to the expressions for all dual conformal invariant master integrals up to three loops, we have performed a complete calculation in dimensional regularization of the planar two-loop master integrals. The results, given in appendix A, can be used in future studies of light-by-light scattering at NLO, keeping the exact mass dependence. For previous studies in various mass regimes, see refs. [50,51].
The integrals we computed may also be helpful for tackling problems on the 2013 Les Houches QCD wishlist 13 , in particular when one wants to keep the exact dependence on quark masses. Indeed, the three-loop family we computed contain several integrals relevant 12 In that context, we wish to mention that differential equations involving integrals at different loop orders have a natural explanation in terms of supersymmetry [44,45]. Also, we would expect the differential equations of [15,16] to appear naturally as part of the system of differential equations described here.
13 https://phystev.in2p3.fr/wiki/2013:groups:sm:nnlowishlist for gg → H at NNNLO. Moreover, the process gg → gH at NNLO involves four-point functions similar to the ones computed here, albeit with one additional mass scale.
Here we report on the computation of all master integrals of the integral family (2.3), see Fig. 1(b), in D dimensions. Contrary to the main text, in this appendix we keep the full dependence on . These integrals may be used in future applications for light-by-light scattering in QED or QCD, and we also hope that they provide a first useful step towards similar but more complicated integrals needed e.g. for gg → Hg at NNLO.

A.1 Choice of integral basis
We apply the method of differential equations as explained in the main text. Using the methods of ref. [7], we arrived at the following convenient basis choice, where the normalization factor c is The basis integrals are depicted (qualitatively) in Fig. 5. Integrals f 1 to f 13 are s-channel triangle integrals, while the subsequent four integrals are t-channel triangle integrals. These integrals were computed previously, to some order in , in ref. [18]. The form of the basis we have chosen has the advantage that all basis elements have uniform weight and therefore are simpler compared to that reference. Some related integrals were also considered in ref. [52]. We wish to comment that, contrary to the D = 4 case, where we could employ the systematic algorithm presented in section 4, in the D-dimensional case finding the basis where 1.1 is satisfied required a fair amount of intuition, and trial and error. Also, the resultingÃ-matrix is no long block-triangular (nor even triangular).

A.2 Differential equations and analytic solution
With the above basis choice, we find the differential equations in the form of eq. (3.6), with A being a 29 × 29 matrix. It is given in electronic form as an ancillary file with the arXiv submission of this paper. As was discussed in the main text, the differential equations can be immediately solved, to any order in , in terms of Chen iterated integrals, cf. eq. (3.12). The relevant boundary condition is The alphabet, i.e. the set of entries of the matrix appearing in the differential equations turns out to be given by  with a new square root, in addition to the letters present in (4.11). We see that the full D-dimensional solution has a more complicated structure than just the four-dimensional dual conformal integrals.
Remarkably, we find that the only new singularity introduced by (A.32) w.r.t. the four-dimensional case is 1 + u + v = 0.

B Generation of IBP identities using the embedding formalism
Here we describe the generation of integration-by-parts identities among dual conformal integrals. We first briefly review the embedding formalism for arbitrary integrals in Ddimensions, following closely ref. [53], and then we set D = 4 and restrict to conformal integrals.

B.1 Loop integrals in the embedding formalism
To each point x µ in D-dimensional spacetime we assign a (D+2)-dimensional vector X a defined up to rescaling: The projective identification naturally removes one degree of freedom from the (D + 2)vector X. The second degree of freedom is removed by the fact that the vector (B.1) is null with respect to the (D + 2)-dimensional metric Thus eq. (B.1) defines a one-to-one map between null projective (D+2)-vectors and points in (compactified) D-dimensional spacetime. Furthermore, the dot product of the representative (B.1) evaluates to (XY ) = −(x − y) 2 , effectively linearizing the propagators.
To handle massive propagators in momentum space, we relax the null condition and modify X − in eq. (B.1) such that X 2 = m 2 . As a concrete example, the inverse propagators defined in eq. (2.1) map to the following (D+2)-vectors: 3) To make formulas invariant under choices of representative (B.1), it is useful to introduce an additional "infinity point" Assigning the vector Y 1 ∝ (k µ 1 , −k 2 1 , 1) to the loop integration variable, in this notation we have 1 All propagator factors are now linear in the integration variable Y 1 . The quadratic nature of the propagators has effectively been transferred to the integration measure: the (D + 2)-vector Y 1 is to be integrated over the null cone Y 2 1 = 0, subject to the GL(1) "gauge symmetry" Y 1 αY 1 . A natural integration measure for this space is . (B.6) Here the notation 1/vol(GL(1)) means to insert appropriate gauge-fixing factors and Faddeev-Popov determinant, which should integrate to unity over each GL(1) gauge orbit. The notation is borrowed from ref. [54]. A particularly simple gauge choice is 1/vol(GL(1)) → δ((Y I) − 1); the Faddeev-Popov determinant is then unity. 14 The general solution to the δ-function constraint in this gauge takes the form (B.1), and in this gauge it is easy to see that the above measure is proportional to the usual one: The left-hand-side is manifestly gauge-invariant (the GL(1) weight of the denominator cancels that of the measure in eq. (B.6)), so the identity does not depend on the gauge choice. For example we could gauge-fix using 1/vol(GL(1)) → δ(((Y ), I ) − 1) for any I , and although the identity would be more difficult to see, it would still be true since it is clear that nothing can depend on I .
Combining the ingredients (B.5) and (B.7) we can write the integrals in our one-loop family as The special simplification in the conformal case = 0 and i a i = 4 is manifest: the dependence on I drops out.
Generalization to higher loop planar integrals is straightforward: internal propagators become, for example, 1/(Y 1 Y 2 ). This will be sufficient for our purposes. It is unclear whether the formalism could be used for non-planar momentum space integrals, although we note that non-planar coordinate-space integrals would pose no problem.

B.2 Integration by parts identities: one-loop case
To derive integration-by-parts identities within the embedding formalism, we start from the usual fundamental identity [22], This expression makes sense (e.g. is invariant under the GL(1) gauge symmetry) provided that the vector is homogeneous: . We stress that the derivative does not need to act on the gauge fixing factor "1/vol(GL(1))" for this identity to be true.
To avoid taking derivatives of δ( 1 2 Y 2 ), however, which would have no interpretations in terms of Feynman integrals, one must require orthogonality: Y i V i (Y ) = 0. This is trivially 14 The triviality of the Faddeen-Popov determinant can be verified easily by integrating over a GL (1) orbit: dα α δ((αY I) − 1) = 1.
solved by V a ∝ Y a , but, due to the homogeneity property, the right-hand-side of (B.9) can be seen to become zero trivially in this case, so one obtains no useful identity in this case. We conclude that V a (Y ) is a homogeneous (D+2)-vector orthogonal to Y a and defined modulo Y a . This effectively reduces its number of vector components to D, as expected. This agrees with other treatments of vector fields in the embedding formalism, see [55,56] and references therein.
Let us now illustrate the rules with a simple example, involving a one-loop bubble: The left-hand side is an allowed IBP vector because the replacement ∂/∂Y a → Y a yields zero, as required. This justifies our omission of the explicit δ( 1 2 Y 2 ) factor, part of the definition of the measure (B.6). The right-hand side (in which we have used the symmetry G a 1 ,a 2 ,a 3 ,a 4 = G a 3 ,a 2 ,a 1 ,a 4 ) can be easily verified to be indeed a valid integral identity.
To illustrate nontrivial identities among conformal integrals, we consider the triangle integrals: where we have set = 0, and suppressed bubble and tadpole topologies from the righthand-sides. One can see that combining these identities allows the triangles to be removed. Further reducing the bubbles, one obtains in this way the following reductions: where the master integrals g 1 and g 2 , defined in eq. (4.8), have at most two propagators. These identities, which hold only when D = 4, explain the absence of conformal triangle integrals in the four-dimensional basis, cf. eq. (3.2).

B.3 Integration by parts identities: multi-loop case
The procedure we have just described generalizes trivially to the multi-loop case. However, new subtleties arise in the limit D → 4.
The first subtlety is that we must carefully avoid divergent integrals. Ultraviolet convergence is already guaranteed for dual conformal integrals, to which we will restrict our attention when discussing D = 4. Infrared divergences, however, could appear as a result of squaring the (massless!) internal propagators.
It is important to stress that, when working without a regulator, one should consider only integrals which are not only finite but convergent, e.g. separately finite in each integration region [54,57].
As a simple sufficient condition for convergence, we have imposed the following constraints on the integrals in our basis: • At two-loops, we require the internal propagator to occur with at most two powers, 1/(Y 1 Y 2 ) 2 . Squared internal propagators, when they appear, must multiply numerators which vanish when Y 1 = Y 2 .
• At three-loops, we impose the preceding conditions separately for each internal propagator. Convergence in the region Y 1 → Y 2 → Y 3 is then ensured by requiring that no more than one squared propagators appears.
These conditions define a restricted set of integrals, within which one can safely operate without regulator in D = 4. A second subtlety when D = 4 is related to contact terms. These appear through the familiar identity which holds when D is exactly equal to four. As a simple way to derive the corresponding identity in the embedding formalism, we may the following result for logarithmically divergent divergent integrals: (B.14) Imagine, now, that we have an IBP vector that contains a squared internal propagator (but with a numerator such that it generates only convergent conformal integrals when D = 4). When we evaluate the IBP identity for = 0, we get all the terms we get naively from D = 4, plus some -dependent terms. Actually the only dependence in the formalism comes from the derivative acting on (Y, I) 2 , producing We notice the / cancelation, which means that this term survives the → 0 limit. This additional term can be interpreted as a contact term as in eq. (B.13): where the δ-function is normalized according to the measure (B.7): The "regular terms" are all those obtained by the naive differentiation procedure in D = 4. The presence of the infinity point in this formula is surprising at first sight, but it has a simple interpretation. To see this, we can consider, as an alternative to dimensional regularization, regulating (B.9) by excising a small neighborhood of Y 1 = Y 2 . The contact term will then arise as a boundary term in the integration by part identity, coming from a small three-sphere surrounding Y 1 = Y 2 . One can verify that this boundary term is independent of the shape of the (small) boundary if and only if V (4)a (Y, Y ) ∝ Y a . This is precisely the condition for (B.16) to be independent of I! Only IBP vectors satisfying this condition have unambiguous meaning without a regulator.
This completes our description of the class of integration-by-part identities we have considered, which generates valid identities among dual conformal integrals defined in D = 4. For the reader's convenience, we recapitulate the full list of constraints: • The identities are generated by vectors ∂ ∂Y a k V a k as in eq. (B.9), where V k is orthogonal to Y k and homogeneous with respect to the loop variables.
• The vector V k does not contain any cubed internal propagator. If it contains a squared internal propagator 1/(Y 1 , Y 2 ) 2 , its numerator becomes proportional to Y a 1 when Y 1 = Y 2 . Furthermore the resulting identities should involve integrals which obey the convergence properties listed at the beginning of this subsection.
Generating identities satisfying these rules is now a simple linear algebra problem. In practice, we used a computer to generate all possible vectors subject to prescribed bounds on the various powers which can appear in the denominator (and hence in the numerator, by the DCI constraint). Within this set the constraints become a set of linear conditions, which can be solved by Gaussian elimination (using, for example, Mathematica's NullSpace[] routine).
We stress that we have not tried to generate the most general possible IBP identities. The restricted class we considered turned out to be sufficient for reducing all two-and threeloop conformal integrals we considered to a minimal basis!

B.4 Derivatives with respect to external parameters
For our restriction to conformal integrals to make sense, it is important that we be able to express all derivative of conformal integrals, in terms of conformal integrals.
This property is actually rather trivial within the embedding formalism. Here we record the relevant formula. The first step is to write the required differential operators, see e.g. (3.1), in terms of the X i variables defined in eq. (B.3). The correct vectors are easily found if one notes that they should preserve the normalization and on-shell constraints X 2 i = m 2 and (X i X i+1 ) = 0, which are satisfied by these expressions. For example, we have The derivative of an integral is obtained immediately. For example, we have a 2 ,a 3 ,a 4 + G a 1 +1,a 2 −1,a 3 ,a 4 + t s G a 1 +1,a 2 ,a 3 −1,a 4 + G a 1 +1,a 2 ,a 3 ,a 4 −1 .
It is apparent that conformal integrals are mapped to conformal integrals.

B.5 Simple implementation of an integral reduction scheme
Well known algorithms and public computer packages allow for the generation and management of integration-by-parts identities, see e.g. [23][24][25][26]. To deal with the tables of identities generated manually as just described, we employed a private implementation of the most basic ideas behind the public packages. We describe it here. The task is to convert raw lists of (mostly useless) identities to useful reductions, which express "complicated" integrals in terms of "simpler" ones. The central tool in this task is Gaussian elimination, as also used in the Laporta algorithm [58]. Specifically, to add new identities to the system, we first reduce them using existing ones, and then rank the surviving integrals in terms of their "complexity". Gaussian elimination is then used to express the most complicated remaining integrals in terms of simpler ones, as much as possible. In our implementation, we kept all identities: the total number of identities stored in memory at a given time equals the total number of linearly independent identities generated so far.
Of critical importance for the success of this strategy is the choice of "complexity" function. We obtained the best results when we tuned it such that the algorithm strived to reduce, in the following order of priority: the number of loops; the number of internal propagators; the powers of each internal propagator; the number of massive propagators; the powers of each propagator. (For conformal integrals, the order of the numerator is not an independent number.) For indication, generating a sufficient set of reductions to compute all derivatives of the two-loop ladder and its sub-topologies, turned out to require approximately 500 identities and about 30 seconds of CPU time on a single laptop computer (endowed with a 2.6GHz processor and 8Go RAM). The three-loop ladder required approximately 2000 identities, of considerably longer complexity, but still only on the order of one hour CPU time. The "tennis court" topology, on the other hand, required over 20000 identities and several days of CPU time. We attribute the relative difficulty in this case to the rather conservative convergence conditions we imposed on our IBP vectors; hopefully this situation could be improved by relaxing these conditions, as discussed further in the conclusion section.

C One-and two-loop box integrals in terms of Goncharov polylogarithms
It turns out that, for the alphabet that occurs in D = 4 at both two-and three-loops, one can perform a change of variables that linearizes it. That is, the three square roots √ 1 + u, √ 1 + v and √ 1 + u + v are simultaneously removable. While we did not use this in the main text, we wish to mention it here for convenience of the reader. This allows for a straightforward expression of all D = 4 integrals in terms of Goncharov polylogarithms. Specifically we are going to discuss the one-and two-loop box integrals at = 0, i.e. g 6 and g 10 .
The change of variable is the following: (C.1) We take w, z real and 0 < w < z < 1 , (C.2) which corresponds to the Euclidean region s/m 2 < 0, t/m 2 < 0, where all functions are real-valued. In terms of those variables, the alphabet needed up to two loops, eq. (4.11), can be written as {z , 1 ± z , w , 1 ± w , w − z , 1 − wz , 1 + w − z + wz , 1 − w + z + wz} . Note that all letters appearing in the alphabet (C.3) are multilinear in w, z. This property implies that we can always write the corresponding integral functions in terms of multiple polylogarithms, as we will see presently. The latter can be defined iteratively as follows, G(a 1 , . . . a n ; z) = For a 1 = 0, we have G( 0 n ; x) = 1/n! log n (x). A subset of these functions for indices 0, ±1 are called harmonic polylogarithms (HPL) [35], and will be denoted by H(a 1 , . . . a n ; z). 15 There are many ways of doing this which lead to different integral representations. A straightforward possibility is to choose the integration contour C of eq. (3.12) to consist of two segments that are straight lines along the w and the z direction, respectively. In formulas, denoting base point and argument of the function by (w i , z i ) and (w f , z f ), respectively, we parametrize the two segments as follows,

D Differential equation at three loops
In this appendix we record the full form of the differential equation at three-loops, for the DCI integrals in four dimensions. Out of the 48 master integrals, we show results only for the 38 of them which are related to the integrals g 37 , g 38 , defined in eqs. (4.16) and (4.17), which control the 3-loop light-by-light amplitude in N = 4. This set includes the 10 twoloop integrals g 1 . . . g 10 as well as 4 other related by s ↔ t symmetry, plus 24 genuine 3-loop integrals which are defined in an ancillary that file also contains the A-matrix reproduced here.
Introducing the abbreviations, our result for the complete three-loop hierarchy can be written: g 12 g 13 g 15 g 16 g 17 g 18 g 19 g 20 g 21