aainstitutetext: Niels Bohr International Academy, University of Copenhagen
Blegdamsvej 17, 2100 København, Denmark

# On Epsilon Factorized Differential Equations for Elliptic Feynman Integrals

Hjalte Frellesvig<EMAIL_ADDRESS>
###### Abstract

In this paper we develop and demonstrate a method to obtain epsilon factorized
differential equations for elliptic Feynman integrals. This method works by
choosing an integral basis with the property that the period matrix obtained
by integrating the basis over a complete set of integration cycles is
diagonal. The method is a generalization of a similar method known to work for
polylogarithmic Feynman integrals. We demonstrate the method explicitly for a
number of Feynman integral families with an elliptic highest sector.

##  1 Introduction

The method of differential equations Barucchi:1973zm ; KOTIKOV1991158 ;
Gehrmann:1999as  is the primary way phenomenologically relevant Feynman
integrals are computed. This method works by taking the derivative with
respect to a kinematic variable (or to several) of the members of a minimal
basis of Feynman integrals for the problem, which we will refer to as master
integrals. The result of taking that derivative may then be mapped back to the
original master integrals using IBP identities Chetyrkin:1981qh ;
Laporta:2001dd , as implemented in for instance the computer codes FIRE
Smirnov:2008iw , Kira Maierhofer:2017gsa , and others Anastasiou:2004vj ;
vonManteuffel:2012np ; Lee:2012cn . This will result in a system of coupled
linear differential equations for the master integrals, which may be solved
with traditional methods. A breakthrough was made with the realization
Henn:2013pwa  that for a large class of Feynman integrals in dimensional
regularization (with $d=d_{0}-2\epsilon$), the solution of the system gets
much simpler, in some cases trivialized, by choosing a special basis with the
property that the differential equation system is epsilon factorized. Denoting
the basis of master integrals as $J$, this means that the differential
equation system in the kinematic variable $x$ may be written as

$\displaystyle\frac{\partial}{\partial x}J$ $\displaystyle=\epsilon A^{(x)}J$
(1)

where $A^{(x)}$ is a matrix dependent on the kinematics but independent of
$\epsilon$.

In the cases discussed in ref. Henn:2013pwa , these matrices have the
additional property that the individual $A^{(x)}$-matrices may be unified as
one matrix $A^{(x)}=\partial M/\partial x$ with the property that the entries
of $M$ consist solely of logarithms of algebraic functions of the kinematics.
A basis of master integrals for which the differential equations have these
two properties (epsilon factorized, dlog form) is known as a canonical basis.

It is a well known fact that not all Feynman integrals can be brought to
canonical form. One class (the simplest) of integrals for which a canonical
form is unobtainable are those referred to as elliptic, a term that comes from
the presence of elliptic curves at a cut surface. Such elliptic Feynman
integrals, of which the fully massive sunrise integrals are a paradigmatic
example, have been a subject of intense study in recent years
Broadhurst:1993mw ; Berends:1993ee ; Bauberger:1994nk ; Bauberger:1994by ;
Bauberger:1994hx ; Caffo:1998du ; Laporta:2004rb ; Groote:2005ay ;
Bailey:2008ib ; MullerStach:2011ru ; brown2011multiple ; CaronHuot:2012ab ;
Groote:2012pa ; Adams:2013nia ; Bloch:2013tra ; Remiddi:2013joa ;
Adams:2014vja ; Broedel:2014vla ; Bloch:2014qca ; Adams:2015gva ;
Broedel:2015hia ; Adams:2015ydq ; Bloch:2016izu ; Passarino:2016zcd ;
Primo:2016ebd ; Remiddi:2016gno ; Bonciani:2016qxi ; Broadhurst:2016myo ;
Adams:2017ejb ; vonManteuffel:2017hms ; Ablinger:2017bjx ; Remiddi:2017har ;
Hidding:2017jkk ; Bourjaily:2017bsb ; Broedel:2017kkb ; Broedel:2017siw ;
Broedel:2017jdo ; Lee:2017qql ; Bogner:2017vim ; Groote:2018rpb ;
Adams:2018bsn ; Adams:2018kez ; Adams:2018yfj ; Ablinger:2018zwz ;
Broedel:2018qkq ; Bogner:2019lfa ; Bourjaily:2020hjv ; Frellesvig:2021vdl ;
Bourjaily:2021vyj . The purpose of this paper is to extend the epsilon
factorized property of the differential equations given by eq. (1) to such
elliptic cases. This has been achieved in the literature for a few specific
examples Adams:2018yfj ; Bogner:2019lfa  already, but this paper will present
a more general algorithm for how to obtain such a basis.

For integrals for which a canonical form is obtainable, several approaches
have been developed to finding it Henn:2014qga ; Lee:2014ioa ; Argeri:2014qva
; Gehrmann:2014bfa ; Hoschele:2014qsa ; Meyer:2016slj ; Frellesvig:2017aai ;
WasserMSc ; Prausa:2017ltv ; Gituliar:2017vzm ; Dlapa:2020cwj ; Henn:2020lye ;
Chen:2020uyk ; Dlapa:2021qsl . One popular approach Henn:2014qga ;
Frellesvig:2017aai ; WasserMSc  includes the requirement that the master
integrals are pure Arkani-Hamed:2010pyv , which means that the value of the
integral on each maximal cut (here in the sense of combinations of unitarity
cuts that fix all degrees of freedom) is just a number, as opposed to a
function of the kinematic variables. For elliptic Feynman integrals performing
a maximal cut in that sense is impossible due to the presence of the elliptic
curve. Yet a natural generalization of that cutting operation is an integral
over a cycle, and such integrals can be performed also for the non-zero genus
surfaces appearing for elliptic Feynman integrals (and beyond).

The algorithm proposed in this paper is the natural generalization of that
maximal cut procedure, or more specifically it is to require that each master
integral is non-vanishing on one and only one of the (members of a set of)
basic cycles. This is similar to the prescriptive unitarity approach proposed
in refs. Bourjaily:2017wjl ; Bourjaily:2021vyj  but here performed at the
integral level for the fully dimensionally regulated integrals.

The use of $d$-dimensional unitarity cuts to investigate Feynman integrals and
simplify their relations, has a long history. This was first done in the
context of dimension shift relations Lee:2012te , and later extended to
differential equations Primo:2016ebd , IBPs Bosma:2017ens , and more
Frellesvig:2017aai ; Primo:2017ipr ; Harley:2017qut . In some of those works
Bosma:2017ens ; Primo:2017ipr  the notion of using a diagonal period matrix as
a condition for obtaining an epsilon factorized differential equation, a
notion central to this work, was prefigured, and so was the suggestion to keep
certain epsilon-dependent prefactors before doing the cut analysis
Frellesvig:2017aai . Thus this paper may rightfully be regarded as the newest
entry that series of works.

The dual vector space nature of the set of integrands $\varphi_{i}$ and
integration contours $\gamma_{j}$ Smirnov:2010hn ; Lee:2013hzt ;
Bitoun:2017nre ; Bosma:2017ens , was clarified in recent work on the
connection between Feynman integrals and the mathematical discipline of
intersection theory Mastrolia:2018uzb ; Frellesvig:2019kgj ;
Frellesvig:2019uqt ; Frellesvig:2020qot ; Mizera:2019vvs ; Weinzierl:2020xyy ;
Caron-Huot:2021xqj ; Caron-Huot:2021iev , where the contours and integrands
are viewed as representatives of a twisted de Rahm homology and cohomology
respectively. This change of perspective makes it natural to look at the
period matrix $P_{ij}\\!\sim\\!\int_{\gamma_{j}}\\!\\!\varphi_{i}$ of pairings
between integrands and contours, and the algorithm proposed here is naturally
formulated in terms of that object.

After introducing the notation in higher detail, we will in sec. 2 do two
examples of the algorithm used in non-elliptic cases, while reformulating it
in a way that generalizes to the elliptic case. Then in sec. 3 we will make
the complete formulation of our algorithm, and in sec. 4 do a number of
successful examples of its use for which many of the resulting expressions are
added as an ancillary file. Finally in sec. 5 we will discuss some loose ends
and open questions, and summarize and conclude. In appendix A we will redo one
of the examples with a different approach.

###  1.1 Notation and conventions

We will in this work encounter the complete elliptic integrals of respectively
first, second and third kind. They are defined as

$\displaystyle K(k^{2})$
$\displaystyle=\int_{0}^{1}\\!\frac{\mathrm{d}x}{\sqrt{(1-x^{2})(1-k^{2}x^{2})}}$
(2) $\displaystyle E(k^{2})$
$\displaystyle=\int_{0}^{1}\\!\frac{\sqrt{1-k^{2}x^{2}}\;\mathrm{d}x}{\sqrt{1-x^{2}}}$
(3) $\displaystyle\Pi(n^{2},k^{2})$
$\displaystyle=\int_{0}^{1}\\!\frac{\mathrm{d}x}{(1-n^{2}x^{2})\sqrt{(1-x^{2})(1-k^{2}x^{2})}}$
(4)

where we use the Mathematica convention in which the arguments are $k^{2}$ and
$n^{2}$. The three complete elliptic integrals have the property that any
integral of the form $\int_{\mathcal{C}}\mathcal{R}(x,Y)\mathrm{d}x$ can be
written in terms of them (along with elementary functions), where
$\mathcal{R}$ is any rational function, $Y$ is the square root of a polynomial
in $x$ of degree three or four, and where $\mathcal{C}$ is any closed contour.

We will in this work only discuss two-loop Feynman integrals, which can be
written as

$\displaystyle I_{\\{a\\}}$
$\displaystyle=\frac{1}{(2\pi)^{d}}\int\frac{\mathrm{d}^{d}k_{1}\mathrm{d}^{d}k_{2}}{\prod_{i}D_{i}^{a_{i}}}$
(5)

where the propagators $D_{i}$ are quadratic functions of the momenta on the
form $k^{2}{-}m^{2}$ and $d$ is the spacetime dimension. As mentioned earlier,
we will consider the integrals in dimensional regularization where
$d=d_{0}-2\epsilon$ with $d_{0}$ being an integer, in this work either $2$ or
$4$.

We will consider our Feynman integrals in Baikov representation, in either its
standard Baikov:1996iu ; Lee:2010wea ; Grozin:2011mt  or loop-by-loop
Frellesvig:2017aai ; Frellesvig:2021vem  version. Thus the integral may be
written

$\displaystyle I_{\\{a\\}}$
$\displaystyle=\mathcal{K}\int_{\mathcal{C}}\frac{u\,\mathrm{d}^{n}x}{\prod_{i}x_{i}^{a_{i}}}$
(6)

where the $x_{i}$ are Baikov variables (that equal the original propagators,
$x_{i}=D_{i}$), the $\alpha_{j}$ are generally irrational powers, and where
$u$ is a multivalued function that will have the general form

$\displaystyle u=\prod_{j}\mathcal{B}_{j}^{\alpha_{j}}$ (7)

where the $\mathcal{B}_{j}$ are polynomial functions called Baikov polynomials
of the $x_{i}$. The prefactor $\mathcal{K}$ will generally be given in terms
of gamma functions of arguments linear in the spacetime dimension. In the
examples we will look at, $\mathcal{K}$ will be a pure function, by which we
will mean a function with the property that the coefficients in the $\epsilon$
expansion has uniform weight Duhr:2014woa  such that the coefficient of
$\epsilon^{n}$ has weight one higher than that of the coefficient of
$\epsilon^{n-1}$ and only numerical prefactors. We will not make any attempts
at generalizing these notions to elliptic objects.

We will in this paper make significant use of generalized unitarity cuts, in
particular maximal cuts which refers to cuts of all genuine propagators (and
thus differs from the notation in the $\mathcal{N}{=}4$ literature, in which a
maximal is defined to fix all degrees of freedom). This is due to the fact
that differential equations Primo:2016ebd , as well as IBP relations
Bosma:2017ens , remain the same after such a cutting procedure, and indeed the
Baikov parametrization is particularly suited for generalized unitarity cuts
at the integral level Frellesvig:2017aai ; Harley:2017qut .

In the examples in this paper, the maximal cut will leave a single integration
to be done, bringing the integrals to a univariate form of eq. (6)

$\displaystyle I_{\\{a\\}}|_{\text{max cut}}\ $
$\displaystyle=\mathcal{K}\int_{\mathcal{C}}u\,\hat{\phi}\,\mathrm{d}z$ (8)

where $z$ is the remaining Baikov variable, and $\hat{\phi}$ is a rational
function of $z$ and the kinematical variables. We will occasionally use the
abbreviation $\phi:=\hat{\phi}\mathrm{d}z$. In general this notation follows
that of the recent work on intersection theory and Feynman integrals such as
refs. Mastrolia:2018uzb ; Frellesvig:2019kgj .

We will also discuss various integration contours. We will name them with the
symbols $\mathcal{C}$ or $\gamma$, and use the notation $\mathcal{C}_{i}$ for
a contour surrounding a pole at $z=i$, and the notation
$\mathcal{C}_{\text{i-j}}$ for a contour surrounding a branch cut between the
points $z=i$ and $z=j$.

###  1.2 Notation - The sunrise integrals

We will in this paper work quite a lot with the sunrise integral with
different mass distributions. Rather than repeating the same definitions over
and over, let us write them here in the most general case:

The three-mass sunrise integral is defined by the complete set of propagators

$\displaystyle D_{1}$ $\displaystyle=k_{2}^{2}-m_{1}^{2}\,,$ $\displaystyle
D_{2}$ $\displaystyle=(k_{1}-k_{2})^{2}-m_{2}^{2}\,,$ $\displaystyle D_{3}$
$\displaystyle=(k_{1}+p)^{2}-m_{3}^{2}\,,$ $\displaystyle D_{4}$
$\displaystyle=(k_{1})^{2}\,,$ $\displaystyle D_{5}$
$\displaystyle=(k_{2}+p)^{2}\,,$ (9)

with $p^{2}=s$, and where only the first three may appear as genuine
propagators.

With these definition, the integrals are given as

$\displaystyle I^{\text{sunrise}}_{a_{1}a_{2}a_{3}a_{4}a_{5}}$
$\displaystyle=\int\frac{\mathrm{d}^{d}k_{1}\mathrm{d}^{d}k_{2}}{(2\pi)^{d}}\frac{D_{4}^{-a_{4}}D_{5}^{-a_{5}}}{D_{1}^{a_{1}}D_{2}^{a_{2}}D_{3}^{a_{3}}}$
(10)

We will always look at the case where $d=2-2\epsilon$, when discussing the
sunrise.

Figure 1: The sunrise integrals

In the loop-by-loop Baikov representation, this integral may be written as

$\displaystyle I^{\text{sunrise}}_{a_{1}a_{2}a_{3}a_{4}0}$
$\displaystyle=\mathcal{K}\int_{\mathcal{C}}\frac{x_{4}^{-a_{4}}\,u\,\mathrm{d}^{4}x}{x_{1}^{a_{1}}x_{2}^{a_{2}}x_{3}^{a_{3}}}$
(11)

where

$\displaystyle u$
$\displaystyle=x_{4}^{\epsilon}\times\big{(}(m^{2}_{1}{-}m^{2}_{2}{+}x_{1}{-}x_{2})^{2}{-}2(m^{2}_{1}{+}m^{2}_{2}{+}x_{1}{+}x_{2})x_{4}{+}x_{4}^{2}\big{)}^{-\frac{1}{2}-\epsilon}$
$\displaystyle\;\;\;\;\times\big{(}(m^{2}_{3}{-}s{+}x_{3})^{2}{-}2(m^{2}_{3}{+}s{+}x_{3})x_{4}{+}x_{4}^{2}\big{)}^{-\frac{1}{2}-\epsilon}$
(12)

is the product of Baikov polynomials raised to their corresponding powers, and
where

$\displaystyle\mathcal{K}$
$\displaystyle=\frac{2^{4\epsilon}}{\pi\Gamma^{2}(\tfrac{1}{2}-\epsilon)}$
(13)

is a pure function. The choice $a_{5}=0$ made in eq. (11) is done in order to
simplify the loop-by-loop parametrization, which in this case is performed by
first integrating out the $k_{2}$-loop. This limitation of the loop-by-loop
Baikov parametrization can be overcome (for instance using a Passarino-Veltman
Passarino:1978jh  inspired approach as described in section 11.1 of ref.
Frellesvig:2019kgj , and for a recent alternative approach to this see ref.
Chen:2022lzr ,) but that additional complication is not necessary for the
present discussion.

On the maximal cut (i.e. of $D_{1}$, $D_{2}$, $D_{3}$, and with $z=D_{4}$),
$u$ reduces to

