Analytic results for the planar double box integral relevant to top-pair production with a closed top loop

In this article we give the details on the analytic calculation of the master integrals for the planar double box integral relevant to top-pair production with a closed top loop. We show that these integrals can be computed systematically to all order in the dimensional regularisation parameter $\varepsilon$. This is done by transforming the system of differential equations into a form linear in $\varepsilon$, where the $\varepsilon^0$-part is a strictly lower triangular matrix. Explicit results in terms of iterated integrals are presented for the terms relevant to NNLO calculations.

In realistic scattering processes we face in general Feynman integrals, which depend on multiple scales. A prominent example is given by the planar double box integral for tt-production with a closed top loop. This integral enters the next-to-next-to-leading order (NNLO) contribution for the process pp → tt. Until quite recently, it has not been known analytically. The existing NNLO calculation for the process pp → tt treats this integral numerically [47][48][49][50][51]. This integral is clearly a cornerstone and should be understood for further progress on the analytical side. In a recent letter we reported how this integral can be treated analytically [52]. With this longer article we would like to give all the technical details.
Our starting point is the method of differential equations [53][54][55][56][57][58][59][60][61][62][63] for the master integrals. In the modern incarnation of this method one tries to find a basis of master integrals J, such that the system of differential equations is in ε-form [59]: where A does not depend on ε. If such a form is achieved, a solution in terms of iterated integrals is immediate, supplemented by appropriate boundary conditions. This strategy has successfully been applied to many Feynman integrals evaluating to multiple polylogarithms. The difficulty of Feynman integral calculations is therefore reduced to finding the right basis of master integrals.
In the case of multiple polylogarithms the transformation from a Laporta basis I to the basis J is rational in the kinematic variables and several algorithms exist to find such a transformation [62,[64][65][66][67][68][69][70][71][72]. In ref. [26] we showed that an ε-form can even be achieved for the equal-mass sunrise / kite system. This is made possible by enlarging the set of transformations from the Laporta basis I to the basis J from rational functions in the kinematic variables towards rational functions in the kinematic variables, the periods of the elliptic curve and their derivatives. We may slightly relax the form of the differential equation and consider where A (0) is strictly lower-triangular and A (0) and A (1) are independent of ε. This does not spoil the property, that the system of differential equations is easily solved in terms of iterated integrals. Since A (0) is strictly lower triangular, one can easily transform the system to an εform by introducing primitives for the terms occurring in the ε 0 -part. In this paper we give for the system of the double box integral a transformation matrix, which transforms the system from a pre-canonical form to a form linear in ε, as in eq. (2). The transformation is rational in the kinematic variables, the periods of the elliptic curves and their derivatives. A subsequent transformation, which brings the system into ε-form is possible, however this would introduce additional transcendental functions. We choose the basis of master integrals J such that A (0) vanishes if either t = m 2 or s = ∞. In the former case the solution reduces to multiple polylogarithms, in the latter case the solution reduces to iterated integrals of modular forms already encountered in the sunrise / kite system.
The attentive reader might have noticed that we put "elliptic curves" into the plural. The system of differential equations for the double box integral is not governed by a single elliptic curve. We find that there are three elliptic curves involved, originating from different sub-topologies. We show how the elliptic curves can be extracted from the maximal cuts. This fact has consequences: Elliptic multiple polylogarithms are defined as iterated integrals on a punctured elliptic curve [14-17, 22, 24, 35-38, 73-78]. Inherent to this definition is the notion of a single elliptic curve. Since in our problem there are three distinct elliptic curves involved, we do not expect our results to be expressible in terms of elliptic multiple polylogarithms (which by definition are tied to one specific elliptic curve). We express our results as iterated integrals of the integration kernels appearing in eq. (2). Let us however also stress the other side of the medal: Our analysis also shows that nothing worse than elliptic curves appears in the calculation. This paper is organised as follows: In section 2 we introduce our notation and review a few basic facts about Feynman integrals. In section 3 we discuss iterated integrals. There are two special cases of iterated integrals, which we briefly review: multiple polylogarithms and iterated integrals of modular forms. In section 4 we discuss in detail the kinematic variables related to the subset of Feynman integrals, which evaluate to multiple polylogarithms. Section 5 is devoted to elliptic curves. After a discussion of the generic quartic case, we show how to extract the elliptic curves from the maximal cuts. We discuss the three occurring elliptic curves in detail. In section 6 we define the basis of master integrals J. The system of differential equations for this set of master integrals is given in section 7. In section 8 we show how the boundary values are obtained. The results for the master integrals up to order ε 4 are given in section 9. Finally, our conclusions are given in section 10. The article is supplemented by an appendix. In appendix A we show Feynman graphs for all master topologies. In appendix B we give an extra relation, which reduces the number of master integrals in the sector 93 from five to four. In appendix C we collected useful information on the modular forms occurring in the s → ∞ limit. Appendix D lists the full set of bondary constants. Appendix E describes the supplementary electronic file attached to this article, which gives the definition of the master integrals, the differential equation and the results in electronic form.
A sector (or topology) is defined by the set of propagators with positive exponents. We define a sector id (or topology id) by id = 9 ∑ j=1 2 j−1 Θ(ν j ).

Chains and cycles
It will be useful to group the internal propagators P j into chains [79]. Two propagators belong to the same chain, if their momenta differ only by a linear combination of the external momenta. Obviously, each internal line can only belong to one chain. In fig. (2) we have three chains: C (1) = {P 8 , P 3 , P 1 , P 2 } , C (2) = {P 9 , P 6 , P 7 , P 5 } , We define a cycle to be a closed circuit in the diagram. We can denote a cycle by specifying the chains which belong to the cycle. In the two-loop diagram of fig. (2) there are three different cycles, given by C (13) , C (23) , C (12) .
Here we used the notation that C (i j) denotes the cycle consisting of the chains C (i) and C ( j) . A cycle corresponds to a one-loop sub-graph. The discussion carries over to sub-topologies of eq. (5) by deleting the appropriate propagators from the chains C (1) , C (2) and C (3) . If a chain is empty, the two-loop integral factorises into two one-loop integrals. For non-empty chains C (1) , C (2) and C (3) we are interested in a cycle with a minimum number of propagators. A cycle with a minimum number of propagators simplifies the calculation of the maximal cut of a Feynman integral within the loop-by-loop approach. The choice of such a cycle may not be unique, for example for the sunrise topology all three cycles have two propagators. For the auxiliary topology shown in fig. (2) and all non-trivial sub-topologies (i.e. sub-topologies which are not products of one-loop integrals) it is always possible to choose either as a cycle with a minimal number of propagators. This is due to the fact that the chain C (3) contains already the minimal number of propagators for a non-trivial chain.

The Feynman parameter representation
It is useful to discuss the Feynman integral representation for the double box integral. The Feynman parameter integral for ν 8 = ν 9 = 0 reads where the integration is over The differential form ω is given by where the hat indicates that the corresponding term is omitted. The graph polynomials are given by The graph polynomial U reads in expanded form Let us further define the derivatives of the graph polynomial F with respect to s and t by The graph polynomials for sub-topologies are obtained from U and F by setting the Feynman parameters to zero which correspond to propagators not present in the sub-topology.

Dimensional shift relations and differential equations
Let us introduce an operator i + , which raises the power of the propagator i by one, e.g.
The dimensional shift relations for integrals with irreducible numerators (i.e. ν 8 < 0 or ν 9 < 0) can be obtained as follows: One first converts to a basis of master integrals with ν 8 = ν 9 = 0 (and raised propagators), applies the dimensional shift relations to the latter and converts back to the original basis. The differential equations for ν 1 , ..., ν 7 ≥ 0 can be obtained from The right-hand side is given by integrals in (D + 2) dimensions with three propagators raised by an additional unit. Reducing these integrals to the basis in D dimensions gives the differential equation. Let us mention that eq. (23) allows us to write the differential equations in a compact way. However, from a computational point of view it is not the most advantageous representation, since it requires reduction of integrals with large ν.
We use the programs Reduze [82], Kira [83], Fire [84] and LiteRed [85,86] to reduce the integrals to master integrals. These programs are based on integration-by-parts identities [87,88] and implement the Laporta algorithm [89]. There are 44 master integrals. We may choose a basis of master integrals for which the auxiliary propagators P 8 and P 9 are not present, at the expense of having higher powers of the propagators for the propagators P 1 -P 7 . It is therefore possible to label these master integrals by I ν 1 ν 2 ν 3 ν 4 ν 5 ν 6 ν 7 . In some sectors we have more than one master integral. The list of master integrals is shown in table 1. The master topologies are shown in appendix A.
We may change the basis of master integrals. Let us denote by I the vector of the precanonical master integrals. The differential equation for I reads for µ = m The matrix-valued one-form A satisfies the integrability condition Under a change of basis [59] one obtains where the matrix A ′ is related to A by A comment is in order: The statement that there are 44 master integrals (and not 45) is already a non-trivial statement. We first run Reduze, Kira and Fire (the latter in combination with LiteRed [85,86] and without). Taking trivial symmetry relations into account, all programs give with standard settings 45 master integrals. However, the reductions disagree for the three most complicated topologies and at first sight the results of two of the three programs seem to violate the integrability condition eq. (25). All these symptoms are resolved once an additional relations is taken into account. The relation is given in appendix B. This relation reduces the number of master integrals in sector 93 from 5 to 4 (and in turn the total number of master integrals from 45 to 44). Imposing this relation, the results from Reduze, Kira and Fire agree and the integrability condition is satisfied. We have found this relation by comparing the output of the three programs above. All differences are proportional to a single equation. In addition, we verified numerically the first few terms in the ε-expansion of this relation. This extra relation comes from a higher sector (i.e. sector 123). We would like to mention that Reduze is able to find the relation and can be forced to use this relation with the command distribute_external 1 . Let us also mention that MINT [90] and AZURITE [91] are programs, which can be used to count the number of master integrals. MINT analyses critical points, AZURITE is based on syzygy relations. Applied to our problem, MINT reports 4 master integrals for sector 93, AZURITE gives 5. Table 1: Overview of the set of master integrals. The first column denotes the number of propagators, the second column labels consecutively the sectors or topologies, the third column gives the sector id (defined in eq. (10)), the fourth column lists the master integrals in the basis I, the fifth column the corresponding ones in the basis J. The last column denotes the kinematic dependence.

Iterated integrals
Let us first review Chen's definition of iterated integrals [92]: Let M be a n-dimensional manifold and a path with start point x i = γ(0) and end point x f = γ(1). Suppose further that ω 1 , ..., ω k are differential 1-forms on M. Let us write for the pull-backs to the interval [0, 1]. For λ ∈ [0, 1] the k-fold iterated integral of ω 1 , ..., ω k along the path γ is defined by We define the 0-fold iterated integral to be We have Let us now specialise to our case of interest: Without loss of generality we may set µ = m in eq. (5). Then the Feynman integrals I ( j) ν 1 ν 2 ν 3 ν 4 ν 5 ν 6 ν 7 appearing in the Laurent expansion of eq. (9) depend only on two dimensionless ratios, which may be taken as In other words, we may view the integrals I ( j) denote the homogeneous coordinates. We will express I ( j) ν 1 ν 2 ν 3 ν 4 ν 5 ν 6 ν 7 as iterated integrals on P 2 (C).
Note that we are free to choose any convenient coordinates on M. One possibility is given by eq. (34). We will refer to this choice as (s,t)-coordinates.
A second possibility is given by the set (x, y), where x and y are related to s and t by We will refer to this choice as (x, y)-coordinates. The (s,t)-coordinates and the (x, y)-coordinates will be our main coordinate systems, although not the only ones. The (s,t)-coordinates are closest to physics, however in these coordinates we encounter square roots already in very simple sub-topologies. The (x, y)-coordinates rationalise the most prominent square root −s(4m 2 − s). However, this is not the only occurring square root. For example, we also encounter the square root −s(−4m 2 − s). In order to simultaneously rationalise both square roots we use the coordinates (x, y), wherex is defined by However, there is a price to pay: Rationalising square roots will increase the degree of the polynomials in intermediate stages of the calculation. For this reason we work bottom-up and treat each sub-topology in a coordinate system adapted to this sub-topology. On top of this there are several elliptic topologies. Here we use the modular parameter τ of the associated elliptic curve as one variable.

Multiple polylogarithms
Multiple polylogarithms are a special case of iterated integrals. For z k = 0 they are defined by [93][94][95][96] G(z 1 , ..., z k ; y) = y 0 dy 1 y 1 − z 1 The number k is referred to as the depth of the integral representation or the weight of the multiple polylogarithm. Let us introduce the short-hand notation G m 1 ,...,m k (z 1 , ..., z k ; y) = G(0, ..., 0 where all z j for j = 1, ..., k are assumed to be non-zero. This allows us to relate the integral representation of the multiple polylogarithms to the sum representation of the multiple polylogarithms. The sum representation is defined by The number k is referred to as the depth of the sum representation of the multiple polylogarithm, the weight is now given by m 1 + m 2 + ...m k . The relations between the two representations are given by Li m 1 ,...,m k (x 1 , ..., x k ) = (−1) k G m 1 ,...,m k 1 x 1 x 2 , ..., 1 Note that in the integral representation one variable is redundant due to the following scaling relation (recall z k = 0): G(z 1 , ..., z k ; y) = G(xz 1 , ..., xz k ; xy).
One can slightly enlarge the set of multiple polylogarithms and define G(0, ..., 0; y) with k zeros for z 1 to z k to be This permits us to allow trailing zeros in the sequence (z 1 , ..., z k ) by defining the function G with trailing zeros via eq. (44) and eq. (45). Please note that the scaling relation eq. (42) does not hold for multiple polylogarithms with trailing zeros. Using the shuffle product it is possible to remove trailing zeros and to express any multiple polylogarithms as a linear combination of terms which involve multiple polylogarithms without trailing zeros and powers of ln y. It will be convenient to introduce the following notation: For differential one-forms we define G(ω 1 , ..., ω k ; y) recursively through Methods for the numerical evaluation of multiple polylogarithms can be found in [97].

Iterated integrals of modular forms
A second special case of iterated integrals are iterated integrals of modular forms. Let f 1 (τ), f 2 (τ), ..., f k (τ) be modular forms of a congruence subgroup. The most prominent example in our application will be the congruence subgroup Γ 1 (6). Let us further assume that f k (τ) vanishes at the cusp τ = i∞. We define the k-fold iterated integral by The case where f k (τ) does not vanishes at the cusp τ = i∞ is discussed in [21,98] and is similar to trailing zeros in the case of multiple polylogarithms. This is easily seen by changing the variable from τ to q: Modular forms have a Fourier expansion around the cusp τ = i∞: f j (τ) vanishes at τ = i∞ if a j,0 = 0. Using the Fourier expansion and integrating term-by-term one obtains the q-series of the iterated integral of modular forms corresponding to eq. (48): a 1,n 1 n 1 + ... + n k ... a k−1,n k−1 n k−1 + n k a k,n k n k q n 1 +...+n k .
4 The kinematic variables for the multiple polylogarithms A large fraction of the integrals under consideration will only depend on the kinematic variable s, but not on t. These integrals are expressible in terms of multiple polylogarithms. Furthermore one finds that all integrals under consideration are expressible in terms of multiple polylogarithms in the special kinematic configuration t = m 2 . Let us now introduce several variables related to the multiple polylogarithms. They all replace the kinematic variable s and rationalise one or several square roots. The difference lies in the square roots they rationalise. Let us start with the simplest case relevant to integrals with a singular point at s = 4m 2 . The variable x replaces s and is defined by [99][100][101][102] The interval s ∈] − ∞, 0] is mapped to x ∈ [0, 1], with the point s = −∞ being mapped to x = 0 and the point s = 0 being mapped to x = 1. This change of variables rationalises the square root −s(4m 2 − s). In more detail we have We will encounter sub-topologies, which have a singular point at s = −4m 2 . The simplest example is given by sector 57. For these integrals we change from the variable s to a variable x ′ defined through The This change of variables rationalises the square root −s(−4m 2 − s), but not the square root −s(4m 2 − s). We have We further have In order to rationalise simultaneously the two square roots −s(4m 2 − s) and −s(−4m 2 − s) we introduce a variablex through Expressing x ′ in terms ofx yields We have In eq. (59) we defined for later convenience the differential forms ω 0 , ω 4 , ω −4 , ω 0,4 and ω −4,0 . Note that 2xdx From eq. (53), eq. (55) and eq. (59) we may read off the alphabets A, A ′ andÃ for the variables x, x ′ andx, respectively. We have Thus, iterated integrals in the variable x involving the differential forms of eq. (53) may be expressed in terms of the smaller class of harmonic polylogarithms [103,104]. The same holds true for iterated integrals in the variable x ′ involving the differential forms of eq. (55). On the other hand, iterated integrals in the variablex involving the differential forms of eq. (59) have a larger alphabet and are expressed in terms of multiple polylogarithms [93][94][95][96].
In practice, the results for the more complicated integrals are most compactly expressed by introducing the notation of eq. (47). All sub-topologies, which depend only on s, can be expressed as iterated integrals with integration kernels given by the five differential one-forms In addition, for t = m 2 (or equivalently y = 1) all master integrals can be expressed as iterated integrals with these integration kernels. From eq. (59) and eq. (60) it is clear that all iterated integrals in these integration kernels are expressible in terms of multiple polylogarithms.

