An Elliptic Generalization of Multiple Polylogarithms

We introduce a class of functions which constitutes an obvious elliptic generalization of multiple polylogarithms. A subset of these functions appears naturally in the \epsilon-expansion of the imaginary part of the two-loop massive sunrise graph. Building upon the well known properties of multiple polylogarithms, we associate a concept of weight to these functions and show that this weight can be lowered by the action of a suitable differential operator. We then show how properties and relations among these functions can be studied bottom-up starting from lower weights.


Introduction
The generalized polylogarithms [1][2][3][4] (also called Goncharov functions) are of common use in the evaluation of Feynman graph amplitudes, especially in the differential equation approach. As it is well known, however, they are not enough to span the full set of functions required to evaluate two-loop Feynman integrals. The obvious obstruction comes from Feynman graphs that fulfil irreducible second-(or higher-) order differential equations, of which the most notable example is indeed the massive two-loop sunrise graph. In spite of the long efforts and the vast literature produced on the subject [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23], the way to further generalize the polylogarithms to accomodate this case remains topic of discussion, with a fascinating crosstalk between particle physics, mathematics and string theory, see for example [24,25].
In this paper we introduce and discuss an elliptic generalization of multiple polylogarithms, referred to in this Introduction with the name EG [n] for short, (more refined notations will be used in the next Sections), defined starting from an integral representation of the form where R 4 (u, b) is the fourth order polynomial in b is any pair of the 4 roots of R 4 (u, b), namely b 1 = 0, b 2 = 4m 2 , b 3 = (W − m) 2 and b 4 = (W + m) 2 and g [n] (u, b) is a generalized polylogarithm in b of degree n with "alphabet" corresponding to the above 4 roots b i (for simplicity we consider mainly u real and in the range 9m 2 < u < ∞, but the continuation to other values of u is almost obvious). The EG [n] (k, u) are generalizations of the integrals The latter can be expressed in terms of two independent master integrals, say I 0 (u) and I 2 (u), which are simple suitable linear combinations of the I(k, u) and correspond to K(x) and E(x), the complete elliptic integrals of first and second kind respectively (hence the elliptic terminology for the new functions); one has for example I 0 (u) = 2 Moreover, up to an inessential numerical factor, I 0 (u) is the phase space of three particles of equal mass m at energy √ u in d = 2 dimensions, (see Section 7) while the newly introduced EG [n] (k, u) are obvious generalizations of the terms which arise when expanding the d-dimensional 3-body phase space in powers of (d − 2). A subset of the functions EG [n] (k, u) is therefore naturally appearing in the (d − 2)-expansion of the imaginary part of the two-loop massive sunrise graph.
I 0 (u) and I 2 (u) satisfy a homogeneous two by two system of linear differential equations in u (see Eq. (2.39)), which can be written as d du I 0 (u) I 2 (u) = B 0,0 (u) B 0,2 (u) B 2,0 (u) B 2,2 (u) where the matrix elements of the 2 × 2 matrix B(u), given in Eq.s (2.21,2.22), contain rational coefficients and poles at u equal to 0, m 2 , 9m 2 . At variance with the I k (u), the study of the functions EG [n] (k, u), for every value of n, requires the introduction of three master integrals instead of two. We consider then three new functions EG [n] k (u), k = 0, 1, 2, which are again simple linear combinations of the above EG [n] (k, u), (see for instance Eq. (2.18)), and we find that they satisfy an inhomogeneous system of differential equations of the form d du EG [n] 0 (u) EG where the matrix B(u) is the same as in the previous homogeneous equation for the I k (u), while the coefficients R [j] k,k ′ (u) of the inhomogeneous terms consist again in general of rational expressions in u with poles at u equal to 0, m 2 , 9m 2 . As we will see, EG [0] k (u) = I k (u), and one finds in particular that EG 1 (u) is a constant, which explains why for n = 0 only two independent functions are needed instead of three.
Note the presence, in the r.h.s. of the above equations, of functions of the same family EG [n] k (u), but with lower values of the index n. In particular, for any given n, we find a two by two system of coupled differential equations, plus a third, simpler, decoupled linear differential equation. This suggests indeed that we can tentatively associate a weight n to the functions EG [n] k (u) with respect to the action of a three by three matricial operator. The latter, though, clearly factorises into the two by two operator (−B(u)+d/du), which directly lowers the weight of EG [n] 0 (u) and EG [n] 2 (u), and the simple first order differential operator, d/du, which lowers the weight of EG 1 (u) (similarly to what happens with the Goncharov polylogarithms). Such a generalized weight will be called E-weight and the functions EG [n] k=0,1,2 (u) also referred to as Epolylogarithms.
As the pair of functions I k (u), in particular, is annihilated by the operator (−B(u) + d/du), the two functions I k (u) can be considered E-polylogarithms of E-weight equal zero. Similarly, at E-weight zero the third function, say I 1 (u), is a constant and is therefore annihilated by d/du.
As we will see, it can be useful to rewrite the homogeneous first order system for I 0 (u), I 2 (u) as a homogeneous second order differential equation for I 0 (u) only; when that is done, one obtains where D(u, d/du), Eq. (2.26), is a suitable second order differential operator. Acting similarly on the functions EG where D(u, d/du) is the same differential operator appearing in the second order equation satisfied by I 0 (u), while the coefficients r k (u) are also rational expressions in u, with poles at u equal to 0, m 2 , 9m 2 . As we can see, in the r.h.s of the above equation we find E-polylogarithms of weight n − 1 and n − 2. We can therefore also say that a function EG [n] 0 (u) satisfying the above equations is an E-polylogarithm of E-weight n under the action of the scalar second order differential operator D(u, d/du). That confirms, of course, that I 0 (u), being annihilated by D(u, d/du) has E-weight equal 0. Alternatively, one could also derive a different second order differential operator, say D 2 (u, d/du), such that and the discussion would apply in the very same way.
In the course of the paper we will also encounter repeated integrations of rational factors times for instance the function I 0 (u); in this picture, they constitute a simple subset of E-polylogarithms, and we will refer to them for simplicity as E 0 -polylogarithms, see Section 3.
The equations for the EG [n] k (u) can further be solved by using the Euler method, which provides representations for the EG [n] k (u) as suitable integrals involving the solutions of the homogeneous equation, i.e. the function I 0 (u) above with the accompanying function J 0 (u), Eq.(2.31), and the inhomogeneous term, providing interesting relations between the EG [n] k (u) and the repeated integrations of products of the I 0 (u), J 0 (u) and the usual (poly)logarithms of u. An example of such relations is (where we have written ln(b/m 2 ), for simplicity, as ln b ). Let us recall, indeed, that one of the musts of an analytic calculation is to discuss as deeply and explicitly as possible the identities which might hold between the various functions introduced in the calculation. This allows one to write the result in a compact form and to understand whether two apparently different formulas are indeed different or equal.
The rest of the paper is organized as follows: in Section 2 we study the (well known) functions I(k, u) and reduce them to three master integrals using integration-by-parts. We then show that one of the three masters is not linearly independent and show how to derive a two by two system of differential equations for the two masters I 0 (u) and I 2 (u), and their accompanying functions J 0 (u) and J 2 (u). In Section 3 we study a first simple class of functions obtained by repeated integrations of products of one of the I k (u) or J k (u) times rational factors. As the complexity of these functions is decreased by differentiation, these functions can be given a simple concept of weight, similar to that of multiple polylogarithms. We call these functions E 0 -polylogarithms and their weight E 0 -weight. Then in Section 4 we study the first example of E-polylogarithm at E-weight one and show how it can be rewritten as the product of a logarithm and the E-weight zero function I 0 (u). Similar relations are derived for all E-weight one functions in Section 5. We extend then our study to higher weights in Section 6 and find explicit relations to simplify E-polylogarithms at E-weight 2. Finally we use our results to give a compact representation of the imaginary part of the twol-loop massive sunrise up to order ǫ 2 in Section 7. We then draw our conclusions and outlook in Section 8.