$\displaystyle u|_{3\times\text{cut}}$
$\displaystyle=z^{\epsilon}\;\big{(}(m^{2}_{1}{-}m^{2}_{2})^{2}{-}2(m^{2}_{1}{+}m^{2}_{2})z{+}z^{2}\big{)}^{-\frac{1}{2}-\epsilon}\;\big{(}(m^{2}_{3}{-}s)^{2}{-}2(m^{2}_{3}{+}s)z{+}z^{2}\big{)}^{-\frac{1}{2}-\epsilon}$
(14)

meaning that any integral in the sunrise family can be written as

$\displaystyle I^{\text{sunrise}}$
$\displaystyle=\mathcal{K}\int_{\mathcal{C}}u|_{3\times\text{cut}}\,\hat{\phi}\,\mathrm{d}z$
(15)

on that cut, where $\hat{\phi}$ will be some rational function of $z$, the
exact form of which will depend on the values of the propagator powers
$a_{i}$.

We note that $u|_{3\times\text{cut}}$ is a pure function divided by the square
root of a polynomial of degree four, a fact we will be using implicitly in the
following.

##  2 Motivation

In this section we will discuss two simple examples of non-elliptic Feynman
integrals, and show how to find their canonical forms with the traditional
method. This method will then be reformulated in terms of diagonalization of
the period matrix, giving a formulation that generalizes directly to the
elliptic case.

###  2.1 The double box

Let us as an illustration of the procedure discuss the well-known fully
massless double box. That family of Feynman integrals (first computed in
Smirnov:1999gc ; Smirnov:1999wz ) was used as the example in ref. Henn:2013pwa
which introduced the concept of epsilon-factorized differential equations. It
is given by the set of propagators

$\displaystyle D_{1}$ $\displaystyle=k_{1}^{2}\,,$ $\displaystyle D_{2}$
$\displaystyle=(k_{1}{+}p_{1})^{2}\,,$ $\displaystyle D_{3}$
$\displaystyle=(k_{1}{+}p_{1}{+}p_{2})^{2}\,,$ $\displaystyle D_{4}$
$\displaystyle=(k_{2}{+}p_{1}{+}p_{2})^{2}\,,$ $\displaystyle D_{5}$
$\displaystyle=(k_{2}{-}p_{4})^{2}\,,$ $\displaystyle D_{6}$
$\displaystyle=k_{2}^{2}\,,$ $\displaystyle D_{7}$
$\displaystyle=(k_{1}{-}k_{2})^{2}\,,$ $\displaystyle D_{8}$
$\displaystyle=(k_{1}{-}p_{4})^{2}\,,$ $\displaystyle D_{9}$
$\displaystyle=(k_{2}{+}p_{1})^{2}\,,$ (16)

of which only the first seven are allowed to appear as genuine propagators.

The family contains eight master integrals of which two are in the highest
sector, and it is those two that will interest us in the following. In ref.
Henn:2013pwa  it was found that a set of integrals giving epsilon-factorized
differential equations are

$\displaystyle J_{1}=cs^{2}tI^{\text{db}}_{111111100}\qquad\text{and}\qquad
J_{2}=cs^{2}I^{\text{db}}_{1111111\scalebox{0.5}[0.75]{$\,-$}10}$ (17)

where $c=\epsilon^{4}s^{-2\epsilon}$ is a pure prefactor to be ignored in the
following. How can we know a priori the two prefactors $s^{2}t$ and $s^{2}$?

One way is to perform the maximal cut (in the sense of cutting the seven
actual propagators) of the integrals in the family. Using the loop-by-loop
Baikov parametrization, this gives

$\displaystyle I^{\text{db}}_{1111111a0}|_{\text{7-cut}}$
$\displaystyle=\frac{-\mathcal{K}}{s^{2}}\int(z-t)^{-1-2\epsilon}z^{-a-1-\epsilon}(s+z)^{\epsilon}\,\mathrm{d}z$
(18)

where

$\displaystyle\mathcal{K}$
$\displaystyle=\frac{4^{\scalebox{0.5}[0.75]{$\,-$}1+2\epsilon}s^{\scalebox{0.5}[0.75]{$\,-$}2\epsilon}\,t^{\epsilon}\,(s{+}t)^{\epsilon}}{\pi^{3}\Gamma^{2}(\tfrac{1}{2}{-}\epsilon)}$
(19)

is a pure function.

The prescription is now to factorize out the pure part of the integrand, which
here corresponds to taking the limit $\epsilon\rightarrow 0$ of the object
under the integral sign. This (ignoring the already pure $\mathcal{K}$) gives
the integrand

$\displaystyle\hat{\Phi}$ $\displaystyle=\frac{-z^{-a}}{s^{2}z(z-t)}$ (20)

Continuing the cutting procedure shows us that $I^{\text{db}}_{111111100}$ has
two residues. One at $z=0$ evaluating to $1/(s^{2}t)$ and the other at $z=t$
evaluating to $-1/(s^{2}t)$. On the other hand
$I^{\text{db}}_{1111111\scalebox{0.5}[0.75]{$\,-$}10}$ also has two residues,
one at $z=t$ evaluating to $-1/(s^{2})$ and the other at $z=\infty$ evaluating
to $1/(s^{2})$. So we see that the two integrals become pure if
$I^{\text{db}}_{111111100}$ is given a prefactor of $s^{2}t$ and
$I^{\text{db}}_{1111111\scalebox{0.5}[0.75]{$\,-$}10}$ a prefactor of $s^{2}$,
as it is given by eq. (17).

A different way of formulating the procedure is by writing up the period
matrix for the two integrals, again after factorizing out and discarding the
pure part of $u$. We pick the two integration cycles as $\gamma_{1}$ being the
cycle $\mathcal{C}_{0}$ surrounding the $(z{=}0)$-pole and $\gamma_{2}$ being
$\mathcal{C}_{\infty}$ surrounding the $(z{=}\infty)$-pole, with the various
contours shown on fig. 2. Giving the two integrals prefactors $f_{1}$ and
$f_{2}$ we get the period matrix

$\displaystyle
P_{ij}=\int_{\gamma_{j}}\\!\\!f_{i}\hat{\Phi}_{i}\mathrm{d}z\quad\Rightarrow\quad
P=2\pi i\left[\begin{array}[]{cc}\frac{f_{1}}{s^{2}t}&0\\\
0&\frac{f_{2}}{s^{2}}\end{array}\right]$ (23)

and requiring the period matrix to be $2\pi i$ times the unit matrix $I$ fixes
the $f_{i}$ to the values found above.

Figure 2: Contours discussed for the double box

Making this procedure more general and algorithmic, we could also write each
of the canonical integrals we want to find as general linear combinations of
the two precanonical intermediate basis integrals $I^{\text{db}}_{111111100}$
and $I^{\text{db}}_{1111111\scalebox{0.5}[0.75]{$\,-$}10}$, that is

$\displaystyle
J_{1}=f_{11}I^{\text{db}}_{111111100}+f_{12}I^{\text{db}}_{1111111\scalebox{0.5}[0.75]{$\,-$}10}$
$\displaystyle
J_{2}=f_{21}I^{\text{db}}_{111111100}+f_{22}I^{\text{db}}_{1111111\scalebox{0.5}[0.75]{$\,-$}10}$
(24)

This gives the period matrix

$\displaystyle P=2\pi
i\left[\begin{array}[]{cc}\frac{f_{11}}{s^{2}t}&\frac{f_{12}}{s^{2}}\\\\[2.84526pt]
\frac{f_{21}}{s^{2}t}&\frac{f_{22}}{s^{2}}\end{array}\right]$ (27)

and once again imposing $P=2\pi iI$ fixes the $f_{ij}$ coefficients uniquely
to

$\displaystyle f_{11}=s^{2}t\,,\qquad f_{12}=0\,,\qquad f_{21}=0\,,\qquad
f_{22}=s^{2}\,,$ (28)

again corresponding to the prefactors given by eq. (17).

###  2.2 The two-mass non-elliptic sunrise

Let us look at another non-elliptic example that looks more similar to the
elliptic examples we will encounter later. That is the two-mass non-elliptic
sunrise integral (sne), defined as in eqs. (9) and (10) but with $m_{3}^{2}=0$
and $m_{1}^{2}=m_{2}^{2}\equiv m^{2}$.

Using the loop-by-loop Baikov parametrization of this integral as in eq. (11),
we may make a univariate representation of the integral on the maximal cut.
Following our prescription we have to factorize out the pure part from the
integrand in order to analyze its behaviour. The integrand then becomes (from
eq. (15))

$\displaystyle U\,\hat{\phi}\qquad\text{with}\qquad
U=\frac{1}{(z-s)\sqrt{z(z-4m^{2})}}$ (29)

If one were to follow the prescription of cutting the last variable in order
to impose a maximum cut of $\pm 1$, one would (formally) first have to get rid
of the square root through rationalization. This can be done by a variable
change such as

$\displaystyle z\;\rightarrow\;m^{2}\frac{(1+y)^{2}}{y}$ (30)

but let us try to proceed without such procedures. The integration plane can
be seen depicted on fig. 3. One can define a basis of independent cycles as a
cycle $\gamma_{1}$ surrounding the $(z=s)$-pole, and one $\gamma_{2}$
surrounding the pole at infinity. Picking the contour around the branch-cut
from $z=0$ to $z=4m^{2}$ as one of the master contours would be equally valid,
but the above choice gives the nicest result.

Figure 3: Contours discussed for the two-mass non-elliptic sunrise

We will make a choice of intermediate basis integrals (precanonicals) as
$I^{\text{sne}}_{11100}$ and
$I^{\text{sne}}_{111\scalebox{0.5}[0.75]{$\,-$}10}$. This corresponds to the
integrands

$\displaystyle\hat{\phi}_{1}=1\qquad\text{and}\qquad\hat{\phi}_{2}=z$ (31)

With this we may compute the integrals on the cycles, giving

$\displaystyle g_{11}$
$\displaystyle=\int_{\gamma_{1}}\\!U\hat{\phi}_{1}\mathrm{d}z=\frac{2\pi
i}{\sqrt{s(s-4m^{2})}}$ $\displaystyle g_{12}$
$\displaystyle=\int_{\gamma_{2}}\\!U\hat{\phi}_{1}\mathrm{d}z=0$
$\displaystyle g_{21}$
$\displaystyle=\int_{\gamma_{1}}\\!U\hat{\phi}_{2}\mathrm{d}z=\frac{2\pi
is}{\sqrt{s(s-4m^{2})}}$ $\displaystyle g_{22}$
$\displaystyle=\int_{\gamma_{2}}\\!U\hat{\phi}_{2}\mathrm{d}z=-2\pi i$ (32)

Writing the two integrals we are looking for as generic linear combinations

$\displaystyle
J_{1}=f_{11}I^{\text{sne}}_{11100}+f_{12}I^{\text{sne}}_{111\scalebox{0.5}[0.75]{$\,-$}10}\qquad\quad\text{and}\qquad\quad
J_{2}=f_{21}I^{\text{sne}}_{11100}+f_{22}I^{\text{sne}}_{111\scalebox{0.5}[0.75]{$\,-$}10}$
(33)

we obtain the period matrix $P_{ij}=f_{ik}g_{kj}$ as

$\displaystyle P$ $\displaystyle=\left[\begin{array}[]{cc}\frac{2\pi
i(f_{11}+sf_{12})}{\sqrt{s(s-4m^{2})}}&-2\pi if_{12}\\\ \frac{2\pi
i(f_{21}+sf_{22})}{\sqrt{s(s-4m^{2})}}&-2\pi if_{22}\end{array}\right]$ (36)

Solving for $P=2\pi iI$ gives

$\displaystyle f_{11}=\sqrt{s(s{-}4m^{2})}\,,\qquad f_{12}=0\,,\qquad
f_{21}=s\,,\qquad f_{22}=-1\,,$ (37)

corresponding to

$\displaystyle
J_{1}=\sqrt{s(s{-}4m^{2})}I^{\text{sne}}_{11100}\qquad\text{and}\qquad
J_{2}=sI^{\text{sne}}_{11100}-I^{\text{sne}}_{111\scalebox{0.5}[0.75]{$\,-$}10}$
(38)

With this we may compute the epsilon factorized differential equation matrix

$\displaystyle\frac{\mathrm{d}J_{i}}{\mathrm{d}s}=\epsilon
A_{ij}J_{j}\qquad\text{with}\qquad
A=\left[\begin{array}[]{cc}\frac{4(s{-}m^{2})}{(4m^{2}{-}s)s}&\frac{3}{\sqrt{s(s{-}4m^{2})}}\\\
\frac{-2}{\sqrt{s(s{-}4m^{2})}}&\frac{1}{s}\end{array}\right]$ (41)

written disregarding integrals in lower sectors (here only the double-tadpole
$I^{\text{sne}}_{11000}$) so we see that our method works in this case as
well. This success is what motivates us to try the same approach in elliptic
cases.

##  3 The proposed algorithm

Motivated by the examples in the previous section, we are now ready to
formulate our algorithm. Having a set of $\nu$ integrals of the form

$\displaystyle J_{i}$
$\displaystyle=\mathcal{K}\int_{\mathcal{C}}u\hat{\varphi}_{i}\mathrm{d}^{n}x$
(42)

where as discussed in section 1.1 $u$ is a multivalued function, and
$\hat{\varphi}_{i}$ a set of $\nu$ rational (in $x$) functions, our claim is
that if the integrand can be written as

$\displaystyle\mathcal{K}u\hat{\varphi}_{i}$
$\displaystyle=\sigma\hat{\Phi}_{i}$ (43)

where $\sigma$ is a pure function, and where the period matrix $P$ with
$P_{ij}=\int_{\gamma_{j}}\\!\\!\hat{\Phi}_{i}\mathrm{d}^{n}x$ is proportional
to the identity matrix as $P=(2\pi i)^{n}I$, then the integrals $J_{i}$ will
fulfill epsilon-factorized differential equations. In cases where
$\mathcal{K}$ is itself pure to begin with, it can be ignored in the above
discussion.

What we will do in practice is to write the basis integrals as a linear
combination of known intermediate basis integrals $I_{j}$ (often called
precanonical in the polylogarithmic case,) that is $J_{i}=f_{ij}I_{j}$ or
correspondingly as $\hat{\Phi}_{i}=f_{ij}\hat{\Phi}^{\text{int}}_{j}$, and
then imposing $P=(2\pi i)^{n}I$ will impose $\nu^{2}$ constraints on the
$f_{ij}$ fixing them all uniquely.

In the examples in the following section, we will write the
$\hat{\Phi}^{\text{int}}_{\;}$ as
$\hat{\Phi}^{\text{int}}_{i}=\hat{\phi}_{i}/Y$ where $Y$ is the square root of
a (monic) polynomial of degree four, corresponding to the elliptic curve
characterizing the problem.

##  4 Examples

In this section we will look at examples of elliptic Feynman integrals which
we can bring into a form that has epsilon-factorized differential equations
with the above algorithm. We will restrict the discussion to cases where the
ellipticity is present in the highest sector only, and where a maximal cut can
bring the integrals in that sector to a univariate form.

###  4.1 The same-mass elliptic sunrise

As a first example let us look at the most basic of elliptic Feynman
integrals, the same-mass elliptic sunrise (s1m). This is defined as in section
1.2, but with the restriction that the three masses are the same, i.e.
$m_{1}^{2}=m_{2}^{2}=m_{3}^{2}=m^{2}$.

This integral family has three master integrals. One is the double-tadpole
(e.g. $I^{\text{s1m}}_{11000}$) while the last two are in the highest sector
and the only ones we will care about in the following where we disregard
subsectors completely. Picking as intermediate basis integrals
$I^{\text{s1m}}_{11100}$ and $I^{\text{s1m}}_{21100}$ we can write the two
candidate integrals as

$\displaystyle J_{1}$
$\displaystyle=c(f_{11}I^{\text{s1m}}_{11100}+f_{12}I^{\text{s1m}}_{21100})\,+\,\text{lower}$
$\displaystyle J_{2}$
$\displaystyle=c(f_{21}I^{\text{s1m}}_{11100}+f_{22}I^{\text{s1m}}_{21100})\,+\,\text{lower}$
(44)

where the “lower” covers potential contributions from the double-tadpole
sector. The $c$ is a constant prefactor, the value of which will not change
the differential equations and which will be ignored in all of the following.

Our goal now is to find values of the $f_{ij}$ that makes the period matrix of
these two integrals the identity matrix.

Parametrizing this integral with the loop-by-loop Baikov parametrization, we
get from eq. (11)

