aainstitutetext: Theoretical Physics Department, CERN, Geneva,
Switzerlandbbinstitutetext: Rudolf Peierls Centre for Theoretical Physics,
Clarendon Laboratory, Parks Road, Oxford OX1 3PU, UK

# Algorithms and tools for iterated Eisenstein integrals

Claude Duhr b Lorenzo Tancredi<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

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.

###### Keywords:

Elliptic polylogarithms, Eisenstein series, Feynman integrals.

††preprint: CERN-TH-2019-209, OUTP-19-14P

##  1 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 Nielsen ; Kummer ,
their re-discovery in the high-energy physics community Remiddi:1999ew ;
Gehrmann:2000zt  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 Goncharov:2010jf ; GoncharovMixedTate ;
Duhr:2011zq ; Duhr:2012fh , 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.

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
Gehrmann:2001pz ; Gehrmann:2001jv ; Buehler:2011ev ; Vollinga:2004sn ;
Frellesvig:2016ske ; Ablinger:2018sat . 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 re-expressing 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 Gehrmann:2014bfa ; Gehrmann:2015ora ; vonManteuffel:2015msa . 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 Sabry . 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 $t\bar{t}$
production at hadron colliders Czakon:2013goa . 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 high-energy physics community Broadhurst:1987ei ;
Bauberger:1994by ; Bauberger:1994hx ; MullerStach:2011ru ; Laporta:2004rb ;
Caffo:2008aw ; Remiddi:2016gno ; Broedel:2017siw ; Remiddi:2013joa ;
Adams:2015ydq ; Caffo:1998du ; Adams:2014vja ; Adams:2013kgc ; Adams:2013nia ;
Adams:2015gva ; Bloch:2016izu ; Bloch:2013tra ; Kniehl:2005bc ; Bogner:2019lfa
. 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 LevinRacinet ;
BrownLevin  and had later found application in the calculation of one-loop
scattering amplitudes in string theory Broedel:2014vla ; Broedel:2018izr ;
Broedel:2015hia ; Broedel:2017jdo . 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 Broedel:2017kkb ;
Broedel:2017siw ; Broedel:2018iwv , which has made it possible to naturally
compute different previously out-of-reach multi-loop Feynman integrals
Broedel:2018qkq ; Blumlein:2018jgc ; Broedel:2019hyg . 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,111And also of
other similar problems depending on one single independent variable. a special
subclass of eMPLs has been shown to be sufficient, i.e. the so-called iterated
Eisenstein integrals Adams:2018ulb ; Broedel:2018rwm ; ZagierModular ;
Adams:2017ejb ; ManinModular .

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.222For
specific instances of elliptic Feynman integrals fast numerical codes were
achieved, see for example refs. Remiddi:2016gno ; Primo:2017ipr ;
vonManteuffel:2017hms ; Bogner:2017vim ; Frellesvig:2019byn . 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 Adams:2017ejb ;
Broedel:2018iwv , as well as the three-loop banana graph Broedel:2019kmn .
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.

##  2 Modular forms and Eisenstein series

