Elliptic polylogarithms and iterated integrals on elliptic curves I: general formalism

We introduce a class of iterated integrals, defined through a set of linearly independent integration kernels on elliptic curves. As a direct generalisation of multiple polylogarithms, we construct our set of integration kernels ensuring that they have at most simple poles, implying that the iterated integrals have at most logarithmic singularities. We study the properties of our iterated integrals and their relationship to the multiple elliptic polylogarithms from the mathematics literature. On the one hand, we find that our iterated integrals span essentially the same space of functions as the multiple elliptic polylogarithms. On the other, our formulation allows for a more direct use to solve a large variety of problems in high-energy physics. We demonstrate the use of our functions in the evaluation of the Laurent expansion of some hypergeometric functions for values of the indices close to half integers.


Introduction
The discovery of the Higgs boson at the LHC and the absence of clear signs of New Physics close to the electroweak scale provide a strong confirmation of the Standard Model of particles physics (SM) up to the TeV scale. At the same time, the SM fails spectacularly at explaining a wide spectrum of important phenomena, like the overwhelming experimental indications for the existence of Dark Matter, and also at providing a consistent framework to describe all known fundamental interactions, including gravity. Given the increasing precision of the measurements carried out at the LHC, equally precise calculations in the SM become then an essential tool to uncover shortcomings of the SM and discover possible elusive signs of New Physics.
The precision calculation of physical observables in perturbative quantum field theories relies on the computation of Feynman graphs with many external legs and loops. As a result of more than a decade of investigation, a lot has been understood on the mathematical properties of Feynman integrals. These developments, in turn, have led to powerful techniques for the analytical and numerical calculation of loop integrals. In particular, the discovery that a huge class of Feynman graphs can be naturally expressed in terms of multiple polylogarithms [1][2][3][4], has made possible the calculation of higher-order corrections to a large number of processes of crucial importance for the precision physics program carried out at the LHC, but that were previously thought to be out of reach. In spite of having been known to mathematicians for a very long time [5,6], during the last decade multiple polylogarithms have become again an active field of research, as the discovery of their underlying Hopf algebra structure has paved the way towards a thorough understanding of their analytical and numerical properties [7][8][9], in particular allowing to systematize the study of non-trivial functional identities between them.
While multiple polylogarithms have a large range of applicability, they are not the end of the story. Recently, a growing number of integrals appearing in higher-order calculations in the SM, but also in N = 4 super-Yang-Mills theory or open string theory, have been identified, whose analytical calculation requires functions beyond multiple polylogarithms [10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][27][28][29]. Indeed, the appearance of non-polylogarithmic structures in quantum field theory had been noticed already more than fifty years ago in the computation of the two-loop corrections to the electron self-energy in QED [30]. The origin of this new class of functions could be traced back to a particular Feynman graph, the so-called two-loop massive sunrise graph, whose analytical properties have been extensively analyzed from different perspectives in the physics and mathematics literature [13][14][15][16][17][18][19][20][21][31][32][33][34][35]. 1 While many analytical representations for the sunrise integral have been worked out, a way to generalize them to more complicated Feynman integrals remains matter of debate.
From a more mathematical point of view, multiple polylogarithms are obtained by integrating rational functions on a Riemann surface of genus zero -a Riemann sphere. This perspective suggests as natural generalisation of multiple polylogarithms the class of functions obtained by integrating rational functions on Riemann surfaces of higher genus. Indeed, it is by now very well known that the two-loop massive sunrise graph, together with many other examples from quantum field theory and string theory, can be represented by functions on a surface of genus one -a torus or, equivalently, an elliptic curve. Genus-one generalisations of multiple polylogarithms have been studied by mathematicians [37][38][39] and, quite naturally, have been dubbed multiple elliptic polylogarithms. Many important properties of these functions are by now well understood; in particular, a very important result of [39] was to introduce a set of integration kernels defined on the torus with at most simple poles, and to show that they allow to perform the integration of any rational function on the corresponding elliptic curve, e.g. the space of functions is closed under taking primitives.
In spite of these appealing results, the use multiple elliptic polylogarithms in the calculation of physically relevant Feynman integrals has been impeded and complicated by the fact that most of the mathematical literature makes only use of the torus formulation of elliptic curves. While in that formulation the integration kernels used to define elliptic polylogarithms have a very clean geometrical interpretation, their application to solve problems in high energy physics is not straightforward and for practical purposes a different language appears to be preferable. The main goal of this paper is therefore to provide a different formulation of elliptic polylogarithms in terms of simple integration kernels, which will hopefully turn out to be easier to use in the context of high energy physics calculation. As a fundamental result of this paper, we will show that, when our set of integration kernels is mapped to the torus, then they (essentially) span the same space as defined by the multiple elliptic polylogarithms found in the mathematical literature. Of course, we find that multiple polylogarithms are just a subset of multiple elliptic polylogarithms.
While in this paper we will focus mainly on the mathematical details at the basis of the construction of elliptic polylogarithms, we will also provide two examples to demonstrate the flexibility of our setup and its applicability to a variety of different problems. An interesting and very non-trivial mathematical problem is the computation of the Laurent expansion of different kinds of hypergeometric functions, when the small expansion parameter (say ) appears in their indices. In particular, it is well known that if the indices are expanded around half-integer values, at order zero in the expansion one is in general left with complete elliptic integrals of first and second kind. In this case, however, no closed analytical formula in terms of familiar functions is known for the higher-order Laurent coefficients. As a natural application of our formalism, we will show how both in the case of the Gauss' hypergeometric function 2 F 1 and of the Appel function F 1 , the coefficients of their Laurent expansion around half-integer values can easily and algorithmically be expressed in terms of our elliptic polylogarithms. The computations presented here are mainly of mathematical interest. The most prominent example in high-energy physics, the sunrise integral, will be calculated in our elliptic language in a companion paper [40], which we provide as a hands-on guide to the use of our functions in high-energy physics for those who prefer not to have to dwell into all mathematical details reported here.
The outline of the paper is as follows: after reviewing multiple polylogarithms as the result of integrating rational functions on the Riemann sphere in section 2, we present the main result of the paper in section 3: a class of multiple elliptic polylogarithms. In section 4 we describe iterated integrals in the torus formulation of elliptic curves and discuss the relation between the torus and our formulation in section 5. A general algorithm for solving integrals with roots of cubic polynomials based on our formalism is provided in section 6. In section 7 we extend our formulation to elliptic curves defined by a general quartic polynomial, and in section 8 we use our framwork to compute several hypergeometric functions 2 F 1 and Appell F 1 functions. Finally, we draw our conclusions in section 9. In appendix A we discuss some technical aspects about regularisation of iterated integrals omitted in the main text.

Integrating rational functions on the Riemann sphere
In this section we give a short review of multiple polylogarithms and how they arise when integrating rational functions on the Riemann sphere. The material in this section is well known, but we discuss it in detail because it serves as a motivation and an illustration of the ideas that will be introduced in subsequent sections.
Let R(x) be a rational function on the Riemann sphere CP 1 with poles at x = c i , for some complex numbers c i . Our goal is to compute the integral of R(x) along some path on the Riemann sphere. The precise form of the path is immaterial for our purposes, and we are only interested in computing a primitive of R. The value of any definite integral along some path can then be obtained by comparing the value of the primitive at the endpoints of the path (plus taking into account possible windings around poles).
Using partial fractioning, we can reduce the computation of the integral to a linear combination of elementary integrals 2 , The previous relations are only valid for k = 1. Indeed, even if we start from a rational function R, it is not always possible to find a rational primitive, but we need to enlarge our space of functions. In this particular case, the obstruction to finding a rational primitive is related to the existence of the logarithm, The reason for the appearance of an obstruction to finding a rational primitive is tightly linked to the fact that the integrands in eq. (2.2) have simple poles. Indeed, the residue theorem implies that the integral along a closed curve encircling a simple pole does not vanish, and as a consequence the integral defines a multi-valued function. Since rational functions are always single-valued, the integral of a function with a simple pole can in general not be expressed in terms of rational functions alone. Next, we want to extend this analysis to iterated integrals over rational functions on the Riemann sphere. We then need to consider new kinds of obstructions, which arise from integrals that cannot be expressed in terms of rational functions and logarithms alone. Here we consider the simplest obstructions that one can encounter when studying iterated integrals on the Riemann sphere, multiple polylogarithms (MPLs), which are defined recursively by [3] G(c 1 , . . . , c n ; x) = x 0 dt f (c 1 , t) G(c 2 , . . . , c n ; t) , The recursion starts at G(; x) = 1. We assume that the c i are independent of x. In the case where (c 1 , . . . , c n ) = (0, . . . , 0) the integral in eq. (2.3) is divergent, and we define instead MPLs satisfy some well-known properties, which we now quickly review. First, MPLs form a shuffle algebra, which allows one to express the product of two MPLs as a linear combination of the same class of functions, where the sum runs over all shuffles of c 1 and c 1 , i.e., over all permutations of c 1 ∪ c 2 that preserve the relative orderings within c 1 and c 2 . Second, if we consider the vector space of all MPLs with coefficients that are rational functions, then this space is closed under taking primitives. More precisely, if R denotes the field of rational functions with poles at most at points c i ∈ S, let A MPL denote the R-algebra generated by all MPLs with singularities at most for x = c i . A MPL is graded by the weight (the multiplication is given by the shuffle product in eq. (2.6)), with A MPL is closed under both differentiation and integration. In particular, for every f ∈ A MPL , there is a primitive F ∈ A MPL such that ∂ x F = f . The goal of this paper is to explain how MPLs and their properties can be generalised to elliptic curves. We define a class of iterated integrals with logarithmic singularities that serve as a basis for obstructions that one can encounter when integrating a rational function on an elliptic curve, and we show that their properties are very similar to the properties of MPLs.

