††institutetext:  aInstitute for Theoretical Physics, University of
Regensburg, 93040 Regensburg, Germany
bBethe Center for Theoretical Physics, Universitaet Bonn, 53115 Bonn, Germany
cIndian Institute of Science Education and Research Mohali,
Knowledge City, Sector 81, SAS Nagar, Manauli, Punjab 140306, India

# Two-loop non-planar four-point topology with massive internal loop

Taushif Ahmeda, Ekta Chaubeyb, Mandeep Kaurc and Sara Maggiob
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We study a set of two-loop non-planar master integrals needed for the NNLO QCD
corrections to diphoton and dijet production at hadron colliders. The top-
sector topology contains an internal massive fermion loop and is known to
contain elliptic curves. Leveraging the method of differential equations, we
provide a comprehensive discussion for deriving an $\epsilon$-factorized
differential equation related to the most intricate sector within the Feynman
integral family. Despite the dependence on multiple scales and the presence of
two elliptic sectors, we demonstrate how to leverage the properties of their
maximal cuts and the factorization of the Picard-Fuchs operator to deal with
the complexity of the analytic computation. In particular, we construct a
transformation matrix that brings the differential equations into a format
enabling the convenient expression of analytic results in terms of Chen’s
iterated integrals.

##  1 Introduction

Over the past few decades, an enormous amount of progress has been made in
under- standing Feynman integrals and their corresponding functional spaces.
This progress is necessitated by the growing need for precise theoretical
calculations to effectively confront the experimental data collected at the
present and potential future colliders. With the upgraded Large Hadron
Collider (LHC) operating at higher luminosity and the promise of even higher
energy reach for future colliders, the contributions of massive quark loops to
scattering processes have become increasingly significant. The multi-loop
amplitudes governing these scattering processes can generally be expressed in
terms of an independent set of integrals called master integrals (MIs). In the
realm of perturbative Quantum Field Theory (pQFT), obtaining full analytical
results for processes involving heavy particles as mediators have been limited
to initial orders due to the non-availability of the results for involved
complicated integrals. In this article, we address such a scenario by focusing
on the analytical computation of a specific class of two-loop Feynman
integrals crucial for higher-order corrections to two very important processes
at the LHC. Specifically, we present the analytic computation of the two-loop
Feynman integrals associated with a multi-scale non-planar topology, which
appears in NNLO QCD corrections to the diphoton and dijet production
processes.

Owning to the relevance of diphoton and dijet production processes in searches
for physics beyond the Standard Model (SM) ATLAS:2015nsi ; CMS:2016gsl ;
CMS:2018dqv  and in testing the SM predictions, obtaining analytic expressions
for their production cross-sections is crucial. These expressions are
essential for improving prediction accuracy and enabling more accurate
comparisons with the experimental data. However, the presence of non-trivial
massive Feynman integrals with dependence on multiple kinematic scales pose
significant challenges in obtaining the analytic expressions for their
production cross-sections at the two-loop level. The state-of-the-art studies
of analytic computation of two-loop integrals contributing to NNLO QCD
corrections to the dijet and diphoton production processes are as follows. The
analytic results for the two-loop massive planar integrals have been known for
quite some time Caron-Huot:2014lda ; Becchetti:2017abb ; H:2023wfg ;
Becchetti:2023wev . These results are obtained in terms of Chen’s iterated
integrals and, whenever possible, expressed in terms of Goncharov’s
polylogarithms (GPLs). The integrals belonging to a non-planar subtopology
keeping the full dependence on the internal top-quark mass have been computed
in refs. vonManteuffel:2017hms  and are identified to be involving an elliptic
curve. Other subtopologies in the context of Higgs boson have been calculated
in refs. Aglietti:2006tp ; Anastasiou:2006hc . The lack of analytic results
for one of the non-planar topologies, as shown in figure. 1, prevents us from
getting the full two-loop amplitude for dijet and diphoton production in a
closed analytic form. In refs. Xu:2018eos ; Wang:2020nnr , the integrals
belonging to this topology were shown to be elliptic and are evaluated
numerically. Recently, the results for all non-planar integrals contributing
to the two-loop form factors for diphoton production have been computed using
semi-analytical methods Hidding:2020ytt ; Becchetti:2023wev . Despite all
these efforts, full analytic results for the non-planar elliptic integral
family are still missing from the literature.

In this article, we study all the two-loop elliptic sector integrals belonging
to the non-planar topology shown in figure. 1, keeping the full dependence on
the top-quark mass ($m_{t}$) running in the loop with the aim of obtaining
their analytic results. In the first step, we perform the integration-by-parts
(IBP) TKACHOV198165 ; CHETYRKIN1981159  reduction for the given topology to
get the minimal set of scalar integrals called master integrals (MIs) using
KIRA (Maierhofer:2017gsa, ) and FIRE (Smirnov:2014hma, ). Then, we set up the
system of differential equations for these master integrals (Kotikov:1990kg, ;
Kotikov:2003tc, ; Remiddi:1997ny, ; GEHRMANN2000485, ; Argeri:2007up, ) and
aim to bring it into an $\epsilon$-factorized form (Henn:2013pwa, ;
Henn:2014qga, ) by suitable basis choice, where the dimensional regulator
$\epsilon$ is related to space-time dimension through $d=4-2\epsilon$. Since
the subsector integrals are well-documented in the literature, obtaining an
$\epsilon$-factorized system of differential equations for them is
straightforward. Additionally, these sub-sector integrals can be expressed
either in terms of multiple polylogarithms (MPLs) Goncharov:1995tdt ;
Goncharov:2001iea ; Remiddi:1999ew  or elliptic multiple polylogarithms
Broedel:2019hyg . However, for the top-sector integrals, characterized by the
presence of an elliptic curve, the results cannot be expressed in terms of
MPLs. The considered topology is particularly intricate because it not only
depends on multiple kinematic scales but also contains an elliptic
subtopology. Despite the current prominence of such integrals in research,
there are no universal algorithms applicable to every case, and each new
topology requires specific attention (Adams:2013nia, ; Adams:2014vja, ;
Adams:2018kez, ; Duhr:2021fhk, ; Laporta:2004rb, ; Muller-Stach:2012tgz, ;
Bloch:2013tra, ; Adams:2015gva, ; Adams:2015ydq, ; Adams:2016xah, ;
Adams:2017tga, ; Bogner:2017vim, ; Broedel:2017kkb, ; Adams:2018yfj, ;
Broedel:2018iwv, ; Adams:2018bsn, ; Broedel:2018qkq, ; Honemann:2018mrb, ;
Bogner:2019lfa, ; Broedel:2019hyg, ; Broedel:2019kmn, ; Weinzierl:2019pfw, ;
Adams:2017ejb, ; Marzucca:2023gto, ; Giroux:2024yxu, ; Gorges:2023zgv, ;
Delto:2023kqv, ). To cast the system of differential equations into the
$\epsilon$-form for the top-sector involving an elliptic curve, we have
generalised the technique discussed in (Pogel:2022vat, ) for the basis
transformation. Finally, we discuss how to use our choice of basis to express
the analytic results of the top-sector master integrals at each order in
$\epsilon$ through Chen’s iterated integrals (Chen:1977oja, ) over elliptic
one-forms.

The paper is organized as follows: In section 2, we introduce our notation for
the kinematics and integral family, as well as discuss our choice of master
integral basis for the topology under consideration. Section 3 is devoted to a
detailed discussion on establishing the $\epsilon$-form of the differential
equation system for the most intricate elliptic top-sector, utilizing the
maximal cut and factorization of Picard-Fuchs operator information. Finally,
we discuss how to express all the integrals belonging to this Feynman integral
family through Chen’s iterated integrals.

##  2 The setup

We consider a generic four-point scattering process involving massless
external particles with momenta $p_{1},\cdots,p_{4}$, all of them we regard as
incoming. The external momenta obey the following kinematics:
$\sum_{i=1}^{4}p_{i}=0,\;p_{i}^{2}=0\;\text{for}\;i=1,2,3,4$. We define the
following abbreviations to represent the independent set of Mandelstam
variables

$\displaystyle s=(p_{1}+p_{2})^{2},\qquad t=(p_{1}+p_{3})^{2}.$ (1)

Let us define a set of two-loop non-planar Feynman integrals through

$\displaystyle
I_{a_{1},\cdots,a_{9}}=\left(\frac{e^{\epsilon\gamma_{E}}}{i\pi^{\frac{d}{2}}}\right)^{2}\int\prod_{i=1}^{2}d^{d}k_{i}\frac{D_{8}^{a_{8}}D_{9}^{a_{9}}}{D_{1}^{a_{1}}D_{2}^{a_{2}}D_{3}^{a_{3}}D_{4}^{a_{4}}D_{5}^{a_{5}}D_{6}^{a_{6}}D_{7}^{a_{7}}}\,,\quad
a_{j}\in\mathbb{Z}\,.$ (2)