The beginning
In this Section we will start by recalling some known results, which will be generalized in the rest of this paper. To begin with, for 9m 2 < u < ∞ we consider the (real) function Eq. (2.1) corresponds, up to a multiplicative constant, to the imaginary part of the equal mass sunrise amplitude in d = 2 dimensions. For convenience of later use, let us recall that its value in the u → 9m 2 limit is lim where K(x) is the complete elliptic integral of the first kind. As a first step, following previous works [9], let us derive a (second order, homogeneous) differential equation for I 0 (u). To that aim, we define the (related) functions where (b i , b j ), as above, are any two (different) roots of the polynomial R 4 (u, b) and n is an integer, so that Eq.(2.1) is recovered for b i = 4m 2 , b j = (W − m) 2 and n = 0. (A trivial remark: it is sufficient to consider for b i , b j only the pairs of adjacent roots, as any other choice is a linear combination of them). One has the (obvious) identity by explicitly carrying out the b-derivative and by using the replacement The above identity holds for any n (except n ≤ −1 if b i = 0 or b j = 0); for integer positive n, by using it (recursively, when needed) one can express any I(b i , b j , n, u), for n ≥ 3, in terms of the 3 master integrals Consider now the (auxiliary) quantities with k = 0, 1, 2. By writing (again) R 4 (u, b) = R 4 (u, b)/ R 4 (u, b), and using Eq.s(2.9) one finds Q(k, u) = n=0,1,2 c(k, n, u)I(b i , b j , n, u) , (2.11) where the coefficients c(k, n, u) are (simple) polynomials in u. From it one gets at once But we can obtain the u-derivative of the Q(k, u) by differentiating directly the definition Eq.(2.10), obtaining where Eq.s(2.9) were again used, and the d(k, n, u) are also (simple) polynomials in u. By writing, for a given value of k, that the r.h.s. of Eq.(2.12) is equal to the r.h.s. of Eq.(2.13) one obtains a linear (homogeneous) equation expressing the u-derivatives of the three master integrals I(b i , b j , n, u), with n = 0, 1, 2, in terms of the same three master integrals. One can then take three such equations, corresponding to three different values of k, say k = 0, 1, 2 for definiteness, and solve them for the three derivatives. The result can be written as where , . (2.17) Eq.(2.14) is a linear homogeneous system of three first order differential equations for the three master integrals I(b i , b j , k, u), k = 0, 1, 2. Quite in general, a three by three first order system is equivalent to a third order differential equation for one of the three functions, say for instance I(b i , b j , 0, u); but we are looking for a second order equation. It is indeed known (see Appendix C) that one of the equations can be decoupled from the other two. To that aim, we introduce a new basis of master integrals according to the definitions In terms of the functions of the new basis the system splits into a very simple equation involving only a result already noted by A.Sabry in his (1962) paper [26] (see also Appendix C), and in a two by two first order system for the other two functions i.e. the system decouples into a (rather simple!) equation for the function I 1 (b i , b j , u) and a two by two first order homogeneous system for the two functions does not appear anymore. Let us just recall that Eq.(2.19) implies that I 1 (b i , b j , u) is constant, with the value of the constant depending on the actual choice of roots (b i , b j ) (see for instance Eq.s(2.35) below). The two by two system can be recast in the form of a single second order homogeneous differential equation for one of the two functions, say I 0 (b i , b j , u); to that aim, we rewrite the first of the Eq.s(2.20) as where D 1 (u, d/du) is the first order differential operator and then evaluate the u-derivative of that same first equation of (2.20). By expressing in the result the derivative of I 2 (b i , b j , u) through the second of the Eq.s(2.20) and then where D(u, d/du) is the second order differential operator . (2.26) Quite in general, the two by two first order differential system in Eq.(2.20) where I 0 (u) is the function already introduced in Eq.(2.1). For obtaining a second solution of the same equation, let us write where we have changed R 4 (u, b) into −R 4 (u, b) to keep the solution real. It is obvious that all the homogeneous relations valid for the generic functions I(b i , b j , n, u) apply as well to the functions I(0, 4m 2 , n, u), as they are equal to the corresponding I(b i , b j , n, u) of Eq.s(2.7), times an overall imaginary factor i. The function J 0 (u), defined as is therefore a second solution of Eq.(2.25). The u → 9m 2 limit of the above function (see for instance Eq.s(8.12) of [16]) is showing in particular, for comparison with Eq.(2.4), that the two functions I 0 (u) and J 0 (u) are linearly independent. An explicit calculation gives also For completeness, we recall (from [16]) also the values of I 1 (u), J 1 (u) (which are constant) As a further remark, given any solution of the second order equation Eq.(2.28) corresponding to the two by two system (2.27), we can complete the pair of solutions by using Eq.(2.24); if I 0 (u), J 0 (u) are the solutions of the second order equation, the accompanying function I 2 (u), J 2 (u) are then given by so that the two independent pairs of solutions of Eq.(2.27) are given by the two columns of the matrix which therefore satisfy (in matricial form) the equations d du For convenience of later use, let us observe here that we have also, according to Eq.s(2.18) and Eq.s(2.35) Let us repeat here that, as anticipated in the Introduction, due to Eq.s (2.36) I k (u) and J k (u) have E-weight equal to zero.
In the range 9m 2 < u < ∞ the two functions I 0 (u), J 0 (u) are real, outside that range they develop also an imaginary part and become complex; the details of their analytic continuation can be found, although with a slightly different notation, in Appendix B of [16].
We can now look back at the three by three system Eq.(2.14). It has three linearly independent solutions, each solution being a set of three functions, namely the three sets I(0, 4m 2 , n, u), I(4m 2 , (W − m) 2 , n, u) and I((W − m) 2 , (W + m) 2 , n, u) with n = 0, 1, 2. With the change of basis of Eq.s(2.18), the two sets I(0, 4m 2 , n, u) and I(4m 2 , (W − m) 2 , n, u) correspond to the decoupled sets J n (u), I n (u) just discussed.
Concerning the third set, an explicit calculation (based on contour integration arguments in the complex plane, see for instance [16]) gives The transformations Eq.s(2.18) then read Note that the pair K 0 (u), K 2 (u), being a solution of the two by two system (2.27), must be a linear combination of the two already discussed pairs of solutions, I 0 (u), I 2 (u) and J 0 (u), J 2 (u), and if fact they are just equal to J 0 (u) and J 2 (u), but K 1 (u) differs from J 1 (u). I k (u), J k (u) and K k (u) for k = 0, 1, 2 are therefore indeed the entries of a 3 × 3 matrix of homogeneous solutions of the 3 × 3 system of differential equations where we used of K 0 (u) = J 0 (u), K 2 (u) = J 2 (u), I 1 (u) = 0, J 1 (u) = −π/3 and K 1 (u) = 2π/3. Besides the homogeneous equations Eq.(2.43) we will consider also the corresponding inhomogeneous equations, namely in matrix form where N 0 (u), N 1 (u), N 2 (u) are the inhomogeneous terms, supposedly known, and the functions g 0 (u), g 1 (u), g 2 (u) are the unknown. The system is equivalent an inhomogeneous second order equation for g 0 (u) and a first order differential equation for g 1 (u), with g 2 (u) given by where D 1 (u, d/du) is given by Eq. (2.24), and N (u) is related to N 0 (u), N 2 (u) of Eq.(2.45) by the relation The solution of Eq.s(2.45) can be obtained by the Euler-Lagrange method; to that aim, given the decoupled form of the system, one can split the problem in two steps. First, one solves the rather trivial first order inhomogeneous differential equation for g 1 (u) by quadrature obtaining where c 1 is an integration constant. Then, in order to solve the two by two coupled system, one considers the two by two matrix of the two independent solutions already introduced in Eq.(2.37), and its determinant, the Wronskian of the system, defined as From the very definition, it satisfies the equation where use is made of Eq.s (2.22), showing that W s (u) is a constant; an explicit calculation gives indeed [16], The inverse of the matrix (2.51) is therefore The Euler-Lagrange method then gives the solutions of the Eq.s(2.45) in the form where c 0 , c 2 are two more integration constants. Let us note that the above formula can be derived by considering the following first order derivatives where B 00 (u) + B 22 (u) = 0 was used. By quadrature one obtains, up to constants, which satisfies the first order homogeneous equation (2.60) Its solution, with I 0 (u), J 0 (u) given by Eq.s(2.1,2.31), is .
The solution of Eq.(2.46) then reads where the two integration constants c 0 , c 2 , to be fixed by the boundary conditions, are the same as in Eq.(2.56); g 2 (u) is then given by Eq.(2.48).
3 Repeated integrations of I 0 (u), J 0 (u) and rational factors In the previous Section, we have introduced the pairs of functions I 0 (u), I 2 (u) and J 0 (u), J 2 (u), and shown their use in writing the solution Eqs.(2.56,2.62) of the inhomogeneous equations Eqs.(2.45,2.46). As it is easy to imagine, the integration of products of those functions times rational factors appears even in the simplest cases. Therefore, before studying the more general E-polylogarithms, we consider now the properties of such (possibly repeated) integrations, discussing the analogy with the ordinary generalized polylogarithms [1][2][3][4], also called Goncharov functions, of common use in the evaluation of Feynman graph amplitudes. The Goncharov functions can be defined as where the parameters p i vary within a given finite set of values, proper of the problem under study. The repeated integrations arise naturally when solving iteratively the differential equations by the Euler approach (i.e. evaluating first the solution of the homogeneous equation and then accounting of the inhomogeneous term with the variation of the constants method). The superscript n is called the degree (or polylogarithmic weight) of the function. In the context of this paper, we will refer to this weight as G-weight for obvious reasons. By construction, these functions satisfy the relation

