Algorithms and tools for iterated Eisenstein integrals

We present algorithms to work with iterated Eisenstein integrals that have recently appeared in the computation of multi-loop Feynman integrals. These algorithms allow one to analytically continue these integrals to all regions of the parameter space, and to obtain fast converging series representations in each region. We illustrate our approach on the examples of hypergeometric functions that evaluate to iterated Eisenstein integrals as well as the well-known sunrise graph.


Introduction and motivation
During the past ten years, the Large Hadron Collider (LHC) at CERN in Geneva has been providing us with an impressive amount of data on the dynamics of the fundamental constituents of matter at unprecedented energies and with unprecedented precision. Interpreting correctly this data against theoretical predictions in the Standard Model offers us a precious opportunity to discover signs of New Physics lurking at higher energy scales. The Standard Model is based on Quantum Field Theory (QFT), and in QFT the probability amplitude for a scattering process is encoded in so-called scattering amplitudes. Scattering amplitudes are in turn usually computed perturbatively through an expansion in Feynman diagrams, whose calculation can be reduced to the evaluation of increasingly complicated multi-loop scalar Feynman integrals. While alternative approaches exist which attempt to avoid the combinatorial complexity introduced by an expansion in off-shell Feynman diagrams, eventually one always falls back to the calculation of scalar multi-loop Feynman integrals. The latter can therefore be thought as irreducible building blocks for the computation of any physical observable. As dictated by causality and unitarity in QFT, scattering amplitudes are typically complex functions with a highly non-trivial branch-cut structure, which stems from those regions of the phase-space where virtual particles running in the loops are allowed to go on-shell. The details of their analytic structure must be encoded in the special functions used to represent them, which is expected to become prohibitively complicated with the increase of the complexity of the process and the perturbative order considered.
In the past two decades important steps forward have been achieved towards the classification of the special functions that appear in multi-loop calculations. The first crucial step in this direction was the realisation of the importance of so-called multiple polylogarithms (MPLs) for scattering amplitudes with (mainly) massless internal and external particles. MPLs are the natural generalisation of the logarithm, and they can be defined (loosely) as the class of special functions required to perform all iterated integrals of rational functions with poles on the complex plane. While MPLs had been known to the mathematical community for over a century [1,2], their re-discovery in the high-energy physics community [3,4] contributed to initiate a new spurt of interest in these functions. This has ultimately led to a much deeper understanding of their analytic and algebraic properties [5][6][7][8], which is important in pehnomenological applications. Indeed, these developments were crucial to make a multitude of precision calculations in the Standard Model possible and, in this way, to provide a theoretical interpretation of the most recent results obtained at the LHC. Scattering amplitudes are typically only the first step towards a complete physical prediction. For example, at the LHC the calculation of differential distributions for interesting observables requires evaluating complicated multi-dimensional phase-space integrals over the scattering amplitudes, which can in general only be performed numerically by Monte Carlo techniques. The latter typically require being able to evaluate numerically the scattering amplitudes (hundreds of) thousands of times, in different regions of the relevant phase-space, typically crossing multiple branch cuts. For this reason, the importance of being able to evaluate numerically fast and with high-precision the special functions required and to properly handle their branch cuts, cannot be overestimated for realistic physical applications.

JHEP02(2020)105
Thanks to the developments hinted to above, today we can claim that both goals are achievable as long as the scattering amplitudes can be expressed in terms of MPLs only. Rather generic codes have been published which can evaluate numerically MPLs for (in principle) arbitrary complex arguments [9][10][11][12][13][14]. Still, in cases which are complicated enough, a brute force application of these algorithms appears unfeasible, as the numerical evaluation for arbitrary phase-space points across multiple branch-cuts tends to become too slow and inefficient for their use in Monte Carlo integration codes. Fortunately, this problem can be greatly alleviated by properly manipulating the relevant MPLs and reexpressing them in terms of new combinations of MPLs which can be more easily evaluated numerically in the region of interest. To achieve this, one typically uses the knowledge of the algebraic and analytical properties of MPLs in order to express MPLs whose numerical evaluation would require crossing complicated branch cuts, in terms of MPLs which are free of those branch cuts in the region considered and whose numerical values can be obtained by fast converging series expansions. While it cannot be proved that this strategy will always be successful, its potential has been shown in many examples, most notably in the evaluation of the two-loop QCD corrections to the productions of pairs of vector bosons at the LHC [15][16][17]. It is worth stressing that this step was essential to make the a large number of phenomenological studies at the LHC possible.
In spite of these impressive results, this is not the end of the story. As it is by now very well known, starting at the two-loop order MPLs are not enough to express all interesting scattering amplitudes. The first occurrence of classes of iterated integrals beyond MPLs can be traced back to a paper by A. Sabry in 1962, where he discovered new functions of elliptic type in the calculation of the two-loop corrections to the electron propagator in QED with massive electrons [18]. Renewed interest in these functions arose when it was realised that they appear in a large number of calculations relevant for collider physics, most notably in the two-loop correction to tt production at hadron colliders [19]. The simplest occurrence of these new functions was identified in the so-called massive two-loop sunrise graph, whose computation has since then received a lot of attention from different sides of the highenergy physics community [20][21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38]. Indeed, it has been recently realised that a natural generalisation of MPLs to include the new class of functions that appear in the sunrise graph exists. These so-called elliptic multiple polylogarithms (eMPLs) can be obtained by considering iterated integrals of rational functions on a genus one complex surface, i.e. a torus, which is well known to be mathematically equivalent to an elliptic curve. Their description as iterated integrals on a complex torus had been formalised by mathematicians already in 2011 [39,40] and had later found application in the calculation of one-loop scattering amplitudes in string theory [41][42][43][44]. More recently, an alternative formulation of eMPLs naturally defined on an elliptic curved defined by a polynomial equation instead of on a torus has been proposed [27,45,46], which has made it possible to naturally compute different previously out-of-reach multi-loop Feynman integrals [47][48][49]. While eMPLs are a rather large class of functions that can encompass different kinds of physical problems, for the particular case of the two-loop sunrise graph, 1 a special subclass of eMPLs has been shown to be sufficient, i.e. the so-called iterated Eisenstein integrals [50][51][52][53][54].

JHEP02(2020)105
In the past years there has been an impressive progress in understanding the analytical and algebraic properties of these new functions. In spite of that, algorithms for the fast and precise numerical evaluation of iterated Eisenstein integrals and, more in general eMPLs, are still missing. 2 Without this step, their use in state-of-the-art physical calculations will not become possible. This paper aims to start closing this gap, at least for the restricted set of one-parameter Feynman integrals that can be expressed in terms of iterated Eisenstein integrals. Examples of such integrals involve the well-known two-loop equal-mass sunrise and kite integrals [46,53], as well as the three-loop banana graph [59]. Indeed, we will show that for the case of iterated Eisenstein integrals, algorithms similar to those used to treat MPLs can be devised, which allow one to efficiently evaluate complicated combinations of these functions for any values of the external kinematics of the problem considered. Similarly to what is done with MPLs, we will show that given a Feynman integral expressed in terms of iterated Eisenstein integrals, irrespective of the kinematical point we are interested in, it is always possible to rewrite it in terms of functions which can be numerically evaluated by a fast converging series expansion, making their use suitable for realistic physical applications. The rest of the paper is organised as follows: in section 2 we review the definition of modular forms in particular showing how to define a parity-invariant basis for the latter, which is explicitly real or imaginary. In section 3 we define iterated integrals of modular forms for this new parity-invariant basis and show how to use it to express a special class of hypergeometric functions. We then have a short intermezzo in section 4 to discuss under which assumptions it is possible to express iterated Eisenstein integrals in terms of simpler MPLs. In section 5 we then come to the central idea of this paper, which consists in an algorithm to perform modular transformations on iterated Eisenstein integrals in order to rewrite them in terms of combinations of iterated Eisenstein integrals which can be more easily numerically evaluated. The details of the numerical evaluation are discussed in section 6, while explicit examples to classes of hypergeometric functions and to the sunrise integral are provided in sections 7 and 8 respectively. Finally, we draw our conclusions in section 9.