$\displaystyle I^{\text{s1m}}_{a_{1}a_{2}a_{3}a_{4}0}\;$
$\displaystyle=\mathcal{K}\int_{\mathcal{C}}\mathrm{d}^{4}x\frac{x_{4}^{-a_{4}}\;u}{x_{1}^{a_{1}}x_{2}^{a_{2}}x_{3}^{a_{3}}}$
(45)

with

$\displaystyle u$
$\displaystyle=x_{4}^{\epsilon}\times\big{(}(x_{1}{-}x_{2})^{2}{-}2(2m^{2}{+}x_{1}{+}x_{2})x_{4}{+}x_{4}^{2}\big{)}^{-\frac{1}{2}-\epsilon}$
$\displaystyle\;\;\;\;\times\big{(}(m^{2}{-}s{+}x_{3})^{2}{-}2(m^{2}{+}s{+}x_{3})x_{4}{+}x_{4}^{2}\big{)}^{-\frac{1}{2}-\epsilon}$
(46)

as by eq. (12). $\mathcal{K}$ given by eq. (13) is a pure function, and $u$ is
pure multiplied with a factor of $Y^{-1}$ where on the maximal cut (i.e. of
$x_{1}$,$x_{2}$,$x_{3}$ and with $x_{4}\rightarrow z$)

$\displaystyle Y$
$\displaystyle=\sqrt{z\,(z{-}4m^{2})\,\big{(}z^{2}-2(m^{2}{+}s)z+(m^{2}{-}s)^{2}\big{)}}\,,$
(47)

Following the algorithm of sec. 3 we factor out the pure part of $u$, leaving
the expressions

$\displaystyle
I^{\text{s1m}}_{11100}\rightarrow\int_{\mathcal{C}}\frac{\hat{\phi}_{1}\mathrm{d}z}{Y}\qquad\qquad\qquad
I^{\text{s1m}}_{21100}\rightarrow\int_{\mathcal{C}}\frac{\hat{\phi}_{2}\mathrm{d}z}{Y}$
(48)

where the integrands are given by

$\displaystyle\hat{\phi}_{1}=1\qquad\text{and}\qquad\hat{\phi}_{2}=\frac{1{+}2\epsilon}{z{-}4m^{2}}$
(49)

of which the latter has been computed as
$\hat{\phi}_{2}=(\partial_{x_{1}}u)/u|_{\text{cut}}$.

Please note the presence of the $\epsilon$-dependence in $\hat{\phi}_{2}$.
This factor comes out of the algorithm in a natural way, and is essential in
order to get the correct epsilon factorization of the differential equation.
Such $\epsilon$-dependence would be missed by a purely two-dimensional
approach!

Figure 4: The two independent contours for the same-mass elliptic sunrise

$Y^{2}$ has four roots

$\displaystyle r_{\text{i}}$ $\displaystyle=0\;,$ $\displaystyle
r_{\text{ii}}$ $\displaystyle=(\sqrt{m^{2}}-\sqrt{s})^{2}\;,$ $\displaystyle
r_{\text{iii}}$ $\displaystyle=(\sqrt{m^{2}}+\sqrt{s})^{2}\;,$ $\displaystyle
r_{\text{iv}}$ $\displaystyle=4m^{2}\;,$ (50)

so we can look at the two basic cycles $\gamma_{1}=\mathcal{C}_{\text{ii-
iii}}$ and $\gamma_{2}=\mathcal{C}_{\text{i-ii}}$ which may be seen on fig. 4.
In the convergent case

$\displaystyle\int_{\mathcal{C}_{i\text{-}j}}$
$\displaystyle=\;\;2\int_{r_{i}}^{r_{j}}$ (51)

allowing us to compute the integrals

$\displaystyle g_{11}$
$\displaystyle=\int_{\gamma_{1}}\frac{\hat{\phi}_{1}\mathrm{d}z}{Y}\\!\\!$
$\displaystyle=\,\frac{4K(k^{2})}{R}$ $\displaystyle g_{21}$
$\displaystyle=\int_{\gamma_{1}}\frac{\hat{\phi}_{2}\mathrm{d}z}{Y}\\!\\!$
$\displaystyle=\,\frac{(1{+}2\epsilon)}{m^{2}R}\left(K(k^{2})+\frac{(\sqrt{m^{2}}{+}\sqrt{s})^{2}}{4m^{2}-(\sqrt{m^{2}}{+}\sqrt{s})^{2}}E(k^{2})\right)$
$\displaystyle g_{12}$
$\displaystyle=\int_{\gamma_{2}}\frac{\hat{\phi}_{1}\mathrm{d}z}{Y}\\!\\!$
$\displaystyle=\,\frac{-4iK(1{-}k^{2})}{R}$ (52) $\displaystyle g_{22}$
$\displaystyle=\int_{\gamma_{2}}\frac{\hat{\phi}_{2}\mathrm{d}z}{Y}\\!\\!$
$\displaystyle=\,\frac{-2i(1{+}2\epsilon)}{R\big{(}4m^{2}{-}(\sqrt{m^{2}}{+}\sqrt{s})^{2}\big{)}}\left(2K(1{-}k^{2})-\frac{(\sqrt{m^{2}}{+}\sqrt{s})^{2}}{2m^{2}}E(1{-}k^{2})\right)$

with

$\displaystyle R$
$\displaystyle=\sqrt{(3\sqrt{m^{2}}{-}\sqrt{s})(\sqrt{m^{2}}{+}\sqrt{s})^{3}}$
(53)

and

$\displaystyle k^{2}$
$\displaystyle:=\frac{16\sqrt{(m^{2})^{3}s}}{(3\sqrt{m^{2}}{-}\sqrt{s})(\sqrt{m^{2}}{+}\sqrt{s})^{3}}$
$\displaystyle 1{-}k^{2}$
$\displaystyle=\frac{(3\sqrt{m^{2}}{+}\sqrt{s})(\sqrt{m^{2}}{-}\sqrt{s})^{3}}{(3\sqrt{m^{2}}{-}\sqrt{s})(\sqrt{m^{2}}{+}\sqrt{s})^{3}}$
(54)

The period matrix ($P_{ij}=\int_{\gamma_{j}}\phi_{i}/Y$) of the two integrals
of eqs. (44) is then

$\displaystyle P$
$\displaystyle=\left[\begin{array}[]{cc}f_{11}g_{11}+f_{12}g_{21}&f_{11}g_{12}+f_{12}g_{22}\\\\[5.69054pt]
f_{21}g_{11}+f_{22}g_{21}&f_{21}g_{12}+f_{22}g_{22}\end{array}\right]$ (57)

Imposing $M=2\pi iI$ gives a unique solution for the $f_{ij}$ as

$\displaystyle f_{11}$
$\displaystyle=iR\,E(1{-}k^{2})-\frac{4im^{2}R}{(\sqrt{m^{2}}{+}\sqrt{s})^{2}}K(1{-}k^{2})$
$\displaystyle f_{12}$
$\displaystyle=\frac{-4im^{2}(\sqrt{m^{2}}{-}\sqrt{s})(3\sqrt{m^{2}}{+}\sqrt{s})R}{(1{+}2\epsilon)(\sqrt{m^{2}}{+}\sqrt{s})^{2}}K(1{-}k^{2})$
$\displaystyle f_{21}$
$\displaystyle=-R\,E(k^{2})-\frac{(\sqrt{m^{2}}{-}\sqrt{s})(3\sqrt{m^{2}}{+}\sqrt{s})R}{(\sqrt{m^{2}}{+}\sqrt{s})^{2}}K(k^{2})$
(58) $\displaystyle f_{22}$
$\displaystyle=\frac{-4m^{2}(\sqrt{m^{2}}{-}\sqrt{s})(3\sqrt{m^{2}}{+}\sqrt{s})R}{(1{+}2\epsilon)(\sqrt{m^{2}}{+}\sqrt{s})^{2}}K(k^{2})$

With this choice we can compute the system of differential equations for the
two integrals. It is epsilon factorized

$\displaystyle\frac{\mathrm{d}J_{i}}{\mathrm{d}s}=\epsilon
A_{ij}J_{j}\;+\;\text{lower}$ (59)

with

$\displaystyle A_{11}$
$\displaystyle=\frac{-1}{2s}-\frac{(3m^{2}{+}s)^{2}K(k^{2})K(1{-}k^{2})}{\pi
s(3\sqrt{m^{2}}{-}\sqrt{s})(\sqrt{m^{2}}{+}\sqrt{s})^{3}}+\frac{3(\sqrt{m^{2}}{+}\sqrt{s})^{2}\,E(k^{2})E(1{-}k^{2})}{\pi
s(\sqrt{m^{2}}{-}\sqrt{s})(3\sqrt{m^{2}}{+}\sqrt{s})}$
$\displaystyle\;\;\;+\frac{(27m^{2}-24m^{3}\sqrt{s}-18m^{2}s-s^{2})E(k^{2})K(1{-}k^{2})}{\pi
s(m^{2}{-}s)(9m^{2}{-}s)}$ $\displaystyle A_{12}$
$\displaystyle=\frac{-16im^{2}(3m^{3}{-}3m^{2}\sqrt{s}{-}3ms{-}\sqrt{s}^{3})K(1{-}k^{2})^{2}}{\pi\sqrt{s}(m^{2}{-}s)(3m^{2}{-}s)(\sqrt{m^{2}}{-}\sqrt{s})^{2}}+\frac{3i(\sqrt{m^{2}}{+}\sqrt{s})^{2}E(1{-}k^{2})^{2}}{\pi
s(\sqrt{m^{2}}{-}\sqrt{s})(3\sqrt{m^{2}}{+}\sqrt{s})}$
$\displaystyle\;\;\;+\frac{2i(3m^{3}-9m^{2}\sqrt{s}-2ms+\sqrt{s}^{3})K(1{-}k^{2})E(1{-}k^{2})}{\pi
s(m^{2}{-}s)(3\sqrt{m^{2}}{-}\sqrt{s})}$ $\displaystyle A_{21}$
$\displaystyle=\frac{i(3m^{2}{+}s)^{2}K(k)^{2}}{\pi
s(3\sqrt{m^{2}}{-}\sqrt{s})(\sqrt{m^{2}}{+}\sqrt{s})^{3}}+\frac{3i(\sqrt{m^{2}}{+}\sqrt{s})^{2}E(k)^{2}}{\pi
s(3\sqrt{m^{2}}{+}\sqrt{s})(\sqrt{m^{2}}{-}\sqrt{s})}$
$\displaystyle\;\;\;-\frac{4i(3m^{2}{-}s)(3m^{2}{+}s)K(k^{2})E(k^{2})}{\pi
s(m^{2}{-}s)(9m^{2}{-}s)}$ (60) $\displaystyle A_{22}$
$\displaystyle=\frac{2(2\sqrt{m^{2}}+\sqrt{s})}{\sqrt{s}(\sqrt{m^{2}}{-}\sqrt{s})(3\sqrt{m^{2}}{+}\sqrt{s})}-\frac{16m^{2}(3m^{3}{-}3m^{2}\sqrt{s}{-}3ms{-}\sqrt{s}^{3})K(k^{2})K(1{-}k^{2})}{\pi\sqrt{s}(m^{2}{-}s)(9m^{2}{-}s)(\sqrt{m^{2}}{+}\sqrt{s})^{2}}$
$\displaystyle\;\;\;-\frac{3(\sqrt{m^{2}}{+}\sqrt{s})^{2}E(k^{2})E(1{-}k^{2})}{\pi
s(\sqrt{m^{2}}{-}\sqrt{s})(3\sqrt{m^{2}}{+}\sqrt{s})}+\frac{(27m^{4}{-}24m^{3}\sqrt{s}{-}18m^{2}s{-}s^{2})K(k^{2})E(1{-}k^{2})}{\pi
s(m^{2}{-}s)(9m^{2}{-}s)}$

and a similar epsilon-factorized differential equation may be found in the
other variable $\mathrm{d}J/\mathrm{d}m^{2}$. This shows that our algorithm
works for elliptic examples too.

###  4.2 The elliptic nonplanar double triangle

The elliptic nonplanar double triangle (npt) Czakon:2008ii ;
vonManteuffel:2017hms , is defined by the set of propagators

$\displaystyle D_{1}$ $\displaystyle=k_{1}^{2}{-}m^{2},$ $\displaystyle D_{2}$
$\displaystyle=(k_{1}{+}p_{1})^{2}{-}m^{2},$ $\displaystyle D_{3}$
$\displaystyle=(k_{1}{-}k_{2}{-}p_{2})^{2}{-}m^{2},$ $\displaystyle D_{4}$
$\displaystyle=(k_{1}{-}k_{2})^{2}{-}m^{2},$ $\displaystyle D_{5}$
$\displaystyle=k_{2}^{2},$ $\displaystyle D_{6}$
$\displaystyle=(k_{2}{+}p_{1}{+}p_{2})^{2},$ $\displaystyle D_{7}$
$\displaystyle=(k_{2}{+}p_{2})^{2},$ (61)

where only the first six may appear as actual propagators. The kinematics is
such that $p_{1}^{2}=p_{2}^{2}=0$ and $(p_{1}{+}p_{2})^{2}=s$.

Figure 5: The elliptic nonplanar double triangle

This sector contains 11 master integrals. The first nine are in subsectors,
and can be put in canonical form in the traditional sense. The last two are in
the highest, elliptic sector, and it is those that will concern us here.
Choosing as intermediate basis integrals the set $I^{\text{npt}}_{1111110}$
and $I^{\text{npt}}_{2111110}$ we will write the two elliptic master integrals
as linear combination of the intermediate basis:

$\displaystyle J_{1}$
$\displaystyle=c(f_{11}I^{\text{npt}}_{1111110}+f_{12}I^{\text{npt}}_{2111110})\,+\,\text{lower}$
$\displaystyle J_{2}$
$\displaystyle=c(f_{21}I^{\text{npt}}_{1111110}+f_{22}I^{\text{npt}}_{2111110})\,+\,\text{lower}$
(62)

Using the Baikov parametrization, the integrals in this sector may on the
maximal cut (of $x_{1}$-$x_{6}$ with $z=x_{7}$) be written as

$\displaystyle I^{\text{npt}}\;$
$\displaystyle=\mathcal{K}\int_{\mathcal{C}}\\!u\hat{\phi}\,\mathrm{d}z$ (63)

Around $d=4$, we have once again that $\mathcal{K}$ is pure, and $u$ is pure
times a factor of $Y^{-1}$ with

$\displaystyle Y$ $\displaystyle=\sqrt{z\,(z{+}s)\,(z^{2}{+}sz{-}4m^{2}s)}\,,$
(64)

where a factor of $s$ has been absorbed in the $\hat{\phi}$. This justifies
studying the integrals with the $u\rightarrow Y^{-1}$ limit taken, and in that
case we may write the intermediate basis integrals as

$\displaystyle
I^{\text{npt}}_{1111110}\rightarrow\int_{\mathcal{C}}\frac{\hat{\phi}_{1}\mathrm{d}z}{Y}\qquad\qquad\qquad
I^{\text{npt}}_{2111110}\rightarrow\int_{\mathcal{C}}\frac{\hat{\phi}_{2}\mathrm{d}z}{Y}$
(65)

where

$\displaystyle\hat{\phi}_{1}=\frac{1}{s}\qquad\text{and}\qquad\hat{\phi}_{2}=\frac{(1{+}2\epsilon)(z{+}s)}{s(z^{2}{+}sz{-}4m^{2}s)}$
(66)

The four roots of $Y^{2}$ are given as

$\displaystyle
r_{\text{i}}=-\tfrac{1}{2}\sqrt{s}(\sqrt{s}{+}\sqrt{16m^{2}{+}s})\,,\qquad
r_{\text{ii}}=-s\,,\qquad r_{\text{iii}}=0\,,\qquad
r_{\text{iv}}=-\tfrac{1}{2}\sqrt{s}(\sqrt{s}{-}\sqrt{16m^{2}{+}s})\,,$ (67)

and we see that a set of independent integration contours can be chosen as

$\displaystyle\gamma_{1}$ $\displaystyle=\mathcal{C}_{\text{ii-iii}}$
$\displaystyle\gamma_{2}$ $\displaystyle=\mathcal{C}_{\text{i-ii}}$ (68)

similarly to the previous example where the contours are depicted on fig. 4.

We note that where in the previous section we always had
$\int_{\mathcal{C}_{\alpha\text{-}\beta}}\phi=2\int_{\alpha}^{\beta}\phi$,
that is not the case here for $\int_{\mathcal{C}_{\text{i-ii}}}\phi_{2}/Y$ due
to a divergence in $z=r_{\text{i}}$. The results for the integrals are