Elliptic curves
In this section we discuss elliptic curves. We start with a review of the general quartic case in section 5.1. The relevant elliptic curves are extracted from the maximal cuts. This is done in section 5.2. We find three different elliptic curves, which we label E (a) , E (b) and E (c) . These curves are discussed individually in sections 5.3 -5.5.

The general quartic case
Let us consider the elliptic curve where the roots z j may depend on variables x = (x 1 , ..., x n ): We use the notation We set Note that we have We define the modulus and the complementary modulus of the elliptic curve E by Note that there are six possibilities of defining k 2 . Our standard choice for the periods and quasiperiods is These periods satisfy the first-order system of differential equations and the Legendre relation The parameter τ and the nome squared q are defined by We have Let us now consider a path γ : where the variable λ parametrises the path. A specific example is the path γ α : [0, 1] → C n , indexed by α = [α 1 : ... : α n ] ∈ CP n−1 and given explicitly by For a path γ we may view the periods ψ 1 and ψ 2 as functions of the variable λ. We then have where This defines the Picard-Fuchs operator along the path γ: The Wronskian is defined by We have Let us now specify to the case, where the base space is given by the two variables (x, y). Eq. (76) and eq. (78) allows us to obtain the Picard-Fuchs operator and the Wronskian from the roots z 1 , z 2 , z 3 and z 4 for the variation of the elliptic curve along the paths We define the Wronskians W x and W y at the point (x, y) as the derivatives in the directions x and y, respectively. Thus We have We may use the Picard-Fuchs operator to eliminate second derivatives: where the subscript x + y refers to the path with (α 1 , α 2 ) = (1, 1). This leaves us with There is a further relation, since we may exchange any derivative of ψ 1 in favour of φ 1 : This yields Using eq. (86) we may eliminate one derivative, say d dx ψ 1 . This leaves us with as expected, since the first cohomology group of an elliptic curve is two dimensional.