Modular forms and Eisenstein series
The purpose of this section is to review in detail the modular group SL(2, Z) and modular forms. Most of the material in this section is well known in the mathematics literature. We include it nevertheless because modular forms and their properties are the foundations on which subsequent sections are built.

The modular group and its congruence subgroups
The group SL(2, Z) consists of all 2 × 2 matrices with integer entries and unit determinant. It is generated by the three elements

JHEP02(2020)105
In other words, every element of SL(2, Z) can be written as a finite product of S's and/or T 's up to a sign. The decomposition of an element of SL(2, Z) into a product of generators can be performed algorithmically (see appendix A for a review). We are often not interested in the whole modular group SL(2, Z), but only in a certain subgroup thereof. Of particular interest are the so-called congruence subgroups. A subgroup Γ ⊆ SL(2, Z) is called a congruence subgroup of level N if it contains the principle congruence subgroup of level N , Particularly important examples of congruence subgroups are For the rest of this section, Γ will denote some congruence subgroup of SL(2, Z) (which may be SL(2, Z) = Γ(1) itself). Γ naturally acts on the upper half-plane H = {τ ∈ C|Im τ > 0} via Möbius transformations, Note that −1 acts trivially on τ , i.e., (−1) · τ = τ . The action in eq. (2.4) partitions the upper half-plane into distinct orbits. The space Y Γ = H/Γ of orbits is called the modular curve for Γ. It can be useful to consider a connected domain in H whose points correspond to the distinct orbits in Y Γ , called a fundamental domain for Γ. A particularly important example of a fundamental domain is the one for the full modular group SL(2, Z), which is F = {τ ∈ H : −1/2 ≤ Re τ < 1/2 and |τ | > 1} ∪ {τ ∈ H : −1/2 < Re τ ≤ 0 and |τ | = 1} .
(2.5) In other words, for every τ ∈ H there are unique γ ∈ SL(2, Z) and τ 0 ∈ F such that τ = γ · τ 0 . Γ also acts on Q ∪ {i∞}, and so it is natural to consider its action on the extended upper half-plane H = H ∪ Q ∪ {i∞}. It partitions Q ∪ {i∞} into a finite number of distinct equivalence classes, called the cusps of Γ. For every c 1 , c 2 ∈ Q ∪ {i∞}, there is a γ ∈ SL(2, Z) such that c 2 = γ · c 1 , which implies that SL(2, Z) has only one cusp. The matrix γ can be constructed explicitly: first, we note that it is sufficient to consider Next, let c 2 = m/n ∈ Q. It follows from the properties of the greatest common divisor that there are integers p and q such gcd(m, n) = mp+nq. 3 If m and n are relatively prime, gcd(m, n) = 1, we can take γ c 2 = m mr−q n nr+p , with r an arbitrary integer. This matrix obviously satisfies c 2 = m/n = γ c 2 · i∞. Moreover, this matrix has unit determinant, and so it defines an element of SL(2, Z).

Modular forms
We start by considering functions from the extended upper half-plane into the complex numbers that are invariant under Γ. A modular function for Γ is a meromorphic function f : , for all γ ∈ Γ. One can show that there are no non-constant holomorphic modular functions, i.e., every non-constant function that is invariant under Γ has at least one pole. If we want to obtain holomorphic functions, we need to consider more general transformation behaviours under Γ. For every positive integer n, we define This defines for every n a genuine action of Γ on functions on the extended upper half-plane.
In particular, we have A modular form of weight n for Γ is a function f : H → C such that 1. f invariant under the action of Γ in eq. (2.6), i.e., it satisfies (f |nγ )(τ ) = f (τ ).
2. f is holomorphic on H.
We can think of modular forms as holomorphic functions, i.e., functions without poles, that are invariant under the action in eq. (2.6). Modular forms of weight n form a finitedimensional vector space M n (Γ). Moreover, M n (Γ) is a graded algebra, i.e., the product of two modular forms of weight n 1 and n 2 is a modular form of weight n 1 + n 2 . Since there are no non-constant holomorphic modular functions, all modular forms of weight zero are constant. If Γ is a congruence subgroup of level N and f ∈ M n (Γ), then we say that f has level N . Every modular form of level N is invariant under T N = 1 N Having identified modular forms for Γ as invariants under the group action in eq. (2.6), it is natural to ask how a modular form f ∈ M n (Γ) transforms for some γ ∈ SL(2, Z) that is not an element of Γ.
We start by assuming that Γ ⊆ SL(2, Z) is a normal subgroup. 4 Then, for every γ ∈ SL(2, Z) and γ 1 ∈ Γ, there is γ 2 ∈ Γ such that γγ 1 = γ 2 γ. We have where the last step follows from the fact that f is a modular form for Γ, and thus invariant.
We conclude that if f is a modular form for a normal subgroup Γ, then f |nγ is also a modular form for Γ of the same weight, for all γ ∈ SL(2, Z).

JHEP02(2020)105
The previous argument relies crucially on Γ being normal in SL(2, Z). However, not all subgroups of SL(2, Z) are normal, not even all congruence subgroups. In particular, the congruence subgroups Γ 0 (N ) and Γ 1 (N ) are in general not normal. The principle congruence subgroups Γ(N ), however, are always normal. We can therefore generalise the previous statement in a weaker form to general congruence subgroups: if f is a modular form for a congruence subgroup Γ of level N , then f |nγ is a modular form for the principle congruence subgroup Γ(N ) of the same weight and level, for all γ ∈ SL(2, Z). We will see at the end of this section that it is in general not possible to make stronger statements, i.e., f |nγ will in general not be invariant under Γ, but it is only invariant under the smaller subgroup Γ(N ).
While the previous argument shows that a modular form for a congruence subgroup Γ of level N in general transforms into a linear combination of modular forms for Γ(N ), it does not allow us to predict this linear combination. This situation changes if we restrict the discussion to a subset of modular forms, the so-called Eisenstein series, which we will study in more detail in the next section.

Eisenstein series
The vector space M n (Γ) of modular forms of weight n for Γ admits a direct sum decomposition M n (Γ) = E n (Γ) ⊕ S n (Γ) . (2.10) Here S n (Γ) denotes the subspace of cusp forms, i.e., the subspace of modular forms of weight n for Γ that vanish on Q ∪ {i∞}. Its (orthogonal) complement is the Eisenstein subspace E n (Γ). So far mostly Eisenstein series have appeared in the context of Feynman integral calculations. We therefore study the structure of the Eisenstein subspaces in detail in the remainder of this section. We first present an explicit basis for the space of Eisenstein series. There are numerous ways to write down a basis for E n (Γ) (cf., e.g., ref. [60]). Our choice is motivated by the fact that Eisenstein series for Γ(N ) are closely related to the Kronecker series and elliptic polylogarithms [40,43,46]. Since Γ(N ) ⊆ Γ for every congruence subgroup Γ of level N , we only discuss the case of the principle congruence subgroup Γ(N ), and we refer to the literature for the other cases (see, e.g., ref. [60]).
We start by writing down a spanning set for E n (Γ(N )) [46], where n and N are integers greater than unity and r, s are integers defined modulo N . These functions transform like eq. (2.6) and form a spanning set for E n (Γ(N )), except for the case n = 2 and (r, s) = (0, 0) mod N , which is not a modular form. More generally, we have [46],

JHEP02(2020)105
The functions defined in eq. (2.11) admit a q-expansion that can be written down in closed form h The Fourier coefficients for k ≥ 1 are given by (cf., e.g., ref. [60]), (2.14) The constant term for n ≥ 2 is For n = 1 the constant term is given by (2.17) Note that all the Fourier coefficients of h and, for every d divides N and 0 ≤ ρ, σ < N , the distribution identity [46], These relations still hold for n = 1, though in that case there are additional relations which can easily be found by looking for linear relations among the q-series. For n ≥ 3, it is possible to explicitly write down a subset of the h (n) N,r,s that are linearly independent [46].

