Elliptic polylogarithms and Feynman parameter integrals

In this paper we study the calculation of multiloop Feynman integrals that cannot be expressed in terms of multiple polylogarithms. We show in detail how certain types of two- and three-point functions at two loops, which appear in the calculation of higher order corrections in QED, QCD and in the electroweak theory (EW), can naturally be expressed in terms of a recently introduced elliptic generalisation of multiple polylogarithms by direct integration over their Feynman parameter representation. Moreover, we show that in all examples that we considered a basis of pure Feynman integrals can be found.


Introduction
Feynman integrals constitute the building blocks for the study of scattering processes in perturbative quantum field theory (QFT). High precision calculations in QFT require the ability to compute increasingly more complicated Feynman integrals which involve many internal loops and external legs. The success of the collider physics program, highlighted in the last years by the impressive results obtained by the LHC at CERN, has pushed the required precision of theoretical computations to an unprecedented level. In order to keep up with the experimental demands, theoretical predictions of processes involving Feynman integrals with internal masses and with at least two loops and up to five external legs have become mandatory. In spite of the extreme complexity of these calculations, the last two decades have witnessed an impressive advancement in our understanding of perturbative QFT and, as a result, of our ability to keep such calculations under control.
Typically, complicated Feynman integrals are computed by means of two seemingly orthogonal methods. On the one hand, one can attempt their direct integration over some integral representation (for example in terms of Feynman parameters or Mellin-Barnes integrals). On the other hand, one can derive differential equations (DE) satisfied by the Feynman integrals and try to solve them [1][2][3][4]. Understanding the importance of multiple polylogarithms (MPLs) [5][6][7] in high-energy physics [8,9], and the study of their analytical, algebraic [10][11][12] and numerical [13] properties, have been crucial steps to systematise both strategies. For what concerns direct integration techniques, this program culminated in the enunciation of the criterion of linear reducibility [14,15], which allows one to define a (quite general) class of Feynman integrals that can be algorithmically expressed in terms of MPLs by direct integration over their Feynman-parameter representation. A similarly important result in the context of the differential-equation method (even if mathematically less well established) is the concept of canonical basis of master integrals [16]. Canonical master integrals fulfil differential equations which admit solutions in terms of iterated integrals over particularly simple kernels, which in turn can be expressed as total differentials of logarithms. If there exists a parametrisation of the external kinematics in terms of which the arguments of the logarithms are all rational functions, the corresponding iterated integrals can be straightforwardly expressed in terms of MPLs. Indeed, both by direct integration through the linear reducibility criterion, and in the case of canonical differential equations 1 , one ends up with iterated integrals in a set of variables over rational functions in these variables. It is then a well-known fact that the space defined by these integrals is spanned by linear combinations of MPLs and rational functions.
Despite their applicability to large classes of problems in high-energy physics, already at the second loop order MPLs are known not to exhaust the whole space of special functions required for the computation of Feynman integrals. Indeed, as early as 1962 A. Sabry, in an attempt to compute the two-loop corrections to the electron propagator in QED, encountered integrals of complicated algebraic functions which could not be evaluated in terms of polylogarithms [17], but instead required the introduction of elliptic integrals and integrals thereof. It was not until the second decade of the twenty-first century that such integrals came back to the centre of investigation in particle physics, when it was realised that similar mathematical objects were required for the computation of multiloop corrections to processes of crucial importance to the physics programme at the LHC, like the production of tt pairs in NNLO QCD. Since then, the development of techniques to treat integrals beyond MPLs has been a very pressing issue, both for their potential phenomenological impact in collider physics, and also for their conceptual relevance . Thanks to this concerted effort, in the last years significant steps have been taken in extending both strategies (i.e. direct integration and differential equations) to the so-called elliptic case, namely when the natural geometry associated to the Feynman graphs under consideration is related to a family of Riemann surfaces of genus one.
The scope of this paper is to show how the algorithms for the direct integration of Feynman integrals in terms of MPLs can be suitably generalised to the elliptic case by exploiting the properties of the elliptic polylogarithms (eMPLs) defined in refs. [40,52,53]. We recall here that these eMPLs are essentially equivalent to the multiple elliptic polylogarithms defined by Brown and Levin in ref. [25]. In particular, by working out different examples explicitly, we demonstrate how to treat those classes of Feynman integrals which do not fulfil the criterion of linear reducibility and, instead, require dealing with iterated integrals over more general rational functions R(x, y), where y = P (x) defines an elliptic curve 2 . We stress that, in order for our approach to be successful, one must deal with integrals where only one single elliptic curve appears and no other square roots are present. In this sense, this paper constitutes a concrete step towards the generalisation of the machinery developed to integrate Feynman integrals in terms of MPLs to the elliptic case.
As an important by-product of our calculations, we show that in all examples that we have considered, a basis of pure master integrals can be defined, following the definition provided in ref. [52]. The existence of a pure basis of master integrals is conjectured in the polylogarithmic case and the results presented in this paper constitute non-trivial evidence of a possible generalisation of the current conjectures to Feynman integrals beyond MPLs.
It is possible to draw a parallel between our approach and the idea of "elliptic linear reducibility" recently put forward in ref. [43]. There, a generalisation of the criterion of linear reducibility to the elliptic case is attempted. Our approach is different since we work in the framework of a well defined class of functions whose algebraic and analytic properties can be studied rigorously and for which a concept of transcendental weight can be defined.
The possibility of evaluating Feynman integrals in terms of a well-known class of functions with understood algebraic properties is not only important for computational reasons, but can also be of conceptual relevance. Indeed, the notion of transcendental weight, which is an integer number associated to a pure function, can be used as an organisational tool to classify different expressions. For specific theories, such as N = 4 super Yang-Mills, scattering amplitudes at every loop order L are believed to be of strictly maximal weight 2L. This property is obeyed by all known examples and can significantly reduce the space of functions needed to represent an amplitude. This fact has been heavily explored by the amplitudes bootstrap community in refs. [54][55][56][57][58][59][60][61][62] in order to obtain results up to four loops and seven external legs or five loops and six external legs. Although every amplitude in N = 4 super Yang-Mills which evaluates to MPLs is indeed of uniform maximal weight, it is known that more complicated functions (of the elliptic kind and beyond) are inevitable also in this theory starting already at two loops [63,64]. Therefore, extending the notion of transcendental weight to functions beyond MPLs is a crucial step in order to test the conjecture that observables in N = 4 super Yang-Mills evaluate to functions of uniform weight.
Before diving into the computations, one more comment is in order. It is very clear to us that, even at two loops, the class of functions that we are considering will probably not be the end of the story, since either multiple elliptic curves [45,46], or entirely new geometrical objects [38,[64][65][66][67] can appear. Still, with this paper we aim to show that our framework is general and flexible enough to cover many problems of direct physical interest and, therefore, deserves to be developed further. We will show explicitly how different twoand three-point functions at two loops can be integrated in terms of eMPLs, discussing the details of the manipulations required to bring the integrals to the correct form.
The paper is organised as follows. We begin in Section 2 with a review of eMPLs and their properties, in particular how to assign them a concept of (uniform) transcendental weight, in view of their usage in the next sections. We then move to explicit applications. In Section 3 we consider a family of two-loop non-planar three-point Feynman integrals, whose calculation is relevant for tt and γγ production at the LHC. This family contains two elliptic master integrals, which we express in terms of eMPLs by direct integration over the Feynman parameters. In Section 4 we consider a similar family of two-loop three-point functions, relevant for the computation of the electroweak form factor. The latter contains three elliptic Feynman integrals, which we also explicitly integrate in terms of eMPLs. In Section 5 we show that the same ideas can be applied also for the two-loop kite integral with different internal masses. Finally, in Section 6 we draw our conclusions.

Review of Elliptic Polylogarithms
Our goal in this paper is to show explicitly how the notion of elliptic multiple polylogarithms (eMPLs) developed in refs. [25,40,44,52] can be put into action for a wide range of Feynman integrals known not to be expressible in terms of ordinary MPLs. Before presenting the inner workings of this framework in specific examples, in this section we review the necessary concepts in the context of both ordinary and elliptic MPLs, and in particular the eMPLs introduced recently in ref. [52]. The literature on eMPLs is vast, so here we content ourselves with summarising only the most important aspects necessary for the present calculations and refer the interested reader to refs. [40,52] for more detailed discussions.
Multiple polylogarithms are multi-valued functions defined recursively as iterated integrals over kernels which are rational functions with at most simple poles. The most well-known examples are the classic polylogarithms Li n (x), of which the logarithm is a special case, General MPLs are functions of many variables a i denoting the poles of the rational integration kernels, as well as the endpoint of the integration contour, G(a 1 , . . . , a n ; x) = x 0 dt t − a 1 G(a 2 , . . . , a n ; t) , G(; x) = 1 .
They satisfy properties such as homotopy invariance (they do not depend on the details of the integration path and as such are functions only of its endpoint x) and shuffle relations, where Σ(k, l) stands for all order-preserving permutations of {a 1 , . . . , a k }∪{a k+1 , . . . , a k+l }, called shuffles. It is possible to assign notions of length and weight to MPLs. The length of an iterated integral (polylogarithmic or not) is always defined as the number of integrations, thus the length of an MPL G(a 1 , . . . , a n ; x) is n. The notion of weight, however, is more subtle. For MPLs, the weight is the same as the length, but as will become clear once we discuss its elliptic version, this is not the general case. One can also assign a notion of weight for constants which correspond to MPLs evaluated at special arguments. While a constant has length zero (there are no integrals left to perform; see ref. [52] for a detailed discussion), the weight remembers that of the iterated integral it originated from. For example, log(−1) = iπ → Weight(iπ) = 1 , Length(iπ) = 0 , ζ n = Li n (1) → Weight(ζ n ) = n , Length(ζ n ) = 0 .
(2.4) Upon total differentiation, MPLs undergo a length drop, and their differential takes a particularly simple form, dG(a 1 , . . . , a n ; z) = n i=1 G(a 1 , . . . ,â i , . . . , a n ; z) d log where we defined a 0 ≡ 0 and a n+1 ≡ z. Functions whose total differential does not contain any homogeneous term are referred to as unipotent, and this concept will become important in the following discussions. Elliptic generalisations of MPLs are functions which behave like MPLs but accommodate (in addition to the kernels 1/(x − a)) functions which are rational in the variables x and y which define an elliptic curve, i.e. [x, y, 1] ∈ CP 2 where x and y satisfy a polynomial equation y 2 = P n (x) of degree n = 3, 4. For our purposes, we consider only the case with n = 4 since the n = 3 case can be seen a gauge-fixed version of the former and the examples we consider arise naturally as square roots of degree-four polynomials. Therefore, we are interested in iterated integrals of rational functions in the variables (x, y) subject to the constraint The elements of the vector a ≡ (a 1 , a 2 , a 3 , a 4 ) are referred to as the branch points of the elliptic curve. The periods and quasi-periods of the elliptic curve are chosen according to (2.7) and K and E denote the complete elliptic integrals of the first and second kind, respectively, (2.10) The function Φ 4 (x, a) entering the integrand of the quasi-periods is defined as where s n ≡ s n ( a) denotes the n th elementary symmetric polynomial in the branch points. The periods and quasi-periods are not independent and satisfy the Legendre relation, Since an elliptic curve y is given in terms of a square root, it is important to make a choice for the signs of the branches of the square root y which is consistent with the conventions for the periods and quasi-periods in eqs. (2.7) and (2.8). In particular, we find it convenient to order the branch points such that ω 1 ∈ R and ω 2 ∈ iR, whenever possible. Note that this implies that λ defined in eq. (2.9) lies between 0 and 1. Moreover, in order to correctly define eMPLs we need to provide a prescription for how to perform the integrals in the regions of interest, namely on the real line and between branch points. In this paper, we deal with Feynman-parameter integrals whose integrands are ratios of polynomials defined over the real numbers. Therefore, depending on the kinematic regime we wish to consider, the branch points can be real or complex so long as the polynomial y 2 = P 4 (x) is real. There are only three possible configurations for the branch points such that y 2 is real, which we consider in turn below (see fig. 1): (i) All branch points are real.
In this situation, we only need to consider integrations over the real axis. We order the branch points according to a 1 < a 2 < a 3 < a 4 and fix the signs of the branches of the square root as (2.13) (ii) Two branch points are real and two are complex conjugate to each other. The configuration of the branch points that feature in our applications is such that one can always impose the following ordering a = (a 1 , a 2 , a 3 , a 4 ) where a 1 = a * 4 ∈ C , Im(a 1 ) > 0 , a 2 , a 3 ∈ R , a 2 > a 3 , (2.14) where the branch points satisfy Re(a 1 ) = Re(a 4 ) = 1/2 . This choice seems particularly ad hoc, but we can motivate it as follows. In this situation, the polynomial P 4 (x) is negative on the real axis for a 3 < x < a 2 and for x = 1/2 + i t, Im(a 4 ) < t < Im(a 1 ). We recall here that, in order to define the periods we need a prescription for how to compute integrals between different branch points.
In this case, in addition to a prescription to integrate along the real axis, it is enough to supplement it with a prescription on the line x = 1/2 + i t for Im(a 4 ) ≤ t ≤ Im(a 1 ), see fig. 1(ii). With this in mind, we define the elliptic curve in all regions of interest as (iii) All branch points are complex and pairwise complex conjugate.
(2.18) eMPLs were originally defined in refs. [25,68,69] as iterated integrals on a complex torus. This description is related to the one presented here through the relation between elliptic curves defined by the equation y 2 = P 4 (x) and a torus defined as the complex plane quotiented by a two-dimmensional lattice Λ = Z ω 1 + Z ω 2 . The ratio of the two lattice periods τ = ω 2 /ω 1 is called the modular parameter and it is easy to see that the lattice Λ remains invariant under modular transformations SL(2, Z) mixing the two periods. Modular transformations act on τ as Möbius transformations. The way to map an elliptic curve to its equivalent torus description is via the function [40], Here z is a variable on the torus, ℘(z) is the Weierstrass ℘-function,s n ≡ s n (a 2 , a 3 , a 4 ) and the s n are the symmetric polynomials defined below eq. (2.11). The κ-function satisfies a differential equation which is identical to the definition of the elliptic curve in eq. (2.6), namely (c 4 κ (z)) 2 = P 4 (κ(z)) , (2.20) and thus one may identify (x, y) ↔ (κ(z), c 4 κ (z)). The inverse of the κ-function is known as Abel's map, which takes a point (x, y) on the elliptic curve to a point z x on the complex torus, Since elliptic curves are isomorphic to complex tori, eMPLs can be described as iterated integrals over functions related to the torus, and were originally defined as such in refs. [25,68,69]. In this context, eMPLs are defined as iterated integrals given by where the integration kernels are the coefficients in the expansion of the Kronecker-Eisenstein series F (z, α, τ ), and θ 1 (z, τ ) is the odd Jabobi theta function with θ 1 (z, τ ) denoting a derivative with respect to its first argument. The eMPLs (2.22) behave similarly to ordinary MPLs in that they also form a shuffle algebra and are unipotent. Moreover, they are pure according to the definition of ref. [52], namely: A function is called pure if it is unipotent and its total differential involves only pure functions and one-forms with at most logarithmic singularities.
In the calculation of Feynman integrals that evaluate to functions of the elliptic kind, the representation of elliptic polylogarithms in terms of a polynomial equation y 2 = P (x) appears more naturally than the torus picture. Therefore in this paper we use the definition of pure eMPLs on the elliptic curve recently put forward in ref. [52]. They are defined as iterated integrals of kernels that are rational functions on the elliptic curve with at most logarithmic singularities in all variables, In contrast with MPLs, for the elliptic case the requirement that all integrations over rational functions on the elliptic curve close on the same space of functions put together with the requirement that all integrals must have at most logarithmic singularities leads to an infinite tower of independent kernels Ψ n for n ∈ Z. This fact can also be seen from the torus description, where an infinite number of kernels are generated by eq. (2.23).
In particular, the kernels in eq. (2.24) depend on a certain kind of functions which are themselves transcendental, namely and Φ 4 (x, a) given in eq. (2.11). The kernels Ψ n entering the eMPLs in eq. (2.24) are spelled out below for |n| = 0, 1, 2. Higher values of n do not appear in the present applications since the corresponding functions would satisfy higher-order differential equations. Before writing down the expressions for the kernels, we introduce some functions which appear as ingredients. The first is the function Z 4 (x, a) defined in eq. (2.25). Likewise, an important element is the image of the point x = −∞ under Abel's map (2.21), It is possible to represent z * in terms of elliptic integrals. In the situation where all roots are real and ordered, it is given by [52], and for other configuration of the branch points (or complex ones), z * may pick up a minus sign depending on the conventions for the branches of the square root. Finally, the kernels entering the pure eMPLs depend on the function G * ( a), which is simply the image of the point z * under g (1) (see eq. (2.23)), As shown in ref. [52], G * can be integrated explicitly in terms of (incomplete) elliptic integrals of the first and second kind. In the situation where the branch points a are real and ordered according to a 1 < a 2 < a 3 < a 4 one finds In the special case where the point z * is of the form 3 for a and b constants, then G * ( a) admits an even simpler form, namely This follows because eq. (2.27) together with eq. (2.30) imply that α = α(λ). At last, we are now ready to write down the expressions for the kernels. For n = 0, there is only one kernel, For n = 1, we have instead four kernels (with c = ∞) where y c ≡ P 4 (c). Finally, for n = 2, we have (with c = ∞), where Z 4 (x, a) stands for a degree-two polynomial in Z 4 (x, a), (2.35) We conclude this short exposition of pure eMPLs with a comment: much like with ordinary MPLs, one can associate a concept of length and of weight to eMPLs and to quantities which arise from evaluating eMPLs at special points, for example the periods and quasi-periods of the elliptic curve defined in eqs. (2.7) and (2.8). In summary we have the values shown in Table 1.
Using the formalism revised in this section, in the rest of this paper we will show how certain Feynman integrals which evaluate to functions beyond MPLs can be brought to neat expressions in terms of combinations of pure eMPLs (2.24) of uniform weight by direct integration of their Feynman parametrisation. Table 1. Length and weight of eMPLs and related constants.

A non-planar triangle with a massive loop
In this first application, we consider the family of two-loop non-planar three-point functions with a massive loop shown in fig. 2. These integrals involve massless propagators and a massive loop with four propagators with mass m. Two external legs are massless, i.e. p 2 1 = p 2 2 = 0 and we set our kinematic variable q 2 = (p 1 + p 2 ) 2 . The family of two-loop integrals is then where the a i are the exponents of the propagators, and we consider only integrals with a 7 < 0. The propagators are The family of integrals in eq. (3.1) was studied in ref. [37] by means of the differentialequation method, where it was shown that there are two master integrals for the top topology which satisfy a coupled two-dimensional system. 4 We choose as basis integrals where we note that both integrals are finite in D = 4 space-time dimensions. In ref. [37], the solution of the system was presented as an expansion in the dimensional-regularisation parameter = (4 − D)/2 with coefficients being integrals over elliptic integrals of the first and second kind multiplied by ordinary MPLs and rational functions.
In this paper, we apply the framework of eMPLs on the elliptic curve reviewed in Section 2 in order to obtain analytic expressions for these Feynman integrals and study their properties. In particular, we show that once the two master integrals in the top sector are expressed in terms of our class of functions, one can easily transform them to a new basis of master integrals which are explicitly pure, similarly to what one would do if they were standard MPLs. To this end, we approach the problem through direct integration over the Feynman parameter representation of the integrals and show that all integrations up to the penultimate can be performed in terms of ordinary MPLs with algebraic arguments that can in general depend on an elliptic curve. In order to perform the last integration, we employ partial fractioning to rewrite the integrand as kernels of the eMPLs defined in eq. (2.24) (see eqs. (2.32), (2.33) and (2.34)). We now look in detail at the computation of the two master integrals.

First master integral
We start by considering the first master integral M 1 , defined in eq. (3.3). Our task is to find an order of integration of the Feynman parameters such that linear reducibility is achieved for all integrations except the last one, which will in turn require the introduction of eMPLs. The Feynman parameter representation of the first master integral is given by where U and F are the Symanzik polynomials associated with the graph. Since this integral is finite, it can be evaluated directly in D = 4, where the U-polynomial drops out. The next step is to apply the Cheng-Wu theorem [71] in order to find a particular parametrisation of the integral such that as many Feynman parameter integrations as possible can be done in terms of ordinary MPLs. In practice, the Cheng-Wu theorem allows one to exchange the integration domain ∆ of a Feynman integral by using any subset of propagators Σ such that Here we choose to apply the theorem with Σ = {1, 2, 3, 4}. This particular ordering was used first in ref. [43] to write a one-fold integral representation for a similar integral.
In doing so, the integral (3.1) becomes where the polynomial F with and we have factored out (−q 2 ) and encoded the kinematic dependence in the dimensionless variable .
From now on we set m 2 = 1 for simplicity and restore its dependence at the end using dimensional analysis. The integrals over x 6 , x 5 and x 3 can be done easily in terms of MPLs (2.2) as the integrand is linearly reducible in these variables. This results in a two-fold integral, The next integration to be done is over x 2 . We notice immediately that the overall rational pre-factor is quadratic both in x 2 and x 4 . Indeed, as we will see below, upon integration in either variable this will give rise to a square root of a polynomial of degree four in the other one, defining an elliptic curve. On top of this, performing a study of the symbol alphabet of the combination of MPLs in eq. (3.9), we find the following letters, (3.10) The alphabet above involves quadratic letters both in x 2 and x 4 , such that by rewriting the MPLs in the form G(. . . , x 2 ), one would in general be left with MPLs involving multiple square roots involving x 4 . The presence of these additional square roots could prevent us from performing the last integration in terms of eMPLs algorithmically. However, it is easy to realise that after a simple change of variables such that if we perform the change of variables and exchange the order of integration in eq. (3.9) as follows, we expect to be able to perform the integral in x 4 without introducing any additional square roots inx 2 . Note that, while this transformation linearises the overall symbol of the integrand, individual MPLs may still involve quadratic symbol letters that cancel out in the combination. Indeed, in this case it turns out that in order to rewrite the individual MPLs one needs to introduce the square-root valued letters These letters enter in identities of the type . (3.14) Once these identities are inserted back into eq. (3.9), the dependence on r ± cancels out, as expected from eq. (3.11). Note however that, for this particular example, the cancellation of r ± at this stage is not required for the integration algorithm to go through, since r ± do not depend on the remaining integration variablex 2 . By performing these manipulations and integrating over x 4 through the recursive definition of MPLs (2.2) we arrive at a one-fold integral in the variablex 2 given by where y is the square root of a quartic polynomial, as anticipated, and as such defines an elliptic curve with branch points The ordering of the branch points in b is chosen following the prescription in eq. (2.14) for the kinematic region where a > 1/16, so that b − = (b + ) * are complex. This choice is made for convenience, as it is simpler to perform a numerical evaluation of the final result when the branch points and potential poles are not on the real axis. In eq. (3.15) we have also introduced the shorthand notation for symmetric and anti-symmetric combinations of MPLs depending on the elliptic curve, with the variables R ± given by In order to compute the integral in eq. (3.15), our next task is to recast its integrand in terms of eMPLs whose dependence on the integration variable is of the form E 4 (. . . ;x 2 ) multiplied by eMPL kernels (see eqs. (2.32), (2.33) and (2.34)). This way, the last integral in the variablex 2 can be performed trivially using the recursive definition of eMPLs (2.24). This is a bottom-up procedure in the length of the eMPLs: at a given length L = n, it can be achieved through a sequence of four steps: This procedure guarantees that all eMPLs are of the form E 4 (. . . ;x 2 ) times an elementary eMPL kernel and thus can be trivially integrated using eq. (2.24). Note that this is a standard procedure for ordinary MPLs and the only difference here is that we need to identify and use the elliptic kernels instead.
To illustrate the mechanism, let us first consider the MPLs that do not depend on the elliptic curve, i.e. G( n; x) where the endpoint x and the letters n are either constants or depend on the square-root letters defined in eq. (3.13), but not on y. In this situation, a derivative inx 2 plus partial fractioning will lead to an integrand which depends only on polylogarithmic kernels of the form 1 x 2 −c , with c independent ofx 2 . As such, it can simply be integrated back to an MPL of the form G( m,x 2 ) with m independent ofx 2 .
In a similar fashion, for the situation where the MPLs depend on the elliptic curve, the procedure will lead to a combination of eMPL kernels which can then be integrated to E 4 functions. As an example, consider the weight-one function Taking a derivative with respect tox 2 yields The next step is to integrate the above expression back in terms of E 4 functions. To this end, we rewrite the integrand in terms of the kernels defined in eqs. (2.32) and (2.33), In the expression above, the terms proportional E 4 0 0 ;x 2 , b cancel out since for the elliptic curve under consideration G * ( b) and c 4 are related and in particular G * ( b) evaluates to a simple algebraic function, with where we omitted the dependence of E 4 on the branch points b for clarity and restored the factors of m 2 using dimensional analysis. Note that, in order to obtain the expression above, we used the fact that for this region the expression for G * ( b) reduces to the simple algebraic function in eq. (3.23). We stress once more that a Feynman integral which can be put in this form (i.e. one single prefactor which depends on the elliptic periods times a combination of elliptic polylogarithms of uniform transcendental weight) appears to be the natural generalisation of a pure integral from the polylogarithmic to the elliptic case.
We remind the reader that the expression in eq. (3.25) is valid in the region where a > 1/16 for which two of the roots of the elliptic curve are located outside of the real axis, and the branch points are ordered according to the conventions outlined in eq. (2.14), as shown in eq. (3.17). It is easy to see that in the Euclidean region 0 < a < 1/16 the roots are real and can be ordered on the real axis between 0 and 1 (see fig. 3 below). Moreover, the special points r +/− defined in eq. (3.13) also become real and, in particular, one finds 0 < r − < 1. Therefore, in this region the pole r − lies on the integration contour and eq. (3.25) is not a convenient representation. As an example of the analytic continuation of an expression written in terms of eMPLs, in Appendix A we show how to analytically continue eq. (3.25) to the region 0 < a < 1/16.

Second master integral
Performing similar steps for the second master integral, M 2 in eq. (3.3), we can express it in terms of the same class of functions. Contrary to the first master, we find that the second master integral cannot be written as one single (transcendental) prefactor times a pure combination eMPLs of uniform weight. Instead, by direct integration we find where Ω (tt) 2 andM 1 is the very same combination of eMPLs defined in eq. (3.25).M 2 is instead another (independent) pure combination of eMPLs very similar to eq. (3.25) given bỹ where T 2+ (a) and T 2− (a) are pure functions of uniform weight four given by We note here that, according to the prescription provided in ref. [52], the prefactors Ω and Ω (tt) 2 have a natural interpretation in terms of the maximal cuts of the integrals M 1 and M 2 . In order to see why this is the case, let us recall that the two master integrals M 1 and M 2 fulfil a system of two coupled differential equations which, neglecting the subtopologies, read [37] A complete solution of eq. (3.31) can be obtained by considering the matrix of the maximal cuts of the two master integrals evaluated along two independent integration contours [32,38,72,73]. The number of independent contours is always equal to the number of master integrals. We denote a basis of independent contour by C j , j = 1, 2. By indicating the maximal cut of the integral M i along the contour C j as Cut j (M i ) we find This in turn implies that if we define a new basis of master integrals as then by construction the new basis fulfils Cut j (F i ) = δ ij . Let us write Cut 1 (M j ) = Ω j and Cut 2 (M j ) = H j , and decompose the matrix G as follows Comparing eq. (3.34) with eq. (3.30), one can explicitly make the identifications and Ω (tt) 2 correspond to the maximal cuts of the two master integrals along the first integration contour. Clearly, the choice of which integration contour is considered to be the first one and which the second ones is arbitrary. This is reflected in the ambiguity of the splitting in eq. (3.34). Using this insight, we can verify that indeed Ω Before summarising the main points of these calculations, it is worth stressing once again that the result that we obtained, in particular in terms of two independent combinations of pure functions eq. (3.30), is a-priori non trivial and it provides a strong hint towards the generalisation of the idea of a pure basis of master integrals to the elliptic case.

Summary
With the explicit computation of the family of non-planar triangles considered in this section, we have discovered an elegant structure underlying these integrals. In particular, we highlight here the main features of the computation and result: • eMPLs provide a natural language to express Feynman integrals which depend on an elliptic curve. Their recursive definition (2.24) provides an algorithmic method for computing each integration step.
In the examples considered in this section, we applied the Cheng-Wu theorem in order to delay the appearance of the square root defining the elliptic curve to the last integration, but we stress that this is only for technical simplicity, the framework can accommodate a dependence on the elliptic curve at any step. The only caveat is that, as mentioned earlier, the framework is suited for cases where only a single elliptic curve is present, and thus in earlier integrations steps one cannot handle elliptic curves that still depend on other integration variables.
• The results of the two master integrals shown in eqs. (3.24) and (3.26) are naturally organised in terms of pure building blocks of uniform weight. These results can easily be rewritten in a new basis of pure master integrals. However, in contrast with the polylogarithmic case, the change of basis is not algebraic and depends on the periods and quasi-periods of the elliptic curve (see eqs. (3.25) and (3.27)).
• Since the Feynman-parameter integrals we started from are simply integrals over rational functions, the final result expressed in terms of eMPLs should reflect the property that the original integral has no intrinsic dependence on which sign we choose for the square root y. As such, the final integral should be invariant under the parity transformation y → −y. Using the fact that ω 1 , η 1 and Z 4 (x, a) are parity odd, it is easy to see from the definitions of the eMPLs kernels in eqs. (2.32), (2.33), (2.34) that under y → −y the kernels transform as With this, one can see that every term in the final results for the two master integrals are manifestly parity-even as desired.
After learning the details of this computation and the properties of the result, we now turn our attention to a different integral which has a very similar structure: a particular contribution to the electroweak form factor.

Electroweak form factor
In this section we consider another two-loop three-point function which is at first sight very similar to the one considered in the previous section. The topology is shown in fig. 4 and differs from that of fig. 2 in that two additional internal propagators are massless. This family of integrals contributes to non-planar two-loop corrections to the electroweak form factor and was studied in ref. [74], where an approximate solution as series expansions around all singular points was obtained. We define the integral family as where as before a 7 < 0 and the remaining a i > 0 are the powers of the propagators below, As for the previous family of integrals, we define q 2 = (p 1 + p 2 ) 2 = −m 2 /a, where a is a dimensionless ratio. In contrast with the integrals discussed in Section 3, in this case the top-sector is reduced to three master integrals which satisfy a coupled system of three differential equations. We choose the three master integrals as follows where we notice that all three master integrals are finite in D = 4 space-time dimensions.
We proceed similarly to the previous example and integrate all three master integrals explicitly starting from their Feynman parameter representation. As before, in order to be able to express all three master integrals in terms eMPLs, we need to find an ordering of Feynman parameters (or more precisely an application of the Cheng-Wu theorem) which allows us to perform all integrations either in terms of standard MPLs or by introducing at most one square root of a quartic polynomial in one of the integration variables. This square root will define the elliptic curve associated to the problem. By direct inspection we find that the same application of the Cheng-Wu theorem as in the previous section does the trick. The polynomial equation defining the elliptic curve resulting from it, though, appears to be substantially more complicated, at least once its roots are seen as functions of the kinematic invariant a defined in eq. (3.8) above. The elliptic curve is defined by where the roots are defined as Moreover, in the intermediate integration steps we find standard MPLs which have branch cuts in the special points x ∈ {a, −a, 1 + a, 1 − a, r − , 1 − r − } with Similarly to the integrals considered in the previous section, since the integration over the last Feynman parameter varies in the interval (0, 1), it is convenient to work in a kinematical region where there are no explicit poles on their integration contour. Since, in order to implement Feynman's prescription, the dimensionless ratio a becomes a complex number with a small negative imaginary part, we choose to work with Re(a) > 1, such that the all poles lie outside of the integration contour and, in particular, r +/− are complex conjugate to each other. In this situation, two of the branch points of the elliptic curve in eq. (4.4) are real and located between 0 and 1, whereas the two remaining branch points are complex conjugate to each other and have real part equal to 1/2. We order the branch points of the elliptic curve according to the conventions of Section 2 such that eqs. (2.14) and (2.15) are satisfied, namely with d ± given in eq. (4.5). In order to arrive at compact expressions for the integrals in terms of pure eMPLs, we use the fact that the functions c 4 ( d), Z 4 (x, d) and G * ( d) admit particularly simple representations in this case, namely While the evaluation of the three master integrals proceeds at least conceptually along the same lines described in detail in Section 3, the individual manipulations and the final results are more cumbersome, mainly due to the explicit form of the branch points of the elliptic curve. The result for the first master integral in terms of pure eMPLs is given by (4.10) The three master integrals in eq. (4.3) corresponding to the top-dimensional topologies of this family of integrals follow a structure similar to that observed in the previous section for the triangle with a massive loop. Indeed, after computing the second and third master integrals also through direct integration of the Feynman-parametric integral, it is possible to transform them into a basis of functions which are pure combinations of eMPLs of uniform weight,    (4.11) The entries of the matrix in eq. (4.11) are given by  Similarly to the triangle considered in the previous section, the Ω (ew) i above satisfy the differential equations for the maximal cuts of the three master integrals (see eq. (3.36)), namely: The pure functionÑ 2 is a weight-four function which we can decompose into a part that depends on the variable r − as well as a part with dependence only on the variable a and a piece which is purely polylogarithmic. It is given bỹ (4.14) Finally, the functionÑ 3 is more complicated than the previous two and we prefer not to write it here explicitly. Its expression is attached to the ancillary files of the arXiv submission of this paper. As for the previous case of the tt triangle, we see that following eq. (4.11) the three elliptic Feynman master integrals can be re-expressed in a new basis of three pure master integrals,Ñ 1 ,Ñ 2 andÑ 3 . This provides further evidence that it is possible to find a pure basis to represent master integrals which do not evaluate to MPLs.

Kite with three distinct masses
Having computed several three-point functions by performing the integrations over Feynman parameters in terms of pure eMPLs, we now consider an example of a two-point function with more scales, namely the kite integral with three distinct masses shown in fig. 5. The integral we consider is given by , (5.1) in D = 2 − 2 . A simpler version of this integral, when all three internal masses have the same value m 1 = m 2 = m 3 = m, has been computed in the literature in terms of iterated integrals over products of elliptic integrals and polylogarithms [31] or modular forms [35], in terms of elliptic generalisations of polylogarithms [35] and finally, more recently, in terms of the eMPLs considered here [52]. We consider here the more general case with three different internal masses. We encode the kinematic dependence is the three dimensionless ratios We compute the kite integral in the region 0 < p 2 < min(m 2 1 , m 2 2 , m 2 3 ). The branch points are complex and given by a = a − , a * − , a + , a * + , As in the previous applications for three-point functions, the kite integral can be computed in terms of a pure combination of eMPLs of uniform weight three. In order to arrive at the final expressions, we make use of the following relations valid for the kinematic region we consider, The result for the kite integral (5.1) with three distinct masses in terms of eMPLs is given by where we define the shorthand notation (5.8) The part of the result which depends only on ordinary MPLs,K MPL (a i ), in turn is given byK  a 1 a 3 , 0, a 2 ) + G(a 1 a 3 , 0, 1, a 2 ) + G(a 1 a 3 , 0, a 3 , a 2 ) − G (a 1 a 3 , a 1 , 1, a 2 ) − G(a 1 , a 2 ) − G(a 3 , a 2 ) + 4G(0, a 1 , a 2 ) + 2G(0, a 1 a 3 , a 2 ) − 4G(η, 0, a 2 ) + 2G(η, a 1 , a 2 ) − 4G(η, a 3 , a 2 ) + 2G(η, a 1 a 3 , a 2 ) + 4G(a 1 , 0, a 2 ) − 2G(a 1 , a 1 a 3 , a 2 ) + 4G(a 3 , 0, a 2 ) − 4G(a 3 , η, a 2 ) − 2G(a 3 , a 1 a 3 , a 2 ) G(0, a 1 )