$\displaystyle g_{11}$
$\displaystyle=\int_{\gamma_{1}}\\!\\!\frac{\hat{\phi}_{1}\mathrm{d}z}{Y}\\!\\!\\!$
$\displaystyle=\;\;\frac{8K(k^{2})}{s^{3/2}(\sqrt{16m^{2}{+}s}{+}\sqrt{s})}$
$\displaystyle g_{21}$
$\displaystyle=\int_{\gamma_{1}}\\!\\!\frac{\hat{\phi}_{2}\mathrm{d}z}{Y}\\!\\!\\!$
$\displaystyle=\;\;\frac{-8(1{+}2\epsilon)}{s^{3/2}(16m^{2}{+}s)(\sqrt{16m^{2}{+}s}{+}\sqrt{s})}\left(K(k^{2})+\frac{\sqrt{16m^{2}{+}s}{+}\sqrt{s}}{\sqrt{16m^{2}{+}s}{-}\sqrt{s}}\,E(k^{2})\right)$
$\displaystyle g_{12}$
$\displaystyle=\int_{\gamma_{2}}\\!\\!\frac{\hat{\phi}_{1}\mathrm{d}z}{Y}\\!\\!\\!$
$\displaystyle=\;\;\frac{-8iK(1{-}k^{2})}{s^{3/2}(\sqrt{16m^{2}{+}s}{+}\sqrt{s})}$
(69) $\displaystyle g_{22}$
$\displaystyle=\int_{\gamma_{2}}\\!\\!\frac{\hat{\phi}_{2}\mathrm{d}z}{Y}\\!\\!\\!$
$\displaystyle=\;\;\frac{i(1{+}2\epsilon)}{s^{3/2}m^{2}}\left(\frac{K(1{-}k^{2})}{\sqrt{16m^{2}{+}s}}-\frac{\sqrt{16m^{2}{+}s}{+}\sqrt{s}}{2(16m^{2}{+}s)}E(1{-}k^{2})\right)$

where we have

$\displaystyle k^{2}$
$\displaystyle:=\frac{4\sqrt{s}\sqrt{16m^{2}{+}s}}{(\sqrt{16m^{2}{+}s}{+}\sqrt{s})^{2}}\,,$
$\displaystyle 1{-}k^{2}$
$\displaystyle=\frac{(\sqrt{16m^{2}{+}s}{-}\sqrt{s})^{2}}{(\sqrt{16m^{2}{+}s}{+}\sqrt{s})^{2}}\,.$
(70)

We may now compute the period matrix $P_{ij}=f_{ik}g_{kj}$, and imposing
$P=2\pi iI$ allows us to uniquely fix the coefficients $f_{ik}$ of eq. (62) to

$\displaystyle f_{11}$
$\displaystyle=\tfrac{1}{2}is^{3/2}(\sqrt{16m^{2}{+}s}{+}\sqrt{s})E(1{-}k^{2})-is^{3/2}\sqrt{16m^{2}{+}s}K(1{-}k^{2})$
$\displaystyle f_{12}$
$\displaystyle=\frac{-is^{3/2}(16m^{2}{+}s)(\sqrt{16m^{2}{+}s}{-}\sqrt{s})K(1{-}k^{2})}{2(1{+}2\epsilon)}$
$\displaystyle f_{21}$
$\displaystyle=-\tfrac{1}{2}s^{3/2}(\sqrt{16m^{2}{+}s}{+}\sqrt{s})E(k^{2})-\tfrac{1}{2}s^{3/2}(\sqrt{16m^{2}{+}s}{-}\sqrt{s})K(k^{2})$
(71) $\displaystyle f_{22}$
$\displaystyle=\frac{-s^{3/2}(16m^{2}{+}s)(\sqrt{16m^{2}{+}s}{-}\sqrt{s})K(k^{2})}{2(1{+}2\epsilon)}$

With this we may compute the differential equation, which once again is
epsilon factorized and it is given by

$\displaystyle\frac{\mathrm{d}J_{i}}{\mathrm{d}s}=\epsilon
A_{ij}J_{j}\;+\;\text{lower}$ (72)

with

$\displaystyle A_{11}$
$\displaystyle=\frac{8(12m^{2}{+}s)K(k^{2})K(1{-}k^{2})}{\pi\sqrt{s}(16m^{2}{+}s)(\sqrt{16m^{2}{+}s}{+}\sqrt{s})}+\frac{2}{\pi
s}\left(1{-}\frac{8m^{2}}{16m^{2}{+}s}{+}\frac{\sqrt{s}}{\sqrt{16m^{2}{+}s}}\right)E(k^{2})E(1{-}k^{2})$
$\displaystyle\;\;\;+\frac{-4(12m^{2}{+}s)K(k^{2})E(1{-}k^{2})}{\pi
s(16m^{2}{+}s)}+\frac{-2(\sqrt{16m^{2}{+}s}{+}\sqrt{s})E(k^{2})K(1{-}k^{2})}{\pi\sqrt{s}(16m^{2}{+}s)}$
$\displaystyle A_{12}$
$\displaystyle=\frac{-64im^{2}K(1{-}k^{2})^{2}}{\pi\sqrt{s}\sqrt{16m^{2}{+}s}(\sqrt{16m^{2}{+}s}+\sqrt{s})}+\frac{i(\sqrt{16m^{2}{+}s}+\sqrt{s})^{2}E(1{-}k^{2})^{2}}{\pi
s(16m^{2}{+}s)}$ $\displaystyle\;\;\;+\frac{4i}{\pi
s}\left(\frac{\sqrt{s}}{\sqrt{16m^{2}{+}s}}-\frac{8m^{2}}{16m^{2}{+}s}\right)K(1{-}k^{2})E(1{-}k^{2})$
$\displaystyle A_{21}$
$\displaystyle=\frac{i(12m^{2}{+}s)(\sqrt{16m^{2}{+}s}-\sqrt{s})^{2}K(k^{2})^{2}}{4m^{2}\pi
s(16m^{2}{+}s)}+\frac{i(\sqrt{16m^{2}{+}s}-\sqrt{s})^{2}E(k^{2})^{2}}{\pi
s(16m^{2}{+}s)}$ $\displaystyle\;\;\;+\frac{-4iK(k^{2})E(k^{2})}{\pi s}$ (73)
$\displaystyle A_{22}$
$\displaystyle=\frac{2(\sqrt{16m^{2}{+}s}-\sqrt{s})^{2}K(k^{2})K(1{-}k^{2})}{\pi
s\sqrt{16m^{2}{+}s}(\sqrt{16m^{2}{+}s}+\sqrt{s})}+\frac{-2}{s\pi}\left(1{-}\frac{8m^{2}}{16m^{2}{+}s}{-}\frac{\sqrt{s}}{\sqrt{16m^{2}{+}s}}\right)E(k^{2})E(1{-}k^{2})$
$\displaystyle\;\;\;+\frac{16m^{2}K(k^{2})E(1{-}k^{2})}{\pi
s(16m^{2}{+}s)}+\frac{-32m^{2}E(k^{2})K(1{-}k^{2})}{\pi
s\sqrt{16m^{2}{+}s}(\sqrt{16m^{2}{+}s}+\sqrt{s})}$

with a similar and likewise epsilon factorized equation to be found for the
other derivative $\mathrm{d}J/\mathrm{d}m^{2}$.

###  4.3 The two-mass elliptic sunrise

Next let us look at the two-mass elliptic sunrise (s2m). This not very well
studied integral is defined as is section 1.2 with
$m_{1}^{2}=m_{2}^{2}=m_{b}^{2}$ and $m_{3}^{2}=m_{a}^{2}$, making it
intermediate between the same-mass elliptic sunrise of section 4.1 and the
three-mass elliptic sunrise of section 4.4 both in terms of the number of
scales and the number of master integrals. The integral family contains five
master integrals; two double-tadpoles and three in the highest elliptic sector
which is what will concern us here.

On the maximal cut in $d=2$ we get that the integrals in the highest sector
may be written as

$\displaystyle
J_{i}=\int_{\mathcal{C}}\frac{\hat{\phi}\,\mathrm{d}z}{Y}\qquad\text{with}\qquad
Y=\sqrt{z\,(z{-}4m_{b}^{2})\,\big{(}z^{2}{-}2(m^{2}_{a}{+}s)z{+}(m^{2}_{a}{-}s)^{2}\big{)}}$
(74)

We will pick the three intermediate basis integrals $I^{\text{s2m}}_{11100}$,
$I^{\text{s2m}}_{21100}$, and
$I^{\text{s2m}}_{111\scalebox{0.5}[0.75]{$\,-$}10}$, which correspond to the
integrands

$\displaystyle\hat{\phi}_{1}=1\,,\qquad\hat{\phi}_{2}=\frac{1{+}2\epsilon}{z{-}4m_{b}^{2}}\,,\qquad\hat{\phi}_{3}=z\,.$
(75)

Once again we will define the master integrals we are looking for, as generic
linear combinations of the intermediate basis integrals

$\displaystyle
J_{i}=f_{i1}I^{\text{s2m}}_{11100}+f_{i2}I^{\text{s2m}}_{21100}+f_{i3}I^{\text{s2m}}_{111\scalebox{0.5}[0.75]{$\,-$}10}$
(76)

$Y^{2}$ has the four roots

$\displaystyle r_{\text{i}}$ $\displaystyle=0\,,$ $\displaystyle
r_{\text{ii}}$ $\displaystyle=(\sqrt{m_{a}^{2}}-\sqrt{s})^{2}\,,$
$\displaystyle r_{\text{iii}}$
$\displaystyle=(\sqrt{m_{a}^{2}}+\sqrt{s})^{2}\,,$ $\displaystyle
r_{\text{iv}}$ $\displaystyle=4m_{b}^{2}\,,$ (77)

and with this we may define our three master contours

$\displaystyle\gamma_{1}=\mathcal{C}_{\text{ii-
iii}}\,,\qquad\gamma_{2}=\mathcal{C}_{\text{i-ii}}\,,\qquad\gamma_{3}=\mathcal{C}_{\infty}\,.$
(78)

as depicted on fig. 6. We will also define the abbreviations

$\displaystyle\delta_{\pm}$ $\displaystyle=\sqrt{m_{a}^{2}}\pm\sqrt{s}\,,$
$\displaystyle\mu_{\pm}$
$\displaystyle=4m_{b}^{2}{-}(\sqrt{m_{a}^{2}}{\pm}\sqrt{s})^{2}\,,$
$\displaystyle\lambda_{\pm}$
$\displaystyle=m_{b}^{2}{\pm}\sqrt{m_{a}^{2}}\sqrt{s}\,,$ $\displaystyle\rho$
$\displaystyle=2m_{b}^{2}{-}m_{a}^{2}{-}s\,,$ (79)

and

$\displaystyle k^{2}$
$\displaystyle=\frac{16m_{b}^{2}\sqrt{m_{a}^{2}}\sqrt{s}}{\delta_{+}^{2}\,\mu_{-}}\,,$
$\displaystyle 1{-}k^{2}$
$\displaystyle=\frac{\delta_{-}^{2}\,\mu_{+}}{\delta_{+}^{2}\,\mu_{-}}\,,$
$\displaystyle n^{2}$ $\displaystyle=\frac{4\lambda_{-}^{2}}{\rho^{2}}\,,$
$\displaystyle n^{\prime 2}$
$\displaystyle=\frac{\lambda_{-}^{2}}{\lambda_{+}^{2}}\,,$ (80)

where we note the relation $(1{-}k^{2})/n^{\prime 2}=1-k^{2}/n^{2}$. With this
we can compute the needed integrals, as
$g_{ij}=\int_{\gamma_{j}}\frac{\hat{\phi}_{i}\mathrm{d}z}{Y}$:

$\displaystyle g_{11}$
$\displaystyle=\frac{4K(k^{2})}{\delta_{+}\,\sqrt{\mu_{-}}}$ $\displaystyle
g_{21}$
$\displaystyle=\frac{-(1{+}2\epsilon)}{m_{b}^{2}\sqrt{\mu_{-}}}\left(\frac{K(k^{2})}{\delta_{+}}+\frac{\delta_{+}E(k^{2})}{\mu_{+}}\right)$
$\displaystyle g_{12}$
$\displaystyle=\frac{-4iK(1{-}k^{2})}{\delta_{+}\,\sqrt{\mu_{-}}}\qquad\qquad$
$\displaystyle g_{22}$
$\displaystyle=\frac{-i(1{+}2\epsilon)}{\mu_{+}\sqrt{\mu_{-}}}\left(\frac{-4K(1{-}k^{2})}{\delta_{+}}+\frac{\delta_{+}E(1{-}k^{2})}{m_{b}^{2}}\right)$
$\displaystyle g_{31}$
$\displaystyle=\frac{4\delta_{-}^{2}}{\delta_{+}\,\lambda_{-}\,\sqrt{\mu_{-}}}\left(m_{b}^{2}K(k^{2})-\frac{\lambda_{+}\mu_{+}\Pi(n^{2},k^{2})}{2\rho}\right)\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$
$\displaystyle g_{32}$
$\displaystyle=i\pi+\frac{-i4m_{b}^{2}}{\delta_{+}\lambda_{-}\sqrt{\mu_{-}}}\left(\delta_{-}^{2}K(1{-}k^{2})+\frac{2\sqrt{m_{a}^{2}}\,\rho\,\sqrt{s}\,\Pi(n^{\prime
2},1{-}k^{2})}{\lambda_{+}}\right)\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$
(81) $\displaystyle g_{13}$ $\displaystyle=0$ $\displaystyle g_{23}$
$\displaystyle=0\qquad\qquad\qquad\qquad\qquad\quad g_{33}=-2\pi i$

making this the first appearance of the complete elliptic integral of the
third kind $\Pi(n^{2},k^{2})$.

Figure 6: The three independent contours for the two-mass elliptic sunrise

With this we may write down the period matrix $P_{ij}=f_{ik}g_{kj}$ and impose
it to be diagonal $P=2\pi iI$ by solving for the $f_{ik}$. We get

$\displaystyle f_{11}$
$\displaystyle=\frac{-4im_{b}^{2}\sqrt{\mu_{-}}K(1{-}k^{2})}{\delta_{+}}+i\delta_{+}\sqrt{\mu_{-}}E(1{-}k^{2})$
$\displaystyle f_{21}$
$\displaystyle=\frac{-\sqrt{\mu_{-}}\mu_{+}K(k^{2})}{\delta_{+}}-\delta_{+}\sqrt{\mu_{-}}E(k^{2})$
$\displaystyle f_{31}$
$\displaystyle=\frac{\delta_{-}^{2}\lambda_{+}\mu_{p}\Pi(n^{2},k^{2})}{\lambda_{-}\rho\pi}\left(\frac{4\mu_{a}^{2}K(1{-}k^{2})}{\delta_{-}^{2}}-E(1{-}k^{2})\right)-\frac{\sqrt{\mu_{-}}\mu_{+}K(k^{2})}{2\delta_{+}}+\frac{\delta_{-}^{2}m_{b}^{2}}{\lambda_{-}}$
$\displaystyle\;\;\;+\frac{4m_{b}^{2}\sqrt{m_{a}^{2}}\sqrt{s}\rho\Pi(n^{\prime
2},1{-}k^{2})}{\lambda_{-}\lambda_{+}\pi}\left(\frac{\mu_{+}K(k^{2})}{\delta_{+}^{2}}K(k^{2})+E(k^{2})\right)-\tfrac{1}{2}\delta_{+}\sqrt{\mu_{-}}E(k^{2})$
$\displaystyle f_{12}$
$\displaystyle=\frac{-4im_{b}^{2}\sqrt{\mu_{-}}\mu_{p}K(1{-}k^{2})}{\delta_{+}(1{+}2\epsilon)}\qquad\qquad\qquad\qquad
f_{22}=\frac{-4m_{b}^{2}\sqrt{\mu_{-}}\mu_{p}K(k^{2})}{\delta_{+}(1{+}2\epsilon)}$
(82) $\displaystyle f_{32}$
$\displaystyle=\frac{2m_{b}^{2}\mu_{+}}{\delta_{+}(1{+}2\epsilon)}\bigg{(}\frac{2\delta_{-}^{2}\lambda_{+}\mu_{+}\Pi(n^{2},k^{2})K(1{-}k^{2})}{\delta_{+}\lambda_{-}\rho\pi}+\frac{8m_{b}^{2}\rho\sqrt{m_{a}^{2}s}\,\Pi(n^{\prime
2},1{-}k^{2})K(k^{2})}{\delta_{+}\lambda_{-}\lambda_{+}\pi}-\sqrt{\mu_{-}}K(k^{2})\bigg{)}$
$\displaystyle f_{13}$ $\displaystyle=0\qquad\qquad\qquad\qquad\qquad
f_{23}=0\qquad\qquad\qquad\qquad\qquad f_{33}=-1$