2)
i.e. the derivative of a function of G-weight n is a function of the same family but of lower weight n − 1 (times a rational factor). For completeness, we can define also which satisfies, obviously, the equation d du such that a function of G-weight equal to zero is annihilated by the first order differential operator d/du. By following as much as possible Eq.(3.1), we start considering the functions defined by repeated integrations for integer n > 0 as follows where the index k takes the two values k = 0 and k = 2, with the p i taking any of the values of the set {0, m 2 , 9m 2 }. Clearly, for n > 0 these functions behave very similarly to the G-functions under differentiation, such that one would be tempted to associate to them a G-weight in the same way. Nevertheless, as already noted, for n = 0, one defines k (.., u) cannot be naïvely extended to n = 0, defining for instance I and therefore I For this reason, without any claim of rigour or completeness, we call the weight of these functions E 0 -weight, in order to clearly distinguish it from the standard polylogarithmic G-weight, but also from the more general E-weight of E-polylogarithms.
The first of Eq.s(3.5) might also be written, recursively, as from which one has at once, for n > 1, d du I and the procedure can be used recursively, down to n = 2. For n = 1, however, the integrations by parts involve the u derivatives of I k (u), which are non zero and can instead be expressed in terms of the same functions times a combination of the same rational factors, see Eq.s (2.20). The direct, naïve integration-by-parts approach is therefore not sufficient in the case of the very first integration involving I k (v) or J k (v); indeed, one has rather to write the complete system of integration by parts identities obtained by considering the products of all the powers of the rational factors times the functions I k (v) or J k (v), and then to solve the system in terms of the master integrals of the problem. The generic identity has the (obvious) form where X(v) stands for the products of all the possible factors {1, v n , 1/v n , 1/(v − m 2 ) n , 1/(v − 9m 2 ) n } times I k (u) or J k (u), and n is any positive integer.
As a result, it turns out that all the integral of the form where X(v) was just defined above, including both I k (v) and J k (v), can be expressed in terms of the four master integrals which involve only I 0 (v), plus terms in I k (u) (not integrated) generated by the integration by parts.
A few examples (written as relations among indefinite integrals, i.e. valid up to a constant) are For the analytical expression of the master integrals of Eq.(3.9) we refer to Appendix A. So far we have considered repeated integrations associated to the pair of functions I 0 (u), I 2 (u); obviously, the procedure applies as well to the other pair of functions, J 0 (u), J 2 (u), which satisfies the same homogeneous equations as I 0 (u), I 2 (u). While the equations (3.10) (defined up to a constant) remain valid under the exchange of the two pairs of functions, the explicit expression of the four master integrals, corresponding to Eq.s (A.1,A.5,A.7,A.9) is of course different.
As a final remark for this Section, consider a function GI n (u) of the form where G [n] (u) is a function of either G-or E 0 -weight n in the sense defined above, (i.e. obtained by n repeated integrations over rational functions) and I 0 (u) is once more the function of Eq.(2.1) 3 ; let us further recall that I 0 (u), J 0 (u) do not possess definite G-or E 0 -weight, so that G [n] (u) cannot be I 0 (u) or a product of I 0 (u) and J 0 (u). Recalling Eq.s(2.34) an elementary calculation gives where the r