A class of elliptic polylogarithms
The goal of this section is to present the main result of this paper: a generalisation of multiple polylogarithms to elliptic curves. The content of this section is self-contained and provides a summary of our main results. In addition, we present a concise review of the minimal background on elliptic curves needed to understand the construction of elliptic polylogarithms. For a more thorough review of the mathematical background on elliptic curves and for proofs of the results, we refer to subsequent sections and to the mathematical literature (see, e.g., ref. [41]).

A lightning review of elliptic curves and elliptic integrals
Consider a cubic polynomial of the form For concreteness, we assume that the a i are all distinct and real, and they are ordered according to a 1 < a 2 < a 3 . Such a polynomial defines an elliptic curve E as the solution set of the polynomial equation Throughout this section we only consider elliptic curves defined by a cubic polynomial, and we will discuss elliptic curves defined by a quartic polynomial in section 7.
In the following it is convenient to work in projective space CP 2 , and we interpret the polynomial equation in terms of homogeneous coordinates [x, y, 1]. Seen as a curve in projective space, E also contains the point at infinity [0, 1, 0]. The points [a i , 0, 1] always lie on the elliptic curve. In the following we will refer to the points [0, 1, 0] and [a i , 0, 1] as branch points (by abuse of language, we will often refer to the a i and a 4 ≡ ∞ themselves as branch points). In addition, for every x = a i , there are precisely two points [x, ±y, 1] on E.
An elliptic curve defines a compact Riemann surface of genus one. It is natural to ask what is the appropriate generalisation of a rational function to an elliptic curve. A rational function on the elliptic curve E is defined to be a rational function in the two variables (x, y), subject to the constraint y 2 = P 3 (x). In other words, a rational function R on E is an expression of the form where p i and q i are polynomials in x. If we multiply the numerator and the denominator by the conjugate of the denominator, q 1 (x) − q 2 (x) P 3 (x), then we can write R in the alternative form for some rational functions R i . At this point we have to make a comment about the choice of the branch of the square root. In the case where the roots of P 3 are real and ordered according to a 1 < a 2 < a 3 , we let Our goal is to study generalisations of polylogarithms to elliptic curves. We start by discussing the simpler case of a single integral of a rational function on E. Such integrals have been studied extensively in mathematics during the 19 th century under the name of elliptic integrals. We briefly review the computation of elliptic integrals, as some of the concepts will prove useful to define iterated integrals on elliptic curves. Using the decomposition in eq. (3.4), we see that the contribution from R 1 (x) is an ordinary integral of a rational function and can be performed in terms of rational functions and logarithms. We will therefore focus on the contribution from R 2 (x). After partial fractioning, we only need to consider integrals of the form where k is an integer and c is a constant. Using integration by parts, we can reduce every integral of this type to a linear combination of the following integrals: .
These integrals cannot be simplified further, and should be thought of as the analogues of the concept of 'master integrals' familiar from the the physics literature. Via a judicious change of variables, each of these integrals can be evaluated in terms of (incomplete) elliptic integrals of the first, second and third kind 3 (3.8) The elliptic integrals in eq. (3.7), or equivalently in eq. (3.8), provide elementary obstructions to finding a rational primitive when integrating on an elliptic curve. This is in complete analogy with the role played by the logarithm in the case of the Riemann sphere, and we will see at the end of this section that, in a certain sense, we can indeed identify the incomplete elliptic integral of the third kind Π with an elliptic generalisation of the logarithm.
The elliptic integrals in eq. (3.7) allow one to define certain 'invariants' that are attached to an elliptic curve. The periods of E are defined by integrating dx/y between two branch points. We define with and K denotes the complete elliptic integral of the first kind, K(λ) = F(1|λ). If E is nondegenerate (in particular if the roots of P 3 are distinct), then the two periods are linearly independent over the real numbers. Moreover, all other integrals of this type evaluate to integer linear combinations of the ω i , Note that when the roots are real, ω 1 can be chosen as real and positive, and ω 2 has a positive imaginary part. Similarly, the quasi-periods of E are defined by with E(λ) = E(1|λ) and we defined where s n ( a) ≡ s n (a 1 , a 2 , a 3 ) is the elementary symmetric polynomial of degree n in three variables. Just like in the case of the periods, a different choice of the branch points in the integration limits will result in an integer linear combination of η 1 and η 2 . The periods and quasi-periods are not independent, but they are related through the Legendre relation, (3.14) The periods and quasi-periods are strictly speaking not invariants of E, but there may be different values of ω i and η i that correspond to the same elliptic curve (for example, we can perform a global rescaling of the periods without changing the geometry). There is an invariant, called the j-invariant, that uniquely characterises an elliptic curve, Two elliptic curves that have the same j-invariant are isomorphic.