Once again the differential equation matrices are epsilon factorized
$\mathrm{d}J_{i}/\mathrm{d}x=\epsilon A^{(x)}_{ij}J_{j}$. The entries are in
general too large to be written here, one example entry is

$\displaystyle A^{(m_{b}^{2})}_{11}$
$\displaystyle=\frac{-12\delta_{-}^{2}\lambda_{+}\mu_{+}\Pi(n^{2},k^{2})K(1{-}k^{2})}{\delta_{+}^{2}\lambda_{-}\mu_{-}\rho\pi}+\frac{24\delta_{-}\texttwosuperior
m_{b}^{2}K(k^{2})K(1{-}k^{2})}{\delta_{+}^{2}\lambda_{-}\mu_{-}\pi}-\frac{2\delta_{+}^{2}E(k^{2})E(1{-}k^{2})}{m_{b}^{2}\mu_{+}\pi}$
$\displaystyle\quad-\frac{8(2\mu_{+}{+}\mu_{-})E(k^{2})K(1{-}k^{2})}{\mu_{-}\mu_{+}\pi}-\frac{1}{m_{b}^{2}}$
(83)

and for the full expressions see the added file.

###  4.4 The three-mass elliptic sunrise

Next we will look at the generic or three-mass elliptic sunrise (s3m), exactly
as it is defined in section 1.2. The integral family contains seven master
integrals; three double-tadpoles and four in the highest, elliptic sector.

Integrals in that highest sector may, on the maximal cut and in two
dimensions, be written as

$\displaystyle J_{i}=\int_{\mathcal{C}}\frac{\hat{\phi}\,\mathrm{d}z}{Y}$ (84)

with

$\displaystyle
Y=\sqrt{\big{(}z^{2}{-}2(m^{2}_{1}{+}m^{2}_{2})z{+}(m^{2}_{1}{-}m^{2}_{2})^{2}\big{)}\big{(}z^{2}{-}2(m^{2}_{3}{+}s)z{+}(m^{2}_{3}{-}s)^{2}\big{)}}$
(85)

We will pick our master integrals as linear combinations of now four
intermediate basis integrals, as

$\displaystyle
J_{i}=f_{i1}I^{\text{s3m}}_{11100}+f_{i2}I^{\text{s3m}}_{21100}+f_{i3}I^{\text{s3m}}_{111\scalebox{0.5}[0.75]{$\,-$}10}+f_{i4}I^{\text{s3m}}_{1110\scalebox{0.5}[0.75]{$\,-$}1}$
(86)

These four intermediate basis integrals correspond to the integrands

$\displaystyle\hat{\phi}_{1}=1\,,\qquad\hat{\phi}_{2}=\frac{(1{+}2\epsilon)(z{+}m_{2}^{2}{-}m_{1}^{2})}{z^{2}{-}2(m^{2}_{1}{+}m^{2}_{2})z{+}(m^{2}_{1}{-}m^{2}_{2})^{2}}\,,\qquad\hat{\phi}_{3}=z\,,$
$\displaystyle\qquad\;\hat{\phi}_{4}=\frac{(m_{1}^{2}{-}m_{2}^{2})(m_{3}^{2}{-}s)}{2z}+\tfrac{1}{2}(m_{1}^{2}{+}m_{2}^{2}{+}m_{3}^{2}{+}s)-\tfrac{1}{2}z$
(87)

where the latter, which may appear incompatible with the loop-by-loop
parametrization used, has been obtained using the Passarino-Veltman inspired
approach described in section 11.1 of ref. Frellesvig:2019kgj .

Figure 7: The four independent contours for the three-mass elliptic sunrise

These four integrands are linear combination of the set

$\displaystyle\hat{\xi}_{1}=1\,,\qquad\hat{\xi}_{2}=\frac{(z{+}m_{2}^{2}{-}m_{1}^{2})}{z^{2}{-}2(m^{2}_{1}{+}m^{2}_{2})z{+}(m^{2}_{1}{-}m^{2}_{2})^{2}}\,,\qquad\hat{\xi}_{3}=z\,,\qquad\hat{\xi}_{4}=\frac{1}{z}\,,$
(88)

of which we will compute the integrals.

The master integration contours will be chosen as

$\displaystyle\gamma_{1}=\mathcal{C}_{\text{ii-
iii}}\,,\qquad\quad\gamma_{2}=\mathcal{C}_{\text{i-ii}}\,,\qquad\quad\gamma_{3}=\mathcal{C}_{\infty}\,,\qquad\quad\gamma_{4}=\mathcal{C}_{0}\,,$
(89)

as depicted on fig. 7, and the four roots of $Y^{2}$ are

$\displaystyle r_{\text{i}}$ $\displaystyle=\big{(}m_{2}{-}m_{1})^{2}\,,$
$\displaystyle r_{\text{ii}}$ $\displaystyle=(m_{3}{-}\sqrt{s})^{2}\,,$
$\displaystyle r_{\text{iii}}$ $\displaystyle=(m_{3}{+}\sqrt{s})^{2}\,,$
$\displaystyle r_{\text{iv}}$ $\displaystyle=(m_{2}{+}m_{1})^{2}\,.$ (90)

We will now introduce the abbreviations

$\displaystyle\delta_{01}$
$\displaystyle=-m_{1}{+}m_{2}{+}m_{3}{-}\sqrt{s}\;\;$
$\displaystyle\delta_{02}$
$\displaystyle=m_{1}{-}m_{2}{+}m_{3}{-}\sqrt{s}\;\;$
$\displaystyle\delta_{03}$ $\displaystyle=m_{1}{+}m_{2}{-}m_{3}{-}\sqrt{s}\;$
$\displaystyle\delta_{1}$ $\displaystyle=-m_{1}{+}m_{2}{+}m_{3}{+}\sqrt{s}$
$\displaystyle\delta_{2}$ $\displaystyle=m_{1}{-}m_{2}{+}m_{3}{+}\sqrt{s}$
$\displaystyle\delta_{3}$ $\displaystyle=m_{1}{+}m_{2}{-}m_{3}{+}\sqrt{s}$
$\displaystyle\delta_{+}$ $\displaystyle=m_{1}{+}m_{2}{+}m_{3}{+}\sqrt{s}$
$\displaystyle\delta_{0}$ $\displaystyle=m_{1}{+}m_{2}{+}m_{3}{-}\sqrt{s}$
$\displaystyle\beta$ $\displaystyle=m_{1}m_{2}m_{3}\sqrt{s}$
$\displaystyle\lambda_{1\pm}$ $\displaystyle=m_{2}m_{3}\pm m_{1}\sqrt{s}$
$\displaystyle\lambda_{2\pm}$ $\displaystyle=m_{1}m_{3}\pm m_{2}\sqrt{s}$
$\displaystyle\lambda_{3\pm}$ $\displaystyle=m_{1}m_{2}\pm m_{3}\sqrt{s}$ (91)
$\displaystyle\mu_{\pm}$
$\displaystyle=m_{1}m_{2}(m_{3}^{2}{-}s)^{2}{\pm}(m_{2}^{2}{-}m_{1}^{2})^{2}m_{3}\sqrt{s}\;\;\;\;\;\;\rho=(m_{1}^{2}{+}m_{2}^{2})(m_{3}^{2}{-}s)^{2}{-}(m_{2}^{2}{-}m_{1}^{2})^{2}(m_{3}^{2}{+}s)\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$
$\displaystyle q$
$\displaystyle=(m_{3}^{2}{-}m_{2}^{2})^{2}{+}2(m_{3}^{2}{+}m_{2}^{2})(m_{1}^{2}{-}s){+}s^{2}{+}2m_{1}^{2}s{-}3m_{1}^{4}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\
$ $\displaystyle\psi_{\pm}$
$\displaystyle=(m_{1}^{2}{\pm}m_{2}^{2}{-}m_{3}^{2}{-}s)$

and

$\displaystyle k^{2}$
$\displaystyle=\frac{16\beta}{\delta_{0}\delta_{1}\delta_{2}\delta_{3}}$
$\displaystyle n^{2}$ $\displaystyle=\frac{4\lambda_{3-}^{2}}{\psi_{+}^{2}}$
$\displaystyle\tilde{n}^{2}$ $\displaystyle=\frac{4\mu_{-}^{2}}{\rho^{2}}$
$\displaystyle(1{-}k^{2})$
$\displaystyle=\frac{\delta_{+}\delta_{01}\delta_{02}\delta_{03}}{\delta_{0}\delta_{1}\delta_{2}\delta_{3}}$
$\displaystyle n^{\prime 2}$
$\displaystyle=\frac{\lambda_{3-}^{2}}{\lambda_{3+}^{2}}$
$\displaystyle\tilde{n}^{\prime 2}$
$\displaystyle=\frac{\mu_{-}^{2}}{\mu_{+}^{2}}$ (92)

where again we have the relation $(1{-}k^{2})/n^{\prime 2}=1-k^{2}/n^{2}$ and
likewise for $\tilde{n}$. With this we may write the integrals
$h_{ij}=\int_{\gamma_{j}}\xi_{i}/Y$:

$\displaystyle h_{11}$
$\displaystyle=\frac{4K(k^{2})}{\sqrt{\delta_{0}\delta_{1}\delta_{2}\delta_{3}}}\qquad\qquad
h_{21}=\frac{-1}{\sqrt{\delta_{0}\delta_{1}\delta_{2}\delta_{3}}\,m_{1}^{2}}\left(K(k^{2})+\frac{qE(k^{2})}{\delta_{+}\delta_{01}\delta_{02}\delta_{03}}\right)$
$\displaystyle h_{31}$
$\displaystyle=\frac{2}{\sqrt{\delta_{0}\delta_{1}\delta_{2}\delta_{3}}\,\lambda_{3-}}\left(2\lambda_{1-}\lambda_{2-}K(k^{2})-\frac{\delta_{+}\delta_{01}\delta_{02}\delta_{03}\lambda_{3+}\Pi(n^{2},k^{2})}{\psi_{+}}\right)$
$\displaystyle h_{41}$
$\displaystyle=\frac{2}{\sqrt{\delta_{0}\delta_{1}\delta_{2}\delta_{3}}\,\mu_{-}}\left(2\lambda_{1-}\lambda_{2-}K(k^{2})-\frac{\delta_{+}\delta_{01}\delta_{02}\delta_{03}\mu_{+}\Pi(\tilde{n}^{2},k^{2})}{\rho}\right)$
$\displaystyle h_{12}$
$\displaystyle=\frac{-4iK(1{-}k^{2})}{\sqrt{\delta_{0}\delta_{1}\delta_{2}\delta_{3}}}\qquad\qquad
h_{22}=\frac{-i\big{(}qE(1{-}k^{2})+4(2\beta{+}m_{1}^{2}\psi_{-})K(1{-}k^{2})\big{)}}{\sqrt{\delta_{0}\delta_{1}\delta_{2}\delta_{3}}\,\delta_{+}\delta_{01}\delta_{02}\delta_{03}m_{1}^{2}}$
$\displaystyle h_{32}$
$\displaystyle=\frac{-4i}{\sqrt{\delta_{0}\delta_{1}\delta_{2}\delta_{3}}\,\lambda_{3-}}\left(\lambda_{1-}\lambda_{2-}K(1{-}k^{2})+\frac{2\beta\psi_{+}\Pi(n^{\prime
2},1{-}k^{2})}{\lambda_{3+}}\right)+i\pi$ (93) $\displaystyle h_{42}$
$\displaystyle=\frac{-4i}{\sqrt{\delta_{0}\delta_{1}\delta_{2}\delta_{3}}\,\mu_{-}}\left(\lambda_{1-}\lambda_{2-}K(1{-}k^{2})+\frac{2\beta\rho\Pi(\tilde{n}^{\prime
2},1{-}k^{2})}{\mu_{+}}\right)-\frac{i\pi}{(m_{2}^{2}{-}m_{1}^{2})(m_{3}^{2}{-}s)}$
$\displaystyle h_{13}$
$\displaystyle=h_{23}=h_{43}=h_{14}=h_{24}=h_{34}=0\qquad h_{33}=-2\pi i\qquad
h_{44}=\frac{2\pi i}{(m_{2}^{2}{-}m_{1}^{2})(m_{3}^{2}{-}s)}$

We may then form the integrals of the $\hat{\phi}$ integrands
$g_{ij}=\int_{\gamma_{j}}\phi_{i}/Y$ as combinations of the above, and then
form the period matrix $P_{ij}=f_{ik}g_{kj}$. Fixing it to be diagonal $P=2\pi
iI$ will fix all 16 $f_{ij}$. We will not list them all here, two examples are

$\displaystyle f_{11}$
$\displaystyle=i\sqrt{\delta_{0}\delta_{1}\delta_{2}\delta_{3}}\left(\frac{4(2\beta{+}m_{1}^{2}\psi_{-})}{q}K(1{-}k^{2})+E(1{-}k^{2})\right)\,,$
$\displaystyle f_{31}$
$\displaystyle=\left(\frac{\delta_{+}\delta_{01}\delta_{02}\delta_{03}}{q}K(k^{2})+E(k^{2})\right)\left(\frac{4\beta\psi_{+}}{\lambda_{3-}\lambda_{3+}\pi}\Pi(n^{\prime
2},1{-}k^{2})-\frac{\sqrt{\delta_{0}\delta_{1}\delta_{2}\delta_{3}}}{2}\right)$
(94)
$\displaystyle\;\;\;+\frac{-\delta_{+}\delta_{01}\delta_{02}\delta_{03}\lambda_{3+}}{\lambda_{3-}\psi_{+}\pi}\left(\frac{4(2\beta{+}m_{1}^{2}\psi_{-})}{q}K(1{-}k^{2})+E(1{-}k^{2})\right)\Pi(n^{2},k^{2})+\frac{\lambda_{1-}\lambda_{2-}}{\lambda_{3-}}.$

With this choice we get the differential equations in epsilon-factorized form
$\mathrm{d}J_{i}/\mathrm{d}x=\epsilon A^{(x)}_{ij}J_{j}$ for all the kinematic
variables $x$. The matrix entries are too large to be written here.

### The added expressions

In many of the examples done in this section, the result were deemed too large
to be written out in full. For that reason a file expressions.m has been added
to the arXiv version of this paper, which contains all the expressions in
Mathematica format. The expressions are named for instance Ass2m for the
matrix $A^{(s)}_{\text{s2m}}$ or frulesnpt for replacement rules for the
$f_{ij}$ in the case of the non-planar double triangle. To obtain the full
list of expressions defined in the file one might use the Mathematica command
Names["Global‘*"] after reading it in.

##  5 Discussion

###  5.1 Form and integration of the differential equations

One topic that has been avoided so far in this paper, is how to solve the
differential equations once the epsilon-factorized form has been obtained. For
Feynman integrals in the traditional canonical form it is often straight
forward to integrate them up order by order in the epsilon expansion, giving
results in the function class of generalized polylogarithms Goncharov:1998kja
, (at least in the case where all square-roots can be rationalized, which is
however known to not always be possible Brown:2020rda ). Many attempts have
been made to extend this function class to elliptic cases and beyond
brown2011multiple ; Adams:2014vja ; Broedel:2017kkb ; Adams:2017ejb , yet none
of those seem directly applicable to the form of the differential equations
derived in this paper. Further investigations will be needed into how to best
proceed analytically from the differential equations found here.

If, however, one is satisfied with a numerical approach, the development of
techniques to numerically integrate this type of differential equations is
developing quickly Boughezal:2007ny ; Czakon:2007qi ; Aglietti:2007as ;
Lee:2017qql ; Mandal:2018cdj ; Moriello:2019yhu ; Hidding:2020ytt . Not all of
these approaches require the differential equations to be epsilon factorized,
but they all benefit from from it and work faster in that case, so even in the
worst case scenario the techniques developed here can be used to speed up the
numerics.

We do however notice one property that the differential equations found on the
previous pages share with those in traditional canonical form and which might
cause optimism for the prospect of analytical integration: They are free of
higher poles at singular points including the point at infinity. This property
is not obvious, it might appear as if for instance the $A_{11}$-element of the
differential equation matrix found for the same-mass sunrise of eqs. (60)
would have a higher pole in $\sqrt{s}=-m$. But that is not the case. Changing
variable to $y=\sqrt{s}+m$, setting $m=1$, and expanding, that matrix element
becomes