A first example of an E-polylogarithm
Having discussed in detail the properties of the functions I 0 (u), J 0 (u) and of (naïve) iterative integrations over the latter with rational functions, we are now ready to consider the main topic of this paper. Let us start with an explicit example, namely the function EI [1] (where, for simplicity, instead of ln(b/m 2 ) we have written ln b), whose value at u = 9m 2 is as can be easily checked by using the change of variable (2.5). Let us shortly comment the somewhat clumsy notation used; in the name EI [1] 0 (0, u), EI stands for Elliptic integral corresponding to the integration range 4m 2 < b < (W − m) 2 , associated to the functions I k (u), the superscript [1] refers to the weight of the (poly)logarithm G(0, b) = ln b, the arguments (0, u) refer to the "letter" 0 of the (poly)logarithm and (obviously) to the variable u, finally the lower index 0 is the analog of the index k = 0 of I 0 (u) in Eq.(2.1). In this notation, one would have We will work out this example in detail and outline how this generalizes then to higher weights. We can start by deriving a second order differential equation for EI [1] 0 (0, u), by following closely the derivation discussed in Section 2. In analogy with Eq.(2.7) we introduce the auxiliary functions such that clearly EI  Working out the algebra, we are left with an equation, corresponding to Eq.(2.9), whose l.h.s. is the l.h.s. of Eq.(2.9) with the functions I(b i , b j , n, u) replaced by Il(b i , b j , n, u), while the r.h.s. is no longer vanishing, but contains a combination of the I(b i , b j , n, u), due the b-derivative of ln b in Eq.(4.5). That equation can be used for expressing any Il(b i , b j , n, u), with n integer and n > 2, in terms of the three master integrals Il(b i , b j , k, u), k = 0, 1, 2; the homogeneous part of the relations, i.e. the part containing the Il(b i , b j , n, u), has the same coefficients appearing in Eq.(2.9), but in the case of the Il(b i , b j , n, u) there are also inhomogeneous terms, i.e. terms containing not the Il(b i , b j , n, u) but the I(b i , b j , n, u).
We can continue by introducing, in analogy with Eq.(2.10), the auxiliary quantities differentiating them with respect to u etc., we arrive, in analogy to Eq.(2.14), to the following three by three linear system of first order differential equations: where the coeffcients of the homogeneous part, the C nk (u) are the same as in Eq.(2.14), while the C nk (u) are new, similar coefficients (which we do not write here for brevity) multiplying the I(b i , b j , n, u).
Following Eq.s(2.18), we introduce a new basis of master integrals with the definitions In terms of the functions of the new basis and of the I n (b i , b j , u) the system becomes and d du .  nk (u) (which have a similar structure and are not written here again for brevity) are the coefficients of the inhomogeneous terms containing the I n (b i , b j , u). Again, the two by two system can be recast in the form of a single second order homogeneous differential equation for Il 0 (b i , b j , u); the result can be written as