A class of elliptic polylogarithms
We now introduce a class of elliptic generalisations of MPLs, (3.16) with n i ∈ Z and c i ∈ C ≡ C ∪ {∞}, and the recursion starts with E 3 (; x, a) = 1. The vector a = (a 1 , a 2 , a 3 ) encodes the zeroes of the polynomial P 3 , and therefore defines the elliptic curve E. In the following we will always assume that the elliptic curve is fixed, and we will keep the dependence of all quantities on a implicit. The integration kernels ϕ n that appear in eq. (3.16) are chosen such as to satisfy the following basic properties: 1. The functions ϕ n are non-trivial, i.e., they cannot be written as total derivatives of a rational function on E (because otherwise the integration would be trivial). As such, they will be tightly related to the irreducible integrands encountered in eq. (3.7).
2. The kernels are linearly independent in the sense that there is no linear combination that evaluates to a total derivative.
3. Since elliptic polylogarithms should have at most logarithmic singularities (but no poles), each ϕ n can have at most simple poles.
The explicit form of the integration kernels is discussed in the remainder of this section, together with some of the main properties of the iterated integrals E 3 . For n = 0, we define The right-hand side of eq. (3.17) is independent of c, and we will often set c = 0 in ϕ 0 . The integral of ϕ 0 is closely related to the incomplete elliptic integral of the first kind F in eq. (3.8). ϕ 0 defines a rational function on E that is free of poles. Indeed, ϕ 0 has no poles for any finite value of x. Letting 4 x = 1/u 2 , we see that there is no pole at infinity, The quantity dx/y is often referred to as the holomorphic differential on E. While ϕ 0 is free of poles, the functions ϕ ±1 (c, x) have a simple pole at x = c, where we define y c ≡ P 3 (c). Let us make some comments about these functions. First, using integration by parts we see that any integral involving ϕ −1 with c = a i can be reduced to simpler integrals. Hence, the functions ϕ −1 (a i , x) are not part of our basis of integration kernels. Second, we see that the function ϕ 1 is independent of the branch points a i , and it agrees with the integration kernel f in eq. (2.4). In other words, MPLs are a subset of the elliptic polylogarithms, E 3 This is important in physics applications, where it often happens that the differential equations satisfied by an elliptic Feynman integral involve ordinary MPLs in the inhomogeneous term. Equation (3.7) contains an additional integral which is not covered by the integration kernels defined so far, namely the integral over x dx/y. We have (using again x = 1/u 2 ) The integrand in eq. (3.21) has a double pole at infinity (or equivalently at u = 0), so it violates our criterion of a basis of integration kernels with at most simple poles. We can, however, obtain an independent function with a simple pole by computing the primitive in eq. (3.21). In addition, we have some freedom in how to define the basis element that corresponds to a double pole at infinity, because we can add any function without poles to the integrand, and in particular any multiple of the holomorphic differential. Here, we define The term proportional to the quasi-period η 1 is purely conventional at this point, and its role will only become clear in section 5. We define a primitive of Φ 3 by The function Z 3 is regular everywhere, except at infinity where it has a simple pole, We therefore include the function Z 3 into our basis of integration kernels with at most simple poles, To summarise, we see that for each c ∈ C \ {a 1 , a 2 , a 3 } there are two basis elements ϕ ±1 (c, x) with a simple pole at x = c. If c = a i , then only ϕ 1 (a i , x) appears. One of the main differences between polylogarithmic functions on curves of genus zero and one is that, while in the case of the Riemann sphere it is possible to find a purely rational basis of integration kernels with at most simple poles, this is no longer true when working on an elliptic curve, and we need to include the transcendental function Z 3 .
The fact that Z 3 is not rational has further consequences. In the case of ordinary MPLs, we can always reduce any product or power of the integration kernels f in eq. (2.4) to a linear combination of integrals of the same kernels. For example, any product can be linearised using partial fractioning, and all higher powers can be recast in the form of a total derivative in x, Similar relations, however, do not exist for the function Z 3 , and so we need to complement our basis by all possible products and powers that involve this function. Let us start by discussing the case of higher powers of Z 3 . Since Z 3 has a simple pole at infinity, Z n 3 has a pole of order n, and so it violates the condition that our integration kernels should have at most simple poles. We therefore need to include appropriate subtraction terms that allow us to remove the poles of higher order. While there is some arbitrariness in the precise definition of the subtraction terms, a suitable choice will be constructed in section 5. More precisely, based on the results of refs. [29,39,42], we define in section 5 a family of polynomials Z n of degree n such that Z (n) The coefficients appearing in Z n are themselves polynomials in (x, y). The explicit form of the polynomials Z n is rather involved, so we content ourselves to only present the explicit result for Z (3.28) Using the asymptotic behaviour of Z 3 in eq. (3.24), it is easy to check that the previous expressions remain finite as x → ∞ (Note that eq. (3.5) implies y < 0 for large x). We define the integration kernels ϕ n (∞, x), with n ≥ 2, by The previous equation remains valid for n = 1 if we define Z 3 (x) = Z 3 (x). We expect that in applications to two-loop Feynman integrals only Z (n) 3 for small values of n show up, because the transcendental weight of a two-loop amplitude in four dimensions is constrained to be less or equal to four 6 .
Next, we need to consider products between ϕ ±1 (c, x) and powers of Z 3 . We define (3.30) Let us make some comments about these definitions. Since Z (n) 3 is free of poles for n ≥ 2, the integration kernels in eq. (3.30) have at most simple poles. The term proportional to Z 3 (x) Z (n−1) 3 (x) is included to remove the pole at infinity of dx/(x − c). We use integration by parts to show that we do not need to consider kernels of the form ϕ −n (a i , x), because they can always be reduced to other classes of integrals.
To summarise, the integration kernels that define the elliptic polylogarithms in eq. (3.16) are given by the compact formulas (valid for n > 1), (3.31) The functions ϕ n (c, x) form the complete set of integration kernels needed to define elliptic polylogarithms. For each c ∈ C\{a 1 , a 2 , a 3 } there are two infinite families ϕ ±n (c, x), n ∈ N, with a simple pole at x = c. If c = a i , 1 ≤ i ≤ 4 (with a 4 = ∞), the family with negative index is reducible, and we only need to consider the single infinite family ϕ n (a i , x). While it is manifest from eq. (3.31) that the ϕ ±n have at most simple poles, we have not yet shown that they satisfy the remaining two conditions spelled out at the beginning of this section. In section 5 we show that the integration kernels ϕ ±n are in one-to-one correspondence with the integration kernels that define the elliptic polylogarithms considered in refs. [29,39,42]. The latter are shown to define a complete and independent set of integration kernels in ref. [39], and so the same holds true for the functions defined in eq. (3.31).
Before proceeding further, we stress that some of the elliptic polylogarithms obtained by integrating once over the simplest kernels in eq (3.31), can be written in terms of the elliptic integrals in eq. (3.8). We find, for 0 < x < a 1 < a 2 < a 3 < c, Similar relations can be obtained for other regions. Note that the signs in the previous formula are connected to our prescription for the branches of the square root in eq. (3.5). We see from the previous formula that the elliptic integral of the third kind Π plays a role very similar to to the logarithm in the case of the Riemann sphere. We cannot give a closed analytical expression for the integral of ϕ 1 (∞, x), because it requires integrating first over Z 3 (x), cf. eq. (3.31). We can, however, provide an analytic expression for Z 3 (x) itself in terms of incomplete elliptic integrals. It is more convenient to do it in the region 0 < a 1 < a 2 < x < a 3 , where we have Let us discuss some of the properties of the iterated integrals defined by eq. (3.31). We have presented in detail the case of an elliptic curve defined by a cubic polynomial, cf. eq. (3.2). The results of this section are not specific to cubic polynomials, and in section 7 we define a similar set of functions for elliptic curves defined by a quartic polynomial. While the analytic form of some of the kernels, in particular of Φ 3 and Z 3 , are different, the overall picture and the counting of the independent integration kernels are the same: for each c ∈ C \ {a 1 , . . . , a 4 } there is a double infinite family of integration kernels, while for each branch point there is a single infinite family.
The functions E 3 in eq. (3.16) satisfy all the basic properties of iterated integrals. In particular, they form a shuffle algebra and similarly for d. At this point we have to make a comment about potential divergencies at x = 0. Since the lower integration boundary in eq. (3.16) is t = 0, the integral potentially diverges whenever c k = 0. We need to regularise this divergence (e.g., by introducing a suitable subtraction term), and there is a certain degree of arbitrariness in this procedure. However, we would like to do this in a way that preserves the shuffle algebra structure. In the case of ordinary MPLs, this is achieved by a modified integration contour, see eq. (2.5). In appendix A we present a regularisation for E 3 that preserves the shuffle algebra structure.
Our paper is not the first to consider elliptic generalisations of MPLs [12,15,17,18,22,29,39,42,43] or iterated integrals over integration kernels that involve square roots of polynomials of degree greater than two [44] have been considered. In particular, in refs. [29,39,42] elliptic polylogarithms have been defined as iterated integrals on a punctured torus, and a complete basis for the corresponding integration kernels on the torus has been defined. Moreover, it was proven in ref. [39] that every iterated integral on the punctured torus can be expressed in terms of these functions. In the remainder of this paper we study the relationship between the elliptic polylogarithms of refs. [29,39,42] and the iterated integrals defined in eq. (3.16). One of our main results is that the iterated integrals E 3 are equivalent to the elliptic polylogarithms of refs. [29,39,42] 7 , and we can cast every iterated integral E 3 as a linear combination of the elliptic polylogarithms of refs. [29,39,42], and vice-versa. A concrete algorithm how to perform this translation is presented in section 5.
Let us conclude this section by discussing another similarity between ordinary and elliptic polylogarithms, which is at the same time one of the main results of this paper. At the end of section 2 we have seen that the algebra A MPL generated by all MPLs with rational functions as coefficients is closed under both differentiation and integration. There is a similar result for the elliptic polylogarithms. Indeed, in ref. [39] it was shown that every iterated integral on the punctured torus can be expressed in terms of elliptic poylogarithms (and rational functions). Since our iterated integrals define essentially the same class of functions as the elliptic polylogarithms of refs. [29,39,42], we conclude that one can evaluate all iterated integrals on an elliptic curve E in terms of the iterated integrals E 3 . More precisely, let R 3 ≡ C(x, y)/ y 2 = P 3 (x) denote the field of rational functions of the elliptic curve E. We consider the algebra A 3 over R 3 generated by Z 3 (x) and all elliptic polylogarithms Seen as a vector space over R 3 , the algebra A 3 admits the presentation We call the quantity l = m + k the total length. The algebra A 3 shares many of the properties of A MPL of section 2. First, it is easy to check that A 3 is filtered by the total length 8 , where elements of A 3,l are linear combinations of the form m, n, c m+| n|≤l Second, A 3 is closed under differentiation with respect to x. If the coefficients in the linear combination (3.37) are constants, then differentiation lowers a non-zero total length by one unit, because it lowers the length of an elliptic polylogarithm, and the derivative of Z 3 is a rational function. In particular, we see that if Z (m) x has total length m + | n| − 1. Finally, while the closure under differentiation is immediate, it is less obvious to see that A 3 is also closed under integration. More precisely, for every . We provide an explicit and constructive proof of this statement in section 6, where we present an algorithm to explicitly compute the primitive. This algorithm is in fact a generalisation of the classical algorithm to evaluate integrals of rational functions on elliptic curves reviewed at the beginning of this section, and it extends this classical algorithm to include elliptic polylogarithms and the function Z 3 . Said differently, the classical algorithm allows one to find a primitive in the space A 3,0 of functions of total length zero (which are just rational functions), while our extension generalises it to functions of arbitrary total length. 0 Figure 1: The lattice Λ spanned by the two periods ω 1 and ω 2 . The grey-shaded area is the fundamental parallelogram.

Elliptic curves and iterated integrals on a torus
The aim of this section is to provide the necessary mathematical background to understand how the iterated integrals E 3 defined in the previous section are connected to the elliptic polylogarithms that appear in the mathematics literature [39,42]. The material in this section is not new and is in principle well known. We include it nonetheless because we feel that some of these topics are rarely discussed in the Feynman integral literature.

Elliptic functions
Consider two complex numbers ω 1 and ω 2 that are linearly independent over the real numbers. For concreteness, we will assume that ω 1 is real and positive while ω 2 is purely imaginary with a positive imaginary part. We can define a lattice (see fig. 1) Such a lattice is a discrete additive subgroup of C. The torus associated to the lattice Λ is defined as the quotient C/Λ of C by the lattice. In other words, we identify two complex numbers whenever they differ by an element from the lattice. The torus is then obtained by identifying opposite sites of the fundamental parallelgram {r ω 1 + s ω 2 : 0 ≤ r, s < 1}. We now study functions on the torus C/Λ. In order to be well-defined, any function on the torus must be invariant under translations by the periods ω i , i.e., it must be a periodic function, f (z+ω i ) = f (z), i = 1, 2. An elliptic function is a meromorphic periodic function. The singularity structure of an elliptic function is very constrained. In particular, every non-constant elliptic function must have at least two poles on the torus (counted with multiplicity, e.g., a double-pole counts as two poles), and the number of zeroes must equal the number of poles (again, counted with multiplicities). The prototypical example of an elliptic function is the Weierstrass ℘ function, (4.2) We will always keep implicit the dependence of the Weierstrass ℘ function on the periods. The Weierstrass ℘ function is by construction periodic, and it has a double pole at every lattice point z ∈ Λ. In addition, it is an even function, ℘(−z) = ℘(z). The derivative of a periodic function is still periodic, and so the derivative of an elliptic function is itself elliptic. Since ℘ is even and has a double pole, its derivative must define an odd function with a triple pole (all definitions are understood modulo translations by the lattice). Hence, ℘ must have three zeros, and these are located precisely at the half-periods and so ℘ (ω i /2) = 0. The Weierstrass ℘ function and its derivative are not only the prototypical examples of functions that are both meromorphic and periodic, but they play a fundamental role in the theory of elliptic functions. In fact, they are sufficient to recover all elliptic functions. More precisely, one can show that every elliptic function can be written as a rational function in ℘ and ℘ . Note that the set of all elliptic functions forms a field, and so we can identify the field of elliptic functions with the field of rational functions in (℘, ℘ ).
Finally, let us note that the Weierstrass ℘ function and its derivative are not independent, but they are coupled via a non-linear differential equation where g i and e i are constants that depend on the two periods ω 1 and ω 2 . By differentiation we see that all the higher derivatives of ℘ are polynomials in (℘, ℘ ), in agreement with the fact that every elliptic function is a rational function in (℘, ℘ ). For example, we have