$\displaystyle A^{\text{s1m};(y)}_{11}$
$\displaystyle=\frac{1}{\pi}\bigg{(}\big{(}6{-}2\pi
i{+}6\log(y/4)\big{)}\frac{1}{y}+\tfrac{1}{4}\big{(}15{-}7\pi
i{+}(33{-}6i\pi)\log(y/4){+}18\log^{2}(y/4)\big{)}$
$\displaystyle+\tfrac{1}{64}\big{(}273{-}39i\pi{+}(309{-}60i\pi)\log(y/4){+}180\log^{2}(y/4)\big{)}y+\mathcal{O}(y^{2})\bigg{)}$
(95)

and we see that no higher pole in $y$ is present, the candidate having
canceled with a term from the expansion of the elliptic integrals.

###  5.2 Freedom in basis choice

There is some amount of freedom in the algorithm discussed on the previous
pages. This freedom corresponds to performing row operation on the period
matrix. For instance one might define the master contours $\gamma_{i}$ as some
linear combination of the ones used in this paper. In all the elliptic
examples we chose $\gamma_{1}=\mathcal{C}_{\text{ii-iii}}$ and
$\gamma_{2}=\mathcal{C}_{\text{i-ii}}$, but one might equally well have chosen
for instance $\gamma_{2}=\mathcal{C}_{\text{iii-iv}}$ which would yield a
slightly different set of integrals in the end, and also there is nothing
preventing the choice of some $\mathbb{C}$-linear combination of contours,
such as $\gamma=\tfrac{5}{7}\mathcal{C}_{\text{ii-
iii}}+\tfrac{8i}{13}\mathcal{C}_{\text{i-ii}}$. Yet these kinds of
redefinitions seem to provide no simplification of the integrals or their
differential equations. Likewise one might take linear combinations of the
resulting integrals, since if the set $J$ has epsilon factorized differential
equations so will a $\mathbb{C}$-linear combination $\tilde{J}=BJ$ with
$B_{ij}\in\mathbb{C}$. For example for the case of the two-mass elliptic
sunrise, it seems from eqs. (82) as if it would provide a slight
simplification to perform the replacement
$J_{3}^{\text{s2m}}\rightarrow\tilde{J}_{3}^{\text{s2m}}=J_{3}^{\text{s2m}}-\tfrac{1}{2}J_{2}^{\text{s2m}}$.

Yet as mentioned in the introduction, there are previous examples in the
literature in which elliptic Feynman integrals have been put into a form that
allows for epsilon factorized differential equations, and which are not equal,
or in a form related through the above mentioned transformations, to the
expressions found here. In particular the same-mass sunrise of section 4.1 has
been discussed in ref. Adams:2018yfj , and the three-mass sunrise of section
4.4 has been discussed in ref. Bogner:2019lfa . It must be admitted that the
integrals found in those references are “nicer” than those found with the
algorithm outlined here. In particular the matrix of $f_{ij}$ coefficients are
in both cases found to be something more akin to a lower triangular form and
the differential equation matrix is nicer as well as we will see. Focusing on
the same-mass sunrise, the form found in ref. Adams:2018yfj  has $f_{12}=0$
which means that the master integrals there may be written

$\displaystyle\tilde{J}^{\text{s1m}}_{1}$
$\displaystyle=cf_{11}I^{\text{s1m}}_{11100}$
$\displaystyle\tilde{J}^{\text{s1m}}_{2}$
$\displaystyle=c(f_{21}I^{\text{s1m}}_{11100}+f_{22}I^{\text{s1m}}_{21100})$
(96)

where the values are found to

$\displaystyle f_{11}$ $\displaystyle=\frac{-\epsilon\pi R}{4K(k^{2})}$
$\displaystyle f_{21}$
$\displaystyle=\frac{RE(k^{2})}{4\pi}+\frac{\big{(}3(m^{2}{-}s)(9m^{2}{-}s)+2\epsilon(45m^{4}{-}30m^{2}s{+}s^{2})\big{)}K(k^{2})}{12\pi
R}$ (97) $\displaystyle f_{22}$
$\displaystyle=\frac{m^{2}R(\sqrt{m^{2}}{-}\sqrt{s})(3\sqrt{m^{2}}{+}\sqrt{s})K(k^{2})}{\pi(\sqrt{m^{2}}{+}\sqrt{s})^{2}}$

with the definitions of $k$ and $R$ being as given in sec. 4.1. We see that
these expressions are quite different from what was found in eqs. (58). Aside
from having three coefficients rather than four, eqs. (97) has simpler
kinematics dependence, and only complete elliptic integrals of $k^{2}$ appear
where eqs. (58) also contained complete elliptic integrals of argument
$(1{-}k^{2})$. On the other hand the $\pi$ and $\epsilon$ dependence is
simpler in eqs. (58).

From this we may compute the differential equations. They may be written as
$\mathrm{d}\tilde{J}/\mathrm{d}s=\epsilon A\tilde{J}^{\text{s1m}}$ with

$\displaystyle A_{11}$
$\displaystyle=\frac{9m^{4}{+}10m^{2}s{-}3s^{2}}{2s(9m^{2}{-}s)(m^{2}{-}s)}\qquad\qquad\qquad\qquad\qquad\quad\;\;A_{12}=\frac{3\pi^{2}(\sqrt{m^{2}}{+}\sqrt{s})^{2}}{4s(\sqrt{m^{2}}{-}\sqrt{s})(3\sqrt{m^{2}}{+}\sqrt{s})K(k^{2})^{2}}$
$\displaystyle A_{21}$
$\displaystyle=\frac{(3m^{2}{+}s)^{4}K(k^{2})^{2}}{3\pi^{2}s(\sqrt{m^{2}}{-}\sqrt{s})(3\sqrt{m^{2}}{-}\sqrt{s})^{2}(\sqrt{m^{2}}{+}\sqrt{s})^{4}(3\sqrt{m^{2}}{+}\sqrt{s})}\qquad\;A_{22}=\frac{9m^{4}{+}10m^{2}s{-}3s^{3}}{2s(9m^{2}{-}s)(m^{2}{-}s)}$
(98)

which is simpler than the expressions in eqs. (60) by a substantial amount, we
see for instance that the only elliptic integral appearing is $K(k^{2})$, and
that the diagonal entries are free of even this. On the other hand the
$K(k^{2})$ appears both in numerator and denominator, where in eqs. (60) the
elliptic integrals appeared in the numerator only.

We may compute the period matrix of this example. The result is $P_{ij}$ with

$\displaystyle P_{11}$
$\displaystyle=-\pi\epsilon\qquad\qquad\qquad\qquad\qquad\qquad
P_{12}=\frac{i\pi\epsilon K(1{-}k^{2})}{K(k^{2})}$ $\displaystyle P_{21}$
$\displaystyle=\frac{2\epsilon
K(k^{2})}{\pi}\left(\frac{2(3m^{2}{-}s)(3m^{2}{+}s)K(k^{2})}{3R^{2}}-E(k^{2})\right)$
(99) $\displaystyle P_{22}$ $\displaystyle=\frac{-i}{2}-\frac{2i\epsilon
K(k^{2})}{\pi}\left(\frac{(3\sqrt{m^{2}}{+}\sqrt{s})(3m^{3}{-}9m^{2}\sqrt{s}{-}3ms{+}s^{3/2})}{R^{2}}K(1{-}k^{2})+E(1{-}k^{2})\right)$

and in addition it is worth noticing that the determinant of this matrix is
given as

$\displaystyle\text{det}(P)$
$\displaystyle=\tfrac{1}{2}i\pi\epsilon(1{+}2\epsilon)$ (100)

We see that $P_{11}$ and $\text{det}(P)$ both are constants, in the sense that
they are free of any kinematic dependence including through elliptic
integrals. This may be used as a guideline for obtaining an
$\epsilon$-factorized form alternative to the one discussed in this paper. We
do an example of this in Appendix A, but it is not clear if this can be
generalized to elliptic sectors with more than two master integrals.

It would be great if some principle could be found to a priori generate the
transformation between the forms discussed in this paper, and forms similar to
the one discussed above.

###  5.3 The number of cycles

In principle there is a one-to-one correspondence between the number of master
integrals and the number of independent cycles Lee:2013hzt ; Bitoun:2017nre ;
Frellesvig:2019uqt . This correspondence is utilized by the Lee-Pomeransky
criterion Lee:2013hzt , which states that the number of master integrals $\nu$
is given as

$\displaystyle\nu$ $\displaystyle=\text{number of solutions to
``}\omega=0\text{''}\qquad\text{where}\qquad\omega:=\mathrm{d}\text{log}(u)$
(101)

We can test this for the examples discussed in secs. 2 and 4 of this paper.
The result is that the numbers agree in all the cases, except for the same-
mass sunrise of sec. 4.1 and the nonplanar double triangle of sec. 4.2, which
are the two cases in which an elliptic sector has two master integrals. For
the case of the same-mass sunrise, this “miscounting” is well known, and is
discussed in further detail in refs. Frellesvig:2019kgj ; Weinzierl:2020xyy
where it is shown to be caused by peculiarities in the interplay between the
loop-by-loop Baikov parametrization and the maximal cut.

Focusing on the same-mass sunrise it might seem as if a valid third master
integral would be $I_{111\scalebox{0.5}[0.75]{$\,-$}10}^{\text{s1m}}$
corresponding to $\hat{\phi}_{3}{=}z$, with a corresponding third master
contour surrounding the pole at infinity $\gamma_{3}{=}\mathcal{C}_{\infty}$
making the set of cycles look as in the case of two-mass elliptic sunrise of
fig. 6. Yet we know from IBP relations that

$\displaystyle I_{111\scalebox{0.5}[0.75]{$\,-$}10}^{\text{s1m}}$
$\displaystyle=\frac{3m^{2}{+}s}{3}I_{11100}^{\text{s1m}}\,+\,\text{lower}$
(102)

and it might be instructive to see how this relation is realized over the
basic cycles. Using the $\hat{\phi}_{i}$ and $\gamma_{j}$ from sec. 4.1, and
again defining $g_{ij}=\int_{\gamma_{j}}\hat{\phi}_{i}\mathrm{d}z/Y$, we get

$\displaystyle g_{31}$
$\displaystyle=\int_{\gamma_{1}}\\!\frac{\hat{\phi}_{3}\mathrm{d}z}{Y}=2(\sqrt{m^{2}}{-}\sqrt{s})\left(2\sqrt{m^{2}}K(k^{2})-(3\sqrt{m^{2}}{+}\sqrt{s})\Pi(n^{2},k^{2})\right)\\!\big{/}R$
$\displaystyle g_{32}$
$\displaystyle=\int_{\gamma_{2}}\\!\frac{\hat{\phi}_{3}\mathrm{d}z}{Y}=-4i\left(\sqrt{m^{2}}(\sqrt{m^{2}}{-}\sqrt{s})K(1{-}k^{2})+2\sqrt{m^{2}}\sqrt{s}\Pi(n^{\prime
2},1{-}k^{2})\right)\\!\big{/}R+i\pi$ (103) $\displaystyle g_{13}$
$\displaystyle=\int_{\gamma_{3}}\\!\frac{\hat{\phi}_{1}\mathrm{d}z}{Y}=0\qquad\qquad
g_{23}=\int_{\gamma_{3}}\\!\frac{\hat{\phi}_{2}\mathrm{d}z}{Y}=0\qquad\qquad
g_{33}=\int_{\gamma_{3}}\\!\frac{\hat{\phi}_{3}\mathrm{d}z}{Y}=-2\pi i$

where $k^{2}$ and $R$ are as defined in sec. 4.1, and where

$\displaystyle
n^{2}=\frac{4m^{2}}{(\sqrt{m^{2}}{+}\sqrt{s})^{2}}\qquad\quad\text{and}\qquad\quad
n^{\prime
2}=\frac{(\sqrt{m^{2}}{-}\sqrt{s})^{2}}{(\sqrt{m^{2}}{+}\sqrt{s})^{2}}$ (104)

If eq. (102) were to hold on the contour $\gamma_{j}$ it would mean that
$g_{3j}=((3m^{2}{+}s)/3)g_{1j}$, so let us investigate the status of that
relation on each contour.

Starting by $\gamma_{1}$, it might seem unlikely for the relation to hold
since $g_{31}$ contains $\Pi(n^{2},k^{2})$, which $g_{11}$ does not. But it
does hold exactly, due to the relation

$\displaystyle\Pi\\!\left(\frac{4}{(1{+}x)^{2}},\frac{-16x}{(x{-}3)(1{+}x)^{3}}\right)$
$\displaystyle=\frac{2x}{3(x{-}1)}K\\!\left(\frac{-16x}{(x{-}3)(1{+}x)^{3}}\right)$
(105)

(valid for $x\in[-1,1]$). On the other hand for contour $\gamma_{3}$ the
relation clearly does not hold since $g_{13}=0$, while $g_{33}=-2\pi i$ which
is then the amount with which the relation is broken. At the last contour
$\gamma_{2}$ we get using a similar relation for a special point of
$\Pi(n^{2},k^{2})$, that

$\displaystyle g_{32}-\frac{3m^{2}{+}s}{3}g_{12}$ $\displaystyle=\frac{2\pi
i}{3}$ (106)

so we see that when the peculiarities of the loop-by-loop Baikov
parametrization ruins the counting of master contours, relations such as eq.
(102) break only by a factor proportional to $2\pi i$.

For the nonplanar double triangle of sec. 4.2 the behaviour is extremely
similar. There the extra IBP-derived integral relation relating an integral
with a pole at infinity to the original two, is

$\displaystyle I_{111111\scalebox{0.5}[0.75]{$\,-$}1}^{\text{npt}}$
$\displaystyle=\frac{-s}{2}I_{1111110}^{\text{npt}}\,+\,\text{lower}$ (107)

and that relation holds exactly on contour $\gamma_{1}$ and is broken with
factors proportional to $2\pi i$ on $\gamma_{2}$ and a potential $\gamma_{3}$
surrounding a pole at infinity. The main difference is that the
$\Pi(n^{2},k^{2})$ special value relation that makes this explicit is even
simpler in that case.

###  5.4 Freedom in intermediate basis choice

One might consider if the intermediate bases chosen in the examples discussed
in secs. 2 and 4 have some special property that allows for the algorithm to
go through that would not be shared with any arbitrary basis choice. The
answer to this is mostly no, any set of valid master integrals might be chosen
as the intermediate basis. Let us illustrate that with an example, that of the
two-mass elliptic sunrise of sec. 4.3. Instead of the choice of
($I^{\text{s2m}}_{11100}$, $I^{\text{s2m}}_{21100}$,
$I^{\text{s2m}}_{111\scalebox{0.5}[0.75]{$\,-$}10}$) made in sec. 4.3, one
might instead pick for instance the set ($I^{\text{s2m}}_{11100}$,
$I^{\text{s2m}}_{111\scalebox{0.5}[0.75]{$\,-$}10}$,
$I^{\text{s2m}}_{111\scalebox{0.5}[0.75]{$\,-$}20}$) which is the choice
naturally made by Kira.

With this one can go through the same procedure as before, define

$\displaystyle\tilde{J}_{i}=\tilde{f}_{i1}I^{\text{s2m}}_{11100}+\tilde{f}_{i2}I^{\text{s2m}}_{111\scalebox{0.5}[0.75]{$\,-$}10}+\tilde{f}_{i3}I^{\text{s2m}}_{111\scalebox{0.5}[0.75]{$\,-$}20}$
(108)

compute the period matrix, and impose it to be proportional to $I$ thereby
fixing the $\tilde{f}_{ij}$. The resulting $\tilde{J}_{i}$ equal those from
sec. 4.3 up to integrals in lower sectors as well as terms of order
$\epsilon$:

$\displaystyle\tilde{J}_{i}=J_{i}+\mathcal{O}(\epsilon)+\text{lower}$ (109)

Additionally the new integrals have epsilon factorized differential equations
$\mathrm{d}\tilde{J}_{i}/\mathrm{d}x=\epsilon\tilde{A}^{(x)}_{ij}\tilde{J}_{j}$
in each of the three kinematic variables $x\in\\{m_{a}^{2},m_{b}^{2},s\\}$.
The expressions for $\tilde{f}$ and the $\tilde{A}^{(x)}$ can be found in the
file added to the arXiv version of this paper.