u) is given by
The differential operators D(u, d/du), D 1 (u, d/du) in the two above equations are of course the same as those defined in Eq.s(2.26,2.24).
We can now specialize the formulas to the case b 1 = 4m 2 , b 2 = (W − m) 2 with W 2 = u > 9m 2 , i.e., in the notation of Eq.(4.1), By recalling also Eq.s(2.29) and (2.36), we find finally that Eq.(4.10) becomes where, according to the definitions Eq.(4.3), we can write in the r.h.s., instead of I k (u), EI  (4.14) and W (u) is the wronskian given in Eq. (2.61), whose value we remind here Substituting explicitly the value of the Wronskian and the result at weight zero we are left with EI [1] 0 (0, u) = c (1) where we introduced the compact notation We need therefore to understand integrals of the form u dv 1 ; v n ; 1 v n ; Not all these integrals are linearly independent, as we will we show now by using integration by parts identities. In order to see this, let us define the other function 18) such that, in the notation of (2.20), By using Eq. (2.54) it is easy to see that such that, by choosing to re-express F 0,2 (u, v)I 0 (v) in terms of I 2 (v)F 0,0 (u, v), we see we should generate all integration by parts identities of the form where the X j (u) are appropriate boundary terms; note that, for simplicity, we write the IBPs as relations among primitives, i.e. without specifying the lower integration boundary. This means that all relations we provide here are given up to boundary terms. By proceeding similarly to the general algorithm described in [27], we generate a large number of identities for different numerical values of the powers n and solve the system of equations. We find in this way that all integrals can be expressed in terms of 6 master integrals, which we choose as follows plus simpler terms, i.e. terms which do not require integrating over the functions F 0,0 (u, v) and F 0,2 (u, v).
In particular, we find that one of the integrals in Eq. (4.15) can be re-expressed as linear combination of the other three as follows where we see the appearance of a simpler integral, which reminds of the shuffle identities for polylogarithms. We stress again, that these relations are given up to boundary terms. By using this identity in Eq. (4.15) we find at once EI [1] where the second line is obtained fixing properly the boundary conditions. It is very interesting to notice that all occurrences of integrals over elliptic integrals have cancelled out leaving space to a simple product of a logarithm and an elliptic integral.