A parity-invariant spanning set
In the previous section we have presented a spanning set of Eisenstein series for Γ(N ) and we have described all the relations among these functions. This is important in applications, because it allows us to construct explicit bases for E n (Γ(N )). The bases obtained in this way have the disadvantage that the resulting functions will in general be complex-valued. In applications, however, it can be useful to work with basis elements that are manifestly real, at least in certain regions of the parameter space. In this section we present a spanning set of functions for E n (Γ(N )) that are explicitly real when τ is purely imaginary. We start by analysing how h We see that when τ is purely imaginary, complex conjugation amounts to changing the sign of s. Based on this observation, we define, where the second equality follows from the reflection identity in eq. (2.18). These functions also admit a representation in terms of lattice sums similar to eq. (2.11), a n,N,r,s (τ ) = 1 2 (2.25) 6 We are grateful to Nils Matthes for suggesting this counterexample.

JHEP02(2020)105
Equation (2.22) implies that the functions a n,N,r,s (τ ) and b n,N,r,s (τ ) are real whenever τ is purely imaginary. Since they are linear combinations of the Eisenstein series h (n) N,r,s (τ ), they also form a spanning set for E n (Γ(N )).We can easily determine a basis of a n,N,r,s (τ ) and b n,N,r,s (τ ). In particular, for n ≥ 3 a basis for E n (Γ(N )) is The set B 2 N,n depends on whether the weight n is even or odd. If n is even, we have B 2 N,n = a n,N,k,0 , a n,N,k,N/2 , a n,N,0,k , a n,N,N/2,k : 0 < k ≤ N/2 and gcd(N, k) = 1 , (2.28) while for n odd, we have Note that when N is odd, functions with index N/2 are absent.

Definition and basic properties
In this section we define our main objects of interest, namely iterated integrals of modular forms and iterated Eisenstein integrals. Consider a set of modular forms h j (τ ). We define differential forms ω j = dτ 2πi h j (τ ) and consider the iterated integrals [54,61], The recursion starts with I(; τ 0 , τ ) = 1. The number of integrations k is called the length of the iterated integral. We can also define a notion of weight [47]. Assume that all modular forms have level N and that they admit q-expansions of the form where each Fourier coefficient a n (h j ) is a period of weight equal to n j (cf. for example the Fourier coefficients of the Eisenstein series in eq. (2.14), which are all algebraic multiples of powers of π). We define the weight of I(h 1 , . . . , h k ; τ 0 , τ ) as −k + k j=1 n j . The integrals in eq. (3.1) have all the standard properties of iterated integrals, cf. e.g., ref. [62]. In particular, they form a shuffle algebra, 3)

JHEP02(2020)105
where the sum runs over all shuffles of (h 1 , . . . , h k ) and (h k+1 , . . . , h l ), i.e., over all permutations of their union that preserve the ordering within each set. The integrals satisfy the path composition and reversal formulas, It is easy to see that I(h 1 , . . . , h k ; τ 0 , τ ) is linear in its first k arguments, where the c i are constants. Finally, if all the modular forms h j are real on the imaginary axis, then I(h 1 , . . . , h k ; τ 0 , τ ) is real whenever τ and τ 0 are purely imaginary. The path composition formula in eq. (3.4) allows us to consider only a fixed base point τ 0 , and from now on we will always choose τ 0 = i∞. We will use the shorthand notation I(h 1 , . . . , h k ; τ ) = I(h 1 , . . . , h k ; i∞, τ ). However, care is needed to interpret these integrals, because strictly speaking they are divergent and require regularisation. In order to see that, we change variables from τ to q N . We can write Hence, if the constant term in the q-expansion is non zero, ω j has a pole at q N = 0, or equivalently τ = i∞. As a result, the integral I(h 1 , . . . , h k ; τ ) requires regularisation whenever a 0 (h k ) = 0. In ref. [61] it was shown how to replace logarithmically-divergent iterated integrals of modular forms by suitably regularised versions, such that all algebraic properties (e.g., the shuffle algebra) are preserved. This is achieved by interpreting eq. (3.1) as a linear combination of absolutely convergent integrals multiplied by powers of τ [61], where the map R is defined by [61] where denotes the shuffle product, and we defined JHEP02(2020)105

The elliptic 2 F 1 function
In this section we illustrate the concepts from the previous section on a concrete example, namely the family of integrals defined by where n i ∈ Z and a, b, c are complex numbers. These integrals are closely related to a family of hypergeometric 2 F 1 functions, (3.11) The integral in eq. (3.10) defines a meromorphic function of , and we interpret the integral as a Laurent series in . In the following we assume that z is real. As long as z < 1 the integral is real order by order in , but it develops an imaginary part when z > 1. We assume that the branches of all functions are determined by assigning a small positive imaginary part to z. The dependence of the Laurent coefficients on a, b, c is polynomial. For simplicity, we only discuss the case a = b = c = 1. All conclusions remain true in the general case.
Using integration-by-parts identities, one can show that all integrals in the family in eq. (3.10) can be expressed as a linear combination of two independent master integrals, which can be chosen as The square root in the integrand defines an elliptic curve via the polynomial equation The master integrals were evaluated in refs. [45,46,51] orderby-order in in terms of elliptic polylogarithms and iterated Eisenstein integrals for the congruence subgroup Γ(2). The result can be cast in the form: where S(z) denotes the matrix . (3.14) Here K(z) and E(z) denote the complete elliptic integrals of the first and second kind, These integrals are real and well defined for z < 1. Here we only focus on the region 0 < z < 1 and we defer a discussion of the analytic continuation to the other regions to JHEP02(2020)105 section 7. In this region the functions U i (τ ) are power series in whose coefficients can be expressed in terms of iterated Eisenstein integrals. The variable τ ∈ H on which the iterated Eisenstein integrals depend is related to the variable z in eq. (3.10) via Note that τ is purely imaginary for 0 < z < 1. The functions K(z) and i K(1 − z) are in fact periods of the elliptic curve defined by the polynomial equation The first few orders in the Laurent expansion of the functions U i (τ ) read and U 2,0 (τ ) = 0 , Let us make some comments about these results. First, since τ is purely imaginary for 0 < z < 1, we see that U 1 and U 2 are manifestly real and imaginary respectively (the explicit factor of i in U 2 is cancelled by another explicit factor of i in eq. (3.14)). Second, the functions U i are pure functions of uniform weight [63] in the sense of ref. [47]. We emphasise that this statement depends on our choice for the periods of the elliptic curve: in eq. (3.14) we have singled out the periods K(z) and i K(1 − z). We could of course have made a different choice of periods, and any two choices are related by an SL(2, Z) transformation, such that τ transform as in eq. (2.4). A different choice would have led to a different form for the matrix S(z) and the pure functions U i (τ ). We will come back to this point in section 7 when we will discuss the analytic continuation of the integrals outside the range 0 < z < 1. Here we only mention that, if we denote by S(γ, z) the matrix corresponding JHEP02(2020)105 to choosing another basis of periods (cf. eq. (3.20)), then this matrix will be related to S(z) = S(1, z) through the relation Note that the matrices M (γ, τ ) respect the group law of SL(2, Z),

Expressing iterated Eisenstein integrals in terms of multiple polylogarithms
In this section we discuss and clarify the relationship between iterated Eisenstein integrals and another class of iterated integrals that often show up in Feynman integral computations, namely multiple polylogarithms. We start by briefly reviewing multiple polylogarithms and their main properties, before we discuss their connection to iterated Eisenstein integrals in subsequent sections.

Multiple polylogarithms
Multiple polylogarithms (MPLs) are a generalisation of the ordinary logarithm and the classical polylogarithm function. They are defined by the iterated integrals [6,64], G(a 1 , . . . , a n ; z) = z 0 dt t − a 1 G(a 2 , . . . , a n ; z) , where the a i and z are (constant) complex numbers and the recursion starts from G(; z) = 1.
In the case a n = 0, the integral in eq. (4.1) is divergent, and we define instead, The number of integrations n is called the weight of the MPL. MPLs contain the ordinary logarithm and the classical polylogarithms as special cases, e.g., In order to understand the relation between MPLs and (certain) iterated Eisenstein integrals, it is instructive to understand the geometrical background leading to the definition of MPLs. Consider the space of rational functions in one variable x with poles at most at some fixed positions x = a i ∈ Σ (we assume Σ finite). It is easy to see that not every rational function of this type has a primitive within this space of functions, but we need to enlarge the space by including logarithms of the form G(a i , x) = log(1 − x/a i ). Indeed,