There are however cases in which not any choice of intermediate basis is
equally convenient. This is true in the cases that have the “miscounting” of
independent cycles discussed in sec. 5.3, such as the same-mass elliptic
sunrise. There we saw that a counting of independent cycles seems to include a
third master contour surrounding a pole at infinity $\mathcal{C}_{\infty}$,
and that the relation reducing integrals with support there to integrals with
support only on the original two $\mathcal{C}_{\text{ii-iii}}$ and
$\mathcal{C}_{\text{i-ii}}$ is realized only up to factors proportional to
$2\pi i$. For the same-mass elliptic sunrise, ($I^{\text{s1m}}_{11100}$,
$I^{\text{s1m}}_{111\scalebox{0.5}[0.75]{$\,-$}20}$) which is chosen as the
default by Kira would be a valid alternative to the intermediate basis
discussed in section 4.1. But since the integrand of
$I^{\text{s1m}}_{111\scalebox{0.5}[0.75]{$\,-$}20}$ has support on
$\mathcal{C}_{\infty}$, such a factor proportional to $2\pi i$ has to be
removed by hand before a basis with epsilon factorized differential equations
is obtained.

###  5.5 Further discussion, open questions, and conclusions

In all the examples in this paper, we only discussed the highest sector for
the integral families we looked at, which allowed us to only approach the
system on its maximal cut. Yet for practical applications one would need the
whole system of differential equations to be in epsilon-factorized form, not
just the highest sector. Presumably the approach discussed here could be
generalized away from the maximal cut, by introducing a multivariate notion of
master contours, and indeed this is the approach taken in the prescriptive
unitarity scheme Bourjaily:2017wjl ; Bourjaily:2021vyj . Yet to bring
subsectors into epsilon-factorized form, I do not believe this is needed. It
will be enough to fit coefficients of lower sector integrals using an ansatz
with its free coefficients fixed from imposing the epsilon-factorized property
of the differential equations, starting from the highest subsectors going
down, along the lines of the method proposed in ref. Gehrmann:2014bfa  for the
non-elliptic case.

Likewise the discussion here was limited to the case where all subsectors were
free of elliptic contributions and could be put into canonical form in the
traditional sense. For integrals for which that is not the case (such as the
kite integrals that will be obtained by promoting $D_{4}$ and $D_{5}$ to
genuine propagators for the elliptic sunrise integrals) presumably a similar
approach can be taken, but that is a matter for further study.

The canonical form proposed in ref. Henn:2013pwa  not only implied an epsilon-
factorized form for the differential equations of eq. (1), but further imply
that the distinct differential equations for each kinematic variable are
joined in one differential of dlog-form $\mathrm{d}J=\epsilon(\mathrm{d}M)J$
where $M$ is a matrix containing only logarithms of algebraic functions of the
kinematic variables. An obvious question is if the differential equations
found here can likewise be unified into one differential form
$\mathrm{d}J=\epsilon(\mathrm{d}M)J$, but where $M$ obviously no longer will
contain solely logarithms. In the polylogarithmic case, the arguments of those
logarithms are the symbol letters Goncharov:2010jf ; Duhr:2012fh  of the
problem, so an obvious open question is if the entries of $M$ in the elliptic
case, if it even exists, will be related to the entries of the elliptic
generalization of the coaction Broedel:2018iwv  in a similar straight forward
fashion. For the polylogarithmic case, the fact that the algorithm discussed
in this paper works in the first place, may be seen as following from the
existence of an iterated dlog representation for the whole integral Arkani-
Hamed:2010pyv ; WasserMSc ; Herrmann:2019upk ; Dlapa:2021qsl , and the
question of how such a representation might generalize to the elliptic case is
presumably also linked to the above.

The answers to such questions will help clarify what the “correct”
generalization of terms such as pure Broedel:2018qkq , canonical, weight, etc.
that have proven extremely helpful in the understanding of polylogarithmic
Feynman integrals, will be to the elliptic case and beyond, if such
generalizations can even be made in a useful way. I did not attempt to make
any such generalizations in this work, but I hope that my approach contributes
with insights that may help bring us to a point where such a step might be
taken.

The algorithm presented in this paper is able to systematically bring
differential equations for elliptic Feynman integrals into a form for which
the differential equations are epsilon factorized. We demonstrated this for a
number of examples of varying complexity, and it is the hope that this
development will be one step on the way towards bringing the understanding of
elliptic Feynman integrals to the same level as the polylogarithmic case.

#### Acknowledgments

The author would like to thank Pierpaolo Mastrolia and the rest of the
intersection group in Padova, as well as Cristian Vergu and the rest of the
NBIA amplitudes group, for many helpful discussions. Particular thanks go to
Pierpaolo Mastrolia, Manoj Mandal, Matt von Hippel, and Stefan Weinzierl, for
carefully reading through the manuscript at the draft stage and providing
helpful feedback.

This project has received funding from the European Union’s Horizon 2020
research and innovation program under the Marie Skłodowska-Curie grant
agreement No. 847523 ’INTERACTIONS’. The work has been partially supported by
a Carlsberg Foundation Reintegration Fellowship.

##  Appendix A A different form for the elliptic nonplanar double triangle

In the example discussed in section 5.2, we saw that a valid form for the
equal mass sunrise integral Adams:2018yfj  resulting in an
$\epsilon$-factorized differential equation could be made without the diagonal
period matrix otherwise discussed in this paper. That example had $f_{12}=0$
and $P_{11}$ and $\text{det}(P)$ both being constants. Using these properties
as a guideline, we were able to find an epsilon factorized basis for the
nonplanar double triangle of section 4.2, which has two master integrals in
the elliptic sector similarly to the equal mass sunrise.

We write the general form

$\displaystyle\tilde{J}^{\text{npt}}_{1}$
$\displaystyle=cf_{11}I^{\text{npt}}_{1111110}$
$\displaystyle\tilde{J}^{\text{npt}}_{2}$
$\displaystyle=c(f_{21}I^{\text{npt}}_{1111110}+f_{22}I^{\text{npt}}_{2111110})$
(110)

where we have used $f_{12}=0$. Imposing $P_{11}$ and $\text{det}(P)$ being
constants puts two constraints on the remaining $f_{ij}$ while the third may
be fitted by requiring the differential equation to be $\epsilon$-factorized.
We find the values

$\displaystyle f_{11}$
$\displaystyle=\frac{\epsilon\,s^{3/2}(\sqrt{16m^{2}{+}s}{+}\sqrt{s})}{4K(k^{2})}$
$\displaystyle f_{21}$
$\displaystyle=-s^{3/2}\left((\sqrt{16m^{2}{+}s}{+}\sqrt{s})E(k^{2})+(\sqrt{16m^{2}{+}s}{-}\sqrt{s})\left(1+\epsilon\frac{24m^{2}{+}s}{m^{2}}\right)K(k^{2})\right)$
(111) $\displaystyle f_{22}$
$\displaystyle=-s^{3/2}(16m^{2}{+}s)(\sqrt{16m^{2}{+}s}{-}\sqrt{s})K(k^{2})$

which corresponds to the period matrix (given the two contours defined in eqs.
(68))

$\displaystyle P_{11}$ $\displaystyle=2\epsilon\qquad\qquad\qquad\qquad\qquad
P_{12}=-2i\epsilon\frac{K(1{-}k^{2})}{K(k^{2})}$ $\displaystyle P_{21}$
$\displaystyle=2\epsilon
K(k^{2})\left(8E(k^{2})-\frac{(16m^{2}{+}s)(\sqrt{16m^{2}{+}s}{-}\sqrt{s})}{m^{2}(\sqrt{16m^{2}{+}s}{+}\sqrt{s})}K(k^{2})\right)$
(112) $\displaystyle P_{22}$ $\displaystyle=4i\pi+2i\epsilon
K(k^{2})\left(8E(1{-}k^{2})+\frac{\big{(}(8m^{2}{+}s)\sqrt{16m^{2}{+}s}-(24m^{2}{+}s)\sqrt{s}\big{)}}{m^{2}(\sqrt{s}{+}\sqrt{16m^{2}{+}s})}K(1{-}k^{2})\right)$

with the property

$\displaystyle\text{det}(P)$ $\displaystyle=8i\pi\epsilon(1{+}2\epsilon)$
(113)

and we see the two requirements imposed above, being apparent.

The corresponding $s$-differential equation
$\mathrm{d}\tilde{J}^{\text{npt}}_{i}/\mathrm{d}s=\epsilon\tilde{A}_{ij}\tilde{J}^{\text{npt}}_{j}$
has

$\displaystyle\tilde{A}_{11}$
$\displaystyle=\frac{-(8m^{2}{+}s)}{s(16m^{2}{+}s)}\qquad\qquad\qquad\qquad\;\;\tilde{A}_{12}=\frac{(\sqrt{16m^{2}{+}s}{+}\sqrt{s})^{2}}{8\sqrt{2}s(16m^{2}{+}s)K(k^{2})^{2}}$
$\displaystyle\tilde{A}_{21}$
$\displaystyle=\frac{8\sqrt{2}(8m^{2}{+}s)^{2}K(k^{2})^{2}}{s(16m^{2}{+}s)(\sqrt{16m^{2}{+}s}{+}\sqrt{s})^{2}}\qquad\qquad\qquad\tilde{A}_{22}=\frac{-(8m^{2}{+}s)}{s(16m^{2}{+}s)}$
(114)

with a similar differential equation in the other variable $m^{2}$. This form
of the differential equation resembles that given in eq. (98), and as such it
will likely be suitable for the integration approach proposed in ref.
Adams:2018yfj .

It would be interesting to generalize these considerations to a more general
elliptic case, or preferably to find a principle that allows for reliably
generating a transformation that can map between the forms discussed in the
bulk of this paper and forms similar to the above, but to what extend that is
possible is a question for the future.