Derivation of all the relations at weight one
One might wonder whether the relation above is an accident or if, instead, such relations are more general. It is not difficult to repeat the same exercise (i.e. deriving a second order differential equation, solving it, using integration by parts and fixing the boundary conditions) for all the other weight-one possibilities. Nevertheless, we find it more illuminating to follow a different (but of course equivalent) approach.
As we have seen, the operator D(u, d/du) can be conveniently used to effectively reduce the weight of the E-polylogarithms associated with the functions I 0 (u) and J 0 (u) 4 . Following the example of generalized polylogarithms, we can therefore imagine to study the E-polylogarithms bottom-up, starting from weight one, and applying at each step the operator D(u, d/du) to reduce the complexity to the previous weight, which can be considered as understood.
In order to see how this works, let us look again at the example above. The function EI [1] 0 (0, u) is an E-polylogarithm of weight one. From the discussion at the end of Section 3, and in particular from Eq.(3.11), it is easy to see that similarly also the six functions are E-polylogarithms of weight one. It is then natural to consider the following linear combination where c j are constants. We can now apply the operator D(u, d/du) on the function A(u) and fix the coefficients c j such that By applying the operator D(u, d/du) on each of the terms, we produce terms of weight zero, i.e. combinations of rational functions and I 0 (u), J 0 (u), I 2 (u), J 2 (u). By collecting for the independent terms and requiring the coefficients to be zero we find, as expected This implies of course that and therefore, by Euler variation of the constants EI [1] whereĉ j , j = 1, 2 are two numerical constants. By imposing the boundary conditions at u = 9m 2 (according to Eq.(2.32) J 0 (u) has a logarithm singularity at that point, while all the other terms are regular, so that c 2 = 0), we immediately findĉ 1 =ĉ 2 = 0 , reproducing in this way the result in Eq. (4.23).
In order to complete the exercise, we should remember that at order one we have two more functions to compute, namely EI [1] 1 (0, u) and EI [1] 2 (0, u). Clearly we see that, once EI [1] 0 (0, u) is known, then Eqs.(4.8, 4.11) allow us to compute EI [1] 1 (0, u) and EI [1] 2 (0, u). In particular, EI [1] 2 (0, u) can be obtained from EI [1] 0 (0, u) by simple differentiation, while EI [1] 1 (0, u) fulfils a first order differential equation which can be solved by quadrature. Let us then proceed and compute them. From Eq. (4.11) we find immediately EI [1] and working out the (straightforward) derivatives one finds easily Finally, let us consider EI As lim u→9m 2 EI so that EI [1] 1 (0, u) is given by the quadrature formula We are not able to simplify this expression further as we saw that the two integrals are linearly independent from each other, see Eq. (3.9). We can nevertheless use Eq.s(A.1-A.8), where S(u, b) and U (u, b) are defined, obtaining We can now integrate by parts in b the last term of the above equation; by using the definition of EI [1] 1 (0, u) Eq.(5.9) and the second of Eq.(A.8) one finds the identity (5.14)

The relations at weight one
Clearly, the procedure outlined above to compute EI [1] k (0, u), with k = 0, 1, 2, can be easily repeated for all other weight-one functions EI [1] k (p i , u), and for those involving the function J 0 (u), EJ [1] k (p i , u). We proceed as follows 1-First we use the second order differential operator D(u, d/du) to determine relations between the functions EI [1] 0 (p i , u), EJ [1] 0 (p i , u) and the simpler products of logarithms with I 0 (u) and J 0 (u) functions, Eq. (5.2). Surprisingly, at this order this allows us to rewrite all the functions of this form, where p i is on the the zeros in b of R 4 (u, b), as linear combinations of products of I 0 (u) or J 0 (u) and logarithms.
2-With this results at hand, we obtain the corresponding ones for the EI [1] 2 (p i , u) and EJ [1] 2 (p i , u) by differentiation.
3-Finally, we obtain an expression for the functions EI [1] 1 (p i , u) and EJ [1] 1 (p i , u) by integrating by quadrature their first order differential equation.
We list here explicitly all the relations we find for the functions EI [1] 0 (p i , u) and EJ [1] 0 (p i , u)); for clarity we use the notation in terms of the b integration. We find: Note that if u > 9m 2 all the appearing quantities are real, but the identities are of course valid in general if the proper analytic continuation is taken. For ease of typing, once more, we wrote ln(b − 4m 2 ), ln(u − 9m 2 ) etc. instead of ln (b − 4m 2 )/m 2 , ln (u − 9m 2 )/m 2 . We do not report here all the corresponding relations for the EI [1] 2 (p i ; u), EJ [1] 2 (p i ; u) and EI [1] 1 (p i ; u), EJ [1] 1 (p i ; u) for brevity, but it should be clear that they follow the same pattern as the ones for EI [1] 2 (0, u) and EI [1] 1 (0, u), derived respectively in Eqs. (4.11, 5.12).
Summarizing, the action of the differential operator D(u, d/du) on the E-polylogarithms of weight one associated to the functions I 0 (u) and J 0 (u) allows to reduce their weight and to determine algorithmically surprising (and somewhat unexpected) relations between E-polylogarithms and products of simple logarithms and the functions I 0 (u) and J 0 (u).