JHEP02(2020)105
we cannot find any rational function R(x) such that its derivative equals 1/(x − a i ). This 'obstruction' to finding a rational primitive is what naturally leads to the study of the (algebraic) de Rham cohomology groups H 1 dR (X). In our case the space X is the Riemann sphere with the points a i ∈ Σ removed, X = CP 1 \ Σ, and the 'obstruction' to finding a rational primitive can be cast in the form of the statement that H 1 dR (CP 1 \ Σ) is generated by the logarithmic one-forms dx/(x − a i ), The logarithmic one-forms are precisely the integration kernels that appear in the definition of the MPLs in eq. (4.1). Loosely speaking, the first de Rham cohomology group of a Riemann surface allows us to classify those meromorphic functions on CP 1 \ Σ that we cannot trivially integrate in terms of meromorphic functions, but that require the introduction of new classes of functions. In the next section we will see that this cohomological intepretation allows us to express certain classes of iterated Eisenstein integrals in terms of MPLs.

Iterated Eisenstein integrals and MPLs
In section 2 we have defined the modular curve Y Γ = H/Γ. It can be shown that Y Γ is always a compact Riemann surface with a finite number of punctures, i.e., with a finite number of points removed. The punctures correspond to the cusps of Γ. A compact Riemann surface is characterised by its genus g, which, loosely speaking, counts the number of handles of the surface. Every compact Riemann surface of genus g = 0 is isomorphic to the Riemann sphere CP 1 . In the previous section we have described the first de Rham cohomology group of a Riemann sphere with a finite number of points removed, and we have argued that the iterated integrals on that space are naturally connected to MPLs. It is therefore natural to expect that in cases where the modular curve Y Γ has genus zero, there is a connection between iterated integrals of modular forms and MPLs. In this section we briefly review this connection. The results of this section are well known in the mathematics literature. We review them here in a way appropriate for a physics audience. Every meromorphic differential one-form on a Riemann surface, for example, every one-form of the form dz R(z), where z is a local coordinate on the Riemann surface and R(z) is a meromorphic function, belongs to one of the following three classes: • Differentials of the first kind are globally holomorphic, i.e., they have no poles anywhere on the Riemann surface.
• Differentials of the second kind have poles with vanishing residues.
• Differentials of the third kind have poles with non-vanishing residues.
The first de Rham cohomology group of a Riemann surface of genus g is generated by precisely g differentials of the first and g differentials of the second kind, as well as one differential of the third kind for every puncture. In particular, for Riemann surfaces of genus zero the first cohomology group is entirely generated by differentials of the third kind, in agreement with eq. (4.4).

JHEP02(2020)105
It turns out that there is a connection between modular forms of weight two for Γ and the first de Rham cohomology group of Y Γ . To see this, we start by noting that In other words, dτ transforms like a modular form of weight −2. Hence, every modular form h(τ ) of weight two defines a differential dτ h(τ ) that is invariant under Γ, and therefore furnishes a well-defined differential on Y Γ . In addition, there is a one-to-one map between the differentials of the first kind on Y Γ (rather, the classes they define in the first de Rham cohomology group) and the cusp forms of weight two. In particular, this means that the genus of Y Γ is equal to the number of cusp forms of weight two for Γ. Similarly, there is a one-to-one map between the differentials of the third kind of Y Γ and modular forms of weight two from the Eisenstein subspace. Assume now that we are given an iterated integral of Eisenstein series of weight two for some congruence subgroup Γ for which Y Γ has genus zero, i.e., Y Γ is the Riemann sphere with a finite number of points removed. This happens precisely when dim S 2 (Γ) = 0. For the congruence subgroups in eqs. (2.2) and (2.3) it is well known for which level the modular curves have genus zero: In these cases there is a modular function f : H → Y Γ for Γ, called the Hauptmodul, such that every modular function for Γ is a rational function of f (τ ). Hauptmodule are known for many of the modular curves considered above, see, e.g., refs. [65,66]. If h(τ ) is an Eisenstein series of weight two for Γ, dτ 2πi h(τ ) is invariant under Γ, and f maps it to a rational differential form on Y Γ . Since we know all the generators of the first de Rham cohomology group of Y Γ from eq. (4.4), we can change variables to t = f (τ ) and write dτ h(τ ) as a linear combination of logarithmic differential forms dt/(t − t i ), where the t i = f (τ i ) run over the images of the cusps of Γ under the Hauptmodul f . The conclusion of this analysis can be summarised as follows: If Γ is a congruence subgroup such that Y Γ has genus zero, and f : H → Y Γ is a Hauptmodul for Γ, then every iterated Eisenstein integral of the form I(h 1 , . . . , h k ; τ ) with h j ∈ E 2 (Γ) can be written as a linear combination of uniform weight k of MPLs evaluated at t = f (τ ).
Let us illustrate this explicitly on the iterated Eisenstein integrals that appear in eqs. (3.18) and (3.19). The relevant modular curve Y (2) has genus zero, and a Hauptmodul is the modular lambda function,

JHEP02(2020)105
where η(τ ) denotes the Dedekind eta function, The modular lambda function is the inverse of eq. (3.16) for 0 < z < 1, i.e., z = λ(τ ). Γ(2) has three cusps, which may be represented by τ = 0, τ = 1 and τ = i∞. The cusps are mapped by the Hauptmodul to the following points on the Riemann sphere, Hence, we can identify the modular curve Y (2) with the Riemann sphere with the points 0, 1 and ∞ removed. We therefore expect that all iterated integrals of Eisenstein series of weight two and level two (except for a 2,2,0,0 ) can be expressed in terms of MPLs in z with singularities only when z ∈ {0, 1, ∞}. The corresponding set of functions is a well-known subset of MPLs known as harmonic polylogarithms [3]. In the remainder of this section we work this out in detail on the iterated integrals of Eisenstein series of weight two that appear in eqs. (3.18) and (3.19). Letting z = λ(τ ), we can match q-expansions to show that dτ 2πi a 2,2,1,0 (τ ) = 1 12 (4.9) Inserting these relations into the integrands of the iterated Eisenstein integrals and using the fact that the lower integration boundary τ = i∞ is mapped to z = 0, we can easily perform all integrations in terms of MPLs. We find  The constants proportional to log 2 require some explanation, because they cannot be recovered by inserting eq. (4.9) into the integrand and performing the z integration over the range [0, z] in a naive way. In fact, this term has its origins in a different choice of regularisation scheme to regulate the logarithmic divergences at τ = i∞ and z = 0, cf. eqs. (3.7) and (4.2).
To understand this, we first note that if we expand eq. (4.6) into a q-series, we find and so G(0; z) = log λ(τ ) = 4 log 2 + log q 2 + O(q 2 ) . We see that terms proportional to (powers of) logarithms of log 2 arise from the fact that z ∼ 16 q 2 as q 2 → 0. These factors can always be recovered from the following recipe: 1. Insert eq. (4.9) into the integrands of the iterated Eisenstein integrals and perform the integrals in the naive way.