We represent the loop momenta through $k_{i}$, Euler-Mascheroni constant by
$\gamma_{E}$ and the inverse propagators by $D_{j}=q_{j}^{2}-m_{j}^{2}+i0^{+}$
with the parameters $q_{j}$ and $m_{j}$ respectively denoting the momentum and
mass. In this article, we focus on a set of Feynman integrals associated with
the topology shown in figure 1.

Figure 1: The non-planar topology where the red and black lines respectively
denote massive and massless propagators.

The corresponding integral family is constructed out of the following set of
inverse propagators

$\displaystyle\big{\\{}$ $\displaystyle
k_{1}^{2},(k_{1}-p_{1})^{2},(k_{1}-p_{1}-p_{2})^{2},k_{2}^{2}-m_{t}^{2},(k_{2}-p_{1}-p_{2}-p_{3})^{2}-m_{t}^{2},$
$\displaystyle(k_{1}-k_{2})^{2}-m_{t}^{2},(k_{1}-k_{2}+p_{3})^{2}-m_{t}^{2},(k_{1}-p_{1}-p_{2}-p_{3})^{2},(k_{2}+p_{1})^{2}-m_{t}^{2}\big{\\}}\,,$
(3)

where $m_{t}$ denotes the mass of top quark appearing in the loop.

###  2.1 Choice of master integrals

By using the integration-by-parts TKACHOV198165 ; CHETYRKIN1981159  identities
for the given topology with the help of KIRA (Maierhofer:2017gsa, ) and FIRE
(Smirnov:2014hma, ), we obtain a basis of 36 master integrals denoted by
$\vec{I}$. We then set up a system of differential equations for these master
integrals with respect to the Mandelstam variables defined in equation (1).
Utilizing IBP relations, the differential equations take the form:

$d\vec{I}=A\vec{I},$ (4)

where $A$ is a $36\times 36$ matrix-valued one-form depending on the
Mandelstam variables and the space-time dimensions ($d=4-2\epsilon$).

Since the basis choice for the master integrals is not unique, we can perform
a basis transformation

$\vec{{J}}=U\vec{I},$ (5)

which leads to a transformed system of differential equations

$d\vec{{J}}=\tilde{A}\vec{{J}},$ (6)

where,

$\tilde{A}=UA\;U^{-1}-UdU^{-1}.$ (7)

If we can find a transformation matrix $U$ such that the transformed
differential equation attains the form

$d\vec{{J}}=\epsilon\tilde{A}\vec{{J}},$ (8)

such that the matrix-valued one-form $\tilde{A}$ does not have any
$\epsilon$-dependence, then the analytic result of $\vec{{J}}$ can be
systematically expressed in terms of Chen’s iterated integrals order by order
in $\epsilon$. Additionally, each entry at a given $\epsilon$ order will have
a uniform transcendental weight. We may slightly relax this condition of
complete $\epsilon$-factorization and have the form of $\tilde{A}$ as

$\tilde{A}=\tilde{A}^{(0)}+\epsilon\tilde{A}^{(1)},$ (9)

where $\tilde{A}^{(0)}$ is strictly lower triangular and both
$\tilde{A}^{(0)}$ and $\tilde{A}^{(1)}$ are independent of $\epsilon$
Adams:2018kez . This form of the transformed matrix $\tilde{A}$ maintains the
property that we can still express the results for $\vec{J}$ in terms of
iterated integrals at every order in $\epsilon$, although there might be a
mixing of functions with different length at a particular order in $\epsilon$.
Here length is defined as the number of integrations that needs to be
performed iteratively. With the modern technology of using differential
equations for computing Feynman integrals, often the bottleneck is the
construction of a transformation matrix that brings our initial “not so good
set of differential equations” to a better form, such as the canonical form in
equation (8) or the special linear form in equation (9). In this work, we
highlight how to obtain both these forms for the top-sector.

Instead of working with basis $\vec{I}$ in equation (4), guided by the works
in (Caron-Huot:2014lda, ; Becchetti:2017abb, ; H:2023wfg, ;
vonManteuffel:2017hms, ), we choose the following definitions of the master
integrals to set up our differential equations, where we have used
$m_{t}^{2}=1$ for convenience:

$\displaystyle J_{1}\;$
$\displaystyle=\;\epsilon^{2}\;\;\textbf{D}^{-}I_{000010100},$ $\displaystyle
J_{2}\;$ $\displaystyle=\;\epsilon^{2}\;\;\textbf{D}^{-}I_{101000100},$
$\displaystyle J_{3}\;$
$\displaystyle=\;\epsilon^{2}\;\sqrt{(-4+s)\;s}\left[I_{002201000}+\frac{1}{2}I_{001202000}\right],$
$\displaystyle J_{4}\;$
$\displaystyle=\;\epsilon^{2}\;\sqrt{s+t}\sqrt{4+s+t}\left[I_{020021000}+\frac{1}{2}I_{010022000}\right],$
$\displaystyle J_{5}\;$
$\displaystyle=\;\epsilon^{2}\;\sqrt{(4-t)\;t}\left[I_{020200100}+\frac{1}{2}I_{010200200}\right],$
$\displaystyle J_{6}\;$
$\displaystyle=\;\frac{1}{2}\epsilon^{2}s\left[I_{001202000}+2(I_{002201000}+\textbf{D}^{-}I_{001101100})\right],$
$\displaystyle J_{7}\;$ $\displaystyle=\;\epsilon^{3}\;s\;I_{002101100},$
$\displaystyle J_{8}\;$
$\displaystyle=\;\epsilon^{2}\;\;(s+t)\;\left[\frac{1}{2}I_{010022000}+I_{020021000}+\textbf{D}^{-}I_{010011100}\right],$
$\displaystyle J_{9}\;$ $\displaystyle=\;\epsilon^{3}\;(s+t)\;I_{020011100},$
$\displaystyle J_{10}\;$
$\displaystyle=\;\frac{1}{2}\;\epsilon^{2}\;t\;\left[(I_{010200200}+2(I_{020200100}+\;\textbf{D}^{-}I_{010101100})\right],$
$\displaystyle J_{11}\;$ $\displaystyle=\;\epsilon^{3}\;t\;I_{020101100},$
$\displaystyle J_{12}\;$ $\displaystyle=\;\epsilon^{4}\;s\;I_{001111100},$
$\displaystyle J_{13}\;$ $\displaystyle=\;\epsilon^{4}\;s\;I_{101011100},$
$\displaystyle J_{14}\;$ $\displaystyle=\;\epsilon^{3}\;s\;I_{101012000},$
$\displaystyle J_{15}\;$
$\displaystyle=\;\epsilon^{2}\;\sqrt{-s(4+s)}\;\left[-s\;I_{201012000}+\epsilon\frac{\textbf{D}^{-}I_{000010100}}{2+4\epsilon}\right],$
$\displaystyle J_{16}\;$ $\displaystyle=\;\epsilon^{4}\;(s+t)\;I_{011101100},$
$\displaystyle J_{17}\;$
$\displaystyle=-\epsilon^{3}\;\sqrt{s\;t\;(-4\;s+s\;t-4\;t)}\;I_{011201100},$
$\displaystyle J_{18}\;$ $\displaystyle=\;\epsilon^{4}\;t\;I_{011111000},$
$\displaystyle J_{19}\;$
$\displaystyle=-\epsilon^{3}\;\sqrt{s\;(s+t)\;(s^{2}-4\;t+s\;t)}\;I_{011112000},$
$\displaystyle J_{20}\;$ $\displaystyle=\;\epsilon^{4}\;s\;I_{010111100},$
$\displaystyle J_{21}\;$
$\displaystyle=\;\epsilon^{3}\;\sqrt{t\;(s+t)\;(t^{2}-4\;s+s\;t)}\;I_{020111100},$
$\displaystyle J_{22}\;$
$\displaystyle=\;\epsilon^{3}\;s\;(I_{010121100}+I_{010211100}),$
$\displaystyle J_{23}\;$
$\displaystyle=\;\epsilon^{4}\;s\;(s+t)\;I_{111011100},$ $\displaystyle
J_{24}\;$ $\displaystyle=\;\epsilon^{4}\;s\;t\;I_{111101100},$ $\displaystyle
J_{25}\;$ $\displaystyle=\;\epsilon^{3}\;s\;\sqrt{(-4+t)\;t}\;I_{111100200},$
$\displaystyle J_{26}\;$
$\displaystyle=\;\epsilon^{3}\;s\;((-1+2\epsilon)\;I_{111100100}-(-4+t)\;I_{111100200}),$
$\displaystyle J_{27}\;$
$\displaystyle=\;\epsilon^{4}\;\sqrt{s}\;\sqrt{t}\;\sqrt{s+t}\;I_{011111100},$
$\displaystyle J_{28}\;$
$\displaystyle=\;\frac{1}{4}\;\epsilon^{3}\;\left[s\;(s+t)\;(I_{011112000}+4I_{011112100})+s\;t\;I_{011201100}-t\;(s+t)\;I_{020111100}\right],$
$\displaystyle J_{29}\;$
$\displaystyle=\;\epsilon^{3}\;s\;\sqrt{s+t}\;\sqrt{4+s+t}\;I_{111012000},$
$\displaystyle J_{30}\;$
$\displaystyle=\;\epsilon^{3}\;s\left[(-1+2\epsilon)\;I_{111011000}+(4+s+t)\;I_{111012000}\right],$
$\displaystyle J_{31}\;$
$\displaystyle=-\frac{\pi\;s^{2}\;\epsilon^{4}}{2\;\text{K}\big{(}\frac{-16}{s}\big{)}}I_{101111100},$
$\displaystyle J_{32}\;$
$\displaystyle=\frac{\epsilon^{3}s}{\pi}\bigg{[}(1+2\epsilon)\;\bigg{(}s\;\text{E}\big{(}-\frac{16}{s}\big{)}-(16+s)\;\text{K}\big{(}-\frac{16}{s}\big{)}\bigg{)}I_{1,0,1,1,1,1,1,0,0}$
$\displaystyle\;-8\;(16+s)\;\text{K}\big{(}-\frac{16}{s}\big{)}I_{1,0,1,2,1,1,1,0,0}\bigg{]}$
$\displaystyle J_{33}\;$ $\displaystyle=\;\epsilon^{4}\;I_{111111100},$
$\displaystyle J_{34}\;$ $\displaystyle=\;\epsilon^{4}\;I_{11111110-2},$
$\displaystyle J_{35}\;$ $\displaystyle=\;\epsilon^{4}\;I_{11111110-1},$
$\displaystyle J_{36}\;$ $\displaystyle=\;\epsilon^{4}\;I_{1111111-1-1}\,.$
(10)

Here $\textbf{D}^{-}$ denotes the dimension shift operator which shifts the
space-time dimension $d$ by two as

$\displaystyle\textbf{D}^{-}I_{a_{1},\cdots,a_{9}}(d)=I_{a_{1},\cdots,a_{9}}(d-2),$
(11)

whereas $K(\cdot)$ and $E(\cdot)$ are respectively the complete elliptic
integral of the first and second kind,

$K(k^{2})=\int_{0}^{1}\frac{dt}{\sqrt{(1-t^{2})(1-k^{2}t^{2})}},\qquad
E(k^{2})=\int_{0}^{1}\sqrt{\frac{1-k^{2}t^{2}}{1-t^{2}}}dt.$ (12)

With this basis choice, the differential equations readily assume a canonical
form, as in equation (8) for the first 31 integrals, and for $J_{32}$ it takes
on a special linear form, as in equation (9), that can be effortlessly
transformed into the canonical one through a simple transformation involving
the primitive, as outlined in Adams:2017tga . Moreover, the integrals
$J_{33}$-$J_{36}$ in the top sector, featuring all seven propagators, along
with $J_{31}$-$J_{32}$, constitute elliptic sectors. In the following
sections, we explain the main ingredients of our work, tackling the analytic
complexity of the top sector due to the elliptic nature of the Feynman
integral.

##  3 The elliptic top-sector

As previously mentioned, all sub-sector integrals associated with the non-
planar family 1 have been known analytically for quite some time.
Nevertheless, writing the analytic results for the full non-planar topology is
challenging due to the presence of elliptic curves, which extend beyond the
well-understood class of polylogarithmic functions. The presence of multiple
kinematic scales adds further to the analytic complexity. This makes the
analytic computation of the entire integral family considerably more
challenging compared to many other multi-scale two-loop Feynman integrals.
Therefore, in the following sections, we focus primarily on the top sector,
i.e. integrals $J_{33}$-$J_{36}$ and provide a highly detailed explanation of
the steps involved in transforming the corresponding differential equations in
a way that facilitates the straightforward articulation of analytic results.
The ideas presented here can be easily generalized to other multi-variate
elliptic Feynman integrals involving multiple elliptic sectors.

###  3.1 The study of maximal cuts

Maximal cuts are one of the primary tools that can be utilized to analyse the
algebraic complexity of multi-loop, multi-scale Feynman integrals before
solving them explicitly BAIKOV1997347 ; Primo:2017ipr ; Frellesvig:2017aai .
In this context, ‘cutting’ a propagator means forcing the propagating particle
to be on-shell. This can be achieved for the propagators raised to power one
by simply replacing each of them with a Dirac $\delta$-function. A maximal cut
of a diagram corresponds to the simultaneous cutting of all its propagators.
Mathematically, a maximal cut corresponds to taking the $n$-fold residue of
the corresponding integrand in the complex plane, where $n$ is the number of
propagators. Furthermore, it is known that the maximal cuts of any given set
of master integrals are the solutions for the corresponding homogeneous system
of differential equations Primo:2016ebd ; Adams:2018kez . Therefore, the
maximal cut information of a given Feynman integral is useful to obtain the
integral representation for the homogeneous solution of the differential
equation satisfied by that integral. The computation of the maximal cut can be
easily carried out in the so-called Baikov representation Baikov:1996iu  using
the loop-by-loop approach Frellesvig:2017aai .

The maximal cut for the top-sector integral $I_{111111100}$, as shown in
figure. 2(b), is given by

$\frac{16}{\pi^{4}}\int_{\mathcal{C}}\frac{dP}{s\;(s+t+P)\;\sqrt{P}\;\sqrt{s+P}\;\sqrt{-4m_{t}^{2}s+sP+P^{2}}}+O(\epsilon)$
(13)

in $d$-dimensions, where the leading term in $\epsilon$ corresponds to the
4-dimensional part. The integration is over the contour $\mathcal{C}$ in the
remaining Baikov variable $P$, after incorporating all the delta distributions
that arise from putting all the propagators on-shell. The contour can be
chosen to lie between any pair of roots of the polynomial in the denominator.

(a) (a)

(b) (b)

Figure 2: The non-planar sub-topology and the top-topology where the red and
black lines respectively denote massive and massless propagators.

On the other hand, the maximal cut of the elliptic sub-sector integral
$I_{101111100}$, as shown in figure. 2(a), is given by

$\frac{8}{\pi^{3}}\int_{\mathcal{C}}\frac{dP}{s\;\sqrt{P}\;\sqrt{s+P}\;\sqrt{-4m_{t}^{2}s+sP+P^{2}}}+O(\epsilon).$
(14)

Looking at the equations (13) and (14), it becomes evident that both of them
contain the same quartic square root, while the $t$-dependence factors out in
equation (13). A square root $\sqrt{Q(x)}$, where $Q(x)$ is a quartic
polynomial with distinct zeroes, defines an elliptic curve. The occurrence of
the same quartic square roots strongly suggests the appearance of a single
elliptic curve in sectors $J_{31}$-$J_{36}$, instead of two different elliptic
curves, a fact that has been overlooked in earlier studies Becchetti:2023wev ;
Xu:2018eos . This information is useful for the construction of an
$\epsilon$-factorised differential equation for the top-sector. Particularly,
we find that the $t$ dependence completely decouples after using the
factorization properties of the corresponding Picard-Fuchs operator Muller-
Stach:2012tgj ; Adams:2017tga  from the sector that is still coupled at the
order $\epsilon^{0}$. We discuss it in detail in the next sections.

###  3.2 Transformation matrix to obtain an $\epsilon$-form for the top-sector
on the maximal cut

In this section, we first discuss the steps to obtain the special linear form
(9) for the differential equations on the maximal cut corresponding to the
elliptic top sector integrals $(J_{33},..,J_{36})$ defined in equation (2.1).
Later on, we comment on obtaining the $\epsilon$-factorized canonical form (8)
from the special linear form (9) for these integrals. Following Adams:2017tga
, we first act on the differential equation for this $4\times 4$ block with a
transformation matrix that decouples the elliptic curve at the order
$\epsilon^{0}$. As a side product, this procedure already transforms the first
2$\times$2 block for the partial differential equation with respect to $t$ in
a $\epsilon$-factorized form. In particular, we use the transformation matrix

$U_{0}=\small{\left(\begin{array}[]{cccc}0&-\frac{2s}{(16+s)(s+2t)}&\frac{16t+2s(4+s+t)}{(16+s)(s+2t)}&\frac{2s}{(16+s)(s+2t)}\\\
0&\frac{8s(12+s)}{(16+s)^{2}(s+2t)}&\frac{-4s(112+s(33+2s))-8(112+s(21+s))t}{(16+s)^{2}(s+2t)}&-\frac{8s(12+s)}{(16+s)^{2}(s+2t)}\\\
s\sqrt{t}\sqrt{s+t}\sqrt{s(-4+t)+t^{2}}&0&0&0\\\ 0&0&0&s\end{array}\right)}$
(15)

to get

$\tilde{A}^{(0)}_{s,U_{0}}=\begin{pmatrix}\frac{2}{s}&\frac{1}{s}&0&0\\\
-\frac{4(49+5s)}{s(16+s)}&\frac{-80-7s}{s(16+s)}&0&0\\\
\frac{\sqrt{t}(-112s-8s^{2}+27st+2s^{2}t+26t^{2}+2st^{2})}{\sqrt{s+t}\sqrt{-4s+st+t^{2}}}&\frac{\sqrt{t}(-64s-4s^{2}+16st+s^{2}t+16t^{2}+st^{2})}{2\sqrt{s+t}\sqrt{-4s+st+t^{2}}}&0&0\\\
-\frac{s}{2}&0&0&0\end{pmatrix},$ (16)
$\tilde{A}^{(0)}_{t,U_{0}}=\begin{pmatrix}0&0&0&0\\\ 0&0&0&0\\\
-\frac{s\left(2s^{2}t-9s^{2}+2st^{2}+24st-112s+24t^{2}\right)}{\sqrt{t}\sqrt{s+t}\sqrt{st-4s+t^{2}}}&-\frac{s(s+16)\sqrt{st-4s+t^{2}}}{2\sqrt{t}\sqrt{s+t}}&0&0\\\
0&0&0&0\end{pmatrix},$ (17)

where $\tilde{A}^{(0)}_{s,U_{0}}$ and $\tilde{A}^{(0)}_{t,U_{0}}$ are the
transformed differential equation matrices with respect to $s$ and $t$,
respectively, obtained after transforming with $U_{0}$. The elliptic nature of
the Feynman integrals becomes evident when we look at the upper left $2\times
2$ block in $\tilde{A}^{(0)}_{s,U_{0}}$, which is still coupled at
$\epsilon^{0}$. Indeed, if we compute the Picard-Fuchs operator at the order
$\epsilon^{0}$ denoted by $L_{0}$ for the first element from the transformed
basis, i.e. the combination of master integrals that we obtain after acting
with the transformation matrix $U_{0}$, which for our case is

$\displaystyle I_{1}\;$
$\displaystyle=\frac{2\epsilon^{4}}{(16+s)(s+2t)}\;\big{(}sI_{1,1,1,1,1,1,1,-1,-1}\;-sI_{1,1,1,1,1,1,1,0,-2}\;+$
$\displaystyle\hskip
113.81102pt(4s+s^{2}+8t+st)I_{1,1,1,1,1,1,1,0,-1}\big{)},$ (18)

we obtain

$L_{0}=\partial^{2}_{s}-\left(-\frac{4}{s}-\frac{2}{16+s}\right)\partial_{s}+\frac{6(6+s)}{s^{2}(16+s)}.$
(19)

This second order differential operator is associated to an elliptic curve and
its two solutions can be chosen as

$\displaystyle\psi_{0}(s)=\frac{32\;E(-\frac{s}{16})}{\pi s^{3/2}(s+16)},$
(20)
$\displaystyle\psi_{1}(s)=\frac{32\;\left(E(\frac{s}{16}+1)-K(\frac{s}{16}+1)\right)}{s^{3/2}(s+16)},$
(21)

where $\psi_{0}(s)$ is a holomorphic function multiplied by
$\frac{1}{s^{3/2}}$ and $\psi_{1}(s)$ is the first logarithmic solution. Now
we extend the method of Pogel:2022vat  with minor modifications to multiple
scales, in order to factorize $\epsilon$ from this $2\times 2$ block. In
particular, we look for new master integrals defined in the following way

$\displaystyle M_{1}=\frac{1}{h(s,t)}I_{1},$ $\displaystyle
M_{2}=\frac{g_{1}(s,t)}{\epsilon}\frac{d}{ds}M_{1}+\frac{g_{2}(s,t)}{\epsilon}\frac{d}{dt}M_{1}-f_{1}(s,t)M_{1}.$
(22)

The above is the most general ansatz with four rational functions:
$h(s,t),\;g_{1}(s,t),\;g_{2}(s,t),$ ${\text{and}}\;f_{1}(s,t)$, that one can
write thinking about the Hodge filtration of a Hodge structure Bonisch:2021yfw
. With this ansatz, the differential equation matrix $\tilde{A}$ for the first
2 $\times$ 2 block transforms again as in equation (7). We find the unknown
functions in (22) by requiring that the transformed matrix $\tilde{A}$ has
terms proportional to $\epsilon$ and as simplest as possible $\epsilon^{0}$
terms. This is done by choosing

$\displaystyle h(s,t)=\frac{1}{s^{2}+16s},$ (23) $\displaystyle g_{1}(s,t)=0,$
(24) $\displaystyle g_{2}(s,t)=\frac{16\;s\;(s+2t)^{2}}{16\;t+s\;(8+t)},$ (25)
$\displaystyle
f_{1}(s,t)=\frac{16\;s^{3}\;(s+2t)}{t\;(s+t)\;(16\;t+s\;(8+t))}.$ (26)

So the next transformation matrix attains the form

$U_{1}=\begin{pmatrix}s\;(16+s)&0\\\
-64\;s^{2}\;(12+s)&-16\;s^{2}\;(16+s)\end{pmatrix}.$ (27)

Acting with another $2\times 2$ rotation, where the entries of the matrix are
constrained to eliminate the $\epsilon^{0}$ terms, we $\epsilon$-factorize the
$2\times 2$ block. This transformation is given by

$U_{2}=\begin{pmatrix}\sqrt{s}\;K(\frac{s}{16}+1)&\frac{1}{8\sqrt{s}}\bigg{(}E(\frac{s}{16}+1)-K(\frac{s}{16}+1)\bigg{)}\\\
\frac{2\sqrt{s}}{\pi}\;K(-\frac{s}{16})&-\frac{1}{4\pi\sqrt{s}}{E(-\frac{s}{16})}\end{pmatrix}.$
(28)

Here the appearance of elliptic functions is not surprising given that the
associated geometry is elliptic in nature. These entries can indeed be
rewritten in terms of the two periods $\psi_{0}$ and $\psi_{1}$ defined in
equations (20) and (21) and their derivatives in the following way

$U_{2}=\frac{s^{3}(s+16)^{2}}{16}\begin{pmatrix}s\psi^{\prime}_{1}-2\psi_{1}&-\psi_{1}\\\
-2(s\psi^{\prime}_{0}-2\psi_{0})&2\psi_{0}\end{pmatrix}.$ (29)

After fixing the $2\times 2$ sub-block, we proceed to $\epsilon$-factorize the
remaining $4\times 4$ differential matrix. It is important to note that at
this stage, we already obtain the special linear form for our differential
equations, which enables us to conveniently write the analytic results for the
elliptic top-sector (on the cut) in terms of iterated integrals. However, as
mentioned before, the results at a particular $\epsilon$-order do not have the
same length of iterated integrals. Let us now look at how to obtain an
$\epsilon$-form for the full $4\times 4$ matrix. After applying the rotations
$U_{1}$ and $U_{2}$, $\tilde{A}$ transforms accordingly as in equation (7). We
find it useful to examine the new differential equation matrix,
$\tilde{A}_{U_{2}U_{1}}$, based on its scale dependence, i.e. we analyze
$\tilde{A}_{U_{2}U_{1},s}$ and $\tilde{A}_{U_{2}U_{1},t}$ individually. We
find that $\tilde{A}_{U_{2}U_{1},s}$ has $\epsilon^{0}$ terms only at the
positions (3,1), (3,2), (4,1) and (4,2) whereas $\tilde{A}_{U_{2}U_{1},t}$ has
non-zero components only at positions (3,1) and (3,2). With this information,
we construct another transformation matrix as follows:

$U_{3}=\begin{pmatrix}1&0&0&0\\\ 0&1&0&0\\\ 0&0&1&0\\\
-\frac{1}{2}s^{2}(s+12)&-\frac{1}{8}s^{2}(s+16)&0&1\end{pmatrix},$ (30)

where the components are fixed by requiring that the $\epsilon^{0}$ terms at
positions (4,1) and (4,2) cancel each other out. The last transformation
$U_{4}$, to fully $\epsilon$-factorize the $4\times 4$ block of the
differential equation system for the top-sector, is obtained in a similar way
by constructing another non-trivial transformation for the elements at
position (3,1) and (3,2) through

$U_{4}=\begin{pmatrix}1&0&0&0\\\ 0&1&0&0\\\ f(s,t)&g(s,t)&1&0\\\
0&0&0&1\end{pmatrix}.$ (31)

The definitions of $f(s,t)$ and $g(s,t)$ are given in the appendix 5.2.
Furthermore, in appendix 5.3, we highlight the significance of studying the
marked point on the torus to simplify these functions. This is obtained after
establishing a correspondence between the elliptic curve and a torus.

In conclusion, we obtain both a completely $\epsilon$-factorized as well as a
special linear form (9) for the $4\times 4$ differential equation block
corresponding to the elliptic top sector. We also find that the analytic one-
forms are simpler if we only transform the starting block for the top-sector
(the one obtained with $J_{i}$ in equation (2.1)), from $U_{0}$ till $U_{3}$,
rather than proceeding to the final transformation matrix $U_{4}$. Hence this
is the choice made in our supplied transformation files. We also write the
corresponding differential equation matrix $\tilde{A}_{f(i,j)}$, with
$f\in\\{s,t\\}$, for the elliptic top sector explicitly in appendix 5.1. The
corresponding differential one-forms are expressed in terms of the periods and
their derivatives given in equations (20) and (21). All the elements are
proportional to $\epsilon$, apart from $\tilde{A}_{f(3,1)}$ and
$\tilde{A}_{f(3,2)}$.

#### Analytic results through Chen’s iterated integrals

The primary focus of this article revolves around addressing the algebraic
complexity arising from two elliptic sectors within a multi-variate non-planar
Feynman integral family. As evident from the previous section, we made a
comprehensive investigation of all the elliptic sectors to manage this
complexity. However, as a secondary outcome, we also provide a choice of basis
for the full system of differential equations as ancillary materials with this
article. This basis can be utilised to express all the appearing Feynman
integrals analytically through Chen’s iterated integrals, albeit with non-
uniform length, at each order in $\epsilon$. This choice has been obtained as
follows. Typically, following the bottom-top approach, selecting an
appropriate basis for the top sector on the maximal cut, with all the sub-
sector integrals already in a canonical form, also simplifies the off-diagonal
differential equation contributions to the top-sector. This is true especially
concerning the $\epsilon$ expansion behaviour. However, we observe that this
is not the case here. In particular, we additionally fix the
$\epsilon$-behaviour of these off-diagonal elements, to limit the one-forms to
have only finite expansions in $\epsilon$. We provide a choice of basis that
is arranged in a lower block-triangular form which nevertheless enables one to
express the analytic results though iterated integrals, after solving the
differential equations order by order in $\epsilon$. We would also like to
mention that after obtaining the system of differential equations with the
provided choice of basis for the full system, further transformations
Meyer:2017joq  can be carried out, to fix only the sub-sector contributions to
the top-sector, to obtain a completely $\epsilon$-factorised differential
equations for the full $36\times 36$ system. However, in this case, the more
transformations carried out, the more intricate the structure of the kernels
within the iterated integrals becomes.

##  4 Conclusions and outlook

In this paper, we studied the analytic structure of non-planar Feynman
integrals crucial for computing NNLO QCD corrections in diphoton and dijet
production. We presented an extensive discussion on deriving an
$\epsilon$-factorized form of the differential equations for the elliptic top
sector, leveraging the geometric properties of the Feynman integral and
insights obtained from the maximal cuts. Contrary to the earlier studies and
despite the presence of an elliptic sub-topology, we identify only one
elliptic curve in the entire system instead of two distinct elliptic curves.
With this information, we exploit the properties of the Picard-Fuchs operator
to find an $\epsilon$-factorised differential equation of the top-sector on
the maximal cut. Additionally, we adapted the algorithm outlined in
Pogel:2022vat  for our multi-variate elliptic system, enabling the derivation
of appropriate ansätze for the top sector. At the end, we also briefly
discussed how to obtain analytic results for the elliptic top sector integrals
in terms of Chen’s iterated integrals using our choice of basis.

Looking ahead, although significant progress has been achieved in
understanding the analytic structure of this family of Feynman integrals,
further investigation is clearly warranted. Specifically, it would be
intriguing to obtain a simpler definition of a good function basis for these
integrals, incorporating information from both the maximal and sub-maximal
cuts from the beginning. This approach would result in requiring fewer
transformations to obtain the $\epsilon$-factorized differential equation for
the entire system, thereby yielding less complicated kernels in the analytic
results. Furthermore, the numerical evaluation of these elliptic kernels needs
to be explored.

## Acknowledgements

We would like to thank Claude Duhr and Vasily Sotnikov for comments on the
manuscript. EC would also like to thank Christoph Dlapa for interesting
discussions at the workshop ‘MathemAmplitudes 2023: QFT at the Computational
Frontier’. SM would like to thank Yu Jiao Zhu for insightful discussions. The
work of EC and SM is funded by the ERC grant 101043686 ‘LoCoMotive’. The work
of TA is supported by Deutsche Forschungsgemeinschaft (DFG) through the
Research Unit FOR 2926, Next Generation perturbative QCD for Hadron Structure:
Preparing for the Electron-ion collider, project number 409651613. MK would
like to acknowledge financial support from IISER Mohali for this work.

##  5 Appendix

###  5.1 The differential one-forms for the elliptic top sector

The matrix valued differential one-forms for the top-sector mentioned in
section 3.2 has the following form:

$\displaystyle\tilde{A}_{f(1,1)}=$
$\displaystyle\frac{\epsilon}{512}(2\text{$\psi_{0}$}\left(\text{$\psi_{1}$}^{\prime}(\text{ds}\;w_{22}-\text{dt}\;w_{3})-2\text{$\psi_{1}$}(\text{ds}\;w_{21}+\text{dt}\;w_{1})\right)+$
$\displaystyle\text{$\psi_{0}$}^{\prime}\left(\text{$\psi_{1}$}^{\prime}(\text{ds}\;w_{24}-\text{dt}\;w_{2})+2\text{$\psi_{1}$}(\text{ds}\;w_{23}-\text{dt}\;w_{4})\right)+512\text{dt}\;w_{16})\,,$
$\displaystyle\tilde{A}_{f(1,2)}=$
$\displaystyle\frac{\epsilon}{1024}\left(4\text{ds}\;w_{25}-\text{dt}\;w_{5}\text{$\psi_{1}$}\text{$\psi_{1}$}^{\prime}\right)\,,$
$\displaystyle\tilde{A}_{f(1,3)}=$
$\displaystyle\frac{\epsilon}{8}\left(\text{$\psi_{1}$}^{\prime}(\text{ds}\;w_{31}-\text{dt}\;w_{11})+2\text{$\psi_{1}$}(\text{ds}\;w_{26}-\text{dt}\;w_{6})\right)\,,$
$\displaystyle\tilde{A}_{f(1,4)}=$
$\displaystyle-\frac{\epsilon}{4}\left(\text{$\psi_{1}$}^{\prime}(\text{dt}\;w_{12}-\text{ds}\;w_{32})+2\text{$\psi_{1}$}(\text{ds}\;w_{27}+\text{dt}\;w_{7})\right)\,,$
$\displaystyle\tilde{A}_{f(2,1)}=$
$\displaystyle\frac{\epsilon}{256}\left(\text{dt}\;w_{5}\text{$\psi_{0}$}\text{$\psi_{0}$}^{\prime}-4\text{ds}\;w_{25}\right)\,,$
$\displaystyle\tilde{A}_{f(2,2)}=$
$\displaystyle\frac{\epsilon}{512}(2\text{$\psi_{0}$}\left(2\text{$\psi_{1}$}(\text{ds}\;w_{21}+\text{dt}\;w_{1})+\text{$\psi_{1}$}^{\prime}(\text{dt}\;w_{4}-\text{ds}\;w_{23})\right)$
$\displaystyle+\text{$\psi_{0}$}^{\prime}\left(\text{$\psi_{1}$}^{\prime}(\text{dt}\;w_{2}-\text{ds}\;w_{24})+\text{$\psi_{1}$}(2\text{dt}\;w_{3}-2\text{ds}\;w_{22})\right)+512\text{dt}\;w_{16})\,,$
$\displaystyle\tilde{A}_{f(2,3)}=$
$\displaystyle-\frac{\epsilon}{4}\left(\text{$\psi_{0}$}^{\prime}(\text{ds}\;w_{31}-\text{dt}\;w_{11})+2\text{$\psi_{0}$}(\text{ds}\;w_{26}-\text{dt}\;w_{6})\right)\,,$
$\displaystyle\tilde{A}_{f(2,4)}=$
$\displaystyle\frac{\epsilon}{2}\text{$\psi_{0}$}^{\prime}(\text{dt}\;w_{12}-\text{ds}\;w_{32})+\epsilon\text{$\psi_{0}$}(\text{ds}\;w_{27}+\text{dt}\;w_{7})\,,$
$\displaystyle\tilde{A}_{f(3,1)}=$
$\displaystyle\frac{1}{64}\left(2\text{$\psi_{0}$}(\text{ds}\;(2w_{29}+w_{30}\epsilon)-\text{dt}\;(w_{10}\epsilon+w_{9}))+\text{$\psi_{0}$}^{\prime}(2\text{ds}\;w_{34}+\text{ds}\;w_{35}\epsilon-2\text{dt}\;w_{14}-\text{dt}\;w_{15}\epsilon)\right)\,,$
$\displaystyle\tilde{A}_{f(3,2)}=$
$\displaystyle\frac{1}{128}\left(2\text{$\psi_{1}$}(\text{ds}\;(2w_{29}+w_{30}\epsilon)-\text{dt}\;(w_{10}\epsilon+w_{9}))+\text{$\psi_{1}$}^{\prime}(2\text{ds}\;w_{34}+\text{ds}\;w_{35}\epsilon-2\text{dt}\;w_{14}-\text{dt}\;w_{15}\epsilon)\right)\,,$
$\displaystyle\tilde{A}_{f(3,3)}=$
$\displaystyle\epsilon(4\text{ds}\;w_{36}-4\text{dt}\;w_{17})\,,$
$\displaystyle\tilde{A}_{f(3,4)}=$
$\displaystyle\epsilon(2\text{ds}\;w_{37}-2\text{dt}\;w_{18})\,,$
$\displaystyle\tilde{A}_{f(4,1)}=$
$\displaystyle\frac{\epsilon}{128}\left(\text{$\psi_{0}$}^{\prime}(\text{ds}\;w_{33}+4\text{dt}\;w_{13})+2\text{$\psi_{0}$}(\text{ds}\;w_{28}+4\text{dt}\;w_{8})\right)\,,$
$\displaystyle\tilde{A}_{f(4,2)}=$
$\displaystyle\frac{\epsilon}{256}\left(\text{$\psi_{1}$}^{\prime}(\text{ds}\;w_{33}+4\text{dt}\;w_{13})+2\text{$\psi_{1}$}(\text{ds}\;w_{28}+4\text{dt}\;w_{8})\right)\,,$
$\displaystyle\tilde{A}_{f(4,3)}=$
$\displaystyle\frac{\epsilon}{2}(\text{dt}\;w_{19}-\text{ds}\;w_{38})\,,$
$\displaystyle\tilde{A}_{f(4,4)}=$
$\displaystyle\epsilon(\text{dt}\;w_{20}-\text{ds}\;w_{39})\,,$

where

$\displaystyle w_{1}=\frac{s^{3}(s+8)(s+12)(s+16)}{s+2t}\,,\qquad
w_{2}=\frac{s^{5}(s+16)^{3}}{s+2t}\,,\qquad
w_{3}=\frac{s^{4}(s+8)(s+16)^{2}}{s+2t}\,,$ $\displaystyle
w_{4}=\frac{s^{4}(s+12)(s+16)^{2}}{s+2t}\,,\qquad
w_{5}=\frac{s^{3}(s+16)(2s(s+8)(s+16)+2s(s+12)(s+16))}{s+2t}\,,$
$\displaystyle
w_{6}=\frac{s^{2}(s+16)\left(s^{2}+(s+12)t^{2}+(s+12)st\right)}{\sqrt{t}\sqrt{s+t}(s+2t)^{2}\sqrt{s(t-4)+t^{2}}}\,,\qquad
w_{7}=\frac{s^{3}(s+8)}{(s+2t)^{2}}\,,\qquad w_{8}=s(s+8)\,,$ $\displaystyle
w_{9}=\frac{s\left(2s^{2}(t-5)+2s\left(t^{2}+8t-48\right)+16t^{2}\right)}{\sqrt{t}\sqrt{s+t}\sqrt{s(t-4)+t^{2}}}\,,w_{10}=\frac{s\left(3s^{3}+2s^{2}(6t+22)+12s(t+4)^{2}+96t^{2}\right)}{\sqrt{t}\sqrt{s+t}\sqrt{s(t-4)+t^{2}}},$
$\displaystyle
w_{11}=\frac{s^{3}(s+16)^{2}\sqrt{t}\sqrt{s+t}}{(s+2t)^{2}\sqrt{s(t-4)+t^{2}}}\,,\qquad
w_{12}=\frac{s^{4}(s+16)}{(s+2t)^{2}}\,,\qquad w_{13}=s^{2}(s+16)\,,$
$\displaystyle
w_{14}=\frac{s^{2}(s+16)\sqrt{st-4s+t^{2}}}{\sqrt{t}\sqrt{s+t}}\,,\qquad
w_{15}=\frac{s^{2}(s+16)\left(3s^{2}+12st+16s+12t^{2}\right)}{\sqrt{t}\sqrt{s+t}\sqrt{s(t-4)+t^{2}}}\,,$
$\displaystyle w_{16}=\frac{s^{2}}{t(s+t)(s+2t)}\,,\qquad
w_{17}=\frac{s^{3}(t-1)+4s^{2}t^{2}+6st^{3}+3t^{4}}{t(s+t)(s+2t)\left(s(t-4)+t^{2}\right)}\,,$
$\displaystyle
w_{18}=\frac{3s^{2}+8st+8t^{2}}{\sqrt{t}\sqrt{s+t}(s+2t)\sqrt{s(t-4)+t^{2}}}\,,\qquad
w_{19}=\frac{s+2t}{\sqrt{t}\sqrt{s+t}\sqrt{s(t-4)+t^{2}}}\,,\quad
w_{20}=\frac{s+2t}{t(s+t)}\,,$ $\displaystyle
w_{21}=\frac{s^{2}((384-s(s(s+16)+24))t+2s(s+10)(5s+48))}{s+2t}\,,$
$\displaystyle
w_{22}=\frac{s^{3}(s+16)\left(s^{3}(t-11)+s^{2}(t-14)(t+8)+4s(t-40)t-64t^{2}\right)}{(s+t)(s+2t)}\,,$
$\displaystyle
w_{23}=\frac{s^{3}(s+16)\left(-s^{2}(7s+64)+(s+4)(s+8)t^{2}+(s(s+4)-48)st\right)}{(s+t)(s+2t)}\,,w_{24}=\frac{s^{5}(s+16)^{2}(t-8)}{s+2t},$
$\displaystyle w_{25}=\frac{s^{3}(s+16)(s((s+8)t-9s-88)-16t)}{s+2t}\,,$
$\displaystyle
w_{26}=\frac{s\sqrt{t}\left(-4(s+8)s^{2}+(s(s+18)+96)t^{2}+(s(s+13)+48)st\right)}{\sqrt{s+t}(s+2t)^{2}\sqrt{s(t-4)+t^{2}}}\,,\qquad
w_{27}=\frac{s^{2}(s-(s+6)t)}{(s+2t)^{2}}\,,$ $\displaystyle
w_{28}=(s+8)(s-2t)\,,\qquad
w_{29}=\frac{\sqrt{t}\left((s+10)t^{2}+s(s+11)t-4s(s+12)\right)}{\sqrt{s+t}\sqrt{s(t-4)+t^{2}}}\,,$
$\displaystyle
w_{30}=\frac{\sqrt{t}\left(2(s+12)(3s+16)t^{2}+3s(s(s+4)-64)t+64s(5s+48)\right)}{(s+16)\sqrt{s+t}\sqrt{s(t-4)+t^{2}}}\,,$
$\displaystyle
w_{31}=\frac{s^{2}(s+16)\sqrt{t}\sqrt{s+t}(s(t-4)+8t)}{(s+2t)^{2}\sqrt{s(t-4)+t^{2}}}\,,\qquad
w_{32}=\frac{s^{3}(s+16)t}{(s+2t)^{2}}\,,\qquad w_{33}=s(s+16)(s-2t)\,,$
$\displaystyle
w_{34}=\frac{s(s+16)\sqrt{t}\sqrt{s(t-4)+t^{2}}}{\sqrt{s+t}}\,,\qquad
w_{35}=\frac{s\sqrt{t}\left(3s^{2}t+2(3s+32)t^{2}+256s\right)}{\sqrt{s+t}\sqrt{s(t-4)+t^{2}}}\,,$
$\displaystyle
w_{36}=\frac{t\left(3s^{2}+s(t+2)t+t^{3}\right)}{s(s+t)(s+2t)\left(s(t-4)+t^{2}\right)}\,,\qquad
w_{37}=\frac{\sqrt{t}\left(-s^{2}(t-16)+2st(t+16)+64t^{2}\right)}{s(s+16)\sqrt{s+t}(s+2t)\sqrt{s(t-4)+t^{2}}}\,,$
$\displaystyle w_{38}=\frac{t^{3/2}}{s\sqrt{s+t}\sqrt{s(t-4)+t^{2}}}\,,\qquad
w_{39}=\frac{s+2t}{s(s+t)}\,.$

In equation (5.1), we have suppressed the explicit dependence of the periods
and their derivatives on $s$.

###  5.2 The definitions of the entries in $U_{4}$

In this section, we first provide the definitions of $f(s,t)$ and $g(s,t)$
mentioned in section 3.2, that appears in the last transformation matrix
$U_{4}$ needed to transform the differential equations corresponding to the
top-sector to an $\epsilon$-factorised form. Later we also perform a study of
the maximal cut again to suggest a change of variables that simplifies these
functions.

$\displaystyle f(s,t)$
$\displaystyle=\frac{s}{\sqrt{t}\sqrt{s+t}}\Bigg{\\{}2(s+12)(s+t)\sqrt{s(t-4)+t^{2}}$
$\displaystyle-\frac{2\sqrt{2}\sqrt{s}\sqrt{t}(5s+64)}{\sqrt{\frac{\big{(}s+\sqrt{s}\sqrt{s+16}+8\big{)}}{s+t}}}\;\;F\bigg{(}\sin^{-1}\bigg{(}\frac{\sqrt{1+\frac{(t-8)\sqrt{s}}{t\sqrt{s+16}}}}{\sqrt{2}}\bigg{)}|\frac{2\sqrt{s}\sqrt{s+16}}{s+\sqrt{s}\sqrt{s+16}+8}\bigg{)}$
$\displaystyle-\frac{\sqrt{2}(s+12)\sqrt{s}(s+t)}{\sqrt{\frac{s+t}{\left(s+\sqrt{s}\sqrt{s+16}+8\right)t}}}\;\;E\bigg{(}\sin^{-1}\bigg{(}\frac{\sqrt{1+\frac{(t-8)\sqrt{s}}{t\sqrt{s+16}}}}{\sqrt{2}}\bigg{)}|\frac{2\sqrt{s}\sqrt{s+16}}{s+\sqrt{s}\sqrt{s+16}+8}\bigg{)}\Bigg{\\}},$
(33) $\displaystyle g(s,t)$
$\displaystyle=-\frac{s(s+16)}{4\sqrt{t}\sqrt{s+t}}\Bigg{\\{}\big{(}-2(s+t)\sqrt{{s(t-4)+t^{2}}}$
$\displaystyle+\frac{8\sqrt{2}\sqrt{s}\sqrt{t}}{\sqrt{\frac{\big{(}s+\sqrt{s}\sqrt{s+16}+8\big{)}}{s+t}}}\;\;F\bigg{(}\sin^{-1}\bigg{(}\frac{\sqrt{1+\frac{(t-8)\sqrt{s}}{t\sqrt{s+16}}}}{\sqrt{2}}\bigg{)}|\frac{2\sqrt{s}\sqrt{s+16}}{s+\sqrt{s}\sqrt{s+16}+8}\bigg{)}$
$\displaystyle+\frac{\sqrt{2}\sqrt{s}(s+t)}{{\sqrt{\frac{s+t}{\big{(}s+\sqrt{s}\sqrt{s+16}+8\big{)}t}}}}\;\;E\bigg{(}\sin^{-1}\bigg{(}\frac{\sqrt{1+\frac{(t-8)\sqrt{s}}{t\sqrt{s+16}}}}{\sqrt{2}}\bigg{)}|\frac{2\sqrt{s}\sqrt{s+16}}{s+\sqrt{s}\sqrt{s+16}+8}\bigg{)}\Bigg{\\}},$
(34)

where $F(\cdot|\cdot)$ and $E(\cdot|\cdot)$ are respectively the elliptic
integral of the first and second kind:

$F(\phi|k^{2})=\int_{0}^{\sin{\phi}}\frac{dt}{\sqrt{(1-t^{2})(1-k^{2}t^{2})}},\qquad
E(\phi|k^{2})=\int_{0}^{\sin{\phi}}\sqrt{\frac{1-k^{2}t^{2}}{1-t^{2}}}dt\,.$
(35)

###  5.3 The study of the marked point on the torus

Below we study the correspondence between the elliptic curve and a torus to
find a possible change of variables that enables us to express the functions
$f(s,t)$ and $g(s,t)$ in a more compact notation. This might be useful when
aiming to express the results through eMPLs in future works. From the maximal
cut in equation (13) we can see that there is a pole in $P$ at $s+t$. A pole
in the maximal cut corresponds to a marked point in the torus. The map between
an elliptic curve and the corresponding torus is given by Abel’s map and as
done in Giroux:2022wav ; Delto:2023kqv  we can compute how the pole is marked
on the torus

$Z_{x_{p}}=\frac{F(\sin^{-1}{u_{x_{p}}}|k)}{2K(k)},$ (36)

where

$k^{2}=\frac{r_{23}\;r_{14}}{r_{13}r_{24}},\qquad r_{ij}=r_{i}-r_{j}$ (37)

and

$u_{x_{p}}=\sqrt{\frac{(x_{p}-r_{1})(r_{2}-r_{4})}{(x_{p}-r_{2})(r_{1}-r_{4})}}\,.$
(38)

In the above equation, $r_{i}$ are the roots of the elliptic curve and $x_{p}$
is the pole on the maximal cut.

With the roots of our elliptic curve

$r_{1}=\frac{1}{2}\left(-s-\sqrt{s}\sqrt{s+16}\right),\quad r_{2}=-s,\quad
r_{3}=0,\quad r_{4}=\frac{1}{2}\left(\sqrt{s}\sqrt{s+16}-s\right),\quad
x_{p}=-s-t,$ (39)

we get the following $k^{2}$ and $u_{x_{p}}$

$k^{2}=\frac{2\sqrt{s}\sqrt{s+16}}{s+\sqrt{s}\sqrt{s+16}+8},\qquad
u_{x_{p}}=\frac{\sqrt{1+\frac{(t-8)\sqrt{s}}{t\sqrt{s+16}}}}{\sqrt{2}},$ (40)

and the following marked point on the torus

$Z_{x_{p}}=\frac{F\bigg{(}\sin^{-1}\bigg{(}\frac{\sqrt{1+\frac{(t-8)\sqrt{s}}{t\sqrt{s+16}}}}{\sqrt{2}}\bigg{)}|\frac{2\sqrt{s}\sqrt{s+16}}{s+\sqrt{s}\sqrt{s+16}+8}\bigg{)}}{2K\bigg{(}\frac{2\sqrt{s}\sqrt{s+16}}{s+\sqrt{s}\sqrt{s+16}+8}\bigg{)}}.$
(41)

We notice that this is the same integral that appears in the expression for
$f(s,t)$ and $g(s,t)$, so we rewrite them using these variables and notice
that their expression simplifies

$\displaystyle f(s,t)=\frac{2s(s+12)\sqrt{s+t}\sqrt{s(t-4)+t^{2}}}{\sqrt{t}}$
$\displaystyle+\frac{\sqrt{2}s^{3/2}}{\sqrt{s+\sqrt{s}\sqrt{s+16}+8}}\left(4(5s+64)K(k)Z_{x_{p}}+(s+12)\left(s+\sqrt{s}\sqrt{s+16}+8\right)E\left(\left.\sin^{-1}(u_{x_{p}})\right|k\right)\right),$
(42) $\displaystyle
g(s,t)=\frac{(s+16)s\sqrt{s+t}\sqrt{s(t-4)+t^{2}}}{2\sqrt{t}}$
$\displaystyle+\frac{\sqrt{2}s^{3/2}(s+16)}{\sqrt{s+\sqrt{s}\sqrt{s+16}+8}}\left(4K(k)Z_{x_{p}}+\frac{1}{4}\left(s+\sqrt{s}\sqrt{s+16}+8\right)E\left(\left.\sin^{-1}(u_{x_{p}})\right|k\right)\right).$
(43)

## References

  * (1) ATLAS collaboration, _Search for new phenomena in dijet mass and angular distributions from $pp$ collisions at $\sqrt{s}=$ 13 TeV with the ATLAS detector_, _Phys. Lett. B_ 754 (2016) 302 [1512.01530].
  * (2) CMS collaboration, _Search for dijet resonances in proton–proton collisions at $\sqrt{s}$ = 13 TeV and constraints on dark matter and other models_, _Phys. Lett. B_ 769 (2017) 520 [1611.03568].
  * (3) CMS collaboration, _Search for physics beyond the standard model in high-mass diphoton events from proton-proton collisions at $\sqrt{s}=$ 13 TeV_, _Phys. Rev. D_ 98 (2018) 092001 [1809.00327].
  * (4) S. Caron-Huot and J.M. Henn, _Iterative structure of finite loop integrals_ , _JHEP_ 06 (2014) 114 [1404.2922].
  * (5) M. Becchetti and R. Bonciani, _Two-Loop Master Integrals for the Planar QCD Massive Corrections to Di-photon and Di-jet Hadro-production_ , _JHEP_ 01 (2018) 048 [1712.02537].
  * (6) A.A. H., E. Chaubey and H.-S. Shao, _Two-loop massive QCD and QED helicity amplitudes for light-by-light scattering_ , 2312.16966.
  * (7) M. Becchetti, R. Bonciani, L. Cieri, F. Coro and F. Ripani, _Two-loop form factors for diphoton production in quark annihilation channel with heavy quark mass dependence_ , _JHEP_ 12 (2023) 105 [2308.11412].
  * (8) A. von Manteuffel and L. Tancredi, _A non-planar two-loop three-point function beyond multiple polylogarithms_ , _JHEP_ 06 (2017) 127 [1701.05905].
  * (9) U. Aglietti, R. Bonciani, G. Degrassi and A. Vicini, _Analytic Results for Virtual QCD Corrections to Higgs Production and Decay_ , _JHEP_ 01 (2007) 021 [hep-ph/0611266].
  * (10) C. Anastasiou, S. Beerli, S. Bucherer, A. Daleo and Z. Kunszt, _Two-loop amplitudes and master integrals for the production of a Higgs boson via a massive quark and a scalar-quark loop_ , _JHEP_ 01 (2007) 082 [hep-ph/0611236].
  * (11) X. Xu and L.L. Yang, _Towards a new approximation for pair-production and associated-production of the Higgs boson_ , _JHEP_ 01 (2019) 211 [1810.12002].
  * (12) G. Wang, Y. Wang, X. Xu, Y. Xu and L.L. Yang, _Efficient computation of two-loop amplitudes for Higgs boson pair production_ , _Phys. Rev. D_ 104 (2021) L051901 [2010.15649].
  * (13) M. Hidding, _DiffExp, a Mathematica package for computing Feynman integrals in terms of one-dimensional series expansions_ , _Comput. Phys. Commun._ 269 (2021) 108125 [2006.05510].
  * (14) F. Tkachov, _A theorem on analytical calculability of 4-loop renormalization group functions_ , _Physics Letters B_ 100 (1981) 65.
  * (15) K.G. Chetyrkin and F.V. Tkachov, _Integration by parts: The algorithm to calculate $\beta$-functions in 4 loops_, _Nucl. Phys. B_ 192 (1981) 159.
  * (16) P. Maierhöfer, J. Usovitsch and P. Uwer, _Kira—A Feynman integral reduction program_ , _Comput. Phys. Commun._ 230 (2018) 99 [1705.05610].
  * (17) A.V. Smirnov, _FIRE5: a C++ implementation of Feynman Integral REduction_ , _Comput. Phys. Commun._ 189 (2015) 182 [1408.2372].
  * (18) A.V. Kotikov, _Differential equations method: New technique for massive Feynman diagrams calculation_ , _Phys. Lett. B_ 254 (1991) 158.
  * (19) A.V. Kotikov, _Some methods to evaluate complicated Feynman integrals_ , _Nucl. Instrum. Meth. A_ 502 (2003) 615 [hep-ph/0303059].
  * (20) E. Remiddi, _Differential equations for Feynman graph amplitudes_ , _Nuovo Cim. A_ 110 (1997) 1435 [hep-th/9711188].
  * (21) T. Gehrmann and E. Remiddi, _Differential equations for two-loop four-point functions_ , _Nuclear Physics B_ 580 (2000) 485.
  * (22) M. Argeri and P. Mastrolia, _Feynman Diagrams and Differential Equations_ , _Int. J. Mod. Phys. A_ 22 (2007) 4375 [0707.4037].
  * (23) J.M. Henn, _Multiloop integrals in dimensional regularization made simple_ , _Phys. Rev. Lett._ 110 (2013) 251601 [1304.1806].
  * (24) J.M. Henn, _Lectures on differential equations for Feynman integrals_ , _J. Phys. A_ 48 (2015) 153001 [1412.2296].
  * (25) A.B. Goncharov, _Polylogarithms in Arithmetic and Geometry_ , in _International Congress of Mathematicians_ , 1995, DOI.
  * (26) A.B. Goncharov, _Multiple polylogarithms and mixed Tate motives_ , math/0103059.
  * (27) E. Remiddi and J.A.M. Vermaseren, _Harmonic polylogarithms_ , _Int. J. Mod. Phys._ A15 (2000) 725 [hep-ph/9905237].
  * (28) J. Broedel, C. Duhr, F. Dulat, B. Penante and L. Tancredi, _Elliptic polylogarithms and Feynman parameter integrals_ , _JHEP_ 05 (2019) 120 [1902.09971].
  * (29) L. Adams, C. Bogner and S. Weinzierl, _The two-loop sunrise graph with arbitrary masses_ , _J. Math. Phys._ 54 (2013) 052303 [1302.7004].
  * (30) 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) 102301 [1405.5640].
  * (31) 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_ , _JHEP_ 10 (2018) 206 [1806.04981].
  * (32) C. Duhr, V.A. Smirnov and L. Tancredi, _Analytic results for two-loop planar master integrals for Bhabha scattering_ , _JHEP_ 09 (2021) 120 [2108.03828].
  * (33) S. Laporta and E. Remiddi, _Analytic treatment of the two loop equal mass sunrise graph_ , _Nucl. Phys. B_ 704 (2005) 349 [hep-ph/0406160].
  * (34) S. Muller-Stach, S. Weinzierl and R. Zayadeh, _From motives to differential equations for loop integrals_ , _PoS_ LL2012 (2012) 005 [1209.3714].
  * (35) S. Bloch and P. Vanhove, _The elliptic dilogarithm for the sunset graph_ , _J. Number Theor._ 148 (2015) 328 [1309.5865].
  * (36) 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) 072303 [1504.03255].
  * (37) 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) 032304 [1512.05630].
  * (38) L. Adams, C. Bogner, A. Schweitzer and S. Weinzierl, _The kite integral to all orders in terms of elliptic polylogarithms_ , _J. Math. Phys._ 57 (2016) 122302 [1607.01571].
  * (39) L. Adams, E. Chaubey and S. Weinzierl, _Simplifying Differential Equations for Multiscale Feynman Integrals beyond Multiple Polylogarithms_ , _Phys. Rev. Lett._ 118 (2017) 141602 [1702.04279].
  * (40) 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 [1705.08952].
  * (41) 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 [1712.07089].
  * (42) L. Adams and S. Weinzierl, _The $\varepsilon$-form of the differential equations for Feynman integrals in the elliptic case_, _Phys. Lett. B_ 781 (2018) 270 [1802.05020].
  * (43) 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 [1803.10256].
  * (44) 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) 142001 [1804.11144].
  * (45) J. Broedel, C. Duhr, F. Dulat, B. Penante and L. Tancredi, _Elliptic Feynman integrals and pure functions_ , _JHEP_ 01 (2019) 023 [1809.10698].
  * (46) I. Hönemann, K. Tempest and S. Weinzierl, _Electron self-energy in QED at two loops revisited_ , _Phys. Rev. D_ 98 (2018) 113008 [1811.09308].
  * (47) 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 [1907.01251].
  * (48) J. Broedel, C. Duhr, F. Dulat, R. Marzucca, B. Penante and L. Tancredi, _An analytic solution for the equal-mass banana graph_ , _JHEP_ 09 (2019) 112 [1907.03787].
  * (49) S. Weinzierl, _Simple differential equations for Feynman integrals associated to elliptic curves_ , in _14th International Symposium on Radiative Corrections: Application of Quantum Field Theory to Phenomenology_ , 12, 2019, DOI [1912.02578].
  * (50) L. Adams and S. Weinzierl, _Feynman integrals and iterated integrals of modular forms_ , _Commun. Num. Theor. Phys._ 12 (2018) 193 [1704.08895].
  * (51) R. Marzucca, A.J. McLeod, B. Page, S. Pögel and S. Weinzierl, _Genus Drop in Hyperelliptic Feynman Integrals_ , 2307.11497.
  * (52) M. Giroux, A. Pokraka, F. Porkert and Y. Sohnle, _The soaring kite: a tale of two punctured tori_ , 2401.14307.
  * (53) L. Görges, C. Nega, L. Tancredi and F.J. Wagner, _On a procedure to derive $\epsilon$-factorised differential equations beyond polylogarithms_, _JHEP_ 07 (2023) 206 [2305.14090].
  * (54) M. Delto, C. Duhr, L. Tancredi and Y.J. Zhu, _Two-loop QED corrections to the scattering of four massive leptons_ , 2311.06385.
  * (55) S. Pögel, X. Wang and S. Weinzierl, _Bananas of equal mass: any loop, any order in the dimensional regularisation parameter_ , _JHEP_ 04 (2023) 117 [2212.08908].
  * (56) K.-T. Chen, _Iterated path integrals_ , _Bull. Am. Math. Soc._ 83 (1977) 831.
  * (57) P. Baikov, _Explicit solutions of the multi-loop integral recurrence relations and its application_ , _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ 389 (1997) 347.
  * (58) A. Primo and L. Tancredi, _Maximal cuts and differential equations for Feynman integrals. An application to the three-loop massive banana graph_ , _Nucl. Phys. B_ 921 (2017) 316 [1704.05465].
  * (59) H. Frellesvig and C.G. Papadopoulos, _Cuts of Feynman Integrals in Baikov representation_ , _JHEP_ 04 (2017) 083 [1701.07356].
  * (60) A. Primo and L. Tancredi, _On the maximal cut of Feynman integrals and the solution of their differential equations_ , _Nucl. Phys. B_ 916 (2017) 94 [1610.08397].
  * (61) P.A. Baikov, _Explicit solutions of the multiloop integral recurrence relations and its application_ , _Nucl. Instrum. Meth. A_ 389 (1997) 347 [hep-ph/9611449].
  * (62) S. Müller-Stach, S. Weinzierl and R. Zayadeh, _Picard-Fuchs equations for Feynman integrals_ , _Commun. Math. Phys._ 326 (2014) 237 [1212.4389].
  * (63) K. Bönisch, C. Duhr, F. Fischbach, A. Klemm and C. Nega, _Feynman integrals in dimensional regularization and extensions of Calabi-Yau motives_ , _JHEP_ 09 (2022) 156 [2108.05310].
  * (64) C. Meyer, _Algorithmic transformation of multi-loop master integrals to a canonical basis with CANONICA_ , _Comput. Phys. Commun._ 222 (2018) 295 [1705.06252].
  * (65) M. Giroux and A. Pokraka, _Loop-by-loop differential equations for dual (elliptic) Feynman integrals_ , _JHEP_ 03 (2023) 155 [2210.09898].