Maximal cuts
In order to identify the elliptic curves associated to a Feynman integrals we use maximal cuts in the Baikov representation [105][106][107][108][109][110][111]. For the maximal cut of an integral we use the loop-byloop approach [109]. Let us first assume that all propagators occur to the power one, although we allow irreducible numerators. We first consider a one-loop sub-graph with a minimal number of propagators. This is equivalent to the statement that the dimension of the sub-space spanned by the external momenta for this sub-graph is minimal. Let us assume that this sub-graph contains (e + 1) propagators P 1 , ..., P e+1 , therefore the sub-space spanned by the external momenta for this sub-graph has dimension e. We change the integration variables for this sub-graph according to where the momenta p 1 , ..., p e denote the linearly independent external momenta for this subgraph, the Gram determinant (in Minkowski space) is defined by and u denotes an (irrelevant) phase (|u| = 1). The maximal cuts are solutions of the homogeneous differential equations [112]. Any such solution remains a solution upon multiplication with a non-zero constant. Therefore the phase u is not relevant. We then repeat this procedure for the second loop, replacing p 1 , ..., p e by the set of independent external momenta for the full graph. For an integral of the form where N(k 1 , k 2 ) is a polynomial in k 1 and k 2 , a maximal cut is given by where the integration measure is re-written according to eq. (88) and the integration is over a (yet to be) specified contour in the variables P j not eliminated by the delta distributions. The exact definition of the integration contour is not relevant for the extraction of the elliptic curve from the maximal cut. We aim for a one-dimensional integral representation for the maximal cut with a constant in the numerator and a square root of a quartic polynomial in the denominator. This defines the elliptic curve. The possible choices for the integration contour are then given by an integration between any pair of roots of the quartic polynomial. In practice we label/order the roots and define the periods by eq. (69). For integrals with ν j > 1 we may compute a maximal cut by first converting to a basis with ν j = 1 and possibly irreducible numerators and then computing the maximal cut in this basis.
Alternatively, we may interpret the delta distribution δ(P j ) as a contour integration along a small circle around P j = 0. We may therefore compute the maximal cut from residues.
Let us look at a few examples. We start with the equal mass sunrise integral (sector 73) in two space-time dimensions. Starting with the sub-loop C 1 first we obtain For the sunrise integral we could equally well start with the sub-loop C 2 . Doing so we find The two representations are related by P ′ = P − 2m 2 . Let us now look at the maximal cut of the double box integral (sector 127), this time in four space-time dimensions. We have One recognises in eq. (92), eq. (93) and eq. (94) the typical period integrals of an elliptic curve. We note that the integrand of eq. (94) differs from the one of eq. (93). The difference is given by the additional term This term vanishes in the limit s → ∞.
Let us now look at the maximal cut in the sector 79. We find with P = P 8 Up to the prefactor, this is the same maximal cut integral as in eq. (94). Therefore the sectors 79 and 127 are associated to the same elliptic curve.
Our next example is the maximal cut in the sector 121. Here we find with P = P 9 + 2m 2 This corresponds to an elliptic curve different from the one found in sectors 79 and 127. In the limit s → ∞ the maximal cut integral reduces again up to a prefactor to the one of eq. (93). The most complicated example is the maximal cut in sector 93. For this sector we find first within the loop-by-loop approach a two-fold integral representation in P 2 and P 8 for the maximal cut. The integrand has a single pole at P 2 = 0. Choosing as a contour for the P 2 -integration a small circle around this pole leads (with P = P 8 ) to We recognise again the elliptic curve of sector 79 and 127. Our last example is the maximal cut in the sector 123. Here we find with P = P 9 + 2m 2 This does not involve an elliptic curve and corresponds to genus zero.