JHEP02(2020)105
2. Use the shuffle algebra properties to express every G(a 1 , . . . , a n ; z) with a n = 0 in terms of MPLs with a n = 0 and powers of G(0; z).
Using this procedure, we can cast the functions U i from section 3.2 in a simpler form involving only a handful of functions which cannot be expressed in terms of MPLs, For z < 0, we obtain, with Note that for z < 0, τ 1 is purely imaginary and Im τ 1 > 0, and the matrix γ 1 lies in SL(2, Z). Similarly, for z > 1, we find, and Again, we find that τ 3 is purely imaginary with positive imaginary part for z > 1, and that the matrix γ 3 lies in SL(2, Z). Hence, if we define γ 2 = ( 1 0 0 1 ) and τ 2 = i K(1 − z)/ K(z) for 0 < z < 1, we see that for each of the three regions there is γ i ∈ SL(2, Z) and τ i ∈ H such that τ = γ i · τ i . We can also define a map that assigns to all real values of z the corresponding value of τ , In the remainder of this section we describe how we can express iterated Eisenstein integrals evaluated at τ = γ i · τ i in terms of iterated Eisenstein integrals in τ i . The motivations for doing this are twofold. First, when z > 1, the integrals develop an imaginary part, and it can be desirable to make the imaginary parts explicit and only work with iterated integrals that are manifestly real. Since the Eisenstein series in eq. (2.24) have the property of being real when evaluated at purely imaginary arguments, this will automatically be the case if we work with integrals evaluated at τ i . Second, as we will see in section 6, we can speed up the numerical evaluation of the iterated Eisenstein integrals in different regions using modular transformation on its arguments.

Modular transformations of iterated Eisenstein integrals
The goal of this section is to describe a procedure to express iterated Eisenstein integrals of the form I(h 1 , . . . , h k ; γ · τ ), with γ ∈ SL(2, Z), in terms of iterated Eisenstein integrals evaluated at τ . The general strategy is simple: we start from the regularised version of the integral in eq. (3.7), perform the change of variables τ → γ · τ in the integrand, and integrate back. Since SL(2, Z) is generated by the three elements S, T and −1 in eq. (2.1), it is sufficient to discuss the cases γ = S and γ = T , because a generic γ ∈ SL(2, Z) can always be written as a (finite) product of S's and T 's, up to a sign.
Before we discuss the procedure in detail, we stress the importance of working with the regularised integrals in eq. (3.7): treating the divergent integrals naively can lead to wrong results! Indeed, consider the integral Next, write τ = 1 +τ = ( 1 1 0 1 ) ·τ . Using the definition of the iterated Eisenstein integrals, we can write, This integral diverges, and it needs to be interpreted as a regularised integral as discussed in section 3.1. As we will illustrate now, as a consequence of the regularisation, we cannot JHEP02(2020)105 do naive changes of variables on the integral. We can change variables according to τ → 1 + τ . The integrand is invariant under this change of variables, because a 4,1,0,0 (1 + τ ) = a 4,1,0,0 (τ ). Naively, the lower integration boundary i∞ does not change under translations, and we are led to believe that I(a 4,1,0,0 ; 1 +τ ) is identical to I(a 4,1,0,0 ;τ ). It is easy to see from eq. (5.7) that this cannot be correct: since e 2πi(τ +1) = e 2πiτ = q 1 , we see that the power series part in eq. (5.7) remains unchanged. The first term, however, does change, and we find, We see that the constant in eq. (5.9) is a consequence of the regularisation. Indeed, if we start from the regularised version in eq. (3.7), we find The second term is an absolutely convergent integral and gives rise to the power series part in eq. (5.7). We can naively perform the change of variables and conclude that it is invariant under translations. The first term immediately produces a constant offset under translations, and so we recover eq. (5.9).
We can apply exactly the same steps to more general modular transformations: we always start from the regularised version in eq. (3.7) and perform a change of variables in the absolutely convergent integrals. Before we proceed, however, we need to restrict the class of iterated integrals that we want to consider. First, we note that we can always use the shuffle algebra to transform all iterated integrals into a new basis of integrals where none of the ω j correspond to a modular form of weight zero (except for explicit powers of τ ), but instead we consider differential forms of the form dτ h j (τ ) τ m j , for some positive integer m j and h j a modular form of weigh n j . An example will clarify this: From here on, and for the rest of this paper, we will restrict ourselves to iterated integrals where, after transforming the integrals to this new basis, only differential forms dτ h j (τ ) τ m j with 0 ≤ m j < n j − 1 and n j > 1 appear. The reason for this restriction will become clear shortly. The same restriction has already appeared in ref. [54]. We emphasise that this restriction is not a problem for the applications that we have in mind, because in all known physics applications only iterated Eisenstein integrals for this restricted class appear.
Let us now consider an integral from this restricted class, and let us work out the change of variables explicitly. We write down formulas valid for arbitrary γ ∈ SL(2, Z), though we keep in mind that we can recover every modular transformation via a suitable JHEP02(2020)105 sequence of S and T transformations. If γ = a b c d ∈ SL(2, Z), we find, where ω γ j denotes the transform of ω j under γ, In the previous section we have seen that if h j (τ ) is a modular form of weight n j for Γ(N ), then so is h j (γ · τ ). In particular, if h j = h (n j ) N,r,s , eq. (2.12) implies, N,rd+sb,rc+sa (τ ) . (5.14) The transformation behaviour of the Eisenstein series a n,N,r,s and b n,N,r,s can easily be recovered from their definition in terms of h We can then simply integrate back in terms of (regulated) iterated integrals of modular forms to obtain the desired result. We stress that at this point it is crucial that we work with the restricted class of integrals: indeed, if we change to the basis of integrals where the differential forms have the form dτ h j (τ ) τ m j , then these differential forms transform as We see that no negative powers appear precisely when the condition 0 ≤ m j ≤ n j − 2 is fulfilled. This condition also automatically enforces n j ≥ 2, which are precisely the conditions that define our restricted class of integrals. We note that this condition is preserved under modular transformations [54] (which can be seen, for example, from the fact that the highest power of τ in eq. (5.16) never exceeds n j − 2 for 0 ≤ m j ≤ n j − 2). In other words, we have described an algorithm that allows us to derive modular transformations of iterated integrals of modular forms from this restricted class, and this restricted class is closed under modular transformations. Second, we see from eq. (5.12) that the lower integration boundaries are constant rational values − d c = γ −1 ·i∞. If we split the path of integration so that all lower integration limits are the infinite cusp, then we are naturally led to consider iterated Eisenstein integrals evaluated at constant rational values. If we decompose γ into a product of S's and T 's, we only obtain constant iterated Eisenstein integrals evaluated at τ = 0, because T · i∞ = i∞ and S · i∞ = 0 . (5.17)

JHEP02(2020)105
Iterated Eisenstein integrals evaluated at τ = 0 are related to a class of transcendental numbers that has recently appeared in the mathematics literature. We will review this class of transcendental constants in the remainder of this section. For modular forms h j , 1 ≤ j ≤ k, of weight n j ≥ 2 and integers m j such that 1 ≤ m j < n j , the associated multiple modular value (MMV) is defined by the integral [54,61]: with m = m 1 + . . . + m k . We note that the restrictions on m j and n j that appear in the definition of MMVs are precisely the conditions that define our restricted class of integrals. The integral in eq. (5.18) may diverge, and we have to interpret it again as a regularised version. We define the weight of an MMV as n 1 + . . . + n k and the length as m. MMVs of small weight and lengths associated to modular forms of level one have been studied in refs. [61,67]. It is easy to see that MMVs form a shuffle algebra, where n is the weight of the MMV, and The reflection identity can easily be obtained by working out the transformation of the corresponding iterated integrals under the modular transformation S : τ → −1/τ . MMVs are (conjecturally transcendental) constants. In applications it is important to obtain their numerical value. Since MMVs are special cases of iterated integrals of modular forms, namely those evaluated at τ = 0, we can try to proceed in the same way as for iterated integrals and evaluate them numerically using q-expansions. In this case, however, the expansion parameter is q N = 1, which leads to badly converging expansions. In the next section we describe a way to accelerate the convergence of the series.

Asymptotic expansions and numerical evaluation of iterated Eisenstein integrals
In this section we discuss how to efficiently evaluate Eisenstein series and their iterated integrals. A key ingredient is the behaviour of iterated Eisenstein integrals under modular JHEP02(2020)105 transformations discussed in section 5. Our strategy is to derive fast converging q-series close to every cusp. We therefore start by discussing how to obtain asymptotic expansions of Eisenstein series and their iterated integrals in general, before we turn to the numerical evaluation in subsequent sections.