E-polylogarithms at weight two and beyond
The detailed study of the E-polylogarithms at weight one revealed surprising identities between the latter and products of complete elliptic integrals and simple logarithms. We would like now to use similar methods to investigate these functions at higher weights. We could of course repeat the same derivation above, say, for the functions derive a second order differential equation satisfied by the latter, and solve it by Eulers variation of constants.
In order to have a better grasp of the general structure, nevertheless, it is useful to study the more general class of functions defined by It is very easy to repeat the same procedure described above and show that all these functions can be expressed in terms of three independent master integrals, say I ǫ (0, u) , I ǫ (1, u) , I ǫ (2, u) .
We can then perform the usual change of basis I 0 (ǫ, u) = I ǫ (0, u) 2) derive a system of differential equations satisfied by these functions, and turn it into a second order differential equation for I 0 (ǫ, u), together with a first order differential equation for I 1 (ǫ, u). Note that in our notation we have The second order differential equation reads together with the equation for I 1 (ǫ, u) We see that there is a residual coupling (suppressed by two powers of ǫ) between I 0 (ǫ, u) and I 1 (ǫ, u). By expanding left-and right-hand-side of Eqs.(6.4, 6.6) and collecting for the terms proportional to ǫ 2 we are left with the following equations D u, d du EI [2] 0 (0, 0, u) = 0 (0, u) , (6.8) while the results at previous orders read Substituting all results explicitly and partial fractioning in v we find EI [2] 0 (0, 0, u) = c 1 I 0 (u) + c 2 2 J 0 (u) and, since lim u→9m 2 EI [2] 1 (0, 0, u) = 0, First of all, let us try to simplify Eq. (6.11). Integrating by parts the first term in dv we get at once EI [2] 1 (0, 0, u) = where in the last line we renamed t → v. Recalling the analytical result for EI [1] 1 (0, u) Eq. (5.12), we see that we have EI [2] 1 (0, 0, u) = 4 3 ln (u − m 2 ) EI [1] 1 (0, u) , indeed, formally similar to the weight-one results for the functions EI [1] 0 (0, u) and EI [1] 2 (0, u) Let us move now to Eq. (6.10) for EI [2] 0 (0, 0, u). At variance with order one, here we need to consider a more general class of integrals u dv 1 ; v n ; 1 v n ; . Following the same logic as at weight one, we generate integration by parts identities of the form and solve the system of equations. Again we work with primitives, up to boundary terms, i.e. the functions X k (u) depend only on the variable u. We find now that for every choice of logarithm, there are again 6 master integrals, which we can choose once more as . More explicitly, once again we find that one of the integrals in Eq. (6.10) can be expressed as linear combination of the others as follows Using this identity in Eq. (6.10) we see that the highest weight do cancel, similarly to the previous order, and we are left with EI [2] 2 J 0 (u) where in the last line we fixed the boundary conditions finding c 2 = 0 . The result in Eq. (6.16) shows interesting features. Indeed, differently from the weight-one case, not all integrals over the functions F 0,0 (u, v) have disappeared. Nevertheless, we see that the piece of highest transcendental weight, i.e. the one involving integrals over F 0,0 (u, v) and logarithms in this case, can indeed be eliminated in favour of a simpler term which contains a logarithm squared multiplied by I 0 (u). The remaining integrals are simpler, as they do not contain any logarithms.