From the torus to the elliptic curve: the Weierstrass model
In section 3 we have considered elliptic curves given by a cubic equation of the form (3.2). In general, there can be several cubic polynomials that define the same elliptic curve E. In particular, it can be shown that via a judicious change of variables every elliptic curve can be represented as the solution set of a cubic equation of the form An equation of this form is called a Weierstrass equation of the elliptic curve. We see that, upon identifying (℘, ℘ ) with (x, y), eq. (4.4) has precisely the form of the Weierstrass equation, and so there is a strong connection between the Weierstrass ℘ function and elliptic curves. The purpose of this section is to show that the torus and the Weierstrass ℘ function provide a natural parametrisation of any elliptic curve E.
Consider an elliptic curve E with periods ω 1 and ω 2 given by the Weierstrass equation (4.6). Comparing eq. (4.6) to the differential equation (4.4) satisfied by the Weierstrass ℘ function, we see that for every point z on the torus C/Λ, the point [℘(z), ℘ (z), 1] lies on the elliptic curve E. We may thus ask the converse question: given a point [x, y, 1] on E, can we find a point z on the torus such that [x, y, 1] = [℘(z), ℘ (z), 1]? It turns out that the answer to this question is positive. If y = 0, we know that ℘ vanishes on the half-periods, and so we have [e i , 0, 1] = [℘(ω i /2), 0, 1]. If y = 0, then the function f (z) ≡ ℘(z) − x is an elliptic function with a double pole at z = 0 on the torus. Since every elliptic function must have the same number of zeroes and poles, f must have two zeroes (or a double zero) on the torus, and so the equation ℘(z) = x has always two solutions, which differ by a sign because ℘ is even. Hence, there is a one-to-one mapping between the points of the torus C/Λ and the elliptic curve E, and the torus provides a canonical way to parametrise any elliptic curve. The map is explicitly given by Since ℘ has a pole at the origin, we see that the point z = 0 on the torus is mapped to the point at infinity on the elliptic curve. The inverse map from the curve E to the torus can also be given explicitly, but before we do so, we discuss in more detail what the map in eq. (4.7) implies for the structure of functions and integrals on the elliptic curve and the torus.
We start by analysing what elliptic functions correspond to under the map in eq. (4.7). We know that every elliptic function is a rational function in (℘, ℘ ). If R is a rational function in two variables, then under the map in eq. (4.7) the elliptic function R(℘(z), ℘ (z)) is mapped to the rational function R(x, y) on the elliptic curve. In other words, under the map in eq. (4.7) the field of elliptic functions is mapped to the field of rational functions on E.
Next, we review some standard material on differential forms on an elliptic curve. An abelian differential on E is a differential one-form of the form dx R(x, y), where R is a rational function on E. It is customary to consider three different types of abelian differentials. An abelian differential of the first kind is holomorphic everywhere on E, and so in particular it has no poles on E. An abelian differential of the second kind is meromorphic, i.e., it is allowed to have poles on E, but the residue at every pole must vanish. Finally, a meromorphic differential with non-vanishing residues is called an abelian differential of the third kind.
Let us now discuss what happens to abelian differentials under the correspondence (4.7). First, it follows from section 3.1 that we can always reduce any abelian differential to a linear combination of the differentials in eq. (3.7). Let us analyse each of these differentials in turn. First, the differential dx/y has no pole, and thus corresponds to a differential of the first kind. Under the map (4.7), dx/y corresponds to the standard holomorphic differential dz on the torus. Indeed, letting x = ℘(z) and using eq. (4.4), we find Since the differential x dx/y in eq. (3.21) gives rise to a double pole without residue at infinity, it defines a differential of the second kind. On the torus it corresponds to the differential ℘(z) dz = dz/z 2 + O(z 0 ). Finally, the differential dx/(y(x − c)) has a simple pole at x = c with residue y c ≡ P 3 (c), and so it defines a differential of the third kind. In general, it is easy to see that an abelian differential corresponds to a differential of the form f (z) dz on the torus, where f is an elliptic function. This in turn puts very strong restrictions on the pole structure of any abelian differential. Since an elliptic function must have at least two poles, it is not possible to find an abelian differential with only one simple pole. This is tightly connected to the fact that the primitive of an elliptic function is not periodic, because the choice of a lower integration boundary breaks the invariance under translations by a period. Instead, the primitive of an elliptic function defines a quasi-periodic function, i.e., a function F such that where C i is a constant that may depend on the period ω i , but it is independent of z. If f is an elliptic function, then any primitive of f has the form F (z) = z z 0 dz f (z ), for some fixed point z 0 . The function F then satisfies eq. (4.9) with The prototypical example of a quasi-periodic function is the Weierstrass zeta function, defined as a primitive of the Weierstrass ℘ function (up to a sign). For z inside the fundamental parallelogram it is given by the integral The Weierstrass zeta function is not invariant under lattice translations, but it transforms according to eq. (4.9), where η i ≡ ζ(ω i /2) are the quasi-periods defined in eq. (3.12). We see from eq. (4.10) that the Weierstrass zeta function has a simple pole at z = 0 (and hence at every point of the lattice Λ). As we will see, it is the primary building block to define differential one-forms with at most simple poles on the torus. As anticipated, this building block is no longer a periodic function, but it is only quasi-periodic. Let us conclude this review with a discussion of the inverse of the map in eq. (4.7). Consider a point [a, b, 1] on the elliptic curve, and for simplicity we assume b > 0. Then we can find a point z a on the torus such that (℘(z a ), ℘ (z a )) = (a, b). The value of z a is given by Abel's map,

Elliptic polylogarithms
In this section we review the construction of elliptic polylogarithms as iterated integrals of differential one-forms with (at most) simple poles on the torus, following closely refs. [29,39,42]. We define a class of iterated integrals by where z i are complex numbers (assumed constant) and n i ∈ N are positive integers. Note that in some cases the integral may require regularisation at z = 0. In the following, we ignore this technical subtlety and refer to the literature on how to consistently define a regularised version of elliptic polylogarithms (cf., e.g., ref. [29]). The integration kernels are defined through a generating series known as the Eisenstein-Kronecker series, where the quantities E j and G j in the exponential denote the Eisenstein functions and series respectively, 9 The Eisenstein series G j can be cast in the form of polynomials in the parameters g 2 and g 3 appearing in the Weierstrass equation. From the series definition of E j we see that the Eisenstein functions satisfy the differential equation Moreover, comparing eq. (4.15) to eq. (4.2), we see that E 2 (z) = ℘(z)+G 2 . Hence, for j > 2 the Eisenstein functions can be expressed in terms of the derivatives of the Weierstrass ℘ function, For j = 1, instead, eq. (4.16) implies that E 1 is a primitive of the Weierstrass ℘ function, and so it is connected to the Weierstrass zeta function. More precisely, we have We see that g (1) (z) has a simple pole at every lattice point, and in particular for z = 0, but it is regular everywhere else. However, g (1) (z) is not invariant under translations by ω 2 , while it is for ω 1 , Since the Eisenstein functions for j ≥ 2 are elliptic functions, the exponential form of eq. (4.14) implies that the functions g (n) for n ≥ 2 can be expressed as polynomials of degree n in g (1) whose coefficients are elliptic functions, where G n is a polynomial of degree n in g (1) (z). For example, we have (4.21) Note that the form of these polynomials is very reminiscent of the polynomials in eq. (3.28). It follows from the exponential form of the Eisenstein-Kronecker series in eq. (4.14) that the two leading coefficients in the polynomial G n are very simple constants, where the dots indicate a polynomial in g 1 of degree n − 2 at most. Let us discuss some properties of the functions g (n) . First, it is easy to see that the g (n) are always invariant under translations by ω 1 , but not by ω 2 . Second, the g (n) are functions with definite parity, and satisfy the Fay identity [39], The Fay identity is a generalisation of partial fractioning for the function f appearing in the definition of ordinary MPLs, cf. eq. (3.26). Finally, it can be shown that for n ≥ 2, the g (n) are regular on the whole complex plane. For example, we have g (1) (z) = 1/z + O(z 0 ) and ℘(z) = 1/z 2 + O(z 0 ). Inserting this into eq. (4.21), we see that all the poles in g (2) (z) and g (3) (z) cancel. Let us conclude this section by making a comment on the iterated integrals defined by eq. (4.13) and the elliptic polylogarithms defined in refs. [29,39,42]. We have seen that the g (n) are not periodic, and so they are strictly speaking not well-defined functions on the torus. This is connected to the fact that there is no elliptic function with just a single simple pole, i.e., there is no function that is both meromorphic and periodic with just one simple pole. Instead of giving up periodicity to define the integration kernels g (n) , we could also consider periodic functions that are not meromorphic, i.e., that depend explicitly on the complex conjugatez. This is the approach taken in refs. [29,39,42], where elliptic polylogarithms are defined through the iterated integrals where the functions f (n) are defined by the generating series The functions f (n) have the same properties as the functions g (n) , except that they are invariant under translations by both ω 1 and ω 2 and have an explicit dependence on the complex conjugate variablez. Note that the dependence on the antiholomorphic variable is simple, because it only arises from the non-holomorphic exponential factor in eq. (4.26).
In general, we have to give up either meromorphicity or periodicity in order to define elliptic polylogarithms. While in the mathematics and the string theory literature it is more common and natural to preserve periodicity, we prefer to work with functions that are manifestly meromorphic, at the expense of giving up periodicity. Our choice is motivated by the following considerations: integrands of Feynman integrals present themselves as purely meromorphic objects, so it feels unnatural to introduce an explicit dependence on the antiholomorphic variable. In addition, it is often much easier to work with meromorphic expressions. For example, it is not possible to integrate by parts in any naive way when working with non-meromorphic functions, because in that case f (z,z)dz is not a total derivative, and so Stokes' theorem does not apply. Finally, we emphasise that in many practical applications the distinction between Γ and Γ is immaterial. Indeed, whenever the integration contour is parallel to the real axis (which happens regularly in applications), the non-holomorphic exponential factor in eq. (4.26) is constant, and so Γ and Γ are related in a trivial way. In addition, as we will see in the next section, in many cases of interest we can exchange g (1) with f (1) without loosing any information.