Asymptotic expansions of iterated Eisenstein integrals
In applications it is often useful to expand Feynman integrals into an asymptotic series, e.g., close to a singular point of the differential equation satisfied by the integral. When transformed to the variable τ , the singular points of the (homogeneous) differential equation correspond to the cusps of the modular curve. We are therefore particularly interested in asymptotic expansions close to a cusp. We start by discussing q-expansions of Eisenstein series, and we comment on iterated integrals at the end of the section. We already know that every modular form admits a q-expansion around the infinite cusp at τ = i∞ (or q N = 0). For Eisenstein series we can write down the q-expansion around the infinite cusp in closed form, cf. eq. (2.13). We now discuss how we can obtain the q-expansion around other cusps. We note that this is well-known in the mathematics literature, but we feel that it is important to discuss this topic here for completeness.
Consider c ∈ Q and an Eisenstein series f ∈ E n (Γ(N )). Since the h (n) N,r,s (τ ) form a spanning set for E n (Γ(N )), it is sufficient to consider the case where f is an element of this set. 7 We have seen in section 2 that there is γ c ∈ SL(2, Z) such that c = γ c ·i∞. Hence, expanding h Using this approach, we can obtain asymptotic expansions close to every cusp (e.g., close to every singular point of the differential equation satisfied by the Feynman integral). The strategy for obtaining asymptotic expansions of iterated Eisenstein integrals follows exactly the same lines: we can work out the modular transformation properties of the integrals using the steps described in section 5 to map any c ∈ Q to the infinite cusp, where we can easily obtain asymptotic expansions. Fast converging expansions are also the basis for the numerical evaluation of Eisenstein series and their iterated integrals. We turn to this topic in the next section.

Numerical evaluation
In this section we discuss how to evaluate Eisenstein series and their iterated integrals numerically in an efficient way. We focus only on Eisenstein series, but all statements can easily be transposed to iterated Eisenstein integrals.
Consider an Eisenstein series h (n) N,r,s (τ ). This function admits a q-expansion, and we can easily write down an arbitrary number of terms in the expansion, cf. eq. (2.13). Whenever Im τ > 0, we have |q N | < 1, and the series is convergent. However, the convergence may be rather slow depending on the value of q N . For example, if τ ∼ i∞ (and Re τ is not too large), we have |q N | 1 and the series converges extremely fast, while, when Im τ is not large, the convergence of the series will in general be much slower. We illustrate this in the second column of table 1, where we show the numerical value obtained for h (4) 2,1,1 (τ ) at τ = 1 2 + i 10 after truncating the q-expansion around the infinite cusp after a certain number of terms. Te expansion parameter is e iπτ i 0.7304 . . ..
The convergence of the series can be considerably accelerated by using modular transformations. Indeed, for every τ ∈ H, there is γ ∈ SL(2, Z) such that γ · τ lies in the fundamental domain for SL(2, Z) (cf. eq. (2.5)). If this is the case, the imaginary part of τ has a lower bound, i.e. Im τ ≥ √ 3 2 and, moreover, it is not possible to increase the imaginary part of τ further by modular transformations. While this does not constitute a proof, because of this argument we expect this representation to be provide in most cases a fast converging series representation.
In appendix A we review how this matrix can be constructed. In the case of our example, we find 1 2

JHEP02(2020)105
where the matrix γ 1/2 = ( 1 0 2 1 ) of section 5. Using the transformation in eq. (6.1) and the expansion in eq. (6.2), we obtain a series representation with expansion parameter e iπτ /(1−2τ ) −i 0.000388 . . .. The convergence of this series is extremely fast, as can be seen from the third column in table 1.
Using this strategy, we can derive fast converging series representations for Eisenstein series and their iterated integrals for all points in the upper half-plane. There is, however, one potential issue which we need to address. Our strategy involves using modular transformations in order to map a point in H to the fundamental domain before expansion. In section 5 we have argued that modular transformations introduce new transcendental constants, and we have argued that these constants are MMVs. Thus, if our strategy for efficient numerical evaluation is supposed to work, we need to be able to obtain numerical approximations for the new transcendental constants. Since all the MMVs relevant here are iterated Eisenstein integrals evaluated at τ = 0, we can use the same strategy to obtain numerical values for these MMVs. At this point, however, we need to address a critical point: since the iterated integrals are evaluated at τ = 0, a naive q-expansion would have q N = 1 as expansion parameter, and the q-series converges very badly. We can use modular transformations to accelerate the convergence, but this procedure may introduce new MMVs that we need to evaluate. We thus seem to be stuck in an infinite loop. In the following we discuss how we can circumvent this problem and evaluate all MMVs we need without introducing any new transcendental constants.
The basic idea is very simple. We can split the path of integration [0, i∞] into the segments [0, i] and [i, i∞]. The path [0, i] can be mapped to [i∞, i] using S, because S · 0 = i∞ and S · i = i. Hence, using the path composition and reversal formulae, as well as the algorithm from the previous section, we can express every MMV as a linear combination of iterated integrals of modular forms evaluated at τ = i. The expansion parameter for the q-expansion is then q N = e −2π/N < 1, and we obtain a fast-converging series representation of the MMV. After these steps all paths extend from i∞ to i, and we do not introduce any new transcendental constants. We can use this procedure to evaluate all MMVs we need to high precision. We have checked that we reproduce all results for MMVs of level one of small weight and length of ref. [67]. The precision we can reach is only limited by the number of terms in the expansion one wishes to include. We have also computed all MMVs up to weight six associated to Eisenstein series for Γ(2), as well as all MMVs associated to Eisenstein series for Γ(6) through weight five. In all cases, we were able to obtain at least 100 digits. Using the PSLQ algorithm, we find that all of these MMVs can be expressed in terms of the transcendental constants related to MPLs evaluated at the sixth root of unity, except for two and three new constants of weight four and five respectively. 7 Application to the elliptic 2 F 1 function In the remainder of this paper we illustrate the concepts of the previous sections on several examples. We start in this section by returning to the elliptic 2 F 1 function of section 3.2, before we discuss examples of Feynman integrals that involve iterated Eisenstein integrals in the next section.

JHEP02(2020)105
7.1 The imaginary part of the elliptic 2 F 1 function So far we have only considered the integrals in eq. (3.10) for 0 < z < 1, where the integral is real. We have already mentioned that for z > 1 the integral develops an imaginary part (we assume that all branches are fixed by assigning a small positive imaginary part to z). In this section we show how we can use modular transformations of iterated Eisenstein integrals to make this imaginary part explicit.
The idea we will follow is very simple. We have seen in eq. (5.4) that for z > 1 we can write τ = γ 3 ·τ 3 , with γ 3 ∈ SL(2, Z) and τ 3 purely imaginary. We can then write all integrals in terms of iterated Eisenstein integrals evaluated at τ 3 . If we work with the parity-invariant spanning set of Eisenstein series of section 2.4, all iterated Eisenstein integral are real for z > 1. Hence, we can explicitly isolate the imaginary parts of the integrals.
In the following we only discuss the functions U i (τ ) defined in eq. (4.13). The integrals T 1 and T 2 can be recovered by multiplying with the matrix S(z) in eq. (3.14), which does not involve any iterated Eisenstein integrals. For the functions U i (τ ), we find where the matrix M (γ −1 3 , τ 3 ) was defined in eq. (3.22). The real parts are given bỹ 3) Since τ 3 is purely imaginary for z > 1, the functionsŨ i andṼ i are manifestly real in this region. Therefore this is the correct separation of U i (τ ) into its real and imaginary parts. Let us make some comments about these results. First we note that the functions are still pure functions of uniform weight, i.e., the property of uniform weight is preserved. The functions U i , however, are no longer pure after analytic continuation: indeed, the matrix M (γ 3 , τ 3 ) in eq. (7.1) contains τ 3 in the denominator, which spoils purity. However, all the JHEP02(2020)105 terms that spoil purity are contained in the matrix M (γ 3 , τ 3 ). This can be explained as follows: the separation into the matrix S(z) and the pure functions U i (τ ) in eq. (3.13) relies on a choice, namely the choice of two periods for the elliptic curve y 2 = x(1 − x)(1 − zx) associated to our problem, or equivalently the choice of τ in eq. (3.16). Of course, we could have chosen a different basis of periods, and any two choices are related by an SL(2, Z) transformation. The matrices S(z) obtained from two different choices are related by the matrix M (γ 3 , τ 3 ), cf. eq. (3.21). Since nothing should depend on this choice, the matrix M (γ 3 , τ 3 ) in eq. (3.21) must be cancelled by a similar contribution coming from the modular transformation of the iterated Eisenstein integrals, so that .