The elliptic curve associated to sector 73
From eq. (93) we may read off the elliptic curve for the sunrise integral: The roots of the quartic polynomial are z (a) This curve has the j-invariant Let us now consider a path γ β in (s,t)-space parametrised by For the Wronskian and the Picard-Fuchs operator d Eq. (86) reduces to the trivial equation We have For a path γ α in (x, y)-space we may use eq. (106) to express λ as a power series in q (a) and vice versa. The point (x, y) = (0, 1) corresponds to τ (a) = i∞.
5.4 The elliptic curve associated to sector 79, sector 93 and sector 127 From eq. (94) we obtain the elliptic curve associated to the double box integral: The roots of the quartic polynomial are now The j-invariant is given by For the Wronskian and the Picard-Fuchs operator d 0,γ β we find for the path γ β defined in eq. (103) Eq. (86) yields In (x, y)-space this translates to We have with For a path γ α in (x, y)-space we may use eq. (117) to express λ as a power series in q (b) and vice versa. The point (x, y) = (0, 1) corresponds to τ (b) = i∞. For y = 1 we have

The elliptic curve associated to sector 121
From eq. (97) we obtain the elliptic curve associated to sector 121: The roots of the quartic polynomial are now The j-invariant is given by For the Wronskian and the Picard-Fuchs operator d Eq. (86) yields In (x, y)-space this translates to We have with For a path γ α in (x, y)-space we may use eq. (127) to express λ as a power series in q (c) and vice versa. The point (x, y) = (0, 1) corresponds to τ (c) = i∞.
For y = 1 we have