## References

  * (1) G. Barucchi and G. Ponzano, Differential equations for one-loop generalized feynman integrals, J. Math. Phys. 14 (1973) 396–401.
  * (2) A. Kotikov, Differential equations method. New technique for massive Feynman diagram calculation, Physics Letters B 254 (1991), no. 1 158 – 164.
  * (3) T. Gehrmann and E. Remiddi, Differential equations for two loop four point functions, Nucl. Phys. B580 (2000) 485–518, [hep-ph/9912329].
  * (4) K. G. Chetyrkin and F. V. Tkachov, Integration by Parts: The Algorithm to Calculate beta Functions in 4 Loops, Nucl. Phys. B192 (1981) 159–204.
  * (5) S. Laporta, High precision calculation of multiloop Feynman integrals by difference equations, Int. J. Mod. Phys. A15 (2000) 5087–5159, [hep-ph/0102033].
  * (6) A. V. Smirnov, Algorithm FIRE – Feynman Integral REduction, JHEP 10 (2008) 107, [arXiv:0807.3243].
  * (7) P. Maierhöfer, J. Usovitsch, and P. Uwer, Kira—A Feynman integral reduction program, Comput. Phys. Commun. 230 (2018) 99–112, [arXiv:1705.05610].
  * (8) C. Anastasiou and A. Lazopoulos, Automatic integral reduction for higher order perturbative calculations, JHEP 07 (2004) 046, [hep-ph/0404258].
  * (9) A. von Manteuffel and C. Studerus, Reduze 2 - Distributed Feynman Integral Reduction, arXiv:1201.4330.
  * (10) R. N. Lee, Presenting LiteRed: a tool for the Loop InTEgrals REDuction, arXiv:1212.2685.
  * (11) J. M. Henn, Multiloop integrals in dimensional regularization made simple, Phys. Rev. Lett. 110 (2013) 251601, [arXiv:1304.1806].
  * (12) D. J. Broadhurst, J. Fleischer, and O. V. Tarasov, Two loop two point functions with masses: Asymptotic expansions and Taylor series, in any dimension, Z. Phys. C 60 (1993) 287–302, [hep-ph/9304303].
  * (13) F. A. Berends, M. Buza, M. Bohm, and R. Scharf, Closed expressions for specific massive multiloop selfenergy integrals, Z. Phys. C 63 (1994) 227–234.
  * (14) S. Bauberger, M. Böhm, G. Weiglein, F. A. Berends, and M. Buza, Calculation of two-loop self-energies in the electroweak Standard Model, Nucl. Phys. B Proc. Suppl. 37 (1994), no. 2 95–114, [hep-ph/9406404].
  * (15) S. Bauberger, F. A. Berends, M. Bohm, and M. Buza, Analytical and numerical methods for massive two loop selfenergy diagrams, Nucl. Phys. B 434 (1995) 383–407, [hep-ph/9409388].
  * (16) S. Bauberger and M. Bohm, Simple one-dimensional integral representations for two loop selfenergies: The Master diagram, Nucl. Phys. B 445 (1995) 25–48, [hep-ph/9501201].
  * (17) M. Caffo, H. Czyz, S. Laporta, and E. Remiddi, The Master differential equations for the two loop sunrise selfmass amplitudes, Nuovo Cim. A 111 (1998) 365–389, [hep-th/9805118].
  * (18) S. Laporta and E. Remiddi, Analytic treatment of the two loop equal mass sunrise graph, Nucl. Phys. B704 (2005) 349–386, [hep-ph/0406160].
  * (19) S. Groote, J. Korner, and A. Pivovarov, On the evaluation of a certain class of Feynman diagrams in x-space: Sunrise-type topologies at any loop order, Annals Phys. 322 (2007) 2374–2445, [hep-ph/0506286].
  * (20) D. H. Bailey, J. M. Borwein, D. Broadhurst, and M. L. Glasser, Elliptic integral evaluations of Bessel moments, J. Phys. A 41 (2008) 205203, [arXiv:0801.0891].
  * (21) S. Müller-Stach, S. Weinzierl, and R. Zayadeh, A Second-Order Differential Equation for the Two-Loop Sunrise Graph with Arbitrary Masses, Commun. Num. Theor. Phys. 6 (2012) 203–222, [arXiv:1112.4360].
  * (22) F. Brown and A. Levin, Multiple Elliptic Polylogarithms, arXiv:1110.6917.
  * (23) S. Caron-Huot and K. J. Larsen, Uniqueness of two-loop master contours, JHEP 10 (2012) 026, [arXiv:1205.0801].
  * (24) S. Groote, J. Korner, and A. Pivovarov, A Numerical Test of Differential Equations for One- and Two-Loop sunrise Diagrams using Configuration Space Techniques, Eur. Phys. J. C 72 (2012) 2085, [arXiv:1204.0694].
  * (25) L. Adams, C. Bogner, and S. Weinzierl, The two-loop sunrise graph with arbitrary masses, J. Math. Phys. 54 (2013) 052303, [arXiv:1302.7004].
  * (26) S. Bloch and P. Vanhove, The elliptic dilogarithm for the sunset graph, J. Number Theor. 148 (2015) 328–364, [arXiv:1309.5865].
  * (27) E. Remiddi and L. Tancredi, Schouten identities for Feynman graph amplitudes; The Master Integrals for the two-loop massive sunrise graph, Nucl. Phys. B880 (2014) 343–377, [arXiv:1311.3342].
  * (28) L. Adams, C. Bogner, and S. Weinzierl, The two-loop sunrise graph in two space-time dimensions with arbitrary masses in terms of elliptic dilogarithms, J. Math. Phys. 55 (2014), no. 10 102301, [arXiv:1405.5640].
  * (29) J. Broedel, C. R. Mafra, N. Matthes, and O. Schlotterer, Elliptic multiple zeta values and one-loop superstring amplitudes, JHEP 07 (2015) 112, [arXiv:1412.5535].
  * (30) S. Bloch, M. Kerr, and P. Vanhove, A Feynman integral via higher normal functions, Compos. Math. 151 (2015), no. 12 2329–2375, [arXiv:1406.2664].
  * (31) L. Adams, C. Bogner, and S. Weinzierl, The two-loop sunrise integral around four space-time dimensions and generalisations of the Clausen and Glaisher functions towards the elliptic case, J. Math. Phys. 56 (2015), no. 7 072303, [arXiv:1504.03255].
  * (32) J. Broedel, N. Matthes, and O. Schlotterer, Relations between elliptic multiple zeta values and a special derivation algebra, J. Phys. A 49 (2016), no. 15 155203, [arXiv:1507.02254].
  * (33) L. Adams, C. Bogner, and S. Weinzierl, The iterated structure of the all-order result for the two-loop sunrise integral, J. Math. Phys. 57 (2016), no. 3 032304, [arXiv:1512.05630].
  * (34) S. Bloch, M. Kerr, and P. Vanhove, Local mirror symmetry and the sunset Feynman integral, Adv. Theor. Math. Phys. 21 (2017) 1373–1453, [arXiv:1601.08181].
  * (35) G. Passarino, Elliptic Polylogarithms and Basic Hypergeometric Functions, Eur. Phys. J. C 77 (2017), no. 2 77, [arXiv:1610.06207].
  * (36) A. Primo and L. Tancredi, On the maximal cut of Feynman integrals and the solution of their differential equations, Nucl. Phys. B916 (2017) 94–116, [arXiv:1610.08397].
  * (37) E. Remiddi and L. Tancredi, Differential equations and dispersion relations for Feynman amplitudes. The two-loop massive sunrise and the kite integral, Nucl. Phys. B907 (2016) 400–444, [arXiv:1602.01481].
  * (38) R. Bonciani, V. Del Duca, H. Frellesvig, J. M. Henn, F. Moriello, and V. A. Smirnov, Two-loop planar master integrals for Higgs$\to 3$ partons with full heavy-quark mass dependence, JHEP 12 (2016) 096, [arXiv:1609.06685].
  * (39) D. Broadhurst, Feynman integrals, L-series and Kloosterman moments, Commun. Num. Theor. Phys. 10 (2016) 527–569, [arXiv:1604.03057].
  * (40) L. Adams and S. Weinzierl, Feynman integrals and iterated integrals of modular forms, Commun. Num. Theor. Phys. 12 (2018) 193–251, [arXiv:1704.08895].
  * (41) A. von Manteuffel and L. Tancredi, A non-planar two-loop three-point function beyond multiple polylogarithms, JHEP 06 (2017) 127, [arXiv:1701.05905].
  * (42) J. Ablinger, J. Blümlein, A. De Freitas, M. van Hoeij, E. Imamoglu, C. G. Raab, C. S. Radu, and C. Schneider, Iterated Elliptic and Hypergeometric Integrals for Feynman Diagrams, J. Math. Phys. 59 (2018), no. 6 062305, [arXiv:1706.01299].
  * (43) E. Remiddi and L. Tancredi, An Elliptic Generalization of Multiple Polylogarithms, Nucl. Phys. B 925 (2017) 212–251, [arXiv:1709.03622].
  * (44) M. Hidding and F. Moriello, All orders structure and efficient computation of linearly reducible elliptic Feynman integrals, JHEP 01 (2019) 169, [arXiv:1712.04441].
  * (45) J. L. Bourjaily, A. J. McLeod, M. Spradlin, M. von Hippel, and M. Wilhelm, Elliptic Double-Box Integrals: Massless Scattering Amplitudes beyond Polylogarithms, Phys. Rev. Lett. 120 (2018), no. 12 121603, [arXiv:1712.02785].
  * (46) J. Broedel, C. Duhr, F. Dulat, and L. Tancredi, Elliptic polylogarithms and iterated integrals on elliptic curves. Part I: general formalism, JHEP 05 (2018) 093, [arXiv:1712.07089].
  * (47) J. Broedel, C. Duhr, F. Dulat, and L. Tancredi, Elliptic polylogarithms and iterated integrals on elliptic curves II: an application to the sunrise integral, Phys. Rev. D 97 (2018), no. 11 116009, [arXiv:1712.07095].
  * (48) J. Broedel, N. Matthes, G. Richter, and O. Schlotterer, Twisted elliptic multiple zeta values and non-planar one-loop open-string amplitudes, J. Phys. A 51 (2018), no. 28 285401, [arXiv:1704.03449].
  * (49) R. N. Lee, A. V. Smirnov, and V. A. Smirnov, Solving differential equations for Feynman integrals by expansions near singular points, JHEP 03 (2018) 008, [arXiv:1709.07525].
  * (50) C. Bogner, A. Schweitzer, and S. Weinzierl, Analytic continuation and numerical evaluation of the kite integral and the equal mass sunrise integral, Nucl. Phys. B 922 (2017) 528–550, [arXiv:1705.08952].
  * (51) S. Groote and J. Körner, Coordinate space calculation of two- and three-loop sunrise-type diagrams, elliptic functions and truncated Bessel integral identities, Nucl. Phys. B 938 (2019) 416–425, [arXiv:1804.10570].
  * (52) L. Adams, E. Chaubey, and S. Weinzierl, Planar Double Box Integral for Top Pair Production with a Closed Top Loop to all orders in the Dimensional Regularization Parameter, Phys. Rev. Lett. 121 (2018), no. 14 142001, [arXiv:1804.11144].
  * (53) L. Adams, E. Chaubey, and S. Weinzierl, Analytic results for the planar double box integral relevant to top-pair production with a closed top loop, arXiv:1806.04981.
  * (54) L. Adams and S. Weinzierl, The $\varepsilon$-form of the differential equations for Feynman integrals in the elliptic case, Phys. Lett. B781 (2018) 270–278, [arXiv:1802.05020].
  * (55) J. Ablinger, J. Blümlein, P. Marquard, N. Rana, and C. Schneider, Automated Solution of First Order Factorizable Systems of Differential Equations in One Variable, Nucl. Phys. B 939 (2019) 253–291, [arXiv:1810.12261].
  * (56) J. Broedel, C. Duhr, F. Dulat, B. Penante, and L. Tancredi, Elliptic Feynman integrals and pure functions, JHEP 01 (2019) 023, [arXiv:1809.10698].
  * (57) C. Bogner, S. Müller-Stach, and S. Weinzierl, The unequal mass sunrise integral expressed through iterated integrals on $\overline{\mathcal{M}}_{1,3}$, Nucl. Phys. B 954 (2020) 114991, [arXiv:1907.01251].
  * (58) J. L. Bourjaily, N. Kalyanapuram, C. Langer, K. Patatoukos, and M. Spradlin, An Elliptic Yangian-Invariant, ‘Leading Singularity’, arXiv:2012.14438.
  * (59) H. Frellesvig, C. Vergu, M. Volk, and M. von Hippel, Cuts and Isogenies, JHEP 05 (2021) 064, [arXiv:2102.02769].
  * (60) J. L. Bourjaily, N. Kalyanapuram, C. Langer, and K. Patatoukos, Prescriptive Unitarity with Elliptic Leading Singularities, arXiv:2102.02210.
  * (61) J. M. Henn, Lectures on differential equations for Feynman integrals, J. Phys. A48 (2015) 153001, [arXiv:1412.2296].
  * (62) R. N. Lee, Reducing differential equations for multiloop master integrals, JHEP 04 (2015) 108, [arXiv:1411.0911].
  * (63) M. Argeri, S. Di Vita, P. Mastrolia, E. Mirabella, J. Schlenk, U. Schubert, and L. Tancredi, Magnus and Dyson Series for Master Integrals, JHEP 03 (2014) 082, [arXiv:1401.2979].
  * (64) T. Gehrmann, A. von Manteuffel, L. Tancredi, and E. Weihs, The two-loop master integrals for $q\overline{q}\to VV$, JHEP 06 (2014) 032, [arXiv:1404.4853].
  * (65) M. Höschele, J. Hoff, and T. Ueda, Adequate bases of phase space master integrals for gg $\to$ h at NNLO and beyond, JHEP 09 (2014) 116, [arXiv:1407.4049].
  * (66) C. Meyer, Transforming differential equations of multi-loop Feynman integrals into canonical form, JHEP 04 (2017) 006, [arXiv:1611.01087].
  * (67) H. Frellesvig and C. G. Papadopoulos, Cuts of Feynman Integrals in Baikov representation, JHEP 04 (2017) 083, [arXiv:1701.07356].
  * (68) P. Wasser, Analytic properties of Feynman integrals for scattering amplitudes, M.Sc. (2016) [https://publications.ub.uni-mainz.de/theses/frontdoor.php?source opus=100001967].
  * (69) M. Prausa, epsilon: A tool to find a canonical basis of master integrals, Comput. Phys. Commun. 219 (2017) 361–376, [arXiv:1701.00725].
  * (70) O. Gituliar and V. Magerya, Fuchsia: a tool for reducing differential equations for Feynman master integrals to epsilon form, Comput. Phys. Commun. 219 (2017) 329–338, [arXiv:1701.04269].
  * (71) C. Dlapa, J. Henn, and K. Yan, Deriving canonical differential equations for Feynman integrals from a single uniform weight integral, JHEP 05 (2020) 025, [arXiv:2002.02340].
  * (72) J. Henn, B. Mistlberger, V. A. Smirnov, and P. Wasser, Constructing d-log integrands and computing master integrals for three-loop four-particle scattering, JHEP 04 (2020) 167, [arXiv:2002.09492].
  * (73) J. Chen, X. Jiang, X. Xu, and L. L. Yang, Constructing canonical Feynman integrals with intersection theory, Phys. Lett. B 814 (2021) 136085, [arXiv:2008.03045].
  * (74) C. Dlapa, X. Li, and Y. Zhang, Leading singularities in Baikov representation and Feynman integrals with uniform transcendental weight, arXiv:2103.04638.
  * (75) N. Arkani-Hamed, J. L. Bourjaily, F. Cachazo, and J. Trnka, Local Integrals for Planar Scattering Amplitudes, JHEP 06 (2012) 125, [arXiv:1012.6032].
  * (76) J. L. Bourjaily, E. Herrmann, and J. Trnka, Prescriptive Unitarity, JHEP 06 (2017) 059, [arXiv:1704.05460].
  * (77) R. N. Lee and V. A. Smirnov, The Dimensional Recurrence and Analyticity Method for Multicomponent Master Integrals: Using Unitarity Cuts to Construct Homogeneous Solutions, JHEP 12 (2012) 104, [arXiv:1209.0339].
  * (78) J. Bosma, M. Sogaard, and Y. Zhang, Maximal Cuts in Arbitrary Dimension, JHEP 08 (2017) 051, [arXiv:1704.04255].
  * (79) A. Primo and L. Tancredi, Maximal cuts and differential equations for Feynman integrals. An application to the three-loop massive banana graph, Nucl. Phys. B921 (2017) 316–356, [arXiv:1704.05465].
  * (80) M. Harley, F. Moriello, and R. M. Schabinger, Baikov-Lee Representations Of Cut Feynman Integrals, JHEP 06 (2017) 049, [arXiv:1705.03478].
  * (81) A. V. Smirnov and A. V. Petukhov, The Number of Master Integrals is Finite, Lett. Math. Phys. 97 (2011) 37–44, [arXiv:1004.4199].
  * (82) R. N. Lee and A. A. Pomeransky, Critical points and number of master integrals, JHEP 11 (2013) 165, [arXiv:1308.6676].
  * (83) T. Bitoun, C. Bogner, R. P. Klausen, and E. Panzer, Feynman integral relations from parametric annihilators, Lett. Math. Phys. 109 (2019), no. 3 497–564, [arXiv:1712.09215].
  * (84) P. Mastrolia and S. Mizera, Feynman Integrals and Intersection Theory, JHEP 02 (2019) 139, [arXiv:1810.03818].
  * (85) H. Frellesvig, F. Gasparotto, S. Laporta, M. K. Mandal, P. Mastrolia, L. Mattiazzi, and S. Mizera, Decomposition of Feynman Integrals on the Maximal Cut by Intersection Numbers, JHEP 05 (2019) 153, [arXiv:1901.11510].
  * (86) H. Frellesvig, F. Gasparotto, M. K. Mandal, P. Mastrolia, L. Mattiazzi, and S. Mizera, Vector Space of Feynman Integrals and Multivariate Intersection Numbers, Phys. Rev. Lett. 123 (2019), no. 20 201602, [arXiv:1907.02000].
  * (87) H. Frellesvig, F. Gasparotto, S. Laporta, M. K. Mandal, P. Mastrolia, L. Mattiazzi, and S. Mizera, Decomposition of Feynman Integrals by Multivariate Intersection Numbers, JHEP 03 (2021) 027, [arXiv:2008.04823].
  * (88) S. Mizera and A. Pokraka, From Infinity to Four Dimensions: Higher Residue Pairings and Feynman Integrals, arXiv:1910.11852.
  * (89) S. Weinzierl, On the computation of intersection numbers for twisted cocycles, J. Math. Phys. 62 (2021), no. 7 072301, [arXiv:2002.01930].
  * (90) S. Caron-Huot and A. Pokraka, Duals of Feynman Integrals, I: Differential Equations, arXiv:2104.06898.
  * (91) S. Caron-Huot and A. Pokraka, Duals of Feynman Integrals, 2: Generalized Unitarity, arXiv:2112.00055.
  * (92) P. A. Baikov, Explicit solutions of the multiloop integral recurrence relations and its application, Nucl. Instrum. Meth. A389 (1997) 347–349, [hep-ph/9611449].
  * (93) R. N. Lee, Calculating multiloop integrals using dimensional recurrence relation and $D$-analyticity, Nucl. Phys. Proc. Suppl. 205-206 (2010) 135–140, [arXiv:1007.2256].
  * (94) A. G. Grozin, Integration by parts: An Introduction, Int. J. Mod. Phys. A26 (2011) 2807–2854, [arXiv:1104.3993].
  * (95) H. Frellesvig and L. Mattiazzi, On the Application of Intersection Theory to Feynman Integrals: The Univariate Case, in MathemAmplitudes 2019: Intersection Theory and Feynman Integrals, 2, 2021.  arXiv:2102.01576.
  * (96) C. Duhr, Mathematical aspects of scattering amplitudes, in Theoretical Advanced Study Institute in Elementary Particle Physics: Journeys Through the Precision Frontier: Amplitudes for Colliders, pp. 419–476, 2015.  arXiv:1411.7538.
  * (97) J. Bosma, K. J. Larsen, and Y. Zhang, Differential equations for loop integrals in Baikov representation, Phys. Rev. D 97 (2018), no. 10 105014, [arXiv:1712.03760].
  * (98) G. Passarino and M. J. G. Veltman, One Loop Corrections for e+ e- Annihilation Into mu+ mu- in the Weinberg Model, Nucl. Phys. B160 (1979) 151–207.
  * (99) J. Chen, X. Jiang, C. Ma, X. Xu, and L. L. Yang, Baikov representations, intersection theory, and canonical Feynman integrals, arXiv:2202.08127.
  * (100) V. A. Smirnov, Analytical result for dimensionally regularized massless on shell double box, Phys. Lett. B460 (1999) 397–404, [hep-ph/9905323].
  * (101) V. A. Smirnov and O. L. Veretin, Analytical results for dimensionally regularized massless on-shell double boxes with arbitrary indices and numerators, Nucl. Phys. B566 (2000) 469–485, [hep-ph/9907385].
  * (102) M. Czakon and A. Mitov, Inclusive Heavy Flavor Hadroproduction in NLO QCD: The Exact Analytic Result, Nucl. Phys. B 824 (2010) 111–135, [arXiv:0811.4119].
  * (103) A. B. Goncharov, Multiple polylogarithms, cyclotomy and modular complexes, Math. Res. Lett. 5 (1998) 497–516, [arXiv:1105.2076].
  * (104) F. Brown and C. Duhr, A double integral of dlog forms which is not polylogarithmic, in MathemAmplitudes 2019, 6, 2020.  arXiv:2006.09413.
  * (105) R. Boughezal, M. Czakon, and T. Schutzmeier, NNLO fermionic corrections to the charm quark mass dependent matrix elements in $\bar{B}\to X_{s}\gamma$, JHEP 09 (2007) 072, [arXiv:0707.3090].
  * (106) M. Czakon and T. Schutzmeier, Double fermionic contributions to the heavy-quark vacuum polarization, JHEP 07 (2008) 001, [arXiv:0712.2762].
  * (107) U. Aglietti, R. Bonciani, L. Grassi, and E. Remiddi, The Two loop crossed ladder vertex diagram with two massive exchanges, Nucl. Phys. B 789 (2008) 45–83, [arXiv:0705.2616].
  * (108) M. K. Mandal and X. Zhao, Evaluating multi-loop Feynman integrals numerically through differential equations, JHEP 03 (2019) 190, [arXiv:1812.03060].
  * (109) F. Moriello, Generalised power series expansions for the elliptic planar families of Higgs + jet production at two loops, JHEP 01 (2020) 150, [arXiv:1907.13234].
  * (110) M. Hidding, DiffExp, a Mathematica package for computing Feynman integrals in terms of one-dimensional series expansions, Comput. Phys. Commun. 269 (2021) 108125, [arXiv:2006.05510].
  * (111) A. B. Goncharov, M. Spradlin, C. Vergu, and A. Volovich, Classical Polylogarithms for Amplitudes and Wilson Loops, Phys. Rev. Lett. 105 (2010) 151605, [arXiv:1006.5703].
  * (112) C. Duhr, Hopf algebras, coproducts and symbols: an application to Higgs boson amplitudes, JHEP 08 (2012) 043, [arXiv:1203.0454].
  * (113) J. Broedel, C. Duhr, F. Dulat, B. Penante, and L. Tancredi, Elliptic symbol calculus: from elliptic polylogarithms to iterated integrals of Eisenstein series, JHEP 08 (2018) 014, [arXiv:1803.10256].
  * (114) E. Herrmann and J. Parra-Martinez, Logarithmic forms and differential equations for Feynman integrals, JHEP 02 (2020) 099, [arXiv:1909.04777].