Numerical evaluation of the elliptic 2 F 1 function
In the previous section we have seen how to analytically continue the master integrals for the family of elliptic 2 F 1 functions to all real values of z. After analytic continuation all imaginary parts are explicit, but the resulting expressions like in eq. (7.1) are not necessarily in a form where we can evaluate all the functions in a fast and stable way. Indeed, we have seen in section 6.2 that when τ has a small imaginary part, the q-expansions converge very slowly. Alternatively, for every (real) value of z, we can find γ z ∈ SL(2, Z) such that γ −1 z · τ (z) lies in the fundamental domain, where τ (z) was defined in eq. (5.6). A priori, the matrix γ z depends on z. We need only six distinct values of γ z to cover all real values of z: where we defined It is easy to check that the variables τ jα ≡ γ jα · τ (z), j ∈ {1, 2, 3} and α ∈ {a, b}, are purely imaginary and Im τ jα ≥ i. In each of the six regions in eq. (7.5) we obtain a representation of the master integrals T i in eq. (3.12) in terms of iterated Eisenstein integrals evaluated at τ jα , which admit fast converging q-expansions. The corresponding representations of T 1 and T 2 are shown in appendix B. As an illustration of the numerical evaluation, we show in figure 1 a plot of the real and imaginary parts of the first three coefficients in the Laurent expansion of T 1 . Let us conclude this section by commenting on the six regions for z that appear in eq. (7.5). The values {0, 1} may not come as a surprise, because they are related to the singular points of the integral. The values {−1, 1/2, 2} need more explanation, because they do not seem to be directly related to any apparent special point of the integral. In order to understand the appearance of these points, we need to analyse the j-invariant of the elliptic curve y 2 = x(1 − x)(1 − zx) associated to our problem. The periods and τ contain redundant information, but the j-invariant is a complex number that uniquely characterises an elliptic curve. The j-invariant can be computed from τ . It is invariant under SL(2, Z) (i.e., it is a modular function for SL(2, Z)), and it can be expressed as a rational function of Eisenstein series of level one:

JHEP02(2020)105
1,0,0 (τ ) 3 20 h (7.7) Since j(τ ) is a modular function for SL(2, Z), it must be a rational function of the Hauptmodul z = λ(τ ). We have We now see that the poles of the j-invariant are z ∈ {0, 1}, while z ∈ {−1, 1/2, 2} corresponds to j(τ ) = 1! The corresponding values of τ are We quote this here merely as an observation, we will see a similar pattern in the next section in the case of elliptic Feynman integrals, where the different regions will be determined by the points where the j-invariant becomes infinite, zero or one.

Asymptotic expansions
Let us conclude this section with an illustration how we can obtain asymptotic expansions of the elliptic 2 F 1 functions order by order in . In section 6.1 we have shown how we can

JHEP02(2020)105
This expansion is the same as the one obtained from the usual series representation of Gauss' hypergeometric function, (a) n (b) n (c) n z n n! . (7.13) 8 Application to the sunrise integral As a further application of the ideas described in the previous sections, we would like to perform the analytic continuation for the class of Feynman integrals required to compute the sunrise graph. In order to keep the formulas as compact as possible, we will limit ourselves to consider the relevant master integrals in d = 2 dimensions, where the master integrals can be chosen to be finite. We recall that the value in d = 4 dimensions can be obtained by a straightforward application of the well-known formulas for the dimensional shift of Feynman integrals [68,69]. We will start with a short recap on the computation of the sunrise graph, which will also allow us to introduce our notation. We consider the family of two-loop integrals defined by the formula S a 1 ,...,a 5 (p 2 , m 2 ; d) where the integration measure is chosen as Using IBP identities [70,71], one can show that the family of integrals defined in eq. (8.1) can be reduced to three master integrals. We choose one of the master integrals to be the tadpole, which is one in our normalisation, S 2,2,0,0,0 (p 2 , m 2 ; 2 − 2 ) = 1 .
Following refs. [26,59], we choose the other two master integrals as follows where we defined the dimensionless variable t = p 2 /m 2 and rescaled by the appropriate powers of m in order to render S 1 ( ; t) and S 2 ( ; t) dimensionless as → 0. Note that both integrals are finite in d = 2 dimensions. From now on, for simplicity, we will set m = 1, as the dependence of the master integrals on m can be always be recovered from simple dimensional analysis.

JHEP02(2020)105
As it is well known, the two non-trivial master integrals of the two-loop massive sunrise graph satisfy a system of two coupled differential equations. With this choice of basis, the equations in d = 2 space-time dimensions read [26] ∂ t S 1 (t) , (8.6) and where we defined S j (t) = S j (0; t) with j = 1, 2. Equation (8.5) is a coupled 2×2 system of equations whose solution requires one to first find a full solution of the corresponding homogeneous system, i.e., a 2 × 2 matrix W S (t) such that As shown for the first time in ref. [24], and generalised in ref. [72], a particular solution can always be found by studying the maximal cut of the relevant graphs. In our case, we can use the first-order differential operator to write W s (t) in the region 0 < t < 1 as are the two periods of the associated elliptic curve, with t ij (t) = t i (t) − t j (t) and the four branching points are At this point, we can define a new set of master integrals T 1 (t), T 2 (t) as follows , (8.12) which by definition fulfil the purely non-homogeneous system of differential equations . (8.13)

JHEP02(2020)105
It is by now well known that a solution of these differential equations can be found in terms of iterated integrals of Eisenstein series for Γ 1 (6) [53,59]. The modular parameter τ is as usual defined as the ratio of the two periods, and we find that the j-invariant of the family of elliptic curves associated to the maximal cut of the sunrise graph reads It turns out that eq. (8.14) can be inverted to give where η(τ ) was defined in eq. (4.7). The function t(τ ) is a Hauptmodul for Γ 1 (6), and so it is invariant under the corresponding modular transformations. Ψ 1 (t), instead, is a modular form of weight n = 1 for Γ 1 (6). Starting from these observations, one can prove that a basis for M n (Γ 1 (6)) is given by the functions [51] f n,p (τ ) = Ψ 1 (t(τ )) n t(τ ) p , 0 ≤ p ≤ n , (8.17) where we assume f 0,0 (τ ) = 1. Finally we introduce the following notation for iterated integrals of modular forms for Γ 1 (6), With these definitions it is immediate to see that as long as 0 < t < 1, the sunrise graph in d = 2 dimensions can be written in terms of iterated integrals of modular forms for Γ 1 (6): where is the Clausen function. We stress that eq. (8.19) is only valid in the region 0 < t < 1, which corresponds to 0 < s < m 2 . For these values of the external invariants the two master integrals of the sunrise graph S 1 , S 2 are real, which remains true as long as s < 9m 2 (including negative values of s). On the other hand, the integrals have a branch cut starting at s = 9m 2 , and they develop an imaginary part for s > 9m 2 . Note that this is in general not true for T 1 and T 2 separately. Before we discuss how to rewrite this solution for different values of s, we will express it in terms of iterated integrals over the parity-invariant spanning set of modular forms for JHEP02(2020)105 Γ(6) defined in section 2.4. By now it should be clear that this is the natural starting point to be able to analytically continue our expression to the whole phase space. Our goal will be to perform modular transformations in order to remap τ to the fundamental domain, whenever needed. Once that is done, we can evaluate the relevant iterated integrals by very fast converging q-expansions. Equation (8.19) is written in terms of modular forms for Γ 1 (6), which is not a normal subgroup of SL(2, Z) and, as it was discussed in section 2.2, we do not expect this set of iterated integrals to be closed under SL(2, Z) transformations. Γ(6), instead, is a normal subgroup of SL(2, Z) and it is therefore the natural arena where we expect our solution to be defined for general values of s. Moreover, by using the parityinvariant spanning set we are also able to make all real and imaginary parts explicit, at least as long as τ is purely imaginary.
It is straightforward to match eq. (8.19) to the basis defined in eq. (2.24) by matching the corresponding q-expansions. We obtain 8 where we have added the superscript (0, 1) to stress that this solution is defined for 0 < t < 1. By inspecting the result, we see that all terms in eq. (8.21) are purely real or imaginary, which guarantees that the master integrals S 1 and S 2 obtained through eq. (8.12) are real, as expected. We are now ready to focus on the analytic continuation and numerical evaluation of eq. (8.21) for all values of t.