Relating iterated integrals on the torus and on the elliptic curve
In the previous sections we have introduced two different classes of generalisations of polylogarithms to elliptic curves: in eq. (3.16) we have defined the functions E 3 as iterated integrals on an elliptic curve with periods ω 1 and ω 2 , and in eq. (4.13) we have defined the functions Γ as iterated integrals on the torus defined by the lattice Λ = Z ω 1 + Z ω 2 . In this section we show that these two sets of functions are just two different representations of the same space of functions under the map in eq. (4.7) which identifies the torus with the elliptic curve.
Before we can discuss how to relate the functions E 3 and Γ, we need to address the issue that the functions E 3 have been defined for elliptic curves defined by arbitrary cubic polynomials, cf. eq. (3.2), while so far we have only defined a map from the torus to an elliptic curve defined by a Weierstrass equation. While we can always change coordinates and write any cubic equation in Weierstrass form, it is convenient to define a map from the torus to CP 2 that lands us directly on the curve defined by the general cubic equation in eq. (3.2). In the following we define such a map, and we show that this map identifies the functions Γ on the torus with the elliptic polylogarithms E 3 .
Consider the elliptic function Using the properties of the Weierstrass ℘ function, one can check that µ satisfies the nonlinear differential equation where c 3 has been defined in eq. (3.10). It is then easy to see that the map sends the torus to E. Using exactly the same argument as for the Weierstrass equation, we can show that this map is invertible. In particular, under this map the point z = 0 is mapped to the point at infinity on E, while the half-periods ω i /2 map to the remaining branch points a i , The holomorphic differential dx/y pulls back to the standard holomorphic differential on the torus, dx/y = dz/c 3 . Just like in the Weierstrass case, every c = a i has two preimages ±z c on the torus such that µ(±z c ) = c. In the following, we assume without loss of generality that Re µ (z c ) > 0 (otherwise we exchange the roles of +z c and −z c ). In the remainder of this section, we show that under the map in eq. (5.3) the functions g (n) and Γ of refs. [29,39,42] map to the functions Z (n) 3 and E 3 defined in section 3. Let S be a finite set of points in C. We consider iterated integrals of rational functions that have poles at most at points in S. For concreteness, we assume that S contains ∞, and we define S ≡ S \ {∞}, and it does not contain any of the zeroes a i of P 3 , nor the point x = 0 (in order to avoid lengthy discussions about the regularisation of eq. (3.16) -see appendix A). We stress that these assumptions are not essential, but they allow us to avoid having to distinguish too many different special cases. Our goal is to show that under the correspondence in eq. (4.7) the complex vector space spanned by the differential forms dx ϕ n (∞, x) and dx ϕ ±n (c, x), c ∈ S (with n ≥ 0), is identified with the vector space spanned by the one-forms dz g (n) (z) and dz g (n) (z ± z c ). A direct consequence of this result is that the iterated integrals E 3 and Γ define the same class of functions, and so the functions E 3 coincide with the elliptic polylogarithms defined in the mathematics literature (up to the different treatment of periodicity vs. memorphicity, cf. section 4).
Let us start by analysing the differential form dx ϕ −1 (c, x). On the torus it corresponds to α is an elliptic function with simple poles only at z = ±z c with residues ±1. The difference is then free of poles, and thus regular everywhere. Using eq. (4.19), it is easy to check that the difference is periodic as a function of z, despite the fact that g (1) is not. Hence, the difference in eq. (5.6) defines an elliptic function without any poles, and must therefore be constant. The constant is easily determined by using the fact that α(z) vanishes for z = 0, and we find We see that on the torus the differential form dx ϕ −1 (c, x) corresponds to a linear combination with constant complex coefficients of the holomorphic differential dz and the forms dz g (1) (z ± z c ). We can apply exactly the same reasoning to the differential form The only difference with respect to the previous case is that the residues at z = ±z c are both +1, and there is also a simple pole at z = 0 with residue −2, corresponding to the pole at infinity of dx/(x − c). We subtract the poles, and by exactly the same reasoning as before we conclude that the difference must be constant. In order to determine this constant, we observe that both β and eq. (5.9) define an odd function, and so they must vanish at the origin. We then have Let us now turn to dx ϕ 1 (∞, x). If x 0 = ℘(z 0 ), eq. (3.23) gives Hence, we immediately find To summarise, the one-forms dx ϕ ±1 are linear combinations of the holomorphic differential and the one-forms dz g (1) . To complete the proof, we need to analyse what happens for n > 1. In section 3 we have not given the complete definition of the polynomials Z n for n > 1, but we have already noted the similarity for n = 2, 3 between the polynomials Z n in eq. (3.28) and G n in eq. (4.21). In general, we define For n = 1 we recover eq. (5.11). Equation (5.13) immediately implies dx ϕ n (∞, x) = c 3 dx y Z (n) 3 (x) = 4 dz g (n) (z) . (5.14) Next, let us turn to the one forms dx ϕ −n (c, x). Using eqns. (5.7) and (5.13), and applying the Fay identity (4.24), we find We can identify dx ϕ −n (c, x) with the combination dz g (n) (z − z c ) − g (n) (z + z c ) , up to terms that involve less than n powers of g (1) . Applying exactly the same reasoning, we see that the corresponding combination with a plus sign is provided by dx ϕ n (c, x), To summarise, we have shown that every differential form dx ϕ ±n can be written as a linear combination of the forms dz g (k) with complex coefficients that are constants with respect to z. As a consequence every iterated integral E 3 can be written as a linear combination with constant complex coefficients of elliptic polylogarithms Γ, and vice-versa. The functions E 3 are simply an alternative linear basis for the shuffle algebra of elliptic polylogarithms. The change of basis is encoded into the relations (5.7), (5.10), (5.13) and (5.15). For example, if c is not a branch point, eq. (5.7) implies where z 0 is the point on the torus such that ℘(z 0 ) = 0 and ℘ (z 0 ) > 0. Conversely, every elliptic polylogarithm Γ can be written as a linear combination of E 3 functions. For example, we have To conclude, let us make a comment about the connection between the iterated integrals Γ and Γ. The difference between the two sets of iterated integrals only lies in the fact that we have to give up either periodicity or holomorphicity to define the integration kernels (cf. the discussion at the end of section 4). In this section we have shown that we can write the iterated integrals E 3 in a natural way as linear combinations of elliptic polylogarithms Γ. In some cases, however, we can also write them in terms of their periodic and non-holomorphic analogues Γ. Indeed, using eq. (4.19) it is easy to check that the combination of g (1) functions in eq. (5.7) is periodic with respect to both ω 1 and ω 2 , and hence we can express it in terms of the periodic functions f (1) , The non-holomorphic contribution in f (1) cancels in the combination in the right-hand side. Similarly, we have In general, we see that we can write E 3 ( n 1 ... n k c 1 ... c k ; x) as a linear combination of the periodic elliptic polylogarithms Γ whenever |n i | ≤ 1 and c i = ∞. In particular, we can write the function E 3 ( −1 c ; x) in eq. (5.17) as

Iterated integrals on elliptic curves: an algorithmic approach
In the previous section we have shown that the functions E 3 are in fact an alternative basis for the space of elliptic polylogarithms Γ. In this section we consider integrals over elliptic polylogarithms multiplied by rational functions, and we present an algorithm to perform such integrals in terms of a well-defined class of functions. At the same time we prove the claim from section 3 that the algebra A 3 in eq. (3.35) is closed under integration. Since A 3 is filtered by the total length, we can restrict the discussion to integrands with a given total length l. Using partial fractioning, we can reduce the problem to the following four types of integrals, x has total length m + | n| = l. If l = 0, we recover the classical case of integrating a rational function on an elliptic curve (see section 3.1). Not all integrals in a given family are independent, but they are related by integration by parts. We discuss each family in turn.
The integrals in the A-family satisfy the recursion 10 , Note that the total length of ∂ x X is l − 1. The recursion therefore allows us to increase the value of k while at the same time to lower the value of the total length. We can therefore recursively lower the value of l, until we reach l = 0 and the integral is elementary. However, the recursion has a singularity for k = −1, and so we cannot reduce integrals of the type For integrals of type B, we only need to consider the case k > 0. We have the recursion We can lower the value of k, and at the same time lower the total length. The recursion has a singularity for k = 1, and so we are left with integrals of the form B c,1 [X ] that cannot be reduced.
Integrals of type C satisfy the recursion The recursion is non-singular for all integer values of k. It has depth three, and so we can express all integrals from the family C through three representatives of this family. We choose these irreducible integrals to correspond to C k [X ], k ∈ {−1, 0, 1}. Finally, the integrals in family D satisfy the recursion The recursion has depth three. However, for k ≤ 0, the integrals reduce to integrals of type C, and so we choose as many of our basis integrals as possible to lie in the family C. The only obstacle to this is the appearance of a singularity for k = 1 in the recursion. Hence, every integral in the family D can be reduced to the integrals D c,1 [X ], as well as integrals of type C.
The previous discussion only applies if c is not a zero of P 3 . Indeed, if for example c = a 1 , the coefficient of D c,k in eq. (6.5) vanishes, and so we cannot use the recursion to relate D c,k to other integrals. Instead, we have a 12 a 13 ) .
(6.6) 10 We do not include boundary terms in the IBPs, because we are only interested in primitives.
In this case the recursion is non-singular for all integer values of k, and so all integrals of the type D a 1 ,k [X ] can be reduced to integrals of type C, in agreement with the fact that we do not need to consider integration kernels ϕ −1 (a 1 , x).
To summarise, using integration by parts we can reduce all integrals to the following six classes of integrals, It is easy to see that these integrals are in one-to-one correspondence with the differential forms ϕ ±n , and they can be performed using the definition of the iterated integrals E 3 in eq. (3.16). The only case that needs some explanation are the integrals . Equivalently, we may consider the integrals In the previous equation, and until the end of this section, we keep the dependence of all quantities on x implicit. Since Φ 3 has a double pole at infinity, it is not part of our basis of integration kernels, but its primitive Z 3 is. We can integrate by parts and we obtain We have, with m > 1, and so we reproduce the original integral that we started from. Hence, the total length has not been lowered, and so we cannot apply our recursive argument based on the total length. In the following we show how this integral can be evaluated. We only discuss the case where E 3 is absent from the integrand, because the argument is independent of the appearance of E 3 . Using eq. (4.22), we can write where Z (m) R is a polynomial of degree at most m − 1 in Z 3 , and we find Hence, we have The right-hand side only involves polynomials in Z 3 of degree strictly less than m, i.e., it has a lower total length, and so recursively we know how to do these integrals.