Conclusions
In this paper we gave a detailed account of the calculation of different classes of twoloop Feynman integrals which evaluate to elliptic polylogarithms. The main goal has been to show that the direct integration algorithms used to deal with polylogarithmic Feynman integrals can be suitably generalised to include a considerably large class of elliptic Feynman integrals. We considered in detail two different two-loop non-planar three-point functions appearing in the production of tt pairs in QCD and in the electroweak form factor, respectively. Subsequently, we computed the famous two-loop kite integral with three different internal masses. As it is well known, the kite integral inherits its elliptic nature from the two-loop massive sunrise subtopology. The two families of three-point functions on the other hand are genuinely elliptic at the level of the top topology and in particular contain elliptic curves that are not directly related to the sunrise case. In spite of their apparent diagrammtic similarity, these two families of three-point functions differ substantially. In fact, in the tt case the elliptic top-sector can be reduced to two independent master integrals (up to simpler polylogarithmic sub-topologies), while in the electroweak form factor case, the elliptic sector is reduced to three independent master integrals (again, up to subtopologies). By direct integration over their Feynman parameter representation, we have shown that the complete set of master integrals for the three problems above can be consistently expressed in terms of elliptic polylogarithms. We have performed all computations in the Euclidean region, such that all our results are real. In doing so, we described how to define our integrals consistently depending on the location in the complex plane of the four branch points that define the elliptic curve. As a crucial result of our calculations, we showed that for all examples considered it was always possible to organise the master integrals of a given topology into a basis of pure building blocks. We stress here once more that by calling it a basis we mean that we can find a number of independent pure building blocks which is equal to the number of master integrals in the problem under consideration. This supports the conjecture according to which a basis of pure functions should always exist for master integrals that can be expressed in terms of MPLs only, and constitutes a strong hint towards its generalisation to the elliptic case.
Finally, there are strong indications that also many two-loop four-point Feynman integrals can be expressed in terms of the same set of functions and that similar considerations on their transcendentality properties apply. We postpone the details of these calculations to future publications.