Relations for E-polylogarithms at weight two
Having discussed explicitly the case with a ln 2 b, we can now in principle study all other weight-two Epolylogarithms, including possibly those containing di-logarithms Li 2 (f (b)) with branches corresponding to the roots of the polynomial R 4 (u, b). We can do this similarly to weight one, namely writing a general Ansatz and using the second order differential operator D(u, d/du) to fix the coefficients.
As exemplification, let us consider the following weight two E-polylogarithms All these functions can be rewritten in the notation EI [2] 0 (p i , p j ; u), up to analytic continuation. This is achieved by simply rewriting the (products) of logarithms as standard multiple-polylogarithms, for example ln (b − 4m 2 ) = G(4m 2 , b) + ln 4m 2 ± i π , (6.18) depending on the imaginary part given to b. We use here a standard representation in terms of logarithms to keep the formulas as clear as possible.
In order to build an Ansatz that is large enough to match all these functions, we should consider all functions that behave as weight two or one E-polylogarithms under the action of the operator D(u, d/du). First of all, we include the simplest E-polylogarithms, obtained by multiplying I 0 (u) or J 0 (u) by standard multiple polylogarithms.
A(u) = a 0 ln u + a 1 ln (u − m 2 ) + a 2 ln (u − 9m 2 ) I 0 (u) + a 3 ln 2 u + a 4 ln 2 (u − m 2 ) + a 5 ln 2 (u − 9m 2 ) I 0 (u) + a 6 ln u ln (u − m 2 ) + a 7 ln u ln (u − 9m 2 ) + a 8 ln (u − m 2 ) ln (u − 9m 2 ) I 0 (u) where the a j and b j are numerical coefficients. Note that here, for simplicity, we did not include dilogarithms, which in a more general case should also be included. Simple (products of) logarithms seem to be enough as long as we limit ourselves to (products of) logarithms in the functions (6.17). We have verified explicitly that allowing for the presence of a di-logarithm under the integration sign, requires also to enlarge the Ansatz Eq. (6.19) allowing for di-logarithms as well. We do not report these results for brevity. The Ansatz Eq. (6.19) is not complete, as we can see from the explicit result in Eq. (6.16). From the discussion in Section 4, it is clear that, in general, we must include in the Ansatz 6 more functions, i.e. the master integrals in Eq. (4.21). We write therefore can be decoupled in a two by two coupled system, plus a decoupled first order differential equation. These equations can be solved by Euler's variation of constants, providing a representation of these functions as iterated integrals over rational factors and products of complete elliptic integrals.
This allows to tentatively associate to the E-polylogarithms a weight, dubbed E-weight, which turns out to be naturally lowered by the action of the corresponding (matricial or higher order) differential operator. In this way we could study properties and relations among E-polylogarithms bottom-up in their E-weight and show, in particular, that all E-weight one E-polylogarithms can be rewritten as products of standard polylogarithms and complete elliptic integrals. Starting at E-weight equal to two, this is not true anymore and E-polylogarithms introduce genuine new structures. Nevertheless, also at E-weight two, we found interesting relations for the highest transcendental piece of the E-polylogarithms in terms of products of weight-two standard polylogarithms and complete elliptic integrals. Finally, we used these results to provide a compact representation for the order ǫ 2 of the imaginary part of the two-loop massive sunrise graph.
While our study is not definitive, it might open interesting possibilities for the systematic study and simplification of functions appearing in the calculation of multiloop Feynman graphs with many scales and/or massive propagators. Indeed, the analytic calculation of Feynman integrals which fulfil higher order differential equations still remains largely out of reach; a first obstruction was given by the absence of a systematic understanding of the solution of their corresponding higher-order homogeneous equations. Quite recently it was shown that the study of the maximal cut of Feynman integrals provides an efficient tool to determine the missing homogeneous solutions [28][29][30][31][32] and this obstruction was partially lifted. 5 Thanks to these developments, in fact, we are now in the position to systematically write integral representations for the solutions of complicated Feynman integrals; the crucial problem remains therefore that of studying the properties of these functions and of the relations among them, which is one of the most important aspect of an analytic calculation.
The methods described in this paper are, at least in principle, not limited to elliptic generalizations of multiple polylogarithms and can instead be equally well applied to the study of functions which fulfil even higher order differential equations. We hope therefore that they can be of some use for a systematic analysis of the properties of Feynman integrals beyond multiple polylogarithms.

A The analytical calculation of four master integrals
Concerning the analytic expressions of the master integrals given in (3.9) an explicit calculation, obtained by using the integral representation Eq.(2.1) and by exchanging the order of integration, gives . (A.10) (Note that the integrand in the r.h.s. of (A.9) is real, even if some square roots are imaginary when √ b < 4m).
B Another integral representation for I 0 (u) As an extension of the procedure outlined in Section 2, we will derive a second order equation for the integral with u in the range 9m 2 < u < ∞ for definiteness. The integral is convergent, but we cannot follow exactly all the steps of Section 2, as the direct use of Eq.(2.7), for instance, would involve meaningless (non-convergent) integrals like To our knowledge, that result was found in 1962 by A.Sabry [26], albeit in a somewhat different notation, see Eq.(88) of [26], for the particular case b i = 4m 2 , b j = (W − m) 2 , and used to derive Eq.(85) of that paper, which in our notation reads The result was independently reobtained in [9], see the derivation of Eq.(7.7) there (and later repeated in Eq.s(A. 8,9,10) of [16]) by using the relation Eq.(C.2), fully equivalent to Eq.(A.8) of the present paper, was already given in [41], just after Eq.(5.8) there (but unfortunately with typing errors!).
As explained in [16], if in (C.3) the end points of the integration are taken to be a different pair of roots of the polynomial R 4 (u, b), one can have a non vanishing result; indeed, for b 1 = 0 and b 2 = 4m 2 one finds where −R 4 (u, b) was introduced to keep everything real. That feature was overlooked in [9], where however the roots (0, 4m 2 ) were not of interest.