Elliptic curves defined by quartic polynomials
So far we have only discussed elliptic curves that are defined through a cubic equation, cf. eq. (3.2), and we have not yet discussed what happens in the case of elliptic curves defined by a quartic polynomial. The quartic case is particularly relevant for physics, because, e.g., the maximal cut of the sunrise integral leads to an elliptic curve defined by a quartic polynomial (cf., e.g., refs. [14,[31][32][33]). Since every elliptic curve can be represented as the zero set of some cubic equation (e.g., its Weierstrass equation) via a suitable change of variables, the results of the previous sections are in principle sufficient to cover all possible elliptic curves. In practise, however, the change of variables to the cubic form can be rather cumbersome, and it may be preferable to have a formulation where one can directly work with general quartic polynomials. In this section we show that this can be achieved, and we formulate all the results of the previous sections for elliptic curves defined by general quartic polynomials, We assume again that the roots a i are real and ordered according to a 1 < a 2 < a 3 < a 4 , and we choose the branches of the square root as follows,

(7.2)
With this convention the periods can be written as and with λ = cr(a 1 , a 4 , a 3 , a 2 ) = a 14 a 23 a 13 a 24 and c 4 = 1 2 √ a 13 a 24 , (7.4) and the j-invariant is given by eq. (3.15). We stress here that eq. (7.2) implies a negative sign for the square-root in the region a 2 < x ≤ a 3 , which in turn determines the overall sign of ω 2 in eq. (7.3).
Let us now discuss the abelian differentials of the second kind on E. The differential x dx/y, which provided the differential of the second kind in the cubic case (cf. eq. (3.13)), is no longer a good candidate. Indeed, letting u = 1/x, we see that Equation (7.5) reveals that x dx/y has a simple pole at infinity, and hence it defines a differential of the third kind. Instead, a valid differential of the second kind in the quartic case is 11 We emphasise that it is mandatory to consider this particular linear combination. Indeed, x 2 dx/y is by itself not a valid differential of the second kind, because it has non-vanishing residue at infinity, We can add any multiple of the holomorphic differential to eq. (7.6), and we still obtain a differential of the second kind. We find it convenient to define As usual, we will keep the dependence on a implicit. The quasi-periods can then be expressed as elliptic integrals of the differential of the second kind (cf. eq. (3.12)),

From the torus to the elliptic curve
Just like in the case of an elliptic curve defined by a quartic polynomial, we can define a map from the torus to the elliptic curve E, where κ denotes the elliptic function withs n ≡ s n (a 2 , a 3 , a 4 ). Using the differential equation (4.4) satisfied by the Weierstrass ℘ function, it is easy to check that κ satisfies the non-linear differential equation a 4 ) . (7.12) 11 We will show in section 7.4 that, unlike in the cubic case, x 2 dx/y cannot be reduced using integrationby-parts identities, and so it defines a genuine master integral in the quartic case.
Hence the image of the torus under the map in eq. (7.10) is precisely the elliptic curve defined by the quartic equation (7.1). The holomorphic differential dx/y on E corresponds to dz/c 4 on the torus, and the inverse of eq. (7.10) is given by a variant of Abel's map in eq. (4.12) The half periods are naturally mapped by κ to the zeroes of the polynomial P 4 , (7.14) At this point we see a main difference between the cubic and quartic cases: while for a cubic polynomial the point at infinity is always a branch point, this is no longer the case for a quartic polynomial, and so the point at infinity of E is not the image under κ of a half-period. Indeed, we see from eq. (7.14) that κ is regular for every half-period, and it is singular only when the denominator in eq. (7.11) vanishes. Hence, there must be two points ±z * on the torus such that denominator in eq. (7.11) vanishes. The value of z * is determined by Abel's map in eq. (7.13), The point z * plays an important role in understanding the structure of differentials on E that have poles at infinity. We know for example that the differential of the second kind dx Φ 4 has a double pole at infinity, and so it must have double poles at z = ±z * on the torus. In the following we derive a formula which makes this explicit. We obviously have We know that κ(z) has poles at z = ±z * . Using the explicit definition of κ in eq. (7.11), it is easy to show that Since the Weierstrass ℘ function has a double pole at the origin, we conclude that the elliptic function is free of poles and thus constant. The value of the constant is where the value of ℘(z * ) is determined by the requirement that the denominator in eq. (7.11) vanishes. We thus have The previous formula makes explicit the fact that dx Φ 4 has a double pole at infinity on E, and thus double poles at z = ±z * on the torus. We can apply exactly the same reasoning to the differential of the third kind with a simple pole at infinity, and we find We see that the right hand side has simple poles both at z = z * and z = −z * , in agreement with the fact that an elliptic function cannot have a single simple pole.

Elliptic polylogarithms associated to curves defined by a quartic equation
The elliptic polylogarithms attached to an elliptic curve defined by a quartic equation are defined in complete analogy to the cubic case in eq. (3.16), with n i ∈ Z and c i ∈ C, and the recursion starts with E 4 (; x, a) = 1. The branch points are encoded in the vector a = (a 1 , a 2 , a 3 , a 4 ). We will keep the dependence of all quantities on a implicit from now on. The integration kernels ψ n take a form very similar to the kernels that appear in the cubic case, cf. eq. (3.31). We can write for n ≥ 2 Note that the choice of the lower integration boundary again breaks the symmetry among the branch points a i . It is easy to check that Z 4 has a simple pole at infinity, and it is regular everywhere else. The functions Z (n) 4 for n > 1 are defined via the polynomials G n introduced in eq. (4.20), ; p 1 (x), p 2 (x, y) , (7.26) with , The form of the arguments of G n in eq. (7.26) will be motivated in the next section.
As in the cubic case, we can rewrite the integrals over some of the differential forms above in terms of incomplete elliptic integrals of first, second and third kind defined in eq. (3.8). Assuming the ordering 0 < x < a 1 < a 2 < a 3 < a 4 < c, it is easy to see that where the cross ratio function cr is defined in eq. (7.4). Similarly, we find, for a 1 < x < a 2 , Similar relations exist for other regions. The signs in the previous relations are related to our convention for the branches of the square root in eq. (7.2). Before we discuss in more detail the relationship between the iterated integrals E 4 and the elliptic polylogarithms Γ in eq. (4.13), let us compare the integration kernels in eq. (7.23) to their cubic analogues in eq. (3.31). First, just like in the cubic case, eq. (7.22) contains ordinary MPLs as a special case, Second, we see that a main difference between the cubic case in eq. (3.31) and the quartic case in eq. (7.23) is the appearance of the integration kernel ψ −n (∞, x), which is absent in the cubic case. This is due to the fact that the point at infinity is a branch point in the cubic case, while it is not for a quartic polynomial. Conversely, we do not need to consider integration kernels of the form ψ −n (a i , x), because they can always be removed using integration by parts. More generally, the counting of the integration kernels is identical in the cubic and quartic cases: for each c = a i there is a double infinite family of integration kernels, while for each branch point a i the infinite family ψ −n (a i , x) is absent. Finally, another difference to the cubic case is the appearance of terms proportional to δ n2 . These terms remove double poles at infinity that arise when multiplying Z 4 by another function with a simple pole at infinity. The role of these terms will become more transparent in the next section when discussing the connection between the iterated integrals E 4 and the elliptic polylogarithms Γ.
Let us conclude by mentioning that the elliptic polylogarithms E 4 satisfy all the standard properties of iterated integrals. In particular, they form a shuffle algebra. In Appendix A we discuss regularisation procedure to regulate divergence at x = 0 in eq. (7.22) in a way that preserves the shuffle algebra structure.