The purpose of this section is to review in detail the modular group
$\operatorname{SL}(2,\mathbb{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.

###  2.1 The modular group and its congruence subgroups

The group $\operatorname{SL}(2,\mathbb{Z})$ consists of all $2\times 2$
matrices with integer entries and unit determinant. It is generated by the
three elements

$S=S^{-1}=\left(\begin{smallmatrix}0&-1\\\
1&0\end{smallmatrix}\right)\,,\qquad T=\left(\begin{smallmatrix}1&1\\\
0&1\end{smallmatrix}\right)\,,\qquad-\mathds{1}=\left(\begin{smallmatrix}-1&0\\\
0&-1\end{smallmatrix}\right)\,.$ (1)

In other words, every element of $\operatorname{SL}(2,\mathbb{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 $\operatorname{SL}(2,\mathbb{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
$\operatorname{SL}(2,\mathbb{Z})$, but only in a certain subgroup thereof. Of
particular interest are the so-called _congruence subgroups_. A subgroup
$\Gamma\subseteq\operatorname{SL}(2,\mathbb{Z})$ is called a congruence
subgroup of level $N$ if it contains the principle congruence subgroup of
level $N$,

$\Gamma(N)=\left\\{\left(\begin{smallmatrix}a&b\\\
c&d\end{smallmatrix}\right)\in\operatorname{SL}(2,\mathbb{Z}):a,d=1\\!\\!\\!\mod
N\textrm{ and }b,c=0\\!\\!\\!\mod N\right\\}\,.$ (2)

Particularly important examples of congruence subgroups are

$\begin{split}\Gamma_{0}(N)&\,=\left\\{\left(\begin{smallmatrix}a&b\\\
c&d\end{smallmatrix}\right)\in\operatorname{SL}(2,\mathbb{Z}):c=0\\!\\!\\!\mod
N\right\\}\,,\\\ \Gamma_{1}(N)&\,=\left\\{\left(\begin{smallmatrix}a&b\\\
c&d\end{smallmatrix}\right)\in\operatorname{SL}(2,\mathbb{Z}):a,d=1\\!\\!\\!\mod
N\textrm{ and }c=0\\!\\!\\!\mod N\right\\}\,.\end{split}$ (3)

For the rest of this section, $\Gamma$ will denote some congruence subgroup of
$\operatorname{SL}(2,\mathbb{Z})$ (which may be
$\operatorname{SL}(2,\mathbb{Z})=\Gamma(1)$ itself). $\Gamma$ naturally acts
on the upper half-plane $\mathbb{H}=\\{\tau\in\mathbb{C}|\textrm{Im
}\tau>0\\}$ via Möbius transformations,

$\gamma\cdot\tau=\frac{a\tau+b}{c\tau+d}\,,\qquad\gamma=\left(\begin{smallmatrix}a&b\\\
c&d\end{smallmatrix}\right)\in\Gamma\,.$ (4)

Note that $-\mathds{1}$ acts trivially on $\tau$, i.e.,
$(-\mathds{1})\cdot\tau=\tau$. The action in eq. (4) partitions the upper
half-plane into distinct orbits. The space $Y_{\Gamma}={\mathbb{H}}/\Gamma$ of
orbits is called the _modular curve_ for $\Gamma$. It can be useful to
consider a connected domain in ${\mathbb{H}}$ whose points correspond to the
distinct orbits in $Y_{\Gamma}$, called a _fundamental domain_ for $\Gamma$. A
particularly important example of a fundamental domain is the one for the full
modular group $\operatorname{SL}(2,\mathbb{Z})$, which is

$\\!\\!\cal F=\\{\tau\in{\mathbb{H}}:{-{1}/{2}}\leq\textrm{Re
}\tau<{1}/{2}\textrm{~{}and~{}}|\tau|>1\\}\cup\\{\tau\in{\mathbb{H}}:\textrm{Re
}\tau\leq 0\textrm{~{}and~{}}|\tau|=1\\}\,.$ (5)

In other words, for every $\tau\in{\mathbb{H}}$ there are unique
$\gamma\in\operatorname{SL}(2,\mathbb{Z})$ and $\tau_{0}\in\cal F$ such that
$\tau=\gamma\cdot\tau_{0}$.

$\Gamma$ also acts on $\mathbb{Q}\cup\\{i\infty\\}$, and so it is natural to
consider its action on the extended upper half-plane
$\overline{\mathbb{H}}=\mathbb{H}\cup\mathbb{Q}\cup\\{i\infty\\}$. It
partitions $\mathbb{Q}\cup\\{i\infty\\}$ into a finite number of distinct
equivalence classes, called the _cusps_ of $\Gamma$. For every
$c_{1},c_{2}\in\mathbb{Q}\cup\\{i\infty\\}$, there is a
$\gamma\in\operatorname{SL}(2,\mathbb{Z})$ such that $c_{2}=\gamma\cdot
c_{1}$, which implies that $\operatorname{SL}(2,\mathbb{Z})$ has only one
cusp. The matrix $\gamma$ can be constructed explicitly: First, we note that
it is sufficient to consider $c_{1}=i\infty$. Indeed if
$c_{1},c_{2}\in\mathbb{Q}$ such that $c_{k}=\gamma_{c_{k}}\cdot i\infty$, then
$c_{2}=\gamma_{c_{2}}\gamma_{c_{1}}^{-1}\cdot c_{1}$. Next, let
$c_{2}=m/n\in\mathbb{Q}$. It follows from the properties of the greatest
common divisor that there are integers $p$ and $q$ such
$\textrm{gcd}(m,n)=mp+nq$.333The integers $p,q$ can easily be obtained, e.g.,
from Mathematica’s implementation GCD of the greatest common divisor. If $m$
and $n$ are relatively prime, $\textrm{gcd}(m,n)=1$, we can take
$\gamma_{c_{2}}=\left(\begin{smallmatrix}m&mr-q\\\
n&nr+p\end{smallmatrix}\right)$, with $r$ an arbitrary integer. This matrix
obviously satisfies $c_{2}=m/n=\gamma_{c_{2}}\cdot i\infty$. Moreover, this
matrix has unit determinant, and so it defines an element of
$\operatorname{SL}(2,\mathbb{Z})$.

###  2.2 Modular forms

We start by considering functions from the extended upper half-plane into the
complex numbers that are invariant under $\Gamma$. A _modular function_ for
$\Gamma$ is a meromorphic function $f:{\overline{\mathbb{H}}}\to\mathbb{C}$
such that $f(\gamma\cdot\tau)=f(\tau)$, for all $\gamma\in\Gamma$. One can
show that there are no non-constant holomorphic modular functions, i.e., every
non-constant function that is invariant under $\Gamma$ has at least one pole.
If we want to obtain holomorphic functions, we need to consider more general
transformation behaviours under $\Gamma$. For every positive integer $n$, we
define

$(f_{|n\gamma})(\tau)=(c\tau+d)^{-n}\,f(\gamma\cdot\tau)\,,\qquad\gamma=\left(\begin{smallmatrix}a&b\\\
c&d\end{smallmatrix}\right)\in\Gamma\,.$ (6)

This defines for every $n$ a genuine action of $\Gamma$ on functions on the
extended upper half-plane. In particular, we have

$(f_{|n\gamma_{1}})_{|n\gamma_{2}}=f_{|n(\gamma_{1}\gamma_{2})}\,.$ (7)

A _modular form of weight $n$ for $\Gamma$_ is a function
$f:\overline{\mathbb{H}}\to\mathbb{C}$ such that

  1. 1.

$f$ invariant under the action of $\Gamma$ in eq. (6), i.e., it satisfies
$(f_{|n\gamma})(\tau)=f(\tau)$.

  2. 2.

$f$ is holomorphic on $\mathbb{H}$.

  3. 3.

$f_{|n\gamma}$ is holomorphic at $i\infty$ for all
$\gamma\in\operatorname{SL}(2,\mathbb{Z})$.

We can think of modular forms as holomorphic functions, i.e., functions
without poles, that are invariant under the action in eq. (6). Modular forms
of weight $n$ form a finite-dimensional vector space $\cal M_{n}(\Gamma)$.
Moreover, $\cal M_{n}(\Gamma)$ 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 $\Gamma$ is a congruence subgroup of level $N$ and $f\in\cal
M_{n}(\Gamma)$, then we say that $f$ has level $N$. Every modular form of
level $N$ is invariant under $T^{N}=\left(\begin{smallmatrix}1&N\\\
0&1\end{smallmatrix}\right)\in\Gamma(N)$. This matrix acts on
$\overline{\mathbb{H}}$ via translations by $N$, $T^{N}\cdot\tau=\tau+N$, and
so every modular form of level $N$ is a periodic function of period $N$ and
admits a Fourier expansion of the form

$f(\tau)=\sum_{n\geq 0}a_{n}\,q_{N}^{n}\,,\qquad q_{N}=e^{2\pi i\tau/N}\,.$
(8)

Having identified modular forms for $\Gamma$ as invariants under the group
action in eq. (6), it is natural to ask how a modular form $f\in\cal
M_{n}(\Gamma)$ transforms for some $\gamma\in\operatorname{SL}(2,\mathbb{Z})$
that is not an element of $\Gamma$.

We start by assuming that $\Gamma\subseteq\operatorname{SL}(2,\mathbb{Z})$ is
a normal subgroup.444A subgroup $H\subseteq G$ is normal if $gH=Hg$, for all
$g\in G$. Then, for every $\gamma\in\operatorname{SL}(2,\mathbb{Z})$ and
$\gamma_{1}\in\Gamma$, there is $\gamma_{2}\in\Gamma$ such that
$\gamma\gamma_{1}=\gamma_{2}\gamma$. We have

$(f_{|n\gamma})_{|n\gamma_{1}}=f_{|n(\gamma\gamma_{1})}=f_{|n(\gamma_{2}\gamma)}=(f_{|n\gamma_{2}})_{|n\gamma}=f_{|n\gamma}\,,$
(9)

where the last step follows from the fact that $f$ is a modular form for
$\Gamma$, and thus invariant. We conclude that if $f$ is a modular form for a
normal subgroup $\Gamma$, then $f_{|n\gamma}$ is also a modular form for
$\Gamma$ of the same weight, for all
$\gamma\in\operatorname{SL}(2,\mathbb{Z})$.

The previous argument relies crucially on $\Gamma$ being normal in
$\operatorname{SL}(2,\mathbb{Z})$. However, not all subgroups of
$\operatorname{SL}(2,\mathbb{Z})$ are normal, not even all congruence
subgroups. In particular, the congruence subgroups $\Gamma_{0}(N)$ and
$\Gamma_{1}(N)$ are in general not normal. The principle congruence subgroups
$\Gamma(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 $\Gamma$ of level $N$, then
$f_{|n\gamma}$ is a modular form for the principle congruence subgroup
$\Gamma(N)$ of the same weight and level, for all
$\gamma\in\operatorname{SL}(2,\mathbb{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\gamma}$ will in general not be invariant under $\Gamma$, but it is only
invariant under the smaller subgroup $\Gamma(N)$.

While the previous argument shows that a modular form for a congruence
subgroup $\Gamma$ of level $N$ in general transforms into a linear combination
of modular forms for $\Gamma(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.

###  2.3 Eisenstein series

The vector space $\cal M_{n}(\Gamma)$ of modular forms of weight $n$ for
$\Gamma$ admits a direct sum decomposition

$\cal M_{n}(\Gamma)=\cal E_{n}(\Gamma)\oplus\cal S_{n}(\Gamma)\,.$ (10)

Here $\cal S_{n}(\Gamma)$ denotes the subspace of cusp forms, i.e., the
subspace of modular forms of weight $n$ for $\Gamma$ that vanish on
$\mathbb{Q}\cup\\{i\infty\\}$. Its (orthogonal) complement is the Eisenstein
subspace $\cal E_{n}(\Gamma)$. 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 $\cal E_{n}(\Gamma)$ (cf., e.g.,
ref. diamond ). Our choice is motivated by the fact that Eisenstein series for
$\Gamma(N)$ are closely related to the Kronecker series and elliptic
polylogarithms BrownLevin ; Broedel:2015hia ; Broedel:2018iwv . Since
$\Gamma(N)\subseteq\Gamma$ for every congruence subgroup $\Gamma$ of level
$N$, we only discuss the case of the principle congruence subgroup
$\Gamma(N)$, and we refer to the literature for the other cases (see, e.g.,
ref. diamond ).

We start by writing down a spanning set for $\cal E_{n}(\Gamma(N))$
Broedel:2018iwv ,

$h_{N,r,s}^{(n)}(\tau)=-\sum_{\begin{subarray}{c}(a,b)\in\mathbb{Z}^{2}\\\
(a,b)\neq(0,0)\end{subarray}}\frac{e^{2\pi i(bs-ar)/N}}{(a\tau+b)^{n}}\,,$
(11)

where $n$ and $N$ are integers greater than unity and $r,s$ are integers
defined modulo $N$. These functions transform like eq. (6) and form a spanning
set for $\cal E_{n}(\Gamma(N))$, except for the case $n=2$ and
$(r,s)=(0,0)\mod N$, which is not a modular form. More generally, we have
Broedel:2018iwv ,

$h_{N,r,s}^{(n)}\left(\frac{a\tau+b}{c\tau+d}\right)=(c\tau+d)^{n}\,h^{(n)}_{N,rd+sb,rc+sa}(\tau)\,,\qquad\gamma=\left(\begin{smallmatrix}a&b\\\
c&d\end{smallmatrix}\right)\in\operatorname{SL}(2,\mathbb{Z})\,.$ (12)

The functions defined in eq. (11) admit a $q$-expansion that can be written
down in closed form

$h_{N,r,s}^{(n)}(\tau)=\sum_{k\geq 0}C_{N,r,s,k}^{(n)}\,q_{N}^{k}\,.$ (13)

The Fourier coefficients for $k\geq 1$ are given by (cf., e.g., ref. diamond
),

$C_{N,r,s,k}^{(n)}=-\frac{(2\pi
i)^{n}}{N^{n}(n-1)!}\sum_{(c_{1},c_{2})\in(\mathbb{Z}/N\mathbb{Z})^{2}}\sum_{\begin{subarray}{c}m|k\\\
k=mc_{2}\end{subarray}}\left[m^{n-1}\,e^{2\pi
i(rc_{2}-(s-m)c_{1})/N}-(m\leftrightarrow-m)\right]\,.$ (14)

The constant term for $n\geq 2$ is

$C_{N,r,s,0}^{(n)}=\frac{(2\pi i)^{n}}{n!}\,B_{n}\left(\frac{s}{N}\right)\,,$
(15)

where $B_{n}(x)$ are the Bernoulli polynomials,

$\frac{t\,e^{xt}}{e^{t}-1}=\sum_{n=0}^{\infty}\frac{t^{n}}{n!}\,B_{n}(x)\,.$
(16)

For $n=1$ the constant term is given by

$C_{N,r,s,0}^{(1)}=\left\\{\begin{array}[]{ll}2\pi
i\left(\frac{s}{N}-\frac{1}{2}\right)\,,&s\neq 0\\!\\!\\!\\!\mod N\,,\\\
0\,,&(r,s)=(0,0)\\!\\!\\!\\!\mod N\,,\\\ \pi\cot\frac{\pi
r}{N}\,,&\textrm{otherwise}\,.\end{array}\right.$ (17)

Note that all the Fourier coefficients of $h_{N,r,s}^{(n)}(\tau)$ are
algebraic multiples of $\pi^{n}$.

The functions $h_{N,r,s}^{(n)}$ form a spanning set of $\cal
E_{n}(\Gamma(N))$, but they are not linearly independent. It is known that for
$n\geq 3$ the dimension of $\cal E_{n}(\Gamma)$ is equal to the number of
cusps of $\Gamma$.555For $n=2$, the dimension is equal to the number of cusps
minus one, while for $n=1$ the dimension of $\cal E_{n}(\Gamma(N))$ is only
half the number of cusps. The only linear relations among the
$h_{N,r,s}^{(n)}$ for $n\geq 2$ are the reflection identity,

$h_{N,r,s}^{(n)}(\tau)=(-1)^{n}\,h_{N,-r,-s}^{(n)}(\tau)\,,$ (18)

and, for every $d$ divides $N$ and $0\leq\rho,\sigma<N$, the distribution
identity Broedel:2018iwv ,

$\sum_{\frac{1}{N}(r,s)\in\frac{1}{N}(\rho,\sigma)+\Lambda_{N/d}^{F}}h_{N,r,s}^{(n)}(\tau)=\left(\frac{d}{N}\right)^{n-2}\,h^{(n)}_{d,\rho,\sigma}(\tau)\,,$
(19)

with

$\Lambda_{N/d}^{F}=\\{(r,s)\in\mathbb{Z}^{2}:0\leq r,s<N\\}\,.$ (20)

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\geq 3$, it is possible to explicitly write down a subset of
the $h_{N,r,s}^{(n)}$ that are linearly independent Broedel:2018iwv .

Let us conclude this section by giving an example that if $f\in\cal
M_{n}(\Gamma)$ and $\gamma\in\operatorname{SL}(2,\mathbb{Z})$, then
$f_{|n\gamma}$ will in general not be in $\cal M_{n}(\Gamma)$, but only in
$\cal M_{n}(\Gamma(N))$.666We are grateful to Nils Matthes for suggesting this
counterexample. Consider the case $N=2$, and $\Gamma=\Gamma_{1}(2)$. The
vector space $\cal E_{2}(\Gamma(2))$ is two-dimensional, and a basis is
$\\{h^{(2)}_{2,1,0},h^{(2)}_{2,0,1}\\}$. The subspace $\cal
E_{2}(\Gamma_{1}(2))\subset\cal E_{2}(\Gamma(2))$ is one-dimensional with
basis $\\{h^{(2)}_{2,1,0}\\}$. Consider $S=\left(\begin{smallmatrix}0&-1\\\
1&0\end{smallmatrix}\right)\notin\Gamma_{1}(2)$. We have

$h^{(2)}_{2,1,0|6S}=h^{(2)}_{2,0,1}\,,\qquad
h^{(2)}_{2,0,1|6S}=h^{(2)}_{2,1,0}\,,$ (21)

where the $h^{(2)}_{2,r,s|6S}$ refer to $h^{(2)}_{2,r,s}$ after acting with
the modular transformation in eq. (6) with $n=6$ and $\gamma=S$. We see that
the subspace $\cal E_{2}(\Gamma_{1}(2))$ is not left invariant by $S$, because
$h^{(2)}_{2,1,0|6S}\notin\cal E_{2}(\Gamma_{1}(2))$.

###  2.4 A parity-invariant spanning set

In the previous section we have presented a spanning set of Eisenstein series
for $\Gamma(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 $\cal E_{n}(\Gamma(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 $\cal
E_{n}(\Gamma(N))$ that are explicitly real when $\tau$ is purely imaginary.

We start by analysing how $h_{N,r,s}^{(n)}(\tau)$ behaves under complex
conjugation when $\tau$ is purely imaginary. Writing $\tau=it$, with $t$ real
and positive, we have

$h_{N,r,s}^{(n)}(it)^{*}=h_{N,r,-s}^{(n)}(it)\,.$ (22)

Indeed, starting from the lattice sum in eq. (11), we find

$h_{N,r,s}^{(n)}(it)^{*}=-\sum_{\begin{subarray}{c}(a,b)\in\mathbb{Z}^{2}\\\
(a,b)\neq(0,0)\end{subarray}}\frac{e^{-2\pi i(bs-
ar)/N}}{(-ait+b)^{n}}=-\sum_{\begin{subarray}{c}(\tilde{a},\tilde{b})\in\mathbb{Z}^{2}\\\
(\tilde{a},\tilde{b})\neq(0,0)\end{subarray}}\frac{e^{-2\pi
i(\tilde{b}s+\tilde{a}r)/N}}{(\tilde{a}it+\tilde{b})^{n}}=h_{N,r,-s}^{(n)}(it)\,.$
(23)

We see that when $\tau$ is purely imaginary, complex conjugation amounts to
changing the sign of $s$. Based on this observation, we define,

$\begin{split}\small{\texteisfont{a}}^{(n)}_{N,r,s}(\tau)&\,=-\frac{1}{2}\left[h_{N,r,s}^{(n)}(\tau)+h_{N,r,-s}^{(n)}(\tau)\right]=-\frac{1}{2}\left[h_{N,r,s}^{(n)}(\tau)+(-1)^{n}h_{N,-r,s}^{(n)}(\tau)\right]\,,\\\
\small{\texteisfont{b}}^{(n)}_{N,r,s}(\tau)&\,=-\frac{1}{2i}\left[h_{N,r,s}^{(n)}(\tau)-h_{N,r,-s}^{(n)}(\tau)\right]=-\frac{1}{2i}\left[h_{N,r,s}^{(n)}(\tau)-(-1)^{n}h_{N,-r,s}^{(n)}(\tau)\right]\,,\end{split}$
(24)

where the second equality follows from the reflection identity in eq. (18).
These functions also admit a representation in terms of lattice sums similar
to eq. (11),

$\begin{split}\small{\texteisfont{a}}_{n,N,r,s}(\tau)&\,=\frac{1}{2}\sum_{\begin{subarray}{c}(a,b)\in\mathbb{Z}^{2}\\\
(a,b)\neq(0,0)\end{subarray}}\frac{e^{-2\pi iar/N}\cos\frac{2\pi
sb}{N}}{(a\tau+b)^{n}}\,,\\\
\small{\texteisfont{b}}_{n,N,r,s}(\tau)&\,=\frac{1}{2i}\sum_{\begin{subarray}{c}(a,b)\in\mathbb{Z}^{2}\\\
(a,b)\neq(0,0)\end{subarray}}\frac{e^{-2\pi iar/N}\sin\frac{2\pi
sb}{N}}{(a\tau+b)^{n}}\,.\end{split}$ (25)

Equation (22) implies that the functions
$\small{\texteisfont{a}}_{n,N,r,s}(\tau)$ and
$\small{\texteisfont{b}}_{n,N,r,s}(\tau)$ are real whenever $\tau$ is purely
imaginary. Since they are linear combinations of the Eisenstein series
$h_{N,r,s}^{(n)}(\tau)$, they also form a spanning set for $\cal
E_{n}(\Gamma(N))$.We can easily determine a basis of
$\small{\texteisfont{a}}_{n,N,r,s}(\tau)$ and
$\small{\texteisfont{b}}_{n,N,r,s}(\tau)$. In particular, for $n\geq 3$ a
basis for $\cal E_{n}(\Gamma(N))$ is

$B_{N,n}=B_{N,n}^{1}\cup B_{N,n}^{2}\,,$ (26)

with

$\begin{split}B_{N,n}^{1}&\,=\left\\{\small{\texteisfont{a}}_{n,N,r,s},\small{\texteisfont{b}}_{n,N,r,s}:0<r,s<N/2\textrm{
and }\textrm{gcd}(N,r,s)=1\right\\}\,.\end{split}$ (27)

The set $B_{N,n}^{2}$ depends on whether the weight $n$ is even or odd. If $n$
is even, we have

$\\!\\!B_{N,n}^{2}=\left\\{\small{\texteisfont{a}}_{n,N,k,0},\small{\texteisfont{a}}_{n,N,k,N/2},\small{\texteisfont{a}}_{n,N,0,k},\small{\texteisfont{a}}_{n,N,N/2,k}:0<k\leq
N/2\textrm{ and }\textrm{gcd}(N,k)=1\right\\}\,,$ (28)

while for $n$ odd, we have

$\\!\\!B_{N,n}^{2}=\left\\{\small{\texteisfont{a}}_{n,N,k,0},\small{\texteisfont{a}}_{n,N,k,N/2},\small{\texteisfont{b}}_{n,N,0,k},\small{\texteisfont{b}}_{n,N,N/2,k}:0<k\leq
N/2\textrm{ and }\textrm{gcd}(N,k)=1\right\\}\,.$ (29)

Note that when $N$ is odd, functions with index $N/2$ are absent.

##  3 Iterated integrals of modular forms

###  3.1 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}(\tau)$. We define differential forms
$\omega_{j}=\frac{d\tau}{2\pi i}\,h_{j}(\tau)$ and consider the iterated
integrals ManinModular ; Brown:mmv ,

$I(h_{1},\ldots,h_{k};\tau_{0},\tau)=\int_{\tau_{0}}^{\tau}\omega_{k}\ldots\omega_{1}=\int_{\tau_{0}}^{\tau}\frac{d\tau^{\prime}}{2\pi
i}\,h_{1}(\tau^{\prime})\,I(h_{2},\ldots,h_{k};\tau_{0},\tau^{\prime})\,.$
(30)

The recursion starts with $I(;\tau_{0},\tau)=1$. The number of integrations
$k$ is called the _length_ of the iterated integral. We can also define a
notion of weight Broedel:2018qkq . Assume that all modular forms have level
$N$ and that they admit $q$-expansions of the form

$h_{j}(\tau)=\sum_{n\geq 0}a_{n}(h_{j})\,q_{N}^{n}\,,\qquad h_{j}\in\cal
M_{n_{j}}(\Gamma)\,,$ (31)

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. (14), which are all algebraic multiples of powers of $\pi$). We define the
weight of $I(h_{1},\ldots,h_{k};\tau_{0},\tau)$ as $-k+\sum_{j=1}^{k}n_{j}$.

The integrals in eq. (30) have all the standard properties of iterated
integrals, cf. e.g., ref. ChenSymbol . In particular, they form a shuffle
algebra,

$I(h_{1},\ldots,h_{k};\tau_{0},\tau)I(h_{k+1},\ldots,h_{l};\tau_{0},\tau)=\sum_{\sigma\in\Sigma(k,l-k)}I(h_{\sigma_{1}},\ldots,h_{\sigma_{l}};\tau_{0},\tau)\,,$
(32)

where the sum runs over all shuffles of $(h_{1},\ldots,h_{k})$ and
$(h_{k+1},\ldots,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,

$\begin{split}I(h_{1},\ldots,h_{k};\tau_{0},\tau)&\,=\sum_{l=0}^{k}I(h_{l+1},\ldots,h_{k};\tau_{0},\tau_{1})I(h_{1},\ldots,h_{l};\tau_{1},\tau)\,,\\\
I(h_{1},\ldots,h_{k};\tau_{0},\tau)&\,=(-1)^{k}I(h_{k},\ldots,h_{1};\tau,\tau_{0})\,.\end{split}$
(33)

It is easy to see that $I(h_{1},\ldots,h_{k};\tau_{0},\tau)$ is linear in its
first $k$ arguments,

$I(\ldots,c_{1}\,h_{j,1}+c_{2}\,h_{j,2},\ldots;\tau_{0},\tau)=c_{1}\,I(\ldots,h_{j,1},\ldots;\tau_{0},\tau)+c_{2}\,I(\ldots,h_{j,2},\ldots;\tau_{0},\tau)\,,$
(34)

where the $c_{i}$ are constants. Finally, if all the modular forms $h_{j}$ are
real on the imaginary axis, then $I(h_{1},\ldots,h_{k};\tau_{0},\tau)$ is real
whenever $\tau$ and $\tau_{0}$ are purely imaginary.

The path composition formula in eq. (33) allows us to consider only a fixed
base point $\tau_{0}$, and from now on we will always choose
$\tau_{0}=i\infty$. We will use the shorthand notation
$I(h_{1},\ldots,h_{k};\tau)=I(h_{1},\ldots,h_{k};i\infty,\tau)$. 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 $\tau$ to $q_{N}$. We can write

$\omega_{j}=\frac{dq_{N}}{q_{N}}\,\frac{N\,a_{0}(h_{j})}{(2\pi i)^{2}}+\cal
O(q_{N}^{0})\,.$ (35)

Hence, if the constant term in the $q$-expansion is non zero, $\omega_{j}$ has
a pole at $q_{N}=0$, or equivalently $\tau=i\infty$. As a result, the integral
$I(h_{1},\ldots,h_{k};\tau)$ requires regularisation whenever
$a_{0}(h_{k})\neq 0$.

In ref. Brown:mmv  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. (30) as a linear combination of absolutely
convergent integrals multiplied by powers of $\tau$ Brown:mmv ,

$I(h_{1},\ldots,h_{k};\tau)=\sum_{p=0}^{n}\frac{a_{0}(h_{p+1})\ldots
a_{0}(h_{k})}{(k-p)!\,(2\pi
i)^{k-p}}\,\tau^{k-p}\int_{i\infty}^{\tau}R[\omega_{p}\ldots\omega_{1}]\,,$
(36)

where the map $R$ is defined by Brown:mmv

$R[\omega_{p}\ldots\omega_{1}]=\sum_{q=0}^{p}(-1)^{p-q}\,\omega_{q+1}^{\infty}\ldots\omega_{p}^{\infty}{\,\amalg\hskip-3.8pt\amalg\,}\omega_{q}\ldots\omega_{1}\,,$
(37)

where ${\,\amalg\hskip-3.8pt\amalg\,}$ denotes the shuffle product, and we
defined

$\omega_{j}^{\infty}=a_{0}(h_{j})\,\frac{d\tau}{2\pi
i}=\frac{dq_{N}}{q_{N}}\,\frac{N\,a_{0}(h_{j})}{(2\pi i)^{2}}\,.$ (38)

###  3.2 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

$\displaystyle
T(n_{1},n_{2},n_{3};z)=\int_{0}^{1}dx\,x^{-1/2+n_{1}+a\,\epsilon}(1-x)^{-1/2+n_{2}+b\epsilon}(1-z\,x)^{-1/2+n_{3}+c\,\epsilon}\,,$
(39)

where $n_{i}\in\mathbb{Z}$ and $a,b,c$ are complex numbers. These integrals
are closely related to a family of hypergeometric ${}_{2}F_{1}$ functions,

$\begin{split}T(n_{1},n_{2},n_{3};z)&\,=\frac{\Gamma(\frac{1}{2}+n_{1}+a\epsilon)\,\Gamma(\frac{1}{2}+n_{2}+b\epsilon)}{\Gamma(1+n_{1}+n_{2}+(a+b)\epsilon)}\\\
&\,\times{}_{2}F_{1}\left(\frac{1}{2}-n_{3}-c\epsilon,\frac{1}{2}+n_{1}+a\epsilon;1+n_{1}+n_{2}+(a+b);z\right)\,.\end{split}$
(40)

The integral in eq. (39) defines a meromorphic function of $\epsilon$, and we
interpret the integral as a Laurent series in $\epsilon$. In the following we
assume that $z$ is real. As long as $z<1$ the integral is real order by order
in $\epsilon$, 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. (39) can be expressed as a linear combination of two independent
master integrals, which can be chosen as

$T_{1}(z)=T(0,0,0;z)\textrm{~{}~{}and~{}~{}}\quad T_{2}(z)=T(1,0,0;z)\,.$ (41)

The square root in the integrand defines an elliptic curve via the polynomial
equation $y^{2}=x(1-x)(1-zx)$. The master integrals were evaluated in ref.
Broedel:2017kkb ; Broedel:2018iwv ; Broedel:2018rwm  order-by-order in
$\epsilon$ in terms of elliptic polylogarithms and iterated Eisenstein
integrals for the congruence subgroup $\Gamma(2)$. The result can be cast in
the form:

$\left(\begin{array}[]{c}T_{1}(z)\\\
T_{2}(z)\end{array}\right)=S(z)\left(\begin{array}[]{c}U_{1}(\tau)\\\
U_{2}(\tau)\end{array}\right)\,,$ (42)

where $S(z)$ denotes the matrix

$S(z)=\left(\begin{array}[]{cc}2\operatorname{K}(z)&0\\\
\frac{2}{z(1+6\epsilon)}\,\left[\operatorname{E}(z)-\frac{1}{3}(2-z)\operatorname{K}(z)\right]&-\frac{i\pi}{z(1+6\epsilon)\operatorname{K}(z)}\\\
\end{array}\right)\,.$ (43)

Here $\operatorname{K}(z)$ and $\operatorname{E}(z)$ denote the complete
elliptic integrals of the first and second kind,

$\operatorname{K}(z)=\int_{0}^{1}\frac{dx}{\sqrt{(1-x^{2})(1-zx^{2})}}\textrm{~{}~{}and~{}~{}}\operatorname{E}(z)=\int_{0}^{1}{dx}\,{\sqrt{\frac{1-zx^{2}}{1-x^{2}}}}\,.$
(44)

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 Section 7. In this region the functions $U_{i}(\tau)$ are
power series in $\epsilon$ whose coefficients can be expressed in terms of
iterated Eisenstein integrals. The variable $\tau\in\overline{\mathbb{H}}$ on
which the iterated Eisenstein integrals depend is related to the variable $z$
in eq. (39) via

$\tau=i\,\frac{\operatorname{K}(1-z)}{\operatorname{K}(z)}\,.$ (45)

Note that $\tau$ is purely imaginary for $0<z<1$. The functions
$\operatorname{K}(z)$ and $i\operatorname{K}(1-z)$ are in fact periods of the
elliptic curve defined by the polynomial equation $y^{2}=x(1-x)(1-zx)$. The
first few orders in the Laurent expansion of the functions $U_{i}(\tau)$ read

$\begin{split}U_{i}(\tau)=\sum_{k=0}^{\infty}\epsilon^{k}U_{i,k}(\tau)\,,\end{split}$
(46)

with

$\displaystyle U_{1,0}(\tau)$ $\displaystyle\,=1\,,$ $\displaystyle
U_{1,1}(\tau)$
$\displaystyle\,=8\,I(\small{\texteisfont{a}}_{2,2,1,0};\tau)+4\,I(\small{\texteisfont{a}}_{2,2,1,1};\tau)-2\pi^{2}I(1;\tau)-4\log
2\,,$ (47) $\displaystyle U_{1,2}(\tau)$
$\displaystyle\,=\frac{\pi^{2}}{6}+8\log^{2}2-16\pi^{2}\,I(1,\small{\texteisfont{a}}_{2,2,1,0};\tau)-8\pi^{2}\,I(1,\small{\texteisfont{a}}_{2,2,1,1};\tau)+180\,I(1,\small{\texteisfont{a}}_{4,2,0,0};\tau)$
$\displaystyle\,-16\pi^{2}\,I(\small{\texteisfont{a}}_{2,2,1,0},1;\tau)+64\,I(\small{\texteisfont{a}}_{2,2,1,0},\small{\texteisfont{a}}_{2,2,1,0};\tau)+32\,I(\small{\texteisfont{a}}_{2,2,1,0},\small{\texteisfont{a}}_{2,2,1,1};\tau)$
$\displaystyle\,-8\pi^{2}\,I(\small{\texteisfont{a}}_{2,2,1,1},1;\tau)+32\,I(\small{\texteisfont{a}}_{2,2,1,1},\small{\texteisfont{a}}_{2,2,1,0};\tau)+16\,I(\small{\texteisfont{a}}_{2,2,1,1},\small{\texteisfont{a}}_{2,2,1,1};\tau)$
$\displaystyle\,-32\log 2\,I(\small{\texteisfont{a}}_{2,2,1,0};\tau)-16\log
2\,I(\small{\texteisfont{a}}_{2,2,1,1};\tau)+8\pi^{2}\log 2\,I(1;\tau)\,,$

and

$\displaystyle U_{2,0}(\tau)$ $\displaystyle\,=0\,,$ $\displaystyle
U_{2,1}(\tau)$ $\displaystyle\,=i\pi\,,$ (48) $\displaystyle U_{2,2}(\tau)$
$\displaystyle\,=8i\pi I(\small{\texteisfont{a}}_{2,2,1,0};\tau)+4i\pi
I(\small{\texteisfont{a}}_{2,2,1,1};\tau)+\frac{90}{i\pi}I(\small{\texteisfont{a}}_{4,2,0,0};\tau)-4i\pi\log
2\,.$

Let us make some comments about these results. First, since $\tau$ 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. (43)). Second, the functions $U_{i}$ are
pure functions of uniform weight ArkaniHamed:2010gh  in the sense of ref.
Broedel:2018qkq . We emphasise that this statement depends on our choice for
the periods of the elliptic curve: in eq. (43) we have singled out the periods
$\operatorname{K}(z)$ and $i\operatorname{K}(1-z)$. We could of course have
made a different choice of periods, and any two choices are related by an
$\operatorname{SL}(2,\mathbb{Z})$ transformation,

$\left(\begin{array}[]{c}i\operatorname{K}(1-z)\\\
\operatorname{K}(z)\end{array}\right)\longrightarrow\gamma\left(\begin{array}[]{c}i\operatorname{K}(1-z)\\\
\operatorname{K}(z)\end{array}\right)\,,\qquad\gamma=\left(\begin{smallmatrix}a&b\\\
c&d\end{smallmatrix}\right)\in\operatorname{SL}(2,\mathbb{Z})\,,$ (49)

such that $\tau$ transform as in eq. (4). A different choice would have led to
a different form for the matrix $S(z)$ and the pure functions $U_{i}(\tau)$.
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(\gamma,z)$ the matrix corresponding to choosing
another basis of periods (cf. eq. (49)), then this matrix will be related to
$S(z)=S(\mathds{1},z)$ through the relation

$S(\gamma,z)=S(z)\,M(\gamma^{-1},\tau)\,,$ (50)

with

$M(\gamma,\tau)=\left(\begin{array}[]{cc}c\tau+d&0\\\
c&\frac{1}{c\tau+d}\end{array}\right)\,.$ (51)

Note that the matrices $M(\gamma,\tau)$ respect the group law of
$\operatorname{SL}(2,\mathbb{Z})$,

$M(\gamma_{1}\gamma_{2},\tau)=M(\gamma_{2},\tau)M(\gamma_{1},\gamma_{2}\cdot\tau)\textrm{~{}~{}and~{}~{}}M(\gamma^{-1},\tau)=M(\gamma,\gamma^{-1}\cdot\tau)^{-1}\,.$
(52)

##  4 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.

###  4.1 Multiple polylogarithms

Multiple polylogarithms (MPLs) are a generalisation of the ordinary logarithm
and the classical polylogarithm function. They are defined by the iterated
integrals GoncharovMixedTate ; Goncharov:2005sla ,

$G(a_{1},\ldots,a_{n};z)=\int_{0}^{z}\frac{dt}{t-a_{1}}\,G(a_{2},\ldots,a_{n};z)\,,$
(53)

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. (53) is
divergent, and we define instead,

$G(\underbrace{0,\ldots,0}_{n\textrm{ times}};z)=\frac{1}{n!}\log^{n}z\,.$
(54)

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.,

$G(\underbrace{0,\ldots,0}_{n-1\textrm{ times}},1;z)=-\textrm{Li}_{n}(z)\,.$
(55)

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}\in\Sigma$ (we assume $\Sigma$ 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, 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}_{\textrm{dR}}(X)$.
In our case the space $X$ is the Riemann sphere with the points
$a_{i}\in\Sigma$ removed, $X=\mathbb{CP}^{1}\setminus\Sigma$, and the
‘obstruction’ to finding a rational primitive can be cast in the form of the
statement that $H^{1}_{\textrm{dR}}(\mathbb{CP}^{1}\setminus\Sigma)$ is
generated by the logarithmic one-forms $dx/(x-a_{i})$,

$H^{1}_{\textrm{dR}}(\mathbb{CP}^{1}\setminus\Sigma)\simeq\left\langle\frac{dx}{x-a_{i}},a_{i}\in\Sigma\right\rangle_{\mathbb{Q}}\,.$
(56)

The logarithmic one-forms are precisely the integration kernels that appear in
the definition of the MPLs in eq. (53). Loosely speaking, the first de Rham
cohomology group of a Riemann surface allows us to classify those meromorphic
functions on $\mathbb{CP}^{1}\setminus\Sigma$ 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.

###  4.2 Iterated Eisenstein integrals and MPLs

In Section 2 we have defined the modular curve $Y_{\Gamma}=\mathbb{H}/\Gamma$.
It can be shown that $Y_{\Gamma}$ 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 $\Gamma$. 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 $\mathbb{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_{\Gamma}$ 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. (56).

It turns out that there is a connection between modular forms of weight two
for $\Gamma$ and the first de Rham cohomology group of $Y_{\Gamma}$. To see
this, we start by noting that

$d(\gamma\cdot\tau)=\frac{d\tau}{(c\tau+d)^{2}}\,,\qquad\gamma=\left(\begin{smallmatrix}a&b\\\
c&d\end{smallmatrix}\right)\in\operatorname{SL}(2,\mathbb{Z})\,.$ (57)

In other words, $d\tau$ transforms like a modular form of weight $-2$. Hence,
every modular form $h(\tau)$ of weight two defines a differential
$d\tau\,h(\tau)$ that is invariant under $\Gamma$, and therefore furnishes a
well-defined differential on $Y_{\Gamma}$. In addition, there is a one-to-one
map between the differentials of the first kind on $Y_{\Gamma}$ (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_{\Gamma}$ is
equal to the number of cusp forms of weight two for $\Gamma$. Similarly, there
is a one-to-one map between the differentials of the third kind of
$Y_{\Gamma}$ 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 $\Gamma$ for which $Y_{\Gamma}$ has
genus zero, i.e., $Y_{\Gamma}$ is the Riemann sphere with a finite number of
points removed. This happens precisely when $\textrm{dim }\cal
S_{2}(\Gamma)=0$. For the congruence subgroups in eqs. (2) and (3) it is well
known for which level the modular curves have genus zero:

  * •

$Y_{0}(N)=Y_{\Gamma_{0}(N)}$ has genus zero if and only if
$N\in\\{1,\ldots,10,12,13,16,18,25\\}$.

  * •

$Y_{1}(N)=Y_{\Gamma_{1}(N)}$ has genus zero if and only if
$N\in\\{1,\ldots,10,12\\}$.

  * •

$Y(N)=Y_{\Gamma(N)}$ has genus zero if and only if $N\in\\{1,\ldots,5\\}$.

In these cases there is a modular function $f:\mathbb{H}\to Y_{\Gamma}$ for
$\Gamma$, called the _Hauptmodul_ , such that every modular function for
$\Gamma$ is a rational function of $f(\tau)$. Hauptmodule are known for many
of the modular curves considered above, see, e.g., refs. Maier ; Yang . If
$h(\tau)$ is an Eisenstein series of weight two for $\Gamma$,
$\frac{d\tau}{2\pi i}\,h(\tau)$ is invariant under $\Gamma$, and $f$ maps it
to a rational differential form on $Y_{\Gamma}$. Since we know all the
generators of the first de Rham cohomology group of $Y_{\Gamma}$ from eq.
(56), we can change variables to $t=f(\tau)$ and write $d\tau\,h(\tau)$ as a
linear combination of logarithmic differential forms $dt/(t-t_{i})$, where the
$t_{i}=f(\tau_{i})$ run over the images of the cusps of $\Gamma$ under the
Hauptmodul $f$. The conclusion of this analysis can be summarised as follows:

> If $\Gamma$ is a congruence subgroup such that $Y_{\Gamma}$ has genus zero,
> and $f:\mathbb{H}\to Y_{\Gamma}$ is a Hauptmodul for $\Gamma$, then every
> iterated Eisenstein integral of the form $I(h_{1},\ldots,h_{k};\tau)$ with
> $h_{j}\in\cal E_{2}(\Gamma)$ can be written as a linear combination of
> uniform weight $k$ of MPLs evaluated at $t=f(\tau)$.

Let us illustrate this explicitly on the iterated Eisenstein integrals that
appear in eqs. (47) and (48). The relevant modular curve $Y(2)$ has genus
zero, and a Hauptmodul is the modular lambda function,

$\lambda(\tau)=16\,\frac{\eta(\tau/2)^{8}\eta(2\tau)^{16}}{\eta(\tau)^{24}}\,,$
(58)

where $\eta(\tau)$ denotes the Dedekind eta function,

$\eta(\tau)=q^{1/24}\,\prod_{n=1}^{\infty}(1-q^{n})\,,\qquad q=e^{2\pi
i\tau}\,.$ (59)

The modular lambda function is the inverse of eq. (45) for $0<z<1$, i.e.,
$z=\lambda(\tau)$. $\Gamma(2)$ has three cusps, which may be represented by
$\tau=0$, $\tau=1$ and $\tau=i\infty$. The cusps are mapped by the Hauptmodul
to the following points on the Riemann sphere,

$\lambda(0)=1\,,\qquad\lambda(1)=\infty\,,\qquad\lambda(i\infty)=0\,.$ (60)

Hence, we can identify the modular curve $Y(2)$ with the Riemann sphere with
the points $0$, $1$ and $\infty$ removed. We therefore expect that all
iterated integrals of Eisenstein series of weight two and level two (except
for $\small{\texteisfont{a}}_{2,2,0,0}$) can be expressed in terms of MPLs in
$z$ with singularities only when $z\in\\{0,1,\infty\\}$. The corresponding set
of functions is a well-known subset of MPLs known as harmonic polylogarithms
Remiddi:1999ew . 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. (47) and (48). Letting $z=\lambda(\tau)$, we can match $q$-expansions to
show that

$\begin{split}\frac{d\tau}{2\pi
i}\,\small{\texteisfont{a}}_{2,2,1,0}(\tau)&\,=\frac{1}{12}\,\frac{dz}{z-1}-\frac{1}{6}\,\frac{dz}{z}\,,\\\
\frac{d\tau}{2\pi
i}\,\small{\texteisfont{a}}_{2,2,1,1}(\tau)&\,=\frac{1}{12}\,\frac{dz}{z-1}+\frac{1}{12}\,\frac{dz}{z}\,.\end{split}$
(61)

Inserting these relations into the integrands of the iterated Eisenstein
integrals and using the fact that the lower integration boundary
$\tau=i\infty$ is mapped to $z=0$, we can easily perform all integrations in
terms of MPLs. We find

$\begin{split}I(\small{\texteisfont{a}}_{2,2,1,0};\tau)&\,=\frac{1}{12}G(1;z)-\frac{1}{6}G(0;z)+\frac{2}{3}\log
2\,,\\\
I(\small{\texteisfont{a}}_{2,2,1,1};\tau)&\,=\frac{1}{12}G(1;z)+\frac{1}{12}G(0;z)-\frac{2}{3}\log
2\,.\end{split}$ (62)

The constants proportional to $\log 2$ require some explanation, because they
cannot be recovered by inserting eq. (61) 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 $\tau=i\infty$ and $z=0$, cf. eqs. (36) and (54).
To understand this, we first note that if we expand eq. (58) into a
$q$-series, we find

$z=\lambda(\tau)=16\,q_{2}+\cal O(q_{2}^{2})\,,$ (63)

and so

$G(0;z)=\log\lambda(\tau)=4\,\log 2+\log q_{2}+\cal O(q_{2})\,.$ (64)

We see that terms proportional to (powers of) logarithms of $\log 2$ arise
from the fact that $z\sim 16\,q_{2}$ as $q_{2}\to 0$. These factors can always
be recovered from the following recipe:

  1. 1.

Insert eq. (61) into the integrands of the iterated Eisenstein integrals and
perform the integrals in the naive way.

  2. 2.

Use the shuffle algebra properties to express every $G(a_{1},\ldots,a_{n};z)$
with $a_{n}=0$ in terms of MPLs with $a_{n}\neq 0$ and powers of $G(0;z)$.

  3. 3.

Perform the replacement $G(0;z)\to G(0;z)-4\log 2$. This replacement accounts
for the difference in regularisation scheme.

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,

$\begin{split}U_{1}(\tau)&\,=1+\epsilon\,\left[G(1;z)-G(0;z)-2\pi^{2}I(1;\tau)\right]\\\
&\,+\epsilon^{2}\,\Bigg{[}180\,I(1,\small{\texteisfont{a}}_{4,2,0,0};\tau)+2\pi^{2}\left(G(0;z)-G(1;z)\right)I(1;\tau)\\\
&\,+G(0,0;z)-G(0,1;z)-G(1,0;z)+G(1,1;z)+\frac{\pi^{2}}{6}\Bigg{]}+\cal
O(\epsilon^{3})\,,\\\
U_{2}(\tau)&\,=i\pi\,\epsilon+i\pi\,\epsilon^{2}\,\left[G(1;z)-G(0;z)-\frac{90}{\pi^{2}}\,I(\small{\texteisfont{a}}_{4,2,0,0};\tau)\right]+\cal
O(\epsilon^{3})\,.\end{split}$ (65)

##  5 Modular transformations of iterated integrals

###  5.1 Motivation

At the end of the previous section we have given the first three orders in the
$\epsilon$ expansion of the functions $U_{1}$ and $U_{2}$ which describe the
familiy of integrals in eq. (3.2) in terms of MPLs and iterated Eisenstein
integrals. The expressions in eq. (65), however, are only valid in the region
where $0<z<1$. Indeed, in that region the variable $\tau$ defined in eq. (45)
is purely imaginary, because both $\operatorname{K}(z)$ and
$\operatorname{K}(1-z)$ are real. If $z\notin[0,1]$, we need to analytically
continue the elliptic integrals using the formula,

$\operatorname{K}(z+i0)=\frac{1}{\sqrt{z}}\,\left[\operatorname{K}(1/z)+i\,\operatorname{K}(1-1/z)\right]\,,\quad
z>1\,.$ (66)

For $z<0$, we obtain,

$\tau=1+\frac{i\operatorname{K}\left(\frac{1}{1-z}\right)}{\operatorname{K}\left(\frac{z}{z-1}\right)}=\gamma_{1}\cdot\tau_{1}\,,$
(67)

with

$\gamma_{1}=\left(\begin{smallmatrix}1&1\\\
0&1\end{smallmatrix}\right)\textrm{~{}~{}and~{}~{}}\tau_{1}=i\,\frac{\operatorname{K}\left(\frac{1}{1-z}\right)}{\operatorname{K}\left(\frac{z}{z-1}\right)}\,.$
(68)

Note that for $z<0$, $\tau_{1}$ is purely imaginary and
$\textrm{Im}\,\tau_{1}>0$, and the matrix $\gamma_{1}$ lies in
$\operatorname{SL}(2,\mathbb{Z})$. Similarly, for $z>1$, we find,

$\tau=\frac{i\operatorname{K}\left(1-\frac{1}{z}\right)}{\operatorname{K}\left(\frac{1}{z}\right)+i\,\operatorname{K}\left(1-\frac{1}{z}\right)}=\gamma_{3}\cdot\tau_{3}\,,$
(69)

with

$\gamma_{3}=\left(\begin{smallmatrix}0&1\\\
-1&1\end{smallmatrix}\right)\textrm{~{}~{}and~{}~{}}\tau_{3}=i\frac{\operatorname{K}\left(\frac{1}{z}\right)}{\operatorname{K}\left(1-\frac{1}{z}\right)}\,.$
(70)

Again, we find that $\tau_{3}$ is purely imaginary with positive imaginary
part for $z>1$, and that the matrix $\gamma_{3}$ lies in
$\operatorname{SL}(2,\mathbb{Z})$. Hence, if we define
$\gamma_{2}=\left(\begin{smallmatrix}1&0\\\ 0&1\end{smallmatrix}\right)$ and
$\tau_{2}=i\operatorname{K}(1-z)/\operatorname{K}(z)$ for $0<z<1$, we see that
for each of the three regions there is
$\gamma_{i}\in\operatorname{SL}(2,\mathbb{Z})$ and $\tau_{i}\in\mathbb{H}$
such that $\tau=\gamma_{i}\cdot\tau_{i}$. We can also define a map that
assigns to all real values of $z$ the corresponding value of $\tau$,

$\tau(z)=\left\\{\begin{array}[]{ll}\gamma_{1}\cdot\tau_{1}\,,&\textrm{ if
}z<0\,,\\\ \gamma_{2}\cdot\tau_{2}\,,&\textrm{ if }0\leq z\leq 1\,,\\\
\gamma_{3}\cdot\tau_{3}\,,&\textrm{ if }z>1\,.\end{array}\right.$ (71)

In the remainder of this section we describe how we can express iterated
Eisenstein integrals evaluated at $\tau=\gamma_{i}\cdot\tau_{i}$ in terms of
iterated Eisenstein integrals in $\tau_{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. (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 $\tau_{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.

###  5.2 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},\ldots,h_{k};\gamma\cdot\tau)$, with
$\gamma\in\operatorname{SL}(2,\mathbb{Z})$, in terms of iterated Eisenstein
integrals evaluated at $\tau$. The general strategy is simple: we start from
the regularised version of the integral in eq. (36), perform the change of
variables $\tau^{\prime}\to\gamma\cdot\tau^{\prime}$ in the integrand, and
integrate back. Since $\operatorname{SL}(2,\mathbb{Z})$ is generated by the
three elements $S$, $T$ and $-\mathds{1}$ in eq. (1), it is sufficient to
discuss the cases $\gamma=S$ and $\gamma=T$, because a generic
$\gamma\in\operatorname{SL}(2,\mathbb{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. (36): treating the divergent integrals
naively can lead to wrong results! Indeed, consider the integral

$I(\small{\texteisfont{a}}_{4,1,0,0};\tau)=-\frac{\pi^{2}}{180}\,\log
q_{1}-\pi^{2}\,\left[\frac{4}{3}\,q_{1}+6q_{1}^{2}+\cal
O(q_{1}^{3})\right]\,.$ (72)

Next, write $\tau=1+\tilde{\tau}=\left(\begin{smallmatrix}1&1\\\
0&1\end{smallmatrix}\right)\cdot\tilde{\tau}$. Using the definition of the
iterated Eisenstein integrals, we can write,

$I(\small{\texteisfont{a}}_{4,1,0,0};1+\tilde{\tau})=\int_{i\infty}^{1+\tilde{\tau}}\frac{d\tau^{\prime}}{2\pi
i}\,\small{\texteisfont{a}}_{4,1,0,0}(\tau^{\prime})\,.$ (73)

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 do naive changes of variables on
the integral. We can change variables according to $\tau^{\prime}\to
1+\tau^{\prime}$. The integrand is invariant under this change of variables,
because
$\small{\texteisfont{a}}_{4,1,0,0}(1+\tau^{\prime})=\small{\texteisfont{a}}_{4,1,0,0}(\tau^{\prime})$.
Naively, the lower integration boundary $i\infty$ does not change under
translations, and we are led to believe that
$I(\small{\texteisfont{a}}_{4,1,0,0};1+\tilde{\tau})$ is identical to
$I(\small{\texteisfont{a}}_{4,1,0,0};\tilde{\tau})$. It is easy to see from
eq. (72) that this cannot be correct: since $e^{2\pi
i(\tilde{\tau}+1)}=e^{2\pi i\tilde{\tau}}=q_{1}$, we see that the power series
part in eq. (72) remains unchanged. The first term, however, does change, and
we find,

$I(\small{\texteisfont{a}}_{4,1,0,0};1+\tilde{\tau})=I(\small{\texteisfont{a}}_{4,1,0,0};\tilde{\tau})-\frac{i\pi^{3}}{90}\,.$
(74)

We see that the constant in eq. (74) is a consequence of the regularisation.
Indeed, if we start from the regularised version in eq. (36), we find

$I(\small{\texteisfont{a}}_{4,1,0,0};\tau)=-\frac{i\pi^{3}}{90}\,\tau+\int_{i\infty}^{\tau}\frac{d\tau^{\prime}}{2\pi
i}\,\left[\small{\texteisfont{a}}_{4,1,0,0}(\tau^{\prime})-\frac{\pi^{4}}{45}\right]\,.$
(75)

The second term is an absolutely convergent integral and gives rise to the
power series part in eq. (72). 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. (74).

We can apply exactly the same steps to more general modular transformations:
we always start from the regularised version in eq. (36) 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 $\omega_{j}$ correspond to a modular form of weight zero (except for
explicit powers of $\tau$), but instead we consider differential forms of the
form $d\tau\,h_{j}(\tau)\,\tau^{m_{j}}$, for some positive integer $m_{j}$ and
$h_{j}$ a modular form of weigh $n_{j}$. An example will clarify this:

$\begin{split}I(1,h,1;\tau)&\,=I(1;\tau)\,I(h,1;\tau)-2I(h,1,1;\tau)\\\
&\,=\frac{\tau}{2\pi i}\int_{i\infty}^{\tau}\frac{d\tau^{\prime}}{2\pi
i}\,h(\tau^{\prime})\,I(1;\tau^{\prime})-2\int_{i\infty}^{\tau}\frac{d\tau^{\prime}}{2\pi
i}\,h(\tau^{\prime})\,I(1,1;\tau^{\prime})\\\ &\,=\frac{\tau}{2\pi
i}\int_{i\infty}^{\tau}\frac{d\tau^{\prime}}{2\pi
i}\,h(\tau^{\prime})\,\frac{\tau^{\prime}}{2\pi
i}-\int_{i\infty}^{\tau}\frac{d\tau^{\prime}}{2\pi
i}h(\tau^{\prime})\,\left(\frac{\tau^{\prime}}{2\pi
i}\right)^{2}\,.\end{split}$ (76)

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\tau\,h_{j}(\tau)\,\tau^{m_{j}}$ with $0\leq
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.
ManinModular . 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 $\gamma\in\operatorname{SL}(2,\mathbb{Z})$, though we keep in mind
that we can recover every modular transformation via a suitable sequence of
$S$ and $T$ transformations. If $\gamma=\left(\begin{smallmatrix}a&b\\\
c&d\end{smallmatrix}\right)\in\operatorname{SL}(2,\mathbb{Z})$, we find,

$\begin{split}I(h_{1},\ldots,h_{k};\gamma\cdot\tau)&\,=\sum_{p=0}^{n}\frac{a_{0}(h_{p+1})\ldots
a_{0}(h_{k})}{(2\pi
i)^{k-p}(k-p)!}\,(\gamma\cdot\tau)^{k-p}\int_{i\infty}^{\gamma\cdot\tau}R[\omega_{p}\ldots\omega_{1}]\\\
&\,=\sum_{p=0}^{n}\frac{a_{0}(h_{p+1})\ldots a_{0}(h_{k})}{(2\pi
i)^{k-p}(k-p)!}\,\left(\frac{a\tau+b}{c\tau+d}\right)^{k-p}\int_{-d/c}^{\tau}R[\omega_{p}^{\gamma}\ldots\omega_{1}^{\gamma}]\,,\end{split}$
(77)

where $\omega_{j}^{\gamma}$ denotes the transform of $\omega_{j}$ under
$\gamma$,

$\omega_{j}=\frac{d\tau}{2\pi
i}\,h_{j}(\tau)\xrightarrow{\,\,\,\gamma\,\,\,}\omega_{j}^{\gamma}=\frac{d\tau}{2\pi
i}\,{(c\tau+d)^{-2}}\,h_{j}(\gamma\cdot\tau)\,.$ (78)

In the previous section we have seen that if $h_{j}(\tau)$ is a modular form
of weight $n_{j}$ for $\Gamma(N)$, then so is $h_{j}(\gamma\cdot\tau)$. In
particular, if $h_{j}=h_{N,r,s}^{(n_{j})}$, eq. (12) implies,

$\omega_{j}^{\gamma}=\frac{d\tau}{2\pi
i}\,(c\tau+d)^{n_{j}-2}\,h^{(n_{j})}_{N,rd+sb,rc+sa}(\tau)\,.$ (79)

The transformation behaviour of the Eisenstein series
$\small{\texteisfont{a}}_{n,N,r,s}$ and $\small{\texteisfont{b}}_{n,N,r,s}$
can easily be recovered from their definition in terms of $h_{N,r,s}^{(n)}$
and $h_{N,r,-s}^{(n)}$, cf. eq. (24). The additional powers of $\tau$ that are
generated in the right-hand side can be expressed in terms of iterated
integrals of modular forms of weight 0,

$\tau^{p}=(2\pi i)^{p}\,p!\,I(\underbrace{1,\ldots,1}_{p};\tau)\,.$ (80)

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\tau\,h_{j}(\tau)\,\tau^{m_{j}}$, then these differential forms transform as

$d\tau\,h_{j}(\tau)\,\tau^{m_{j}}\xrightarrow{\,\,\,\gamma\,\,\,}d\tau\,h_{j}(\tau)\,(a\tau+b)^{m_{j}}\,(c\tau+d)^{n_{j}-m_{j}-2}\,.$
(81)

We see that no negative powers appear precisely when the condition $0\leq
m_{j}\leq n_{j}-2$ is fulfilled. This condition also automatically enforces
$n_{j}\geq 2$, which are precisely the conditions that define our restricted
class of integrals. We note that this condition is preserved under modular
transformations ManinModular  (which can be seen, for example, from the fact
that the highest power of $\tau$ in eq. (81) never exceeds $n_{j}-2$ for
$0\leq m_{j}\leq 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. (77) that the lower integration boundaries are
constant rational values $-\frac{d}{c}=\gamma^{-1}\cdot i\infty$. 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 $\gamma$ into a product
of $S$’s and $T$’s, we only obtain constant iterated Eisenstein integrals
evaluated at $\tau=0$, because

$T\cdot i\infty=i\infty{\textrm{~{}~{}and~{}~{}}}S\cdot i\infty=0\,.$ (82)

Iterated Eisenstein integrals evaluated at $\tau=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\leq j\leq k$, of weight $n_{j}\geq 2$ and
integers $m_{j}$ such that $1\leq m_{j}<n_{j}$, the associated _multiple
modular value_ (MMV) is defined by the integral ManinModular ; Brown:mmv :

$\begin{split}\Lambda&(h_{k},\ldots,h_{1};m_{k},\ldots,m_{1})=\\\
&\,=(-i)^{m}\,\int_{0}^{i\infty}d\tau_{1}\,h_{1}(\tau_{1})\,\tau_{1}^{m_{1}-1}\int_{0}^{\tau_{1}}\ldots\int_{0}^{\tau_{k-1}}d\tau_{k}\,h_{k}(\tau_{k})\,\tau_{k}^{m_{k}-1}\,,\end{split}$
(83)

with $m=m_{1}+\ldots+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. (83) may
diverge, and we have to interpret it again as a regularised version. We define
the weight of an MMV as $n_{1}+\ldots+n_{k}$ and the length as $m$. MMVs of
small weight and lengths associated to modular forms of level one have been
studied in ref. Brown:mmv ; Brown:mmv2 . It is easy to see that MMVs form a
shuffle algebra,

$\begin{split}\Lambda&(h_{l},\ldots,h_{k+1};m_{l},\ldots,m_{k+1})\,\Lambda(h_{k},\ldots,h_{1};m_{k},\ldots,m_{1})\\\
&\,=\sum_{\sigma\in\Sigma(k,l)}\Lambda(h_{\sigma_{l}},\ldots,h_{\sigma_{1}};m_{\sigma_{k}},\ldots,m_{\sigma_{1}})\,,\end{split}$
(84)

where the sum runs over all shuffles of $k$ and $l-k$ elements. Let us
conclude by discussing an identity relating different MMVs. The
$S$-transformation maps the integration path $[i\infty,0]$ to itself (up to
reversal of the orientation), and so this transformation induces a
transformation on MMVs. We also note that there is a reflection identity for
MMVs associated to Eisenstein series for $\Gamma(N)$,

$\Lambda(h_{N,r_{k},s_{k}}^{(n_{k})},\ldots,h_{N,r_{1},s_{1}}^{(n_{1})};m_{k},\ldots,m_{1})=i^{n}\,\Lambda(h_{N,\tilde{r}_{1},\tilde{s}_{1}}^{(n_{1})},\ldots,h_{N,\tilde{r}_{k},\tilde{s}_{k}}^{(n_{k})};\tilde{m}_{1},\ldots,\tilde{m}_{k})\,,$
(85)

where $n$ is the weight of the MMV, and

$\tilde{m}_{j}=n_{j}-m_{j}\,,\qquad\tilde{r}_{j}=N-s_{j}\,,\qquad\tilde{s}_{j}=r_{j}\,.$
(86)

The reflection identity can easily be obtained by working out the
transformation of the corresponding iterated integrals under the modular
transformation $S:\tau\to-1/\tau$.

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 $\tau=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.

##  6 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 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.

###  6.1 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 $\tau$,
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 $\tau=i\infty$
(or $q_{N}=0$). For Eisenstein series we can write down the $q$-expansion
around the infinite cusp in closed form, cf. eq. (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\in\mathbb{Q}$ and an Eisenstein series $f\in\cal
E_{n}(\Gamma(N))$. Since the $h^{(n)}_{N,r,s}(\tau)$ form a spanning set for
$\cal E_{n}(\Gamma(N))$, it is sufficient to consider the case where $f$ is an
element of this set.777We could also work with the spanning set formed by the
$\small{\texteisfont{a}}_{n,N,r,s}$ and $\small{\texteisfont{b}}_{n,N,r,s}$.
The conversion between the two sets is trivial, cf. eq. (24). We have seen in
Section 2 that there is $\gamma_{c}\in\operatorname{SL}(2,\mathbb{Z})$ such
that $c=\gamma_{c}\cdot i\infty$. Hence, expanding $h^{(n)}_{N,r,s}(\tau)$
around the cusp $\tau=c$ is equivalent to expanding
$h^{(n)}_{N,r,s}(\gamma_{c}\cdot\tau^{\prime})$ around $\tau^{\prime}=i\infty$
using eqs. (12) and (13). Let us illustrate this on an example: Assume we want
to expand $f(\tau)=h^{(4)}_{2,1,1}(\tau)$ into a series close to $\tau=1/2$.
We have $\gamma_{1/2}=\left(\begin{smallmatrix}1&0\\\
2&1\end{smallmatrix}\right)$, and eq. (12) gives

$h^{(4)}_{2,1,1}(\tau)=h^{(4)}_{2,1,1}(\gamma_{1/2}\cdot\tau^{\prime})=(2\tau^{\prime}+1)^{4}\,h^{(4)}_{2,1,1}(\tau^{\prime})\,.$
(87)

Using eq. (13) we can obtain the desired expansion,

$\displaystyle h^{(4)}_{2,1,1}(\tau)$
$\displaystyle\,=(2\tau^{\prime}+1)^{4}\,\frac{\pi^{4}}{360}\left[7+240\,e^{i\pi\tau^{\prime}}-240\,e^{2i\pi\tau^{\prime}}+6720\,e^{3i\pi\tau^{\prime}}+\ldots\right]$
(88)
$\displaystyle\,=(2\tau-1)^{-4}\,\frac{\pi^{4}}{360}\left[7+240\,e^{i\pi\tau/(1-2\tau)}-240\,e^{2i\pi\tau/(1-2\tau)}+6720\,e^{3i\pi\tau/(1-2\tau)}+\ldots\right].$

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\in\mathbb{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.

###  6.2 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,r,s}^{(n)}(\tau)$. This function admits a
$q$-expansion, and we can easily write down an arbitrary number of terms in
the expansion, cf. eq. (13). Whenever $\textrm{Im}\,\tau>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 $\tau\sim
i\infty$ (and $\textrm{Re}\,\tau$ is not too large), we have $|q_{N}|\ll 1$
and the series converges extremely fast. However, when $\textrm{Im}\,\tau$ is
not large, the convergence of the series can be extremely slow. We illustrate
this in the second column of Table 1, where we show the numerical value
obtained for $h_{2,1,1}^{(4)}(\tau)$ at $\tau=\frac{1}{2}+\frac{i}{10}$ after
truncating the $q$-expansion around the infinite cusp after a certain number
of terms. Te expansion parameter is $e^{i\pi\tau}\simeq i\,0.7304\ldots$.

$n_{\textrm{max}}$ | direct | transformed
---|---|---
1 | $36.538+47.432i$ | $1183.791-15.756i$
2 | $36.538+47.432i$ | $1183.797-15.756i$
3 | $36.538-661.090i$ | $1183.797-15.756i$
4 | $18.056-661.090i$ | $1183.797-15.756i$
5 | $18.056+1039.856i$ | $1183.797-15.756i$
10 | $642.471+1470.891i$ | $1183.797-15.756i$
30 | $1172.467+51.758i$ | $1183.797-15.756i$
60 | $1183.787-15.799i$ | $1183.797-15.756i$
70 | $1183.797-15.753i$ | $1183.797-15.756i$
75 | $1183.797-15.757i$ | $1183.797-15.756i$
Table 1: Numerical evaluation of $h_{2,1,1}^{(4)}(\tau)$ at
$\tau=\frac{1}{2}+\frac{i}{10}$ using a truncated $q$-expansion with
$n_{\textrm{max}}$ terms. The column labeled ‘direct’ is obtained by expanding
$h_{2,1,1}^{(4)}(\tau)$ around $\tau=i\infty$, while the numbers in the column
labeled ‘transformed’ are obtained by mapping $\tau$ to the fundamental domain
before expansion.

The convergence of the series can be considerably accelerated by using modular
transformations. For every $\tau\in\mathbb{H}$, there is
$\gamma\in\operatorname{SL}(2,\mathbb{Z})$ such that $\gamma\cdot\tau$ lies in
the fundamental domain for $\operatorname{SL}(2,\mathbb{Z})$ (cf. eq. (5)). In
Appendix A we review how this matrix can be constructed. In the case of our
example, we find

$\frac{1}{2}+\frac{i}{10}=\gamma_{1/2}\cdot\left(-\frac{1}{2}+\frac{5i}{2}\right)\,,$
(89)

where the matrix $\gamma_{1/2}=\left(\begin{smallmatrix}1&0\\\
2&1\end{smallmatrix}\right)$ of Section 5. Using the transformation in eq.
(87) and the expansion in eq. (88), we obtain a series representation with
expansion parameter $e^{i\pi\tau/(1-2\tau)}\simeq-i\,0.000388\ldots$. 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
$\mathbb{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
$\tau=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 $\tau=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\infty]$ into the segments $[0,i]$ and $[i,i\infty]$. The path $[0,i]$
can be mapped to $[i\infty,i]$ using $S$, because $S\cdot 0=i\infty$ and
$S\cdot i=i$. Hence, using the path composition and reversal formulæ, 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
$\tau=i$. The expansion parameter for the $q$-expansion is then
$q_{N}=e^{-2\pi/N}<1$, and we obtain a fast-converging series representation
of the MMV. After these steps all paths extend from $i\infty$ 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.
Brown:mmv2 . 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 $\Gamma(2)$, as well as all
MMVs associated to Eisenstein series for $\Gamma(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.

###  7.1 The imaginary part of the elliptic ${}_{2}F_{1}$ function

So far we have only considered the integrals in eq. (39) 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. (69) that for
$z>1$ we can write $\tau=\gamma_{3}\cdot\tau_{3}$, with
$\gamma_{3}\in\operatorname{SL}(2,\mathbb{Z})$ and $\tau_{3}$ purely
imaginary. We can then write all integrals in terms of iterated Eisenstein
integrals evaluated at $\tau_{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}(\tau)$ defined in eq.
(65). The integrals $T_{1}$ and $T_{2}$ can be recovered by multiplying with
the matrix $S(z)$ in eq. (43), which does not involve any iterated Eisenstein
integrals. For the functions $U_{i}(\tau)$, we find

$\left(\begin{array}[]{c}U_{1}(\tau)\\\
U_{2}(\tau)\end{array}\right)=\left(\begin{array}[]{c}U_{1}(\gamma_{3}\cdot\tau_{3})\\\
U_{2}(\gamma_{3}\cdot\tau_{3})\end{array}\right)=M(\gamma_{3},\tau_{3})^{-1}\,\left(\begin{array}[]{c}\tilde{U}_{1}(\tau_{3})+i\pi\,\tilde{V}_{1}(\tau_{3})\\\
\tilde{U}_{2}(\tau_{3})+i\pi\,\tilde{V}_{2}(\tau_{3})\end{array}\right)\,,$
(90)

where the matrix $M(\gamma_{3}^{-1},\tau_{3})$ was defined in eq. (51). The
real parts are given by

$\displaystyle\tilde{U}_{1}(\tau_{3})$
$\displaystyle\,=1+\epsilon\,\Bigg{[}G(1;{1}/{z})-2\pi^{2}\,I(1;\tau_{3})\Bigg{]}+\epsilon^{2}\,\Bigg{[}180I(1,\small{\texteisfont{a}}_{4,2,0,0};\tau_{3})$
$\displaystyle\,-2\pi^{2}G(1;{1}/{z})\,I(1;\tau_{3})+G(1,1;{1}/{z})-\frac{\pi^{2}}{3}\Bigg{]}+\cal
O(\epsilon^{3})\,,$ (91) $\displaystyle\tilde{U}_{2}(\tau_{3})$
$\displaystyle\,=-1-\epsilon\,G(1;{1}/{z})+\epsilon^{2}\,\Bigg{[}\frac{5\pi^{2}}{3}-180I(\small{\texteisfont{a}}_{4,2,0,0},1;\tau_{3})-G\left(1,1;{1}/{z}\right)\Bigg{]}+\cal
O(\epsilon^{3})\,.$

The imaginary parts are

$\displaystyle\tilde{V}_{1}(\tau_{3})$
$\displaystyle\,=-2I(1;\tau_{3})-2\epsilon\,G(1;{1}/{z})\,I(1;\tau_{3})-2\epsilon^{2}\Bigg{[}180\,I(1,\small{\texteisfont{a}}_{4,2,0,0},1;\tau_{3})$
$\displaystyle\,+G(1,1;{1}/{z})\,I(1;\tau_{3})+\frac{15}{\pi^{2}}\zeta_{3}\Bigg{]}+\cal
O(\epsilon^{3})\,,$ (92) $\displaystyle\tilde{V}_{2}(\tau_{3})$
$\displaystyle\,=\epsilon+\epsilon^{2}\Bigg{[}G(1;{1}/{z})-\frac{90}{\pi^{2}}\,I(\small{\texteisfont{a}}_{4,2,0,0};\tau_{3})\Bigg{]}+\cal
O(\epsilon^{3})\,.$

Since $\tau_{3}$ is purely imaginary for $z>1$, the functions $\tilde{U}_{i}$
and $\tilde{V}_{i}$ are manifestly real in this region. Therefore this is the
correct separation of $U_{i}(\tau)$ 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(\gamma_{3},\tau_{3})$
in eq. (90) contains $\tau_{3}$ in the denominator, which spoils purity.
However, all the terms that spoil purity are contained in the matrix
$M(\gamma_{3},\tau_{3})$. This can be explained as follows: The separation
into the matrix $S(z)$ and the pure functions $U_{i}(\tau)$ in eq. (93) 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
$\tau$ in eq. (45). Of course, we could have chosen a different basis of
periods, and any two choices are related by an
$\operatorname{SL}(2,\mathbb{Z})$ transformation. The matrices $S(z)$ obtained
from two different choices are related by the matrix $M(\gamma_{3},\tau_{3})$,
cf. eq. (50). Since nothing should depend on this choice, the matrix
$M(\gamma_{3},\tau_{3})$ in eq. (50) must be cancelled by a similar
contribution coming from the modular transformation of the iterated Eisenstein
integrals, so that

$\left(\begin{array}[]{c}T_{1}(z)\\\
T_{2}(z)\end{array}\right)=S(z)\left(\begin{array}[]{c}U_{1}(\tau)\\\
U_{2}(\tau)\end{array}\right)=S(\gamma_{3},z)\left(\begin{array}[]{c}\tilde{U}_{1}(\tau_{3})+i\pi\,\tilde{V}_{1}(\tau_{3})\\\
\tilde{U}_{2}(\tau_{3})+i\pi\,\tilde{V}_{2}(\tau_{3})\end{array}\right)\,.$
(93)

This explains the appearance of the matrix $M(\gamma_{3},\tau_{3})$ in eq.
(90).

###  7.2 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. (90) 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 $\tau$ has a small imaginary part, the
$q$-expansions converge very slowly.

Alternatively, for every (real) value of $z$, we can find
$\gamma_{z}\in\operatorname{SL}(2,\mathbb{Z})$ such that
$\gamma_{z}^{-1}\cdot\tau(z)$ lies in the fundamental domain, where $\tau(z)$
was defined in eq. (71). A priori, the matrix $\gamma_{z}$ depends on $z$. We
need only six distinct values of $\gamma_{z}$ to cover all real values of $z$:

$\gamma_{z}=\left\\{\begin{array}[]{ll}\gamma_{1a}\,,&\textrm{ if
}z\leq-1\,,\\\ \gamma_{1b}\,,&\textrm{ if }-1<z<0\,,\\\
\gamma_{2a}\,,&\textrm{ if }0\leq z<1/2\,,\\\ \gamma_{2b}\,,&\textrm{ if
}1/2\leq z\leq 1\,,\\\ \gamma_{3a}\,,&\textrm{ if }1<z\leq 2\,,\\\
\gamma_{3b}\,,&\textrm{ if }z>2\,,\end{array}\right.$ (94)

where we defined

$\begin{split}\gamma_{1a}&\,=\left(\begin{smallmatrix}-1&1\\\
-1&0\end{smallmatrix}\right)\,,\qquad\,\,\,\phantom{\gamma_{1}}\gamma_{2a}=\gamma_{2}=\left(\begin{smallmatrix}1&0\\\
0&1\end{smallmatrix}\right)\,,\qquad\,\,\,\,\gamma_{3a}=\left(\begin{smallmatrix}0&1\\\
-1&1\end{smallmatrix}\right)\,,\\\
\gamma_{1b}&\,=\gamma_{1}=\left(\begin{smallmatrix}1&1\\\
1&0\end{smallmatrix}\right)\,,\qquad\gamma_{2b}=\left(\begin{smallmatrix}0&-1\\\
1&0\end{smallmatrix}\right)\,,\qquad\gamma_{3b}=\gamma_{3}=\left(\begin{smallmatrix}-1&0\\\
-1&-1\end{smallmatrix}\right)\,.\end{split}$ (95)

It is easy to check that the variables
$\tau_{j\alpha}\equiv\gamma_{j\alpha}\cdot\tau(z)$, $j\in\\{1,2,3\\}$ and
$\alpha\in\\{a,b\\}$, are purely imaginary and $\textrm{Im }\tau_{j\alpha}\geq
i$. In each of the six regions in eq. (94) we obtain a representation of the
master integrals $T_{i}$ in eq. (41) in terms of iterated Eisenstein integrals
evaluated at $\tau_{j\alpha}$, 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 Fig. 1 a plot of
the real and imaginary parts of the first three coefficients in the Laurent
expansion of $T_{1}$.

Figure 1: The real and imaginary parts of the first three orders in the
$\epsilon$ expansion of the master integral $T_{1}$.

Let us conclude this section by commenting on the six regions for $z$ that
appear in eq. (94). 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 $\tau$ contain redundant information, but the
$j$-invariant is a complex number that uniquely characterises an elliptic
curve. The $j$-invariant can be computed from $\tau$. It is invariant under
$\operatorname{SL}(2,\mathbb{Z})$ (i.e., it is a modular function for
$\operatorname{SL}(2,\mathbb{Z})$), and it can be expressed as a rational
function of Eisenstein series of level one:

$j(\tau)=\frac{20\,h_{1,0,0}^{(4)}(\tau)^{3}}{20\,h_{1,0,0}^{(4)}(\tau)^{3}+49\,h_{1,0,0}^{(6)}(\tau)^{2}}\,.$
(96)

Since $j(\tau)$ is a modular function for $\operatorname{SL}(2,\mathbb{Z})$,
it must be a rational function of the Hauptmodul $z=\lambda(\tau)$. We have

$j(\tau)=\frac{4\left(1-z+z^{2}\right)^{3}}{27\,(1-z)^{2}\,z^{2}}\,.$ (97)

We now see that the poles of the $j$-invariant are $z\in\\{0,1\\}$, while
$z\in\\{-1,1/2,2\\}$ corresponds to $j(\tau)=1$! The corresponding values of
$\tau$ are

$\lambda(1+i)=-1\,,\qquad\lambda(i)=1/2\,,\qquad\lambda\left(\frac{1}{2}+\frac{i}{2}\right)=2\,.$
(98)

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.

###  7.3 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
$\epsilon$. In Section 6.1 we have shown how we can expand Eisenstein series
and their iterated integrals into a $q$-series close to a cusp. In many
applications, however, one would like to have expansions in terms of the
original variables of the problem (e.g., kinematic variables in the case of
Feynman integrals). In this section we illustrate how we can obtain expansions
in $z$ close to the singular points $z=0$ or $z=1$.

In the following we only discuss in detail the expansion of $T_{1}(z)$ and
$T_{2}(z)$ close to $z=0$. The case $z=1$ can be dealt with in the the same
way. We have $\tau(0)=i\infty$, and so expanding for small values of $z$ is
equivalent to expanding around $\tau(z)\sim i\infty$. We can easily expand all
iterated Eisenstein integrals in eq. (65) into a $q$-series around
$\tau=i\infty$ and then convert the $q$-series into a series in $z$ by
expanding eq. (45). We find

$\tau(z)=\frac{i}{\pi}\,(4\log 2-\log
z)-\frac{i}{2\pi}\,z-\frac{13i}{64\pi}\,z^{2}-\frac{23i}{192\pi}\,z^{3}-\frac{2701i}{32768\pi}\,z^{4}+\cal
O(z^{5})\,.$ (99)

or equivalently,

$q_{2}=e^{i\pi\tau(z)}=\frac{z}{16}+\frac{z^{2}}{32}+\frac{21\,z^{3}}{1024}+\frac{31\,z^{4}}{2048}+\cal
O(z^{5})\,.$ (100)

Inserting eq. (100) into the $q$-series, we find

$\begin{split}T_{1}(z)&\,=\pi\,\left[1+\frac{1}{4}\,z+\frac{9}{64}\,z^{2}+\cal
O(z^{3})\right]\\\ &\,-\epsilon\,\pi\,\Bigg{[}4\log 2+z\left(\frac{1}{2}+\log
2\right)+z^{2}\left(\frac{27}{64}+\frac{9}{16}\log 2\right)+\cal
O(z^{3})\Bigg{]}\\\
&\,+\epsilon^{2}\,\pi\,\Bigg{[}\frac{\pi^{2}}{6}+8\log^{2}2+z\left(\frac{\pi^{2}}{24}+2\log^{2}2+2\log
2\right)\\\
&\,\qquad\qquad+z^{2}\left(\frac{23}{64}+\frac{3\pi^{2}}{128}+\frac{9}{8}\log^{2}2+\frac{27}{16}\log
2\right)+\cal O(z^{3})\Bigg{]}+\cal O(\epsilon^{3})\,,\\\
T_{2}(z)&\,=\pi\,\Bigg{[}\frac{1}{3z}-\frac{1}{12}-\frac{11}{192}\,z-\frac{29}{768}\,z^{2}+\cal
O(z^{3})\Bigg{]}\\\ &\,+\epsilon\,\pi\,\Bigg{[}-\frac{4}{3z}\,\log
2-\frac{1}{6}+\frac{1}{3}\log 2+z\left(\frac{25}{192}+\frac{11}{48}\log
2\right)\\\ &\,\qquad\qquad+z^{2}\left(\frac{283}{2304}+\frac{29}{192}\log
2\right)+\cal O(z^{3})\Bigg{]}\\\
&\,+\epsilon^{2}\,\pi\,\Bigg{[}\frac{1}{z}\,\left(\frac{\pi^{2}}{18}+\frac{8}{3}\log^{2}2\right)-\frac{\pi^{2}}{72}-\frac{2}{3}\log^{2}2+\frac{2}{3}\log
2\\\
&\,\qquad\qquad+z\left(-\frac{13}{192}-\frac{11\pi^{2}}{1152}-\frac{11}{24}\log^{2}2-\frac{25}{48}\log
2\right)\\\
&\,\qquad\qquad+z^{2}\left(-\frac{35}{256}-\frac{29\pi^{2}}{4608}-\frac{29}{96}\log^{2}2-\frac{283}{576}\log
2\right)\Bigg{]}+\cal O(\epsilon^{3})\,.\end{split}$ (101)

This expansion is the same as the one obtained from the usual series
representation of Gauss’ hypergeometric function,

${}_{2}F_{1}(a,b;c;z)=\sum_{n=0}^{\infty}\frac{(a)_{n}(b)_{n}}{(c)_{n}}\,\frac{z^{n}}{n!}\,.$
(102)

##  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 Tarasov:1996br ; Lee:2009dh .

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

$\begin{split}S_{a_{1},\dots,a_{5}}&(p^{2},m^{2};d)\\\
&\,=\int\mathfrak{D}^{d}\ell_{1}\mathfrak{D}^{d}\ell_{2}\frac{(\ell_{1}\cdot
p)^{a_{4}}(\ell_{2}\cdot
p)^{a_{5}}}{[\ell_{1}^{2}-m^{2}]^{a_{1}}[\ell_{2}^{2}-m^{2}]^{a_{2}}[(\ell_{1}-\ell_{2}-p)^{2}-m^{2}]^{a_{3}}}\
,\end{split}$ (103)

where the integration measure is chosen as

$\int\mathfrak{D}^{d}\ell=\frac{1}{\Gamma\left(2-\frac{d}{2}\right)}\int\frac{d^{d}\ell}{i\pi^{d/2}}\,.$
(104)

Using IBP identities Chetyrkin:1981qh ; Tkachov:1981wb , one can show that the
family of integrals defined in eq. (103) 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\epsilon)=1\,.$ (105)

Following ref. Remiddi:2016gno ; Broedel:2019kmn , we choose the other two
master integrals as follows

$\begin{split}\cal
S_{1}(\epsilon;t)&=-m^{2}\,S_{1,1,1,0,0}(p^{2},m^{2};2-2\epsilon)\ ,\\\ \cal
S_{2}(\epsilon;t)&=-\left[\frac{1}{3}(t^{2}-6t+21)-12\epsilon(t-1)\right]m^{2}\,S_{1,1,1,0,0}(p^{2},m^{2};2-2\epsilon)\\\
&\quad-2(t-1)(t-9)m^{4}\,S_{2,1,1,0,0}(p^{2},m^{2};2-2\epsilon),\end{split}$
(106)

where we defined the dimensionless variable $t=p^{2}/m^{2}$ and rescaled by
the appropriate powers of $m$ in order to render $\cal S_{1}(\epsilon;t)$ and
$\cal S_{2}(\epsilon;t)$ dimensionless as $\epsilon\to 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.

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
Remiddi:2016gno

$\begin{split}\partial_{t}\begin{pmatrix}\cal S_{1}(t)\\\ \cal
S_{2}(t)\end{pmatrix}\,=\,B(t)\begin{pmatrix}\cal S_{1}(t)\\\ \cal
S_{2}(t)\end{pmatrix}+\begin{pmatrix}0\\\ 1\end{pmatrix}\ ,\end{split}$ (107)

with

$\begin{split}B(t)\,&=\,\frac{1}{6\,t\,(t-1)(t-9)}\begin{pmatrix}3(3+14t-t^{2})&-9\\\
(t+3)(3+75t-15t^{2}+t^{3})&-3(3+14t-t^{2})\end{pmatrix}\,,\end{split}$ (108)

and where we defined $\cal S_{j}(t)=\cal S_{j}(0;t)$ with $j=1,2$. Equation
(107) is a coupled $2\times 2$ system of equations whose solution requires one
to first find a full solution of the corresponding homogeneous system, i.e., a
$2\times 2$ matrix $\mathcal{W}_{S}(t)$ such that

$\partial_{t}\mathcal{W}_{S}(t)=B(t)\mathcal{W}_{S}(t)\,.$ (109)

As shown for the first time in ref. Laporta:2004rb , and generalised in ref.
Primo:2016ebd , 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

$\mathcal{D}_{t}=\tfrac{1}{3}(3+14t-t^{2})-\tfrac{2}{3}(t-9)(t-1)\partial_{t}\,,$
(110)

to write $\mathcal{W}_{s}(t)$ in the region $0<t<1$ as

$\mathcal{W}_{S}(t)=\begin{pmatrix}\Psi_{1}(t)&\Psi_{2}(t)\\\
\mathcal{D}_{t}\Psi_{1}(t)&\mathcal{D}_{t}\Psi_{1}(t)\\\ \end{pmatrix}\,,$
(111)

where

$\begin{split}\Psi_{1}(t)&=\frac{4}{[(3-\sqrt{t})(1+\sqrt{t})^{3}]^{1/2}}\,{\operatorname{K}}\left(\frac{t_{14}(t)t_{23}(t)}{t_{13}(t)t_{24}(t)}\right)\,,\\\
\Psi_{2}(t)&=\frac{4i}{[(3-\sqrt{t})(1+\sqrt{t})^{3}]^{1/2}}\,{\operatorname{K}}\left(\frac{t_{12}(t)t_{34}(t)}{t_{13}(t)t_{24}(t)}\right)\,,\end{split}$
(112)

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

$t_{1}(t)=-4\,,\quad t_{2}(t)=-(1+\sqrt{t})^{2}\,,\quad
t_{3}(t)=-(1-\sqrt{t})^{2}\,,\quad t_{4}(t)=0\,.$ (113)

At this point, we can define a new set of master integrals $\cal
T_{1}(t),\;\cal T_{2}(t)$ as follows

$\begin{pmatrix}\cal S_{1}(t)\\\ \cal S_{2}(t)\end{pmatrix}=\cal
W_{S}(t)\begin{pmatrix}\cal T_{1}(t)\\\ \cal T_{2}(t)\end{pmatrix}\,,$ (114)

which by definition fulfil the purely non-homogeneous system of differential
equations

$\partial_{t}\left(\begin{matrix}\cal T_{1}(t)\\\ \cal
T_{2}(t)\end{matrix}\right)=\frac{1}{4\pi i}\begin{pmatrix}-\Psi_{2}(t)\\\
\phantom{-}\Psi_{1}(t)\end{pmatrix}\,.$ (115)

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 $\Gamma_{1}(6)$
Adams:2017ejb ; Broedel:2019kmn . The modular parameter $\tau$ is as usual
defined as the ratio of the two periods,

$\tau(t)=\frac{\Psi_{2}(t)}{\Psi_{1}(t)}\,,$ (116)

and we find that the $j$-invariant of the family of elliptic curves associated
to the maximal cut of the sunrise graph reads

$j(t)=\frac{\left(t^{4}-12t^{3}+30t^{2}+228t+9\right)^{3}}{1728(t-9)^{2}(t-1)^{6}t}\,.$
(117)

It turns out that eq. (116) can be inverted to give

$t(\tau)=9\frac{\eta(\tau)^{4}\eta(6\tau)^{8}}{\eta(2\tau)^{8}\eta(3\tau)^{4}}\,,$
(118)

where $\eta(\tau)$ was defined in eq. (59). The function $t(\tau)$ is a
Hauptmodul for $\Gamma_{1}(6)$, and so it is invariant under the corresponding
modular transformations. $\Psi_{1}(t)$, instead, is a modular form of weight
$n=1$ for $\Gamma_{1}(6)$. Starting from these observations, one can prove
that a basis for $\mathcal{M}_{n}(\Gamma_{1}(6))$ is given by the functions
Broedel:2018rwm

$f_{n,p}(\tau)=\Psi_{1}(t(\tau))^{n}\,t(\tau)^{p}\,,\quad 0\leq p\leq n\,,$
(119)

where we assume $f_{0,0}(\tau)=1$. Finally we introduce the following notation
for iterated integrals of modular forms for $\Gamma_{1}(6)$,

$\mathcal{I}\hskip-2.27626pt\left(\begin{smallmatrix}n_{1}&\ldots&n_{k}\\\
p_{1}&&p_{k}\end{smallmatrix};\tau\right)=I(f_{n_{1},p_{1}},\ldots,f_{n_{k},p_{k}};\tau)=\int_{i\infty}^{\tau}d\tau^{\prime}\,f_{n_{1},p_{1}}(\tau^{\prime})\,\mathcal{I}(f_{n_{2},p_{2}},\ldots,f_{n_{k},p_{k}};\tau^{\prime})\,.$
(120)

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 $\Gamma_{1}(6)$:

$\begin{split}\mathcal{T}_{1}(t(\tau))&\,=\frac{{\rm
Cl}_{2}(\pi/3)}{2\pi}-\frac{1}{24\pi^{2}}\,\left[\mathcal{I}\hskip-2.27626pt\left(\begin{smallmatrix}3&0\\\
3&0\end{smallmatrix};\tau\right)-10\,\mathcal{I}\hskip-2.27626pt\left(\begin{smallmatrix}3&0\\\
2&0\end{smallmatrix};\tau\right)+9\,\mathcal{I}\hskip-2.27626pt\left(\begin{smallmatrix}3&0\\\
1&0\end{smallmatrix};\tau\right)\right]\,,\\\
\mathcal{T}_{2}(t(\tau))&\,=\frac{1}{24\pi^{2}}\,\left[\mathcal{I}\hskip-2.27626pt\left(\begin{smallmatrix}3\\\
3\end{smallmatrix};\tau\right)-10\,\mathcal{I}\hskip-2.27626pt\left(\begin{smallmatrix}3\\\
2\end{smallmatrix};\tau\right)+9\,\mathcal{I}\hskip-2.27626pt\left(\begin{smallmatrix}3\\\
1\end{smallmatrix};\tau\right)\right]\,,\end{split}$ (121)

where

${\rm
Cl}_{2}(x)\,=\,\frac{i}{2}(\text{Li}_{2}(e^{-ix})-\text{Li}_{2}(e^{ix}))\ .$
(122)

is the Clausen function. We stress that eq. (121) 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 $\mathcal{S}_{1}$,
$\mathcal{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 $\mathcal{T}_{1}$ and
$\mathcal{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 $\Gamma(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 $\tau$ to the
fundamental domain, whenever needed. Once that is done, we can evaluate the
relevant iterated integrals by very fast converging $q$-expansions. Equation
(121) is written in terms of modular forms for $\Gamma_{1}(6)$, which is not a
normal subgroup of $\operatorname{SL}(2,\mathbb{Z})$ and, as it was discussed
in Section 2.2, we do not expect this set of iterated integrals to be closed
under $\operatorname{SL}(2,\mathbb{Z})$ transformations. $\Gamma(6)$, instead,
is a normal subgroup of $\operatorname{SL}(2,\mathbb{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 parity-invariant spanning set we are
also able to make all real and imaginary parts explicit, at least as long as
$\tau$ is purely imaginary.

It is straightforward to match eq. (121) to the basis defined in eq. (24) by
matching the corresponding $q$-expansions. We obtain888We stress here that,
since $\mathcal{M}_{n}(\Gamma_{1}(6))\subset\mathcal{M}_{n}(\Gamma(6))$, this
is always possible.

$\begin{split}\mathcal{T}^{(0,1)}_{1}(t(\tau))&\,=\frac{{\rm
Cl}_{2}(\pi/3)}{2\pi}-2\left[I(\small{\texteisfont{a}}_{3,6,1,0},1;\tau)-5I(\small{\texteisfont{a}}_{3,6,1,3},1;\tau)+5I(\small{\texteisfont{a}}_{3,6,2,3},1;\tau)\right]\,,\\\
\mathcal{T}^{(0,1)}_{2}(t(\tau))&\,=-\frac{i}{\pi}\left[I(\small{\texteisfont{a}}_{3,6,1,0};\tau)-5I(\small{\texteisfont{a}}_{3,6,1,3};\tau)+5I(\small{\texteisfont{a}}_{3,6,2,3};\tau)\right]\,,\end{split}$
(123)

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.
(123) are purely real or imaginary, which guarantees that the master integrals
$\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ obtained through eq. (114) are real,
as expected. We are now ready to focus on the analytic continuation and
numerical evaluation of eq. (123) for all values of $t$.

###  8.1 Analytic continuation

In the rest of this section we will study the behaviour of the solution in eq.
(123) 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 $\tau$ to the fundamental domain in each
region. Before doing so, we show how the homogeneous solutions in eq. (112),
and therefore $\tau$, can be extended from the region $0<t<1$ to arbitrary
(real) values of $t$. Indeed, the analytical expression for $\tau$ itself
depends on the homogeneous solutions of the differential equation (107), which
have been given explicitly in eq. (112) for $0<t<1$. Since eq. (107) has
regular singular points at $t\in\\{0,1,9,\pm\infty\\}$, in general we expect
eq. (112) 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. (112).

While this is standard material, we report here the main steps for convenience
of the reader. First of all, let us define

$\lambda_{1}=\frac{t_{14}(t)t_{23}(t)}{t_{13}(t)t_{24}(t)}\,,\qquad\lambda_{2}=\frac{t_{13}(t)t_{24}(t)}{t_{12}(t)t_{34}(t)}\,,$
$\lambda_{3}=\frac{1}{2}\left(1-\frac{t^{2}-6t-3}{(t-1)\sqrt{t_{12}(t)t_{13}(t)}}\right)\,,$

where the roots are defined in eq. (113). In order to highlight the local
character of the solutions in eq. (112), we slightly change of notation and
write the local solutions as

$\begin{split}\Phi^{(0,1)}_{1}(t)&=\frac{4}{\sqrt{t_{13}t_{24}}}\,{\operatorname{K}}\left(\lambda_{1}\right)\,,\qquad\Phi_{2}^{(0,1)}(t)=i\,\frac{4}{\sqrt{t_{13}t_{24}}}\,{\operatorname{K}}\left(1-\lambda_{1}\right)\,,\end{split}$
(124)

where we are using now the letter $\Phi$ instead of the letter $\Psi$ and we
added the superscript $(0,1)$ to indicate that this solution is valid only for
$0<t<1$. The reason for this change of notation will become clear soon. In
order to obtain homogeneous solutions valid for all values of $t$, we define
three pairs of auxiliary functions $\Phi_{i}^{(a,b)}(t)$, which are all local
homogeneous solutions of eq. (107) and which are well defined and purely real
(or imaginary) in the corresponding region $a<t<b$. Following ref.
Remiddi:2016gno ; Bogner:2017vim  we can take

$\begin{split}&\Phi_{1}^{(1,9)}(t)=\frac{4}{\sqrt{t_{14}t_{23}}}\,{\operatorname{K}}\left(1/\lambda_{1}\right)\,,\qquad\Phi_{2}^{(1,9)}(t)=i\,\frac{4}{\sqrt{t_{14}t_{23}}}\,{\operatorname{K}}\left(1-1/\lambda_{1}\right)\,.\\\
&\Phi_{1}^{(9,\infty)}(t)=\frac{4}{\sqrt{t_{12}t_{34}}}\,{\operatorname{K}}\left(\lambda_{2}\right)\,,\qquad\Phi_{2}^{(9,\infty)}(t)=i\,\frac{4}{\sqrt{t_{12}t_{34}}}\,{\operatorname{K}}\left(1-\lambda_{2}\right)\,.\\\
&\Phi_{1}^{(-\infty,0)}(t)=\frac{4}{(t_{12}t_{34}t_{13}t_{24})^{1/4}}\,{\operatorname{K}}\left(\lambda_{3}\right)\,,\qquad\Phi_{2}^{(-\infty,0)}(t)=i\,\frac{1}{(t_{12}t_{34}t_{13}t_{24})^{1/4}}\,{\operatorname{K}}\left(1-\lambda_{3}\right)\,.\end{split}$
(125)

These functions form a complete set of local homogeneous solutions from which
we can patch together the global solutions $\Psi_{j}(t)$ with the correct
analyticity properties for all values of $t$:

$\begin{split}&\begin{pmatrix}\Psi_{1}(t)&\Psi_{2}(t)\end{pmatrix}=\begin{pmatrix}\Phi_{1}^{(a,b)}(t)&\Phi_{2}^{(a,b)}(t)\end{pmatrix}M^{(a,b)}\,\quad\mbox{for}\quad
a<t<b\,,\end{split}$ (126)

where $M^{(a,b)}$ are $2\times 2$ matrices with numerical entries. By
definition, we clearly must have

$M^{(0,1)}=\begin{pmatrix}1&0\\\ 0&1\end{pmatrix}\,,$

which allows us to recover the original definition in eq. (112) for $0<t<1$.
The remaining matrices can be found by imposing that the homogeneous solution
is analytic in the neighbourhood of each regular singular point $t=a$ (with
possibly logarithmic singularities at $t=a$). With our choice of functions
$\Phi_{i}^{(a,b)}$ they read

$M^{(1,9)}=\begin{pmatrix}1&0\\\ 3&1\end{pmatrix}\,,\quad
M^{(9,\infty)}=\begin{pmatrix}-2&-1\\\ 3&1\end{pmatrix}\,,\quad
M^{(-\infty,0)}=\begin{pmatrix}1&1/2\\\ 0&2\end{pmatrix}\,.$ (127)

Once the continued homogenous solutions $\Psi_{j}$ have been determined, the
complete matrix of solutions $\mathcal{W}_{S}(t)$ of eq. (111) can be extended
from the region $0<t<1$ to all real values of $t$ by simply using the
piecewise definition of $\Psi_{j}(t)$ of eq. (LABEL:eq:global_SR_solution) in
eq. (111). In the same way, we can extend the definition of $\tau(t)$ in eq.
(116) to all values of $t$.

With our correctly continued homogeneous solutions, we can study when we need
to perform a modular transformation to remap $\tau$ 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
$\tau$ as function of $t$, see fig. 2.

Figure 2: The real and imaginary parts of $\tau$.

It is easy to see that our definition for $\tau$ 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 $\tau$ 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 $\gamma_{t}\in\operatorname{SL}(2,\mathbb{Z})$,
such that $\gamma_{t}^{-1}\cdot\tau(t)$ lies in the fundamental domain.
Defining the four special points

$\begin{split}&r_{1}=-3\,,\quad r_{2}=5-2^{1/3}\,4\,,\quad
r_{3}=3+2\sqrt{3}(1-3^{1/4}\sqrt{2})\,,\\\ &r_{4}=3+2\sqrt{3}\,,\quad
r_{5}=3+2\sqrt{3}(1+3^{1/4}\sqrt{2})\,,\end{split}$ (128)

we find

$\gamma_{t}=\left\\{\begin{array}[]{ll}\gamma_{0a}\,,&\textrm{ if
}t<r_{1}\,,\\\ \gamma_{0b}\,,&\textrm{ if }r_{1}\leq t<r_{2}\,,\\\
\gamma_{1a}\,,&\textrm{ if }r_{2}<t<r_{3}\,,\\\ \gamma_{1b}\,,&\textrm{ if
}r_{3}<t<1\,,\\\ \gamma_{2a}\,,&\textrm{ if }1\leq t<r_{4}\,,\\\
\gamma_{2b}\,,&\textrm{ if }r_{4}\leq t<9\,,\\\ \gamma_{3a}\,,&\textrm{ if
}9\leq t<r_{5}\,,\\\ \gamma_{3b}\,,&\textrm{ if }t\geq
r_{5}\,,\end{array}\right.$ (129)

with

$\begin{split}&\gamma_{0a}\,=\left(\begin{smallmatrix}-1&1\\\
-2&1\end{smallmatrix}\right)\,,\quad\gamma_{0b}\,=\left(\begin{smallmatrix}0&-1\\\
1&-1\end{smallmatrix}\right)\,,\quad\gamma_{1a}\,=\left(\begin{smallmatrix}1&0\\\
0&1\end{smallmatrix}\right)\,,\quad\gamma_{1b}\,=\left(\begin{smallmatrix}0&-1\\\
1&0\end{smallmatrix}\right)\,,\\\
&\gamma_{2a}\,=\left(\begin{smallmatrix}0&-1\\\
1&-3\end{smallmatrix}\right)\,,\quad\gamma_{2b}\,=\left(\begin{smallmatrix}-1&0\\\
-3&-1\end{smallmatrix}\right)\,,\quad\gamma_{3a}\,=\left(\begin{smallmatrix}-1&1\\\
-3&2\end{smallmatrix}\right)\,,\quad\gamma_{3b}\,=\left(\begin{smallmatrix}-1&-1\\\
-2&-3\end{smallmatrix}\right)\,.\end{split}$ (130)

Recalling the expression for the $j$-invariant in eq. (117), we see that the
apparently random special points $t=r_{i}$, $i-1,..,5$ where the
$\gamma_{t}\in\operatorname{SL}(2,\mathbb{Z})$ is discontinuous, are always
points $t_{i}$ such that $j(t_{i})\in\\{0,1,\infty\\}$, 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}\in\\{0,1,9,\infty\\}$, 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...\leq
t<0.0167...$, the relevant solutions (either the ones valid in the range
$(0,1)$ or the ones valid in the range $(-\infty,0)$) can be taken as they
are, without the need of remapping the corresponding $\tau$ to the fundamental
domain. By defining $\tau_{j}=\gamma_{j}^{-1}\cdot\tau$, we can plot the
remapped value in function of $t$ and check that it indeed always lies in the
fundamental domain, see fig. 3.

Figure 3: The real and imaginary parts of $\gamma_{j}^{-1}\cdot\tau^{(a,b)}$.

With the matrices defined in eq. (129), we can now perform the relevant
transformations on the iterated integrals in eq. (123) and obtain new
combinations of integrals evaluated in the fundamental domain. Clearly, there
is nothing to do in region $1a$, as $\gamma_{1a}=\mathds{1}$. The first non-
trivial transformation is in region $1b$, i.e. for $r_{3}<t<1$ where the
transformation induced by $\gamma_{1b}\in\operatorname{SL}(2,\mathbb{Z})$ has
to applied in order to map $\tau$ back to the fundamental domain. We define
$\tau=\gamma_{1b}\cdot\tau_{1b}$, and apply the algorithms described in this
paper to obtain

$\displaystyle\mathcal{T}^{(1b)}_{1}(t(\tau_{1b}))$
$\displaystyle\,=\frac{1}{\pi}\left[I(\small{\texteisfont{b}}_{3,6,0,1};\tau_{1b})-5I(\small{\texteisfont{b}}_{3,6,3,1};\tau_{1b})+5I(\small{\texteisfont{b}}_{3,6,3,2};\tau_{1b})\right]$
(131) $\displaystyle\mathcal{T}^{(1b)}_{2}(t(\tau_{1b}))$
$\displaystyle\,=-\frac{i}{8}\left[\pi-16I(\small{\texteisfont{b}}_{3,6,0,1},1;\tau_{1b})+80I(\small{\texteisfont{b}}_{3,6,3,1},1;\tau_{1b})-80I(\small{\texteisfont{b}}_{3,6,3,2},1;\tau_{1b})\right]\,,$

where now $\tau_{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 $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ as a function of $t$. Since
the result is now always evaluated for $\tau(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$
terms of the expansion and running for a few seconds on a laptop computer we
can get the plots in fig. 4.

(a) The master integral $\mathcal{S}_{1}$.

(b) The master integral $\mathcal{S}_{2}$.

Figure 4: The real and imaginary parts of the two master integrals of the
sunrise graph.

Our numerical results have been checked thoroughly against the known
analytical results in ref. Remiddi:2016gno  and the numerical program SecDec
Borowka:2015mxa ; Borowka:2017idc . 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 Remiddi:2016gno ; Bogner:2017vim  and the
three-loop equal-mass banana graph Bloch:2014qca ; Primo:2017ipr ;
Broedel:2019kmn .

##  9 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 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.

## Acknowledgments

We thank Nils Matthes for interesting discussions and Brenda Penante for
collaboration during the early stages of this project. This research was
supported by the ERC starting grant 637019 “MathAm” and by the Royal Society
through the grant URF\R1\191125.

##  Appendix A The modular group $\operatorname{SL}(2,\mathbb{Z})$

In this appendix we review some useful algorithms related to the modular group
$\operatorname{SL}(2,\mathbb{Z})$. We start by an algorithm to decompose any
element of $\operatorname{SL}(2,\mathbb{Z})$ into a finite product of the
generators $S$, $T$ and $-\mathds{1}$. In a second subsection we discuss, for
some $\tau\in\mathbb{H}$, how to find
$\gamma\in\operatorname{SL}(2,\mathbb{Z})$ such that $\gamma\cdot\tau$ lies in
the fundamental domain for $\operatorname{SL}(2,\mathbb{Z})$.

###  A.1 Decomposing elements into a product of $S$’s and $T$’s

Consider $\gamma=\left(\begin{smallmatrix}\gamma_{11}&\gamma_{12}\\\
\gamma_{21}&\gamma_{22}\end{smallmatrix}\right)\in\operatorname{SL}(2,\mathbb{Z})$.
We start by describing what happens if we multiply $\gamma$ from the left by
$S$ or $T$. If $N$ is any integer, we have

$S\gamma=\left(\begin{smallmatrix}-\gamma_{21}&-\gamma_{22}\\\
\gamma_{11}&\gamma_{12}\end{smallmatrix}\right)\textrm{~{}~{}and~{}~{}}T^{N}\gamma=\left(\begin{smallmatrix}\gamma_{11}+N\gamma_{21}&\gamma_{12}+N\gamma_{22}\\\
\gamma_{21}&\gamma_{22}\end{smallmatrix}\right)\,.$ (132)

We now describe the algorithm to decompose $\gamma$ into a product of
generators.

  1. -

Step 1: If $|\gamma_{11}|<|\gamma_{21}|$, replace $\gamma$ by

$\gamma=\left(\begin{smallmatrix}\gamma_{11}&\gamma_{12}\\\
\gamma_{21}&\gamma_{22}\end{smallmatrix}\right)\longrightarrow
S\gamma=\left(\begin{smallmatrix}-\gamma_{21}&-\gamma_{22}\\\
\gamma_{11}&\gamma_{12}\end{smallmatrix}\right)\,.$ (133)

After this step, we can assume that $|\gamma_{11}|>|\gamma_{21}|$.

  2. -

Step 2: If $\gamma_{21}=0$, go to Step 3. Otherwise, we use Euclidean division
to write

$\gamma_{11}=\gamma_{21}\,q+r\,,\qquad 0\leq r<|\gamma_{21}|\,,$ (134)

and we replace $\gamma$ by

$\gamma\longrightarrow
ST^{-q}\gamma=\left(\begin{smallmatrix}-\gamma_{21}&-\gamma_{22}\\\
\gamma_{11}-q\gamma_{21}&\gamma_{12}-q\gamma_{22}\end{smallmatrix}\right)=\left(\begin{smallmatrix}-\gamma_{21}&-\gamma_{22}\\\
r&\gamma_{12}-q\gamma_{22}\end{smallmatrix}\right)\,.$ (135)

If $r\neq 0$, return to Step 2. At each iteration the value of $r$ is strictly
decreasing, and since $r>0$, it will eventually reach 0.

  3. -

Step 3: The matrix now has the form
$\gamma=\left(\begin{smallmatrix}\gamma_{11}&\gamma_{12}\\\
0&\gamma_{22}\end{smallmatrix}\right)$. Since $\det\gamma=1$, we must have
$\gamma_{11}=\gamma_{22}=\pm 1$. Replace $\gamma$ by

$\gamma\longrightarrow
T^{-\gamma_{11}\gamma_{12}}\gamma=\left(\begin{smallmatrix}\gamma_{11}&0\\\
0&\gamma_{22}\end{smallmatrix}\right)=\left(\begin{smallmatrix}\pm 1&0\\\
0&\pm 1\end{smallmatrix}\right).$ (136)

At the end of this procedure we have constructed a product of $S$’s and $T$’s
which is the inverse of our original $\gamma$ (up to a sign). This easily
fournishes the desired decomposition into the generators $S$ and $T$, using
the fact that $S^{-1}=S$ and $(T^{N})^{-1}=T^{-N}$.

###  A.2 Mapping a point in the upper half-plane to the fundamental domain

In this appendix we describe how to find for every $\tau\in\mathbb{H}$ an
element $\gamma\in\operatorname{SL}(2,\mathbb{Z})$ such that $\gamma\cdot\tau$
lies in the fundamental domain $\cal F$ defined in eq. (5). We will construct
$\gamma$ iteratively using the following algorithm:

  1. -

Step 1: We can easily find an integer $N$ such that
$-\frac{1}{2}\leq\textrm{Re }\tau+N<\frac{1}{2}$. Replace $\tau$ by

$\tau\longrightarrow T^{N}\cdot\tau=\tau+N\,.$ (137)

  2. -

Step 2: If $|\tau|>1$, we are done. Otherwise, replace $\tau$ by

$\tau\longrightarrow S\cdot\tau=-1/\tau\,.$ (138)

If after this replacement $\tau\in\cal F$, we are done. Otherwise we return to
Step 1. At each iteration the imaginary part of $\tau$ is strictly increasing,
and therefore the algorithm will eventually terminate.

##  Appendix B Results for the elliptic ${}_{2}F_{1}$ function

In this appendix we present the results for the two master integrals $T_{1}$
and $T_{2}$ in eq. (41), in each of the six regions defined in eq. (94). In
each region the master integrals can be written as (cf. eq. (93))

$\left(\begin{array}[]{c}T_{1}(z)\\\
T_{2}(z)\end{array}\right)=S_{j\alpha}(z)\left(\begin{array}[]{c}U_{1}^{{j\alpha}}(\tau_{j\alpha})\\\
U_{2}^{{j\alpha}}(\tau_{j\alpha})\end{array}\right)\,,\qquad
j\in\\{1,2,3\\}\,,\,\,\alpha\in\\{a,b\\}\,.$ (139)

In the following we list the matrices $S_{j\alpha}(z)$ and the first few terms
in the $\epsilon$-expansion of the functions
$U_{i}^{{j\alpha}}(\tau_{j\alpha})$ for each region.

$S_{1a}(z)=\left(\begin{array}[]{cc}\frac{2i\operatorname{K}\left(\frac{1}{1-z}\right)}{\sqrt{1-z}}&0\\\
-\frac{2i\sqrt{1-z}\operatorname{E}\left(\frac{1}{1-z}\right)}{z(1+6\epsilon)}-\frac{2i(2z-1)\operatorname{K}\left(\frac{1}{1-z}\right)}{3\sqrt{1-z}z(1+6\epsilon)}&-\frac{\pi\sqrt{1-z}}{z(1+6\epsilon)\operatorname{K}\left(\frac{1}{1-z}\right)}\\\
\end{array}\right)\,.$ (140)

$\begin{split}U_{1}^{1a}(\tau_{1a})&\,=-2i\pi
I(1;\tau_{1a})+i\pi\,\epsilon\left[1-2\log\left(1-\frac{1}{z}\right)I(1;\tau_{1a})\right]\\\
&\,+i\pi\,\epsilon^{2}\Bigg{[}-360I(1,\small{\texteisfont{a}}_{4,2,0,0},1;\tau_{1a})+\left(\frac{4\pi^{2}}{3}-\log^{2}\left(1-\frac{1}{z}\right)\right)I(1;\tau_{1a})\\\
&\,+\log\left(1-\frac{1}{z}\right)-\frac{30}{\pi^{2}}\,\zeta_{3}\Bigg{]}+\cal
O(\epsilon^{3})\,.\end{split}$ (141)

$\begin{split}U_{2}^{1a}(\tau_{1a})&\,=-1-\epsilon\,\log\left(1-\frac{1}{z}\right)+\epsilon^{2}\,\Bigg{[}-180I(\small{\texteisfont{a}}_{4,2,0,0},1;\tau_{1a})-\frac{1}{2}\log^{2}\left(1-\frac{1}{z}\right)\\\
&\,+\frac{2\pi^{2}}{3}\Bigg{]}+\cal O(\epsilon^{3})\,.\end{split}$ (142)

$S_{1b}(z)=\left(\begin{array}[]{cc}\frac{2\operatorname{K}\left(\frac{z}{z-1}\right)}{\sqrt{1-z}}&0\\\
\frac{2\sqrt{1-z}\operatorname{E}\left(\frac{z}{z-1}\right)}{z(1+6\epsilon)}-\frac{2(2-z)\operatorname{K}\left(\frac{z}{z-1}\right)}{3\sqrt{1-z}z(1+6\epsilon)}&-\frac{i\pi\sqrt{1-z}}{z(1+6\epsilon)\operatorname{K}\left(\frac{z}{z-1}\right)}\\\
\end{array}\right)\,.$ (143)

$\begin{split}U_{1}^{1b}(\tau_{1b})&\,=1+\epsilon\,\Bigg{[}\log\left(1-\frac{1}{z}\right)-2\pi^{2}I(1;\tau_{1b})\Bigg{]}+\epsilon^{2}\,\Bigg{[}180I(1,\small{\texteisfont{a}}_{4,2,0,0};\tau_{1b})\\\
&\,-2\pi^{2}\log\left(1-\frac{1}{z}\right)I(1;\tau_{1b})+\frac{1}{2}\log^{2}\left(1-\frac{1}{z}\right)+\frac{\pi^{2}}{6}\Bigg{]}+\cal
O(\epsilon^{3})\,.\end{split}$ (144)

$\begin{split}U_{2}^{1b}(\tau_{1b})&\,=i\pi\,\epsilon+i\pi\,\epsilon^{2}\Bigg{[}\log\left(1-\frac{1}{z}\right)-\frac{90}{\pi^{2}}\,I(\small{\texteisfont{a}}_{4,2,0,0};\tau_{1b})\Bigg{]}+\cal
O(\epsilon^{3})\,.\end{split}$ (145)

$S_{2a}(z)=\left(\begin{array}[]{cc}2\operatorname{K}(z)&0\\\
\frac{2\operatorname{E}(z)}{z(1+6\epsilon)}-\frac{2(2-z)\operatorname{K}(z)}{3z(1+6\epsilon)}&-\frac{i\pi}{z(1+6\epsilon)\operatorname{K}(z)}\\\
\end{array}\right)\,.$ (146)

$\begin{split}U_{1}^{2a}(\tau_{2a})&\,=1+\epsilon\,\Bigg{[}\log\left(\frac{1-z}{z}\right)-2\pi^{2}I(1;\tau_{2a})\Bigg{]}+\epsilon^{2}\,\Bigg{[}180I(1,\small{\texteisfont{a}}_{4,2,0,0};\tau_{2a})\\\
&\,-2\pi^{2}\log\left(\frac{1-z}{z}\right)I(1;\tau_{2a})+\frac{1}{2}\log^{2}\left(\frac{1-z}{z}\right)+\frac{\pi^{2}}{6}\Bigg{]}+\cal
O(\epsilon^{3})\,.\end{split}$ (147)

$\begin{split}U_{2}^{2a}(\tau_{2a})&\,=i\pi\,\epsilon+i\pi\,\epsilon^{2}\,\Bigg{[}\log\left(\frac{1-z}{z}\right)-\frac{90}{\pi^{2}}\,I(\small{\texteisfont{a}}_{4,2,0,0};\tau_{2a})\Bigg{]}+\cal
O(\epsilon^{3})\,.\end{split}$ (148)

$S_{2b}(z)=\left(\begin{array}[]{cc}-2i\operatorname{K}(1-z)&0\\\
\frac{2i\operatorname{E}(1-z)}{z(1+6\epsilon)}-\frac{2i(z+1)\operatorname{K}(1-z)}{3z\
(1+6\epsilon)}&\frac{\pi}{z(1+6\epsilon)\operatorname{K}(1-z)}\\\
\end{array}\right)\,.$ (149)

$\begin{split}U_{1}^{2b}(\tau_{2b})&\,=+2i\pi
I(1;\tau_{2b})+i\pi\,\epsilon\,\Bigg{[}2\log\left(\frac{1-z}{z}\right)I(1;\tau_{2b})-1\Bigg{]}\\\
&\,+i\pi\,\epsilon^{2}\,\Bigg{[}360I(1,\small{\texteisfont{a}}_{4,2,0,0},1;\tau_{2b})+\left(\log^{2}\left(\frac{1-z}{z}\right)-\frac{4i\pi^{3}}{3}\
\right)I(1;\tau_{2b})\\\
&\,-\log\left(\frac{1-z}{z}\right)+\frac{30}{\pi^{2}}\,\zeta_{3}\Bigg{]}+\cal
O(\epsilon^{3})\,.\end{split}$ (150)

$\begin{split}U_{2}^{2b}(\tau_{2b})&\,=1+\epsilon\,\log\left(\frac{1-z}{z}\right)+\epsilon^{2}\,\Bigg{[}180I(\small{\texteisfont{a}}_{4,2,0,0},1;\tau_{2b})+\frac{1}{2}\log^{2}\left(\frac{1-z}{z}\right)\\\
&\,-\frac{2\pi^{2}}{3}\Bigg{]}+\cal O(\epsilon^{3})\,.\end{split}$ (151)

$S_{3a}(z)=\left(\begin{array}[]{cc}\frac{2i\operatorname{K}\left(\frac{z-1}{z}\right)}{\sqrt{z}}&0\\\
-\frac{2i\operatorname{E}\left(\frac{z-1}{z}\right)}{\sqrt{z}(1+6\epsilon)}+\frac{2i(z+1)\operatorname{K}\left(\frac{z-1}{z}\right)}{3z^{3/2}(1+6\epsilon)}&-\frac{\pi}{\sqrt{z}(1+6\epsilon)\operatorname{K}\left(\frac{z-1}{z}\right)}\\\
\end{array}\right)\,.$ (152)

$\displaystyle U_{1}^{3a}(\tau_{3a})$ $\displaystyle\,=1-2i\pi
I(1;\tau_{3a})+\epsilon\,\Bigg{[}\log\left(1-\frac{1}{z}\right)-\left(2\pi^{2}+2i\pi\log\left(1-\frac{1}{z}\right)\right)I(1;\tau_{3a})\Bigg{]}$
$\displaystyle\,+\epsilon^{2}\,\Bigg{[}180I(1,\small{\texteisfont{a}}_{4,2,0,0};\tau_{3a})-360i\pi
I(1,\small{\texteisfont{a}}_{4,2,0,0},1;\tau_{3a})$ (153)
$\displaystyle\,-\left(i\pi\log^{2}\left(1-\frac{1}{z}\right)+2\pi^{2}\log\left(1-\frac{1}{z}\right)+\frac{10i\pi^{3}}{3}\right)I(1;\tau_{3a})+\frac{1}{2}\log^{2}\left(1-\frac{1}{z}\right)$
$\displaystyle\,-\frac{\pi^{2}}{3}-\frac{90i}{3\pi}\,\zeta_{3}\Bigg{]}+\cal
O(\epsilon^{3})\,.$

$\displaystyle U_{2}^{3a}(\tau_{3a})$
$\displaystyle\,=-1+\epsilon\Bigg{[}i\pi-\log\left(1-\frac{1}{z}\right)\Bigg{]}+\epsilon^{2}\,\Bigg{[}-\frac{90i}{\pi}\,I(\small{\texteisfont{a}}_{4,2,0,0};\tau_{3a})-180I(\small{\texteisfont{a}}_{4,2,0,0},1;\tau_{3a})$
$\displaystyle\,-\frac{1}{2}\log^{2}\left(1-\frac{1}{z}\right)+i\pi\log\left(1-\frac{1}{z}\right)+\frac{5\pi^{2}}{3}\Bigg{]}+\cal
O(\epsilon^{3})\,.$ (154)

$S_{3b}(z)=\left(\begin{array}[]{cc}-\frac{2\operatorname{K}\left(\frac{1}{z}\right)}{\sqrt{z}}&0\\\
\frac{2(2z-1)\operatorname{K}\left(\frac{1}{z}\right)}{3z^{3/2}(1+6\epsilon)}-\frac{2\operatorname{E}\left(\frac{1}{z}\right)}{\sqrt{z}(1+6\epsilon)}&\frac{i\pi}{\sqrt{z}(1+6\epsilon)\operatorname{K}\left(\frac{1}{z}\right)}\\\
\end{array}\right)\,.$ (155)

$\displaystyle U_{1}^{3b}(\tau_{3b})$ $\displaystyle\,=-1-2i\pi
I(1;\tau_{3b})+\epsilon\,\Bigg{[}-2i\pi\log\left(1-\frac{1}{z}\right)I(1;\tau_{3b})-\log\left(1-\frac{1}{z}\right)+i\pi\Bigg{]}$
$\displaystyle\,+\epsilon^{2}\,\Bigg{[}-180I(1,\small{\texteisfont{a}}_{4,2,0,0};\tau_{3b})-360i\pi
I(1,\small{\texteisfont{a}}_{4,2,0,0},1;\tau_{3b})$ (156)
$\displaystyle\,+\left(\frac{7i\pi^{3}}{3}-i\pi\log^{2}\left(1-\frac{1}{z}\right)\right)I(1;\tau_{3b})-\frac{1}{2}\log^{2}\left(1-\frac{1}{z}\right)+i\pi\log\left(1-\frac{1}{z}\right)$
$\displaystyle\,+\frac{5\pi^{2}}{6}-\frac{180i}{6\pi}\,\zeta_{3}\Bigg{]}+\cal
O(\epsilon^{3})\,.$

$\begin{split}U_{2}^{3b}(\tau_{3b})&\,=-1-\epsilon\,\log\left(1-\frac{1}{z}\right)\epsilon^{2}\,\Bigg{[}\frac{90i}{\pi}\,I(\small{\texteisfont{a}}_{4,2,0,0};\tau_{3b})-180I(\small{\texteisfont{a}}_{4,2,0,0},1;\tau_{3b})\\\