Modular forms
The integrals which only depend on t, but not on s, are all related to the elliptic curve E (a) . Associated to the curve E (a) are modular forms of Γ 1 (6). The differential one-forms relevant to the integrals dependent on t but not on s are of the form where τ (a) 6 which we substitute for y (or t). Furthermore, f is a modular form of Γ 1 (6) from the set The modular weights are given by 0, 2, 3, 4 and 2, respectively. The non-trivial modular forms are given by In appendix C we collected useful information on the occurring modular forms of Γ 1 (6) and give alternative representations.

The high-energy limit
In the high-energy limit s → ∞ (or equivalently x = 0) the elliptic curves E (b) and E (c) degenerate to the elliptic curve E (a) . In this limit we therefore have only one elliptic curve E (a) . In this limit we may express all master integrals in terms of iterated integrals of modular forms. The resulting set of modular forms is slightly larger than eq. (132). In order to present this set, we first set Relevant to the high-energy limit is the set These are again modular forms of Γ 1 (6) in the variable τ (a) 6 . Additional details are given in appendix C. All integration kernels reduce in the high-energy limit to Q-linear combinations of elements of this set (times (2πi) dτ (a) 6 ).

Master integrals
In this section we define the 44 master integrals in the basis J. They are labelled J 1 − J 30 and J 32 − J 45 . An integral J 31 is missing due to the extra relation given in eq. (196). In the following we set µ = m and D = 4 − 2ε.
The basis J is constructed as follows: The master integrals, which only depend on s do not pose any problems and can in principal be constructed with the algorithms of [66][67][68][69][70][71]. The master integrals, which only depend on t are similar to the kite / sunrise system and can be constructed along the lines of [21,26]. This leaves the master integrals, which depend on s and t. Here we first analyse the diagonal blocks. Combining the information from the maximal cuts with the technique based on the factorisation properties of Picard-Fuchs operators [62] we first obtain an intermediate basis with the desired properties up to sub-topologies. We then use a slightly modified version of the algorithm of Meyer [69,70] for the non-diagonal blocks and fix the sub-topologies.
As usual, A (1) is block-triangular. The system of differential equations simplifies for t = m 2 (corresponding to y = 1) as well as for s = ∞ (corresponding to x = 0). In both limits the matrix A (0) vanishes. In the former case (t = m 2 ) the integration kernels are linear combinations of the one-forms given in eq. (59) and the solution for the master integrals can be expressed in terms of multiple polylogarithms. In the latter case (s = ∞) the integration kernels are of the form where f is a modular form of the congruence subgroup Γ 1 (6) from the set given in eq. (135). In this case the solution for the master integrals can be expressed in terms of iterated integrals of modular forms. For all modular forms from the set given in eq. (135) the modular weight can be inferred from the scaling behaviour under a rescaling of the periods. One has modular weight = scaling power + 2, where the additional 2 is due to the Jacobian obtained by replacing dy by dτ (a)