The relationship between E 4 and Γ
In this section we show that, just like in the cubic case, the elliptic polylogarithms E 4 can be written as linear combinations with constant complex coefficients of the iterated integrals Γ.
We start by showing that the differential one-forms dx ψ ±1 can be written as linear combinations of the one-forms dz g (1) . We use the following identity relating the elliptic function κ in eq. (7.11) and the function g (1) , This relation can be proved in the standard way 12 : Seen as a function of z, the left-hand side has simple poles only for z = ±z c , and the residues at these poles are ±1. The righthand side is periodic, despite individual terms not being periodic. Hence, the difference of the left-and right-hand sides defines an elliptic function without poles, and must thus be constant. The constant is zero, because both sides vanish for z = z * . Using eq. (7.31), we immediately obtain for c = ∞, The corresponding relation for dx ψ −1 (∞, x) is given in eq. (7.21), while dx ψ 1 (∞, x) requires one to know what Z 4 corresponds to when seen as a function on the torus. Using eq. (7.20) and the definition of the Weierstrass zeta function in eq. (4.10), we find dx ψ 1 (∞, x) = −dz g (1) (z − z * ) + g (1) (z + z * ) . (7.33) Let us now analyse the one-forms dx ψ ±n for n > 1. Using eq. (7.31) and (7.33), as well as the explicit definition of κ, we find dx ψ n (∞, x) = dz G n g (1) (z); p 1 (κ(z)), p 2 (κ(z), c 4 κ (z) = dz G n g (1) (z); ℘(z), ℘ (z) = dz g (n) (z) .

(7.34)
Since the remaining one-forms with n > 1 can be obtained by multiplying the functions ψ ±1 by Z (n−1) 4 (x), all other cases can easily be obtained by applying the Fay identity (4.24), just like for the cubic case discussed in section 5. The only difference to the cubic case is the appearance of the terms proportional to Kronecker delta functions for n = 2 in eq. (7.23), which are required to cancel double poles at z = ±z * . For example, we have where in the last step we used eq. (4.21). We see that the double poles contained in the Weierstrass ℘ functions in eq. (4.21) are precisely cancelled by the subtraction term dx = dz κ (z).
To conclude, we see that every one-form dx ψ ±n can be written as a linear combination with constant complex coefficients of one-forms of the form dz g (k) (z ± z c ). Hence, just like in the cubic case, every iterated integral E 4 can be written as a linear combination of elliptic polylogarithms Γ.

Algorithmic integration in the quartic case
In this section we extend the integration algorithm of section 6 to the quartic case. More precisely, let us denote by R 4 ≡ C(x, y)/ y 2 = P 4 (x) the field of rational functions of the elliptic curve E defined by the quartic equation y 2 = P 4 (x), and A 4 is the R 4 -algebra generated by Z 4 (x) and all elliptic polylogarithms E 4 ( n 1 ...n k c 1 ...c k ; x), The total length is defined in the same way as in the cubic case, and it is easy to see that A 4 is filtered by the total length. In the following we present an algorithm which allows us to compute a primitive of every element in A 4 . The algorithm is very similar to the cubic case, so we only highlight here the main similarities and differences. We start by classifying integrals into the four classes defined in eq. (6.1), and for each class we obtain recursion relations using integration by parts. The recursion relations for the types A and B have depth one, and we can reduce every integral in these types to a the integrals A −1 and B c,1 , just like in the cubic case. For integrals of type C, however, the recursion relation has depth four in the quartic case (compared to depth three in the cubic case), and so every integral in this type can be written as a linear combination of C k , k ∈ {−1, 0, 1, 2}. Finally, integrals of type D can be reduced to D c,1 and integrals of type C, and we do not need to consider integrals of type D with c = a i . In the end, we find that every integral can be reduced to the following integrals Comparing these integrals to eq. (6.7), we see that the only difference between the cubic and quartic cases is that we need to include C 2 [X ] into the list of independent integrals. This reflects the fact that x 2 dx/y is related to the abelian differential of the second kind in the quartic case, cf. eq. (7.8).
The independent integrals in eq. (7.37) are in one-to-one correspondance with the integration kernels in eq. (7.23), with the exception of integrals of type C 2 . Equivalently, we may consider the integrals Integrating by parts, we obtain Just like in the cubic case, the integral in the right-hand side does not have lower total length, and so it cannot be done recursively. Indeed, eq. (4.22) and (7.26) imply Inserting eq. (7.40) into eq. (7.39), we obtain an integral of type D with c = a 1 , and so this integral can be reduced to a linear combination of integrals of type C. We find m 2c 4 a 12 a 13 a 14 dx y(x − a 1 ) Z (m) 4 where the dots indicate terms that have lower total length, or they have the same total length and fall into the classes C 0 and C 1 . Finally, inserting eq. (7.41) into eq. (7.40), we have . . , (7.42) where the right-hand side only involves integrals with lower total length or from the class C 0 and C 1 . All the integrals in the right-hand side have lower complexity and can be done recursively. This concludes the proof that every element in A 4 has a primitive that can be computed in an algorithmic way. Just like in the cubic case, this algorithm generalises the classical algorithm to evaluate integrals of rational functions on elliptic curves in terms of elliptic integrals.

Examples
In this section we discuss some examples of how to work with elliptic polylogarithms. The examples are of mostly mathematical nature. For an application to the sunrise integral, we refer to ref. [40].

MPLs depending on square roots of cubic or quartic polynomials
We have already seen that ordinary MPLs are special cases of elliptic polylogarithms, cf. eqns. (3.20) and (7.30). In this section we show that there are other classes of MPLs that can be expressed in terms of elliptic polylogarithms in a natural way.
Consider an MPL of the form G (a 1 (x, y), . . . , a n (x, y); a n+1 (x, y)), where a i (x, y) is a rational function subject to the constraint y 2 = P N (x), N = 3 or 4, In other words, we consider MPLs whose arguments are rational functions on the elliptic curve E defined by the equation y 2 = P N (x). In the following we only discuss the cubic case, N = 3, and the quartic case is similar. Our goal is to show that every MPL of this type can be expressed in a natural way in terms of the elliptic polylogarithms E 3 . The argument proceeds by induction in the weight n of the MPL. We start by discussing the case n = 1. Differentiating with respect to x, we find ∂ x G(a 1 (x, y); a 2 (x, y)) = 1 a 2 (x, y) − a 1 (x, y) a 2 (x, y) − a 2 (x, y) a 1 (x, y) a 1 (x, y) , y) is the (total) derivative with respect to x. The right-hand side of eq. (8.1) is obviously a rational function on E, and so its primitive can be expressed in terms of Z 3 (x) and E 3 ( n c ; x). This concludes the proof that G(a 1 (x, y); a 2 (x, y)) can be expressed in terms of elliptic polylogarithms (and Z 3 ). The case of weight n > 1 then follows by induction. Assume that the claim is true for MPLs up to weight n − 1. We can then differentiate G (a 1 (x, y), . . . , a n (x, y); a n+1 (x, y)) with respect to x, and since differentiation lowers the weight of MPLs, we know by induction that the derivative lies in A 3 . Since every element in A 3 has a primitive in A 3 , we conclude that G (a 1 (x, y), . . . , a n (x, y); a n+1 (x, y)) lies in A 3 , i.e., it can be expressed in terms of elliptic polylogarithms.
In the remainder of this section we illustrate the previous result on some simple examples. Let us consider the following function of weight one, Differentiating with respect to x, we find where b i denote the roots of the cubic polynomial P 3 (x) − 1. It is possible to obtain explicit algebraic expressions for the roots b i in terms of the branch points a i . The expressions are lengthy and not very illuminating for our purposes, so we do not show them here. We only mention the following useful relations, We see from eq. (8.3) that the derivative of f is a rational function on the elliptic curve E, and so f itself can be written in terms of elliptic polylogarithms, where c may depend on the branch points a i , but it is independent of x. The value of c is determined from the fact that the elliptic polylogarithms vanish for x = 0 while f does not. We find Before we turn to examples of higher weight, let us make some comments about the previous formula. We see that have managed to write an ordinary logarithm as a non-trivial linear combination of elliptic polylogarithms (non-trivial in the sense that the elliptic polylogarithms do not simply reduce to ordinary MPLs due to eq. (3.20)). Using eq. (3.32), we know that all the elliptic polylogarithms in the right-hand side of eq. (8.6) can be written as an incomplete ellipitic integral of the third kind, leading to an intriguing relation connecting a logarithm involving the square root of a cubic polynomial and incomplete elliptic integrals. Alternatively, we also know that every E 3 function can be written as a linear combination of iterated integrals Γ on the torus. If z a is such that (µ(z a ), c 4 µ (z a )) = (a, y a ), eq. (5.17) gives where in the last step we used the fact that f (a i ) = 0. Note that we can replace Γ by their periodic analogues Γ in the previous equation, cf. eq. (5.21).
Finally, since f (a i ) = 0, we see that eq. (8.6) implies a linear relation among E 3 −1 b i ; a j , i, j ∈ {1, 2, 3}. Similar relations of this type have appeared in ref. [14], see eq. (7.7) therein, where they have shown up in the context of the two-loop sunrise integral. While these relations were rather mysterious in ref. [14], the analysis of this section reveals their origin: Ordinary logarithms with cubic or quartic roots inside their arguments can always be expressed in terms of elliptic polylogarithms. At some special points where the logarithm vanishes, this implies a linear relation among elliptic polylogarithms evaluated at those points. We stress that these special linear relations do not contradict the linear independence of the integration kernels ϕ −1 (b i , x). Indeed, the functions E 3 Let us now illustrate how the previous story generalises to higher weights. As an example, we consider the functions The decomposition into even and odd parts is not essential to the discussion, but it makes some of the formulas more compact. We only discuss the odd combination f − in detail, because the even case is very similar. Differentiating with respect to x, we find where in the last step we have used eqns. (3.20) and (8.6). As expected, the derivative takes values in the algebra A 3 , and so it admits a primitive inside the same space, i.e., f − can be expressed in terms of elliptic polylogarithms. We find Similarly, we find We see that f ± can always be expressed in terms of elliptic polylogarithms and ordinary MPLs. Using the results of section 5 we can also write the previous relations in terms of the functions Γ or Γ. Let us conclude this section with some comments about the relevance of the results in this section for the computation of Feynman integrals. The results of this section show that it is possible to find combinations of genuine elliptic polylogarithms that evaluate to ordinary MPLs with square-root arguments. These examples show that the distinction between ordinary and elliptic MPLs may not be as clear-cut as sometimes assumed in the physics literature, and some care may be needed when making statements about when a given Feynman integral or amplitude can or cannot be expressed in terms of ordinary MPLs alone. For example, it could be conceivable that an amplitude is expressed as a sum of integrals which individually evaluate to elliptic polylogarithms that cannot be written as ordinary MPLs, but their sum only involves MPLs, through a formula similar to eq. (8.10) or (8.11). where n i and α i are integers, and for concreteness we assume 0 < z < 1. This class of integrals is tightly connected to Gauss' hypergeometric function,

Hypergeometric
We are interested in the Laurent expansion in of the integral. In the case where the exponents in the integrand are integers for = 0, the Laurent coefficients are MPLs that can be computed explicitly using standard techniques. In the case where all the exponents are half-integers, however, the Laurent coefficients are not known in the literature. In the remainder of this section we show that the Laurent coefficients of T can be expressed in terms of elliptic polylogarithms E 3 .
Using integration by parts, one can show that every integral in the family defined by eq. (8.12) can be written as a linear combination of two master integrals. We choose the following basis integrals, Let us discuss the computation of T 1 . After expansion in , the integrand involves powers of logarithms that can be recast in the form of MPLs. Using eq. (3.20), we see that all integrals can be reduced to integrals of the type with c i ∈ {0, 1, λ}. For the first few orders, we find explicitly We see that eq. (8.16) involves at every order in only functions of uniform weight. More precisely, all the terms in the coefficient of k have weight k (we recall that the weight of E 3 ( n 1 ... n k c 1 ... c k ; x) is not just the number of integrations, but it is defined as |n 1 | + . . . + |n k |).
Let us now discuss the computation of the integral T 2 . Using the algorithm described in section 6, we find with where the ε 3 denote regularised values at the upper integration limit. These arise because individual terms may diverge logarithmically at 1, e.g., if δ → 0, we have The finite terms are explicitly given by , .
The logarithmic singularities all cancel in the final expression for T 2 , leaving a finite result. More details on the regularisation can be found in appendix A. Let us make some comments about eq. (8.18). First, we see that the result involves elliptic polylogarithms of the form E 3 ( 2 ... * ... ; x). We thus see the necessity for the integration kernels ϕ ±n with n > 1. Second, let us comment on the weight of the Laurent coefficients in eq. (8.18). Since ω 1 = E 3 ( 0 0 ; 1) has weight zero, we see that all the terms in eq. (8.18) have uniform weight 13 . We have checked that this observation remains true for the first six terms in the Laurent expansion. Furthermore, since the periods ω i have weight zero, the Legendre relation (3.14) implies that we should assign weight one to the quasi-periods η i (and thus also to Z 3 ). We observe that with this assignment of the weight, all the terms in the right-hand side of eq. (8.17) have uniform weight once the overall prefactor is scaled out. We emphasise that in this context it is important that the concept of weight is not associated to the number of integrations!