Analytic continuation
In the rest of this section we will study the behaviour of the solution in eq. (8.21) for different values of t, in order to produce suitable analytical representations for the two integrals, which are valid across the whole phase-space and can be evaluated numerically by fast converging (q-)series expansions.
Our main goal will be to remap τ to the fundamental domain in each region. Before doing so, we show how the homogeneous solutions in eq. (8.10), and therefore τ , can be extended from the region 0 < t < 1 to arbitrary (real) values of t. Indeed, the analytical expression for τ itself depends on the homogeneous solutions of the differential equation (8.5), which have been given explicitly in eq. (8.10) for 0 < t < 1. Since eq. (8.5) has regular singular points at t ∈ {0, 1, 9, ±∞}, in general we expect eq. (8.10) to define local solutions to the differential equations which are only valid in some region, in our case the region 0 < t < 1. Hence, our first goal will be to find an analytic continuation of the local homogeneous solutions to all values of t. This can be achieved by finding suitable auxiliary homogeneous solutions which are well defined in the three remaining regions t < 0, 1 < t < 9 and t > 9, and then appropriately matching them at the branching points with the solution determined in eq. (8.10). Once the continued homogenous solutions Ψ j have been determined, the complete matrix of solutions W S (t) of eq. (8.9) can be extended from the region 0 < t < 1 to all real values of t by simply using the piecewise definition of Ψ j (t) of eq. (8.24) in eq. (8.9). In the same way, we can extend the definition of τ (t) in eq. (8.14) to all values of t.
With our correctly continued homogeneous solutions, we can study when we need to perform a modular transformation to remap τ to the fundamental domain and, in this way, accelerate the convergence of the q-expansion of the corresponding iterated Eisenstein integrals. It is useful to start by plotting τ as function of t, see figure 2. It is easy to see that our definition for τ almost never lies in the fundamental domain, except for a small region around the point t = 0. As we have seen in the case of the 2 F 1 hypergeometric function, one can find a fast-converging representation in the region a < t < b via a modular transformation such that τ is mapped to the fundamental domain. In order to achieve this in general we must use different modular transformations, depending on the value of t. We find in particular that we need to consider 8 different transformations γ t ∈ SL(2, Z), such that γ −1 t · τ (t) lies in the fundamental domain. Defining the four special points we find with (8.28) Recalling the expression for the j-invariant in eq. (8.15), we see that the apparently random special points t = r i , i − 1, . . . , 5 where the γ t ∈ SL(2, Z) is discontinuous, are always points t i such that j(t i ) ∈ {0, 1, ∞}, similar to what was observed for the 2 F 1 function. It is important to stress that these regions do not always overlap with the regions delimited by the regular singular points t i ∈ {0, 1, 9, ∞}, in particular the region 1a lies across the point t = 0. This is not a problem. It just tells us that as long as t varies in the small region r 2 < t < r 3 , i.e. −0.0397 . . . ≤ t < 0.0167 . . ., the relevant solutions (either the ones valid in the range (0, 1) or the ones valid in the range (−∞, 0)) can be taken as they are, without the need of remapping the corresponding τ to the fundamental domain. By defining τ j = γ −1 j · τ , we can plot the remapped value in function of t and check that it indeed always lies in the fundamental domain, see figure 3.
With the matrices defined in eq. (8.27), we can now perform the relevant transformations on the iterated integrals in eq. (8.21) and obtain new combinations of integrals evaluated in the fundamental domain. Clearly, there is nothing to do in region 1a, as γ 1a = 1. The first non-trivial transformation is in region 1b, i.e. for r 3 < t < 1 where the transformation induced by γ 1b ∈ SL(2, Z) has to applied in order to map τ back to the fundamental domain. We define τ = γ 1b · τ 1b , and apply the algorithms described in this paper to obtain where now τ 1b lies in the fundamental domain. The analytic expressions for the master integrals in the remaining regions can be found in appendix C. We can then use these expressions to plot the result for the original master integrals S 1 and S 2 as a function of t. Since the result is now always evaluated for τ (t) in the fundamental domain, we can get very precise numerical evaluations by considering only the first few terms of the corresponding q-series expansions. For example, by keeping 10 JHEP02(2020)105 Re(S1) Im(S1) (a) The master integral S 1 .  terms of the expansion and running for a few seconds on a laptop computer we can get the plots in figure 4. Our numerical results have been checked thoroughly against the known analytical results in ref. [26] and the numerical program SecDec [73,74]. As a concluding remark, we want to point out that all techniques used here can be straightforwardly applied to analytically continue and numerically evaluate other Feynman diagrams that can be expressed in terms of the same class of functions. Obvious examples are given by the so-called Kite integral [26,57] and the three-loop equal-mass banana graph [55,59,75].

Conclusion
In the last years it has started to become clear that multiloop Feynman integrals in dimensional regularisation can be naturally expressed in terms of iterated integrals of sets of differential forms defined on increasingly complicated complex (hyper)-surfaces. The simplest case is represented by multiple polylogarithms, which are defined as iterated integrals of rational functions on the Riemann sphere. The next-to-simplest case is instead constituted by their generalisation to genus-one Riemann surfaces -tori, or equivalently to elliptic curves -and have been dubbed elliptic multiple polylogarithms. A subset of the latter, which has found wide application in a multitude of physically interesting problems which depend on one single ratio of kinematical invariants, are closely related to iterated integrals of Eisenstein series. While a lot is known about the geometrical aspects of the construction of eMPLs and iterated Eisenstein integrals, their use in physically realistic problems requires being able to control the details of their branch cut structure (i.e. their analytic continuation) and to devise strategies to evaluate them numerically with high-precision for any values of their arguments. This problem can be solved in full generality for the simple case of MPLs, but till today still a little was known about how to treat these new classes of functions. In this paper a first step in this direction has been made. We have presented a set of algorithms which allow one to analytically continue iterated integrals of Eisenstein series to any region of the phase-space, in a way that makes it straightforward to evaluate them numerically with very high-precision through fast converging q-series expansions. We have shown the techniques in detail by applying them to a mathematical problem, i.e. the Taylor expansion JHEP02(2020)105 of a class of hypergeometric functions, and to the physical problem of the calculation of the two-loop massive sunrise graph. We expect generalisations of the algorithms described here to be able to address the evaluation of more general eMPLs.
In this appendix we review some useful algorithms related to the modular group SL(2, Z). We start by an algorithm to decompose any element of SL(2, Z) into a finite product of the generators S, T and −1. In a second subsection we discuss, for some τ ∈ H, how to find γ ∈ SL(2, Z) such that γ · τ lies in the fundamental domain for SL(2, Z). After this step, we can assume that |γ 11 | > |γ 21 |.
-Step 3: the matrix now has the form γ = In this appendix we describe how to find for every τ ∈ H an element γ ∈ SL(2, Z) such that γ · τ lies in the fundamental domain F defined in eq. (2.5). We will construct γ iteratively using the following algorithm: -Step 1: we can easily find an integer N such that − 1 2 ≤ Re τ + N < 1 2 . Replace τ by τ −→ T N · τ = τ + N . If after this replacement τ ∈ F , we are done. Otherwise we return to Step 1. At each iteration the imaginary part of τ is strictly increasing, and therefore the algorithm will eventually terminate.