. We may view the entries of
as differential one-forms rational in We observe that each entry of A is homogeneous under a simultaneous rescaling of all periods and their derivatives ψ (r) This allows us to group the entries of A according to the scaling behaviour under a simultaneous rescaling of all periods and their derivatives. We define a m-weight as m-weight = scaling power + 2.
This is an ad hoc definition, which we find useful for bookkeeping. No further properties are implied. In the limit x = 0 the m-weight agrees with the modular weight. Let us note that the requirements that

Integration kernels
Let us now discuss the integration kernels appearing in the matrix A. For our choice of basis J we find 107 Q-independent integration kernels. It is a matter of personal taste, if one considers this number to be large or small. On the one hand, it can be considered a large number as it limits our possibilities to present explicit results on paper: For iterated integrals of depth 4 we face a priori 107 4 combinations. On the other hand, it can be considered a small number, given the fact that we are dealing with a matrix of size 44 × 44. The matrix A can be considered to be sparse. Let us mention explicitly that the number 107 is the number for our choice of master integrals J. Other choices of master integrals (respecting the criteria given at the beginning of this section) may lead to a different number of Q-independent integration kernels.
We group the integration kernels according to their m-weight. We have integration kernels with m-weight 0, 1, 2, 3 and 4. The complexity of the expressions increases with the m-weight. The ε 0 -part A (0) contains only integration kernels of m-weight 3 and 4.
Our naming scheme is as follows: We keep the notation for the integration kernels, which already occurred in the special cases t = m 2 or s = ∞. This concerns Integration kernels appearing in the ε 0 -part A (0) are denoted by where n gives the m-weight, (r) indicates the periods appearing in the integration kernel and j indexes different integration kernels with the same n and (r). We denote integration kernels appearing in the ε 1 -part A (1) generically by The superscript (r) and the second subscript j are optional. The interpretation of the super-and subscripts is as above. For dlog-forms we use the notation These are necessarily of m-weight 2.
Let us now discuss for all m-weights typical examples, the cases of m-weight 0 and 1 are discussed completely. The full list of integration kernels is given in the supplementary electronic file attached to this article. The full list consists of the integration kernels 3,1−3 , η with r, s ∈ {a, b, c} and r = s.