Appell F 1 functions that evaluate to elliptic polylogarithms
In this section we study two different classes of Appell F 1 functions that evaluate to elliptic polylogarithms. The Appell F 1 function admits the integral representation We assume that (a, b 1 , b 2 , c), and thus the exponents in the integrand, depend linearly on . In the following we show that in the case where for = 0 three or more of the exponents are half-integers, the coefficients appearing in the Laurent expansion around = 0 can be expressed in terms of elliptic polylogarithms.

The cubic case
We consider the following family of integrals, A(n 1 , n 3 , n 3 , n 4 ; z 1 , z 2 ) where n i and α i are integers, and for concreteness we assume 0 < z i < 1. These integrals are closely related to the Appell F 1 function defined in eq. (8.21). Using integration by parts, we can write every integral in this family as a linear combination of the following three master integrals with λ i ≡ 1/z i and y 2 = x(x − 1)(x − λ 1 ).
The computation of A 1 and A 2 is completely analogous to the case of the 2 F 1 function. For the first master integral we find, where T 1 is given in eq. (8.16). The second master integral is given by with The third master integral has an additional logarithmic singularity at x = λ 2 , and all the integrations can be performed using the formula

28)
We observe that, just like in the case of the 2 F 1 function, the Laurent coefficients of the master integrals A i have uniform weight.

The quartic case
We now consider the following family of integrals, which corresponds to the case where all the exponents in the integrand in eq. (8.21) evaluate to half-integers for = 0, B(n 1 , n 3 , n 3 , n 4 ; z 1 , z 2 ) (8.30) where n i and α i are integers, and we assume 0 < z 2 < z 1 < 1. Just like in the cubic case, there are three master integrals for this family, which we choose as . Since λ 2 > λ 1 > 1, our convention in eq. (7.2) implies Im y < 0 for 0 < x < 1. All the integrations can easily be done order by order in in terms of elliptic polylogarithms. For the first master integral, we find For the second master integral, we find B 2 (z 1 , z 2 ) = i √ 1 − z 2 2z 1 z 2 4 η 1 ω 1 B 1 (z 1 , z 2 ) + 1 1 + (α 1 + α 2 + α 3 + α 4 ) B 2 (z 1 , z 2 ) , (8.34) with B 2 (z 1 , z 2 ) = − 2πi ω 1 + α 1 E 4 ( 2 0 ; 1) + α 2 ε 4 ( 2 1 ) + α 3 E 4 Finally, the third master, which is related to an abelian differential with a simple pole at infinity, is given by We observe again that all master integrals are uniform in weight order by order in the expansion.

Conclusion
In this paper we have introduced a class of iterated integrals on elliptic curves that have at most logarithmic singularities, and therefore deserve to be called elliptic polylogarithms.
The key idea is that the kernels that define the iterated integrals are obtained by analysing the abelian differentials of the first, second and third kinds on an elliptic curve, requiring the introduction of root-valued integration kernels. This idea by itself is not new, and iterated integrations over kernels involving square roots have been considered before, cf., e.g., refs. [44][45][46]. The main difference between the functions in the literature and our elliptic polylogarithms is that we insist on having integration kernels with at most simple poles, leading to iterated integrals with at most logarithmic singularities on the elliptic curve. For this reason, we do not include abelian differentials of the second kind into our basis of integration kernels (which by definition have poles of higher order), and we are forced to consider integration kernels over transcendental functions, in particular incomplete elliptic integrals of the second kind. The closure of the algebra of integration kernels then forces us to consider two infinite towers of integration kernels for each point on the elliptic curve (one of the two towers is absent for the branch points). A large part of our paper was devoted to studying some of the properties of elliptic polylogarithms. First, we have clarified how the iterated integrals of root-valued elliptic kernels are related to the multiple elliptic polylogarithms considered in the mathematics literature [29,39,42]. We have shown that, under the isomorphism which identifies an elliptic curve with a torus, our integration kernels are mapped to simple linear combinations of the integration kernels introduced in refs. [29,39,42] (up to the distinction that we prefer to work with meromorphic rather than periodic functions). Since the kernels of refs. [29,39,42] are known to be independent, this proves at the same time the independence of our integration kernels. Second, we have provided an explicit algorithm to compute primitives of elliptic polylogarithms multiplied by rational functions. This algorithm generalises to elliptic polylogarithms the classical algorithm to compute primitives of rational functions on an elliptic curve in terms of the elliptic integrals of the first, second and third kinds. Finally, we have applied our results to compute the -expansion of certain classes of 2 F 1 and Appell F 1 function that cannot be expressed in terms of ordinary MPLs. We observe in all cases that it is possible to choose the master integrals in such a way that the results have uniform weight order by order in the -expansion, hinting at an extension of the concept of "pure function" well-known from ordinary polylogarithms. We emphasise that in order for this to be true, it is important that the weight is not identified with the number of integrations.
Our elliptic polylogarithms are very flexible and incorporate large classes of other special functions. We give here a list of functions which can be expressed in terms of our elliptic polylogarithms: • Ordinary MPLs evaluated at arguments that are rational functions on the elliptic curve, cf. eq. (3.20) and section 8.1.
• It was shown [47] that some instances of the ELi functions of refs. [17][18][19] can be expressed in terms of the elliptic polylogarithms of refs. [29,39,42]. As a consequence, these functions can also be written in terms of our elliptic polylogarithms.
• The elliptic generalisations of polylogarithms of ref. [34] are a special case of the functions considered here.
Since many of these functions have appeared in physics computations, we foresee that the elliptic polylogarithms introduced in this paper will play a prominent role in applications to Feynman integrals. In a companion paper [40], we have applied our formalism, in particular the integration algorithm introduced in section 6, to the computation of the sunrise integral with three different masses in two space-time dimensions. It will be exciting to see for which other Feynman integrals our elliptic polylogarithms provide the natural language. We conclude this paper by commenting on some limitations and directions for future research. Throughout this paper, we have assumed that the branch points a i that define the elliptic curve are held constant. In applications, however, it is often the case that the branch points are themselves dynamical variables, and one wishes to differentiate and/or integrate with respect to the branch points. This situation is typical for Feynman integrals, where one usually invokes differential equations to obtain analytic results. The latter lead to iterated integrals in the branch points of the elliptic curve, or equivalently the modular parameter τ that defines the lattice, cf., e.g., [17][18][19][20][21]25]. The functions introduced in this paper are not adapted to such a scenario. It is known, however, that there is a close relationship between the iterated integrals of refs. [29,39,42] and iterated integrals in the modular parameter τ , in particular iterated integrals over modular forms, cf. refs. [42,[48][49][50][51][52]. The latter are known to show up also in Feynman integral computations [21]. It will be important to extend our framework to include differentiation and integration with respect to the branch points, and it will be fascinating to understand in detail the connection between the functions introduced in this paper and the theory of iterated integrals over modular forms. in eq. (2.3)). While easy to implement in practise, this regularisation procedure is rather subtle, because it is a priori not clear that a deformation of the contour will preserve the shuffle algebra properties of MPLs, see eq. (2.6). Indeed, the derivation of eq. (2.6) relies on the fact that both factors in the left-hand side are iterated integrals over the same contour, and so it is not clear that the prescription in eq. (2.5) preserves the shuffle product. In this appendix we review a general procedure to regularise iterated integrals with logarithmic singularities in a way that preserves the shuffle algebra structure, and we show that this procedure reduces to the usual prescription for MPLs in eq. (2.5).

A.1 Shuffle regularisation
Consider a function f (x) that has at most logarithmic singularities. For every point x 0 ∈ C there is a neighbourhood U such that for every x ∈ U we can write where n is a positive integer and the functions f k are holomorphic on U . The regularised value of f at x 0 is defined by (see, e.g., ref. [53]) The map Reg x 0 has the following properties: 1. If f has at most logarithmic singularities, then f 0 is regular at x = x 0 .
2. If f is regular at x = x 0 , then Reg x 0 f (x) = f (x), i.e., the regularised value of a regular quantity does not change.
3. The map Reg x 0 is an algebra homomorphism, i.e., it preserves the multiplication, Reg x 0 (f (x) g(x)) = Reg x 0 f (x) Reg x 0 g(x) . It is usually assumed in the literature that one is talking about the regularised version, and and the map Reg 0 is usually not written explicitly. Equation (A.3) ensures that if we replace all polylogarithms by their values regularised at the origin, the shuffle-algebra structure is preserved. While the previous discussion allows one to define a regularised version of iterated integrals that is compatible with the shuffle algebra structure, it is less well suited for explicit computations, because the regularisation procedure requires one to first work with an arbitrary base-point δ and then to expand around δ = 0 and to remove the logarithmic divergences. In the remainder of this section we show that for ordinary MPLs the effect of the map Reg 0 is equivalent to the well-known prescription in eq. (2.5), and we then propose an extension of this prescription to elliptic polylogarithms.