A An example of analytic continuation
In Section 3.1, we presented a result for the first master integral for the two-loop triangle integral relevant for top production in terms of eMPLs. The result shown in eq. (3.25) was however not well defined in the Euclidean region for which 0 ≤ a ≤ 1 16 as there were poles on the integration contour. In this section, we show how to circumvent this problem by rewriting each E 4 appearing in eq. (3.25) as a well defined combinations of eMPLs without poles on the integration contour. This can be achieved through manipulations similar to how the integrand of eq. (3.15) was converted into eMPLs, i.e. by taking a derivative with respect to the kinematic variables r +/− and integrating back, thus ensuring the endpoint of the contour for every eMPL to be r +/− . In order to fix the imaginary parts of the eMPLs, we fix the imaginary parts of r +/− using a Feynman iε prescription which amounts to giving a small positive imaginary part to q 2 in eq. (3.8). Using eq. (3.13) this gives in the Euclidean region Im(a) > 0 , Im(r − ) > 0 , Im(r + ) < 0 . (A.1) We recall here that for the problem at hand a and r +/− are not independent, see eq. (3.13). Nevertheless, we imagine to treat them as independent variables in what follows.
The expression in eq. (3.25) contains eMPLs of length 3 and 4, which must be recast as eMPLs with endpoints r +/− . Like for ordinary MPLs, this procedure consists of taking a derivative, thus generating lower-length eMPLs and integrating these back in the variable r +/− . Therefore, one needs to work out first the length 1 examples in order to build the higher length iteratively. In the following we show examples of this procedure and in particular how to compute the boundary constants in some intricate cases.