m-weight 0
At m-weight 0 we have three integration kernels. They are given by These are exactly the integration kernels we expect at m-weight 0. We expressed them (compactly) in terms of the variables τ Associated to the elliptic curve E (c) are η (c) Note that in the last expression the square root is rationalised by using the variablex. An alternative form for η which is manifestly rational inx, y and ψ (c) 1 .

m-weight 2
The integration kernels of m-weight 2 are numerous and we only list a few typical cases. The integration kernels for the multiple polylogarithms defined in eq. (59) belong to this class. Furthermore, the modular forms of modular weight 2 clearly belong to this class: The differential one-forms in eq. (155) and eq. (156) are all dlog-forms, depending either on x (or alternatively onx) or y, but not both. There are further dlog-forms, depending on both variables x and y. These are There are six differential one-forms involving ratios of periods, one for each ratio. For example dy.
In addition there are 12 differential one-forms of m-weight 2, which do not belong to any class discussed up to now. An example is given by For our choice of master integrals J we observe that in the integration kernels of m-weight 2 polynomials in denominator occur only as a single power, i.e. there are no higher poles in mweight 2.

m-weight 3
At m-weight 3 we have first of all the modular form of weight 3 from the sunrise sector At m-weight 3 we have integration kernels appearing in the ε 0 -part A (0) , an example is given by In addition, there are integration kernels of m-weight 3 in the ε 1 -part A (1) , an example is given by

m-weight 4
At m-weight 4 we have one modular form of weight 4 from the sunrise sector In addition, we encounter integration kernels appearing in the ε 0 -part A (0) . An example is given by Finally, there are integration kernels of m-weight 4 appearing in the ε 1 -part A (1) . An example is given by

Singularities
As already mentioned, the integration kernels are rational in In the next section we will choose as boundary point the point (x, y) = (0, 1) (or equivalently (s,t) = (∞, m 2 )). This motivates the introduction of the variablẽ Our boundary point is then (x,ỹ) = (0, 0). Of particular interest are the polynomials inx andỹ appearing in the denominator of the integration kernels. There aren't too many. Polynomials, which only depend onx are (compare with eq. (59)) Polynomials, which only depend onỹ are Polynomials, which depend onx andỹ are It is helpful to have corresponding expressions in (s,t)-space: The polynomials Q 10 and Q 11 appear when the expression st + (m 2 − t) 2 is expressed in the variablesx andỹ: .
The polynomial Q 12 is related to .
The polynomials Q 13 and Q 14 are related to the elliptic curve E (b) . We have The polynomial Q 13 enters through eq. (114) or eq. (115), the polynomial Q 14 appears in the Picard-Fuchs operator for ψ The polynomial Q 15 enters through eq. (124) or eq. (125), the polynomials Q 16 and Q 17 appears in the Picard-Fuchs operator for ψ (c) 1 in eq. (122). Note that the polynomials Q 1 , Q 7 , Q 15 and Q 18 vanish for (x,ỹ) = (0, 0).