A.1 Length 1 example
As a warm-up, consider the eMPL E 4 1 r − ; 1 which is simply a logarithm, log(1 − 1/r − ). Taking a derivative with respect to r − we find Then, integrating back we find The integration constant c is then fixed by taking the limit r − → 0. In situations like the example above, there are logarithmic divergences in this limit due to the lower integration boundary at zero. As such, it is necessary to extract divergent logarithms in order to compute the boundary constant. For Im(r − ) > 0 we have for the left-hand-side whereas for the right-hand-side we have thus fixing c = iπ.

A.2 Length 2 example
We now consider a more intricate example of length two and explicit dependence on the elliptic curve, namely E 4 0 1 0 r − ; 1 , in order to illustrate the typical steps involved in the analytic continuation of certain eMPLs. Note that due to the presence of the kernel Ψ 1 (r − , x), the integration on the real axis 0 < x < 1 cannot be performed for 0 < r − < 1 because the integrand has a pole on the contour. This is the case for the tt triangle in the Euclidean region. To avoid this problem, our goal is to rewrite E 4 0 1 0 r − ; 1 as a combination of eMPLs of the form E 4 ( ··· ··· ; r − ) such that there are no poles in the new contour 0 < x < r − . Similarly to the previous length-one example, this is achieved by differentiating the expression with respect to r − and integrating back in terms of eMPLs with a different end point. The first step is where we used the explicit expressions for the eMPLs kernels in eqs. (2.32) and (2.33) and we omit the explicit dependence on the branch points of the elliptic curve, Note that the ordering above differs from that of eq. (3.17) as we are now in a configuration with 0 < a < 1/16 where all roots are real and ordered. Taking the derivative in r − explicitly under the integral sign (recall the definition of Z 4 in eq. (2.25)) and integrating in x 1 and x 2 , we find In order to compute c 1 we consider the limit r − → 0. We know that in this limit Z 4 (r − ) → 0 and thus we are left with only the first term in the integrand (A.10), where we show explicit the Feynman iε prescription according to eq. (A.1). Changing variables to t = x/r − and taking the limit r − → 0, the square roots above become very simple, With this we can obtain the result for the integration constant c 1 , (A.14) We are now ready to plug the result E 4 and thus we get our final expression valid for 0 < a < 1/16, The procedure shown for this particular example can be generalised for arbitrary length and used to transform the complete expression into a result valid in the Euclidean region. For illustrative purposes, in the following section we show the first master integral completely transformed using these techniques.

A.3 Result for top-production first master in the Euclidean region
In the last section, we showed how the eMPLs expression for the first master integral for the top-production triangle given in eq. (3.25) can be turned into an expression valid in the Euclidean region. In this section, we show the result in terms of pure eMPLs, which is valid in the Euclidean region 0 < a < 1/16 and is of uniform weight 4. The result is