Boundary conditions and boundary constants
We integrate the system of differential equations starting from the point (x, y) = (0, 1) (corresponding to s = ∞ and t = m 2 ). In order to do so, we need the boundary constants at this point. The boundary constants may be expressed as a linear combination of transcendental constants. A basis of these transcendental constants up to weight four is given by [113] w = 1 : ln (2), (2), ζ 2 ln 2 (2), ln 4 (2).
The boundary constants for the master integrals, which are products of one-loop integrals are easily computed. For example, the master integral J 1 is a product of two tadpole integrals. The tadpole integral is given by For D = 2 − 2ε, µ = m and ν = 1 we have In addition we need to calculate explicitly the boundary constants for the integrals, which neither depend on s nor on t. There are two such integrals: J 1 (which is also a product of tadpoles) and J 8 (the sunrise integral at the pseudo-threshold). For the latter we proceed as follows: We start from the Feynman parameter integral. We have (180) At each order in ε, the x 4 -integration is easily performed, resulting in multiple polylogarithms G(z 1 , ..., z k ; x 2 ), where the remaining variable x 2 appears in the argument list z 1 , ..., z k . With the methods of [97,114] we convert all polylogarithms to a form, where the parameters z 1 , ..., z k do not depend on x 2 , the simplest example is given by We may then perform the integration over the variable x 2 . The resulting expressions may be simplified with the help of the PSLQ-algorithm [115].
For all other master integrals we obtain the boundary constants from the behaviour at a specific point, where the master integral vanishes or reduces to simpler integrals. The specific points which we consider are (x, y) = (0, 1), (x, y) = (1, 1) and (x, y) = (−1, 1). For the points (x, y) = (1, 1) or (x, y) = (−1, 1) we integrate the system along y = 1 from x = 0 to x = ±1. For y = 1 we only obtain multiple polylogarithms. We evaluate the multiple polylogarithms to high precision [97] and use the PSLQ-algorithm to extract the transcendental constants.
The complete list of boundary constants is given in appendix D.

Integrals, which are expressed in the variablex
The integrals J 22 and J 36 are expressed in terms of multiple polylogarithms in the variablex. We use the notation of eq. (47). The differential one-forms ω 0 , ω 0,4 and ω −4,0 are defined in eq. (59).

Numerical checks
All results have been verified numerically with the help of the program sector_decomposition [116]. The program sector_decomposition allows (as SecDec [117][118][119][120] or FIESTA [121,122]) the numerical evaluation of multi-loop integrals. On the one hand we evaluated all master integrals of the basis I at a few kinematic points numerically with the program sector_decomposition. On the other hand, we evaluated our results in the basis J at the same kinematic points, converted to the basis I and compared the two results. We find good agreement.
The evaluation of the iterated integrals appearing in our results is done as follows: We split the integration path into two pieces: First we integrate in (x,ỹ)-space from (0, 0) to (x, 0), then from (x, 0) to (x,ỹ). The integration along the first part gives only multiple polylogarithms, which can be evaluated to high precision [97]. We use these results as new boundary constants for the integration along the second part. Assuming thatỹ is small, we may expand for the integration along the second part all integration kernels inỹ.
As a reference we give numerical results for the master integrals in the basis J at the kinematic point

Conclusions
In this article we gave a detailed account on the analytic calculation of the master integrals for the planar double box integral relevant to top-pair production with a closed top loop. The planar double box integral depends on two scales and involves several elliptic sub-sectors. We showed that the associated elliptic curves can be extracted from the maximal cuts. We demonstrated that the system involves three inequivalent elliptic curves. To the best of our knowledge, this is the first time an integral involving more than one elliptic curve has been calculated. Elliptic polylogarithms are by definition iterated integrals on a single elliptic curve. Since the system discussed in this article involves three elliptic curves, we do not expect that our results are naturally expressible in terms of elliptic polylogarithms. We showed that the system of differential equations can be transformed to a form linear in ε, where the ε 0 -term is strictly lower-triangular. This system of differential equations is easily solved to any desired order in ε. We expressed our results in terms of iterated integrals and discussed the occurring integration kernels. We believe that the techniques used in this paper are applicable to a wider class of Feynman integrals. for the n-th root of unity. Let us first introduce a set of modular forms for the congruence subgroup Γ 1 (6), which was already encountered in the calculation of the sunrise / kite system [21,26]. We consider the set ψ 1 π , f 1 , f 2 , f 3 , g 2,0 , g 2,1 , g 2,9 .
and ELi n;m (x; y; q) = ∞ ∑ j=1 ∞ ∑ k=1 x j j n y k k m q jk .