aainstitutetext: Physics Department, Torino University and INFN Torino, Via
Pietro Giuria 1, I-10125 Torino, Italybbinstitutetext: Dipartimento di Fisica
e Astronomia, Università di Bologna e INFN, Sezione di Bologna, via Irnerio
46, I-40126 Bologna, Italyccinstitutetext: Physik-Institut, Universität
Zürich, Winterthurerstrasse 190, CH-8057 Zürich, Switzerlandddinstitutetext:
CERN, Theoretical Physics Department, CH-1211 Geneva 23, Switzerland

# Two-loop integrals for $t\bar{t}+$jet production at hadron colliders in the
leading colour approximation

Simon Badger b Matteo Becchetti c Nicolò Giraudo d Simone Zoia
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We compute the differential equations for the two remaining integral
topologies contributing to the leading colour two-loop amplitudes for
$pp\rightarrow t\bar{t}j$. We derive differential equations for the master
integrals by solving the integration-by-parts identities over finite fields.
Of the two systems of differential equations, one is presented in canonical
‘${\rm d}\log$’ form, while the other is found to have an elliptic sector. For
the elliptic topology we identify the relevant elliptic curve, and present the
differential equations in a more general form which depends quadratically on
$\eps$ and contains non-logarithmic one-forms in addition to the canonical
${\rm d}\log$’s. We solve the systems of differential equations numerically
using generalised series expansions with the boundary terms obtained using the
auxiliary mass flow method. A summary of all one-loop and two-loop planar
topologies is presented including the list of alphabet letters for the ‘${\rm
d}\log$’ form systems and high-precision boundary values.

††preprint: CERN-TH-2024-048, ZU-TH 22/24

##  1 Introduction

Precision tests of the Standard Model (SM) in the top quark sector are
essential part of the current physics program at high energy colliders.
Observables involving a top-quark pair have a wide variety of applications,
and can be used to constrain SM parameters and parton distribution functions.
Furthermore, they are also important backgrounds for other searches. Top-quark
pair production in association with at least one hard jet forms a significant
fraction, around 50%, of all top-quark pair events and has been studied
extensively by both ATLAS and CMS experiments CMS:2016oae ; ATLAS:2018acq ;
CMS:2020grm ; CMS:2024ybg .

The normalised differential cross section for the $t{\bar{t}}j$ invariant mass
in top-quark pair production in association with one jet ($pp\to t{\bar{t}}j$)
has been proposed as an observable highly sensitive to the top-quark mass
Alioli:2013mxa . Phenomenological studies and experimental analyses at next-
to-leading order (NLO) in QCD have already demonstrated this as a viable
method for the extraction of the mass with current experimental data
Bevilacqua:2017ipv ; Alioli:2022lqo ; ATLAS:2019guf . While the current
experimental and theoretical uncertainties are comparable, higher order
theoretical predictions will be necessary to make significant improvements as
more data become available.

Precision theoretical predictions for top-quark processes present one of the
most challenging classes of computation. The presence of internal massive
particles creates analytic complexity in the corresponding Feynman integrals,
which means that conventional methods for their numerical evaluation can fail.
How to achieve a fast and stable evaluation of the integrals remains an active
field of research, see Bourjaily:2022bwx  and references therein. Coupling
this with the multiple scales associated with the three-parton final states
makes processes such as $pp\to t{\bar{t}}j$ and $pp\to t{\bar{t}}W/Z/H$
seriously daunting computations.

The NLO QCD corrections have been available since 2007 Dittmaier:2007wz ;
Dittmaier:2008uj . Remarkable theoretical developments since then now enable a
broad range of phenomenological predictions including full off-shell decays
and interfaces with parton shower Melnikov:2010iu ; Alioli:2011as ;
Czakon:2015cla ; Bevilacqua:2015qha ; Bevilacqua:2016jfk , as well as mixed
QCD and EW corrections Gutschow:2018tuk . Improving on the theoretical
precision requires going to next-to-next-to-leading order (NNLO) in QCD and,
while the tree-level and one-loop level amplitude ingredients can be provided
by the current generation of automated tools Campbell:2022qmc , the two-loop
amplitudes remain a bottleneck. In this article, we continue the task of
computing the two-loop QCD corrections to top-quark pair production in
association with a jet that was initiated in refs. Badger:2022mrb ;
Badger:2022hno . The preliminary tasks undertaken so far explored the
kinematic complexity of the helicity amplitudes at one-loop up to
$\mathcal{O}(\eps^{2})$ in the dimensional regularisation parameter $\eps$,
and the computation of master integrals relevant for a single ‘pentagon-box’
two-loop integral family via canonical form differential equations
Kotikov:1990kg ; Remiddi:1997ny ; Henn:2013pwa .

The last few years have seen rapid progress for $2\to 3$ scattering amplitudes
with massless internal propagators, in which a lot has been understood for
both the relevant Feynman integrals and the analytic representation of the
amplitudes. Two new methods have had a major impact on these developments: the
widespread use of finite field modular arithmetic for algebraic manipulations
and rational function reconstruction vonManteuffel:2014ixa ; Peraro:2016wsq ;
Peraro:2019svx , and the construction of special function bases, commonly
referred to as ‘pentagon functions’, suitable for fast and stable numerical
evaluations Gehrmann:2015bfy ; Gehrmann:2018yef ; Chicherin:2020oor ;
Chicherin:2021dyp ; Abreu:2023rco . Thanks to these developments, a plethora
of new results have now been obtained for the one- and two-loop five-particle
Feynman integrals relevant for the cases with no or one external massive
particle Gehrmann:2015bfy ; Papadopoulos:2015jft ; Abreu:2018aqd ;
Chicherin:2018old ; Abreu:2020jxa ; Canko:2020ylt ; Abreu:2021smk ;
Abreu:2023rco , scattering amplitudes Gehrmann:2015bfy ; Abreu:2018aqd ;
Badger:2018enw ; Abreu:2018jgq ; Chicherin:2018yne ; Abreu:2018zmy ;
Badger:2019djh ; Abreu:2019odu ; Hartanto:2019uvl ; Chicherin:2019xeg ;
Abreu:2020cwb ; Chawdhry:2020for ; Caron-Huot:2020vlo ; Abreu:2021asb ;
Badger:2021ega ; Agarwal:2021grm ; Badger:2021imn ; Chawdhry:2021mkw ;
Badger:2021nhg ; Abreu:2021oya ; Agarwal:2021vdh ; Badger:2022ncb ;
Hartanto:2022qhh ; Abreu:2023bdp ; DeLaurentis:2023izi ; Badger:2023mgf ;
DeLaurentis:2023nss ; Agarwal:2023suw , and differential distributions at NNLO
QCD accuracy Chawdhry:2019bji ; Kallweit:2020gcp ; Chawdhry:2021hkp ;
Badger:2021ohm ; Chen:2022ktf ; Buonocore:2022pqq ; Hartanto:2022qhh ;
Hartanto:2022ypo ; Badger:2023mgf ; Mazzitelli:2024ura . The first results for
massless six-particle kinematics have also recently appeared Henn:2021cyv ;
Henn:2024ngj  together with the first steps investigating processes involving
internal massive propagators Badger:2022mrb ; Badger:2022hno ;
FebresCordero:2023gjh ; Buccioni:2023okz ; Agarwal:2024jyq . Approximations of
the two-loop amplitudes for $pp\to t\bar{t}H$ may also be used to make precise
phenomenological predictions Catani:2022mfv ; Wang:2024pmv .

This article presents the completion of the Feynman integral topologies needed
for $pp\to t{\bar{t}}j$ at two loops in the leading colour approximation. The
construction of the master integral bases and their evaluation follow the
methodology used in the previous works Badger:2022mrb ; Badger:2022hno . Of
the two new pentagon-box topologies that appear, one is shown to have some
complicated features that go against the conventional wisdom for what should
appear for a leading colour, planar configuration. In particular, we find that
one sector requires the use of a nested square root in order to rotate the
differential equations into an $\eps$-factorised form (see sec. 3.3). This
feature is perhaps not completely unexpected since it has been observed before
albeit with more complicated kinematics such as $pp\to t{\bar{t}}H$
FebresCordero:2023gjh . The second feature is that one sector contains
elliptic integrals (see sec. 3.5). This is perhaps more surprising, as no sub-
sectors of this topology were previously known to be elliptic, although there
are no principles forbidding it. Both of these features mean that the route
towards a well defined basis of special functions which can can be efficiently
evaluated numerically is unclear. Canonical form differential equations have
been studied in the case of elliptic Feynman integrals although not for such
complicated kinematics Broedel:2017kkb ; Broedel:2018qkq ; Broedel:2018iwv ;
Broedel:2019hyg ; Frellesvig:2023iwr ; Gorges:2023zgv . Therefore in this
article we present a compact representation of the differential equations,
which are not $\eps$-factorised, in terms of a minimal set of independent one-
forms, and show that they can be reliably solved through generalised series
expansions Francesco:2019yqt .

The computation of the differential equations satisfied by the integral bases
requires the solution of systems of Integration-by-Parts (IBPs) identities
Tkachov:1981wb ; Chetyrkin:1981qh , which we generate using LiteRed Lee:2012cn
; Lee:2013mka  and NeatIBP Wu:2023upw . The latter package, in particular,
allows us to obtain optimised systems of IBP relations through the solution of
syzygy equations Gluza:2010ws , this way making their solution substantially
simpler. We solve the IBP relations with the Laporta algorithm Laporta:2001dd
within the FiniteFlow framework Peraro:2019svx ; Peraro:2019okx .

While we do not obtain canonical differential equations for the family
involving elliptic integrals, we do put some effort into making their analytic
properties transparent. In particular, we identify the elliptic curve
associated with the appearing elliptic integrals. Furthermore, we express the
differential equations in terms of a set of linearly independent one-forms. We
arrange the latter to be logarithmic as much as possible, in analogy with the
usual canonical differential equations. Only the one-forms which relate master
integrals involving either the elliptic curve or the nested square root are
non-logarithmic. This makes the expression of the differential equations
compact, and separates out clearly the features in common with the standard
canonical cases from the more complicated ones. We expect that this
information will be precious in view of future work to obtain a canonical form
for the differential equations of this topology as well.

We provide a semi-analytic solution to the differential equations for the
master integrals by means of the multivariate generalisation of the method of
generalised power series expansions Pozzorini:2005ff ; Aglietti:2007as ;
Lee:2017qql ; Lee:2018ojn ; Bonciani:2018uvv ; Fael:2021kyg ; Fael:2022rgm
proposed in ref. Francesco:2019yqt  and implemented in a number of public
Mathematica packages Hidding:2020ytt ; Liu:2022chg ; Armadillo:2022ugh . In
particular, we make use of DiffExp Hidding:2020ytt . The solution to the
differential equations is fully characterised once a set of boundary values is
given. We obtain high-precision numerical boundary values by using the
Mathematica package AMFlow Liu:2022chg , which implements the auxiliary mass
flow method Liu:2017jxz ; Liu:2021wks ; Liu:2022tji . We use the interface to
FiniteFlow (Peraro:2019svx, ) and LiteRed Lee:2012cn ; Lee:2013mka  for the
required IBP reduction.

The paper is structured as follows. In section 2 we define the pentagon-box
topologies under study and describe the computational framework. In section 3
we discuss the construction and the analytic features of the master integral
bases. Section 4 is devoted to the one-form representation of the differential
equations, while in section 5 we discuss the numerical evaluation of the
master integrals. We draw our conclusions and give an outlook on future
developments in section 6. Finally, we identify the elliptic curve underlying
the analytic structure of the Feynman integrals in the elliptic sector in
appendix A, define the one-loop families in appendix B, and describe the
supplementary material in appendix C.

##  2 Notation and definitions

(a) Topology ${\rm PB}_{A}$.

(b) Topology ${\rm PB}_{B}$.

(c) Topology ${\rm PB}_{C}$.

Figure 1: The three pentagon-box topologies contributing to $pp\rightarrow
t\bar{t}j$ in the leading colour limit. Black lines denote massless particles
and red double-lines denote massive particles.

In this article we consider the three two-loop pentagon-box integral
topologies shown in figure 1. In addition, we present an explicitly
‘${\mathrm{d}}\log$’ form for all one-loop pentagon integrals, which we define
in appendix B. The integrals of the pentagon-box topology $F\in\\{{\rm
PB}_{A},{\rm PB}_{B},{\rm PB}_{C}\\}$ have the form

$I_{a_{1},a_{2},a_{3},a_{4},a_{5},a_{6},a_{7},a_{8}}^{\left(F\right),a_{9},a_{10},a_{11}}=\int\mathcal{D}^{d}k_{1}\,\mathcal{D}^{d}k_{2}\,\frac{D_{F,9}^{a_{9}}D_{F,10}^{a_{10}}D_{F,11}^{a_{11}}}{D_{F,1}^{a_{1}}\cdots
D_{F,8}^{a_{8}}}\,,$ (1)

where $a_{1},\cdots,a_{11}\geq 0$ and $d=4-2\eps$. The integration measure is

$\mathcal{D}^{d}k_{i}=\mathrm{e}^{\eps\gamma_{\rm
E}}\dfrac{{\mathrm{d}}^{d}k_{i}}{{\mathrm{i}}\pi^{\frac{d}{2}}}\,.$ (2)

The inverse propagators $D_{F,i}$ are defined in table 1. Note that $D_{F,9}$,
$D_{F,10}$ and $D_{F,11}$ are irreducible scalar products.

| ${\rm PB}_{A}$ | ${\rm PB}_{B}$ | ${\rm PB}_{C}$
---|---|---|---
$D_{F,1}$ | $k_{1}^{2}$ | $k_{1}^{2}-m_{t}^{2}$ | $k_{1}^{2}$
$D_{F,2}$ | $(k_{1}-p_{1})^{2}-m_{t}^{2}$ | $(k_{1}-p_{2})^{2}$ | $(k_{1}-p_{3})^{2}$
$D_{F,3}$ | $(k_{1}-p_{1}-p_{2})^{2}$ | $(k_{1}-p_{2}-p_{3})^{2}$ | $(k_{1}-p_{3}-p_{4})^{2}$
$D_{F,4}$ | $(k_{1}-p_{1}-p_{2}-p_{3})^{2}$ | $(k_{1}-p_{2}-p_{3}-p_{4})^{2}$ | $(k_{1}-p_{3}-p_{4}-p_{5})^{2}$
$D_{F,5}$ | $k_{2}^{2}$ | $k_{2}^{2}-m_{t}^{2}$ | $k_{2}^{2}$
$D_{F,6}$ | $(k_{2}-p_{5})^{2}$ | $(k_{2}-p_{1})^{2}$ | $(k_{2}-p_{2})^{2}-m_{t}^{2}$
$D_{F,7}$ | $(k_{2}-p_{4}-p_{5})^{2}$ | $(k_{2}-p_{1}-p_{5})^{2}$ | $(k_{2}-p_{1}-p_{2})^{2}$
$D_{F,8}$ | $(k_{1}+k_{2})^{2}$ | $(k_{1}+k_{2})^{2}$ | $(k_{1}+k_{2})^{2}$
$D_{F,9}$ | $(k_{1}+p_{5})^{2}$ | $(k_{1}+p_{1})^{2}-m_{t}^{2}$ | $(k_{1}+p_{2})^{2}-m_{t}^{2}$
$D_{F,10}$ | $(k_{2}+p_{1})^{2}-m_{t}^{2}$ | $(k_{2}+p_{2})^{2}-m_{t}^{2}$ | $(k_{2}+p_{3})^{2}$
$D_{F,11}$ | $(k_{2}+p_{1}+p_{2})^{2}$ | $(k_{2}+p_{2}+p_{3})^{2}-m_{t}^{2}$ | $(k_{2}+p_{3}+p_{4})^{2}$
Table 1: Inverse propagators $D_{F,i}$ of the pentagon-box topologies shown in
figure 1.

The external momenta $p_{i}$ are considered as outgoing from the graphs and
all the particles are on-shell, i.e. $p_{1}^{2}=p_{2}^{2}=m_{t}^{2}$ while
$p_{3}^{2}=p_{4}^{2}=p_{5}^{2}=0$. The kinematics of the integrals can be
described in terms of six independent Lorentz invariants. Here we choose the
top-quark mass $m_{t}^{2}$ and the five adjacent scalar products,

$\displaystyle\vec{x}=\\{d_{12},d_{23},d_{34},d_{45},d_{15},m_{t}^{2}\\}\,,$
(3)

where

$d_{ij}=p_{i}\cdot p_{j}\,.$ (4)

The minimal set of master integrals (MIs) is obtained by Integration-by-Parts
(IBP) reduction Chetyrkin:1981qh ; Chetyrkin:1979bj . We generate the systems
of IBPs with the software LiteRed Lee:2012cn ; Lee:2013mka  and NeatIBP
Wu:2023upw , and solve them via the Laporta algorithm Laporta:2000dsw  within
the finite-field framework FiniteFlow (Peraro:2019svx, ). Topology ${\rm
PB}_{A}$ has been considered previously in ref. Badger:2022hno  and has 88
MIs. For the new cases presented here we find 121 MIs for topology ${\rm
PB}_{B}$ and 109 for topology ${\rm PB}_{C}$.111NeatIBP finds additional
symmetry relations with respect to LiteRed, reducing the number of MIs by 2
for each topology ${\rm PB}_{B}$ and ${\rm PB}_{C}$. The missing relations are
a result of both particles 1 and 2 having the same mass.

We find a set of 14 square roots appearing in the differential equations and
associated alphabet for all one-loop and planar two-loop topologies. These are
defined as

$\displaystyle\beta_{12}$
$\displaystyle=\sqrt{\frac{\operatorname{det}G(p_{1},p_{2})}{s_{12}}}=\sqrt{1-\frac{4m_{t}^{2}}{s_{12}}}\,,$
(5) $\displaystyle\beta_{34}$
$\displaystyle=\sqrt{1-\frac{4m_{t}^{2}}{s_{34}}}\,,$ (6)
$\displaystyle\beta_{45}$
$\displaystyle=\sqrt{1-\frac{4m_{t}^{2}}{s_{45}}}\,,$ (7)
$\displaystyle\Delta_{3,1}$
$\displaystyle=\sqrt{-\operatorname{det}G(p_{1},p_{23})}\,,$ (8)
$\displaystyle\Delta_{3,2}$
$\displaystyle=\sqrt{-\operatorname{det}G(p_{2},p_{15})}\,,$ (9)
$\displaystyle\Delta_{3,3}$
$\displaystyle=\sqrt{-\operatorname{det}G(p_{1},p_{25})}\,,$ (10)
$\displaystyle\Delta_{3,4}$
$\displaystyle=\sqrt{-\operatorname{det}G(p_{2},p_{13})}\,,$ (11)
$\displaystyle\Lambda_{1}$
$\displaystyle=\sqrt{\frac{\operatorname{det}Y^{0mm0}_{p_{1}|p_{3}|p_{2}|p_{45}}}{d_{23}^{2}d_{13}^{2}}}\,,$
(12) $\displaystyle\Lambda_{2}$
$\displaystyle=\sqrt{\frac{\operatorname{det}Y^{0mmm}_{p_{2}|p_{3}|p_{4}|p_{51}}}{d_{23}^{2}d_{34}^{2}}}\,,$
(13) $\displaystyle\Lambda_{3}$
$\displaystyle=\sqrt{\frac{\operatorname{det}Y^{m0mm}_{p_{13}|p_{2}|p_{4}|p_{5}}}{d_{24}^{2}d_{45}^{2}}}\,,$
(14) $\displaystyle\Lambda_{4}$
$\displaystyle=\sqrt{\frac{\operatorname{det}Y^{m0mm}_{p_{1}|p_{23}|p_{4}|p_{5}}}{d_{15}^{2}d_{45}^{2}}}\,,$
(15) $\displaystyle\Lambda_{5}$
$\displaystyle=\sqrt{\frac{\operatorname{det}Y^{mmmm}_{p_{12}|p_{3}|p_{4}|p_{5}}}{d_{34}^{2}d_{45}^{2}}}\,,$
(16) $\displaystyle\Lambda_{6}$
$\displaystyle=\sqrt{(s_{15}-s_{23})^{2}+\frac{4s_{34}s_{45}m_{t}^{2}}{s_{12}}}\,,$
(17) $\displaystyle\operatorname{tr}_{5}$
$\displaystyle=4\sqrt{\operatorname{det}G(p_{3},p_{4},p_{5},p_{1})}\,,$ (18)

where $G_{ij}(\vec{v})=v_{i}\cdot v_{j}$ is the Gram matrix,
$p_{ij}=p_{i}+p_{j}$, $s_{ij}=(p_{i}+p_{j})^{2}$, and $m$ denotes $m_{t}$ for
the sake of compactness. $Y^{m_{1}\ldots m_{4}}_{P_{1}|\ldots|P_{4}}$
represents the Cayley matrices associated with one-loop box configurations
appearing at sub-leading colour, which are defined by

$\displaystyle\left[Y_{P_{1}|P_{2}|P_{3}|P_{4}}^{m_{1}m_{2}m_{3}m_{4}}\right]_{ij}=\frac{1}{2}\left[(q_{i}-q_{j})^{2}-m_{i}^{2}-m_{j}^{2}\right]\,,\qquad
q_{i}=\sum_{k=0}^{i-1}P_{k}\,,$ (19)

for $i,j=1,\ldots,4$, with $P_{0}=0$. $\Lambda_{6}$ only appears at two loops.

Note that $\operatorname{tr}_{5}$ is related to the five-point pseudo-scalar
invariant via

$\displaystyle\operatorname{tr}_{5}^{2}={\rm
tr}(\gamma_{5}\not{p}_{3}\not{p}_{4}\not{p}_{5}\not{p}_{1})^{2}\,.$ (20)

We refrain from identifying $\operatorname{tr}_{5}$ with ${\rm
tr}(\gamma_{5}\not{p}_{3}\not{p}_{4}\not{p}_{5}\not{p}_{1})$, as the latter is
a parity-odd object, i.e. it changes sign under space-time parity conjugation
and odd-signature permutations. The parity degree of freedom is required to
describe scattering amplitudes, but not for the computation of Feynman
integrals. We therefore prefer to define $\operatorname{tr}_{5}$ as the
parity-even square root in eq. (18). Furthermore, we will make use of the
short-hand

$\displaystyle{\rm
tr}_{\pm}\left[\not{p}_{i}\not{p}_{j}\not{p}_{k}\not{p}_{l}\right]=\frac{1}{2}\,{\rm
tr}\left[\not{p}_{i}\not{p}_{j}\not{p}_{k}\not{p}_{l}(1\pm\gamma_{5})\right]$
(21)

to make certain expressions more compact. We however remove the parity degree
of freedom in them by replacing the parity-odd trace ${\rm
tr}(\gamma_{5}\not{p}_{3}\not{p}_{4}\not{p}_{5}\not{p}_{1})$ with the parity-
even square root $\operatorname{tr}_{5}$ defined in eq. (18).

In order to write compact expressions for the MIs which satisfy canonical
differential equations, we use choices motivated by the local numerators
introduced in refs. Arkani-Hamed:2010zjl ; Arkani-Hamed:2010pyv ;
Gehrmann:2015bfy ; Badger:2016ozq ; Abreu:2020jxa . In $d=4-2\eps$ dimensions,
we can find numerator insertions that can be written in terms of the $-2\eps$
dimensional components of the loop momenta and are conventionally denoted
$\mu_{ij}$,

$\displaystyle k_{i}$ $\displaystyle=k_{i}^{[4]}+k_{i}^{[-2\eps]},$
$\displaystyle\mu_{ij}$ $\displaystyle=-k_{i}^{[-2\epsilon]}\cdot
k_{j}^{[-2\epsilon]}.$ (22)

We denote these numerator insertions using an additional superscript $[ij]$,
as

$\displaystyle\begin{aligned}
I_{a_{1},a_{2},a_{3},a_{4},a_{5},a_{6},a_{7},a_{8}}^{(F),[ij],a_{9},a_{10},a_{11}}&=\int\mathcal{D}^{d}k_{1}\,\mathcal{D}^{d}k_{2}\,\mu_{ij}\,\frac{D_{F,9}^{a_{9}}D_{F,10}^{a_{10}}D_{F,11}^{a_{11}}}{D_{F,1}^{a_{1}}\cdots
D_{F,8}^{a_{8}}}\,,\\\
I_{a_{1},a_{2},a_{3},a_{4},a_{5},a_{6},a_{7},a_{8}}^{(F),[ij,kl],a_{9},a_{10},a_{11}}&=\int\mathcal{D}^{d}k_{1}\,\mathcal{D}^{d}k_{2}\,\mu_{ij}\,\mu_{kl}\,\frac{D_{F,9}^{a_{9}}D_{F,10}^{a_{10}}D_{F,11}^{a_{11}}}{D_{F,1}^{a_{1}}\cdots
D_{F,8}^{a_{8}}}\,.\end{aligned}$ (23)

##  3 Master integral bases

In this section we discuss the construction of the bases of master integrals
for topology ${\rm PB}_{B}$ and ${\rm PB}_{C}$, and describe their features.
The guiding principle in this construction is the simplification of the
differential equations (DEs) satisfied by the MIs. Let $\vec{\mathcal{I}}_{F}$
be the list of MIs for topology $F$. In general, $\vec{\mathcal{I}}_{F}$
satisfies a system of DEs of the form Barucchi:1973zm ; Kotikov:1990kg ;
Kotikov:1991hm ; Gehrmann:1999as ; Bern:1993kr

$\displaystyle\frac{\partial}{\partial
x_{i}}\,\vec{\mathcal{I}}_{F}(\vec{x},\eps)=A_{x_{i}}^{(F)}\left(\vec{x},\eps\right)\cdot\vec{\mathcal{I}}_{F}(\vec{x},\eps)\,,$
(24)

for every $i=1,\ldots,6$. The matrices $A_{x_{i}}^{(F)}$ are called
_connection matrices_. We rewrite eq. (24) in a more compact form by
introducing the total differential with respect to the kinematic invariants,
${\mathrm{d}}$, as

$\displaystyle{\mathrm{d}}\,\vec{\mathcal{I}}_{F}(\vec{x},\eps)={\mathrm{d}}A^{(F)}\left(\vec{x},\eps\right)\cdot\vec{\mathcal{I}}_{F}(\vec{x},\eps)\,,$
(25)

where ${\mathrm{d}}A^{(F)}$ is the matrix-valued one-form

$\displaystyle{\mathrm{d}}A^{(F)}(\vec{x},\eps)=\sum_{i=1}^{6}A_{x_{i}}^{(F)}\left(\vec{x},\eps\right)\,{\mathrm{d}}x_{i}\,.$
(26)

With a slight abuse of notation, we refer to
${\mathrm{d}}A^{(F)}(\vec{x},\eps)$ as connection matrix as well. The solution
to eq. (25) is enormously simplified if a choice of MIs is found such that the
DEs take the _canonical form_ Henn:2013pwa

${\mathrm{d}}\,\vec{\mathcal{I}}_{F}(\vec{x},\eps)=\eps\,{\mathrm{d}}A^{(F)}(\vec{x})\,\vec{\mathcal{I}}_{F}(\vec{x},\eps)\,,$
(27)

where the dependence of the connection matrices on $\eps$ is factorised, and
${\mathrm{d}}A^{(F)}(\vec{x})$ is a linear combination of logarithmic one-
forms:

${\mathrm{d}}A^{(F)}(\vec{x})=\sum_{i}c^{(F)}_{i}{\mathrm{d}}\log\left(W_{i}(\vec{x})\right)\,.$
(28)

Here, the $c^{(F)}_{i}$ are matrices of rational numbers, and the _letters_
$W_{i}(\vec{x})$ are algebraic functions of the kinematic invariants
$\vec{x}$. Their ensemble, called _alphabet_ , dictates the singularity
structure of the MIs. The factorisation of $\eps$ allows us to express the
solution algorithmically in terms of Chen iterated integrals Chen:1977oja ,
order by order in the Laurent expansion around $\eps=0$. This, in conjunction
with the presence of logarithmic one-forms only, enables the application of a
well-established toolbox of mathematical techniques — most notably the symbol
Goncharov:2010jf  — to write down and manipulate the solution. Building on
this, the method of the so-called _pentagon functions_ Gehrmann:2018yef ;
Chicherin:2020oor ; Badger:2021nhg ; Chicherin:2021dyp ; Badger:2023xtl ;
Abreu:2023rco ; FebresCordero:2023gjh  has proven particularly successful in
the computation of two-loop amplitudes for $2\to 3$ processes.

It is however known that the DEs for Feynman integrals can take more
complicated forms. Indeed, we anticipate that the canonical form in eqs. (27)
and (28) can only be achieved for topologies ${\rm PB}_{A}$ and ${\rm
PB}_{C}$, whereas a generalisation is necessary for topology ${\rm PB}_{B}$.
First of all, it is not proven that one can always factorise the dependence on
$\eps$ in the connection matrices. On top of that, even when $\eps$ is
factorised, one-forms other than ${\mathrm{d}}\log$’s may be necessary (see
e.g. the review Bourjaily:2022bwx  and references therein). For such cases,
the notion of ‘canonical’ DEs is still under debate Broedel:2018qkq ;
Frellesvig:2023iwr . On the one hand, the techniques for bringing the DEs to
an $\eps$-factorised form are much less mature than in the ${\mathrm{d}}\log$
case. Moreover, even when an $\eps$-factorised form is achieved which involves
one-forms more complicated than the ${\mathrm{d}}\log$’s in eq. (28),
manipulating and evaluating the solution efficiently remain challenging. For
these reasons, it is often convenient to resort to more flexible numerical
approaches to solve DEs beyond the ${\mathrm{d}}\log$ case, as opposed to
fully analytic solutions in terms of well understood special functions. The
method of generalised power series expansions Francesco:2019yqt  is proving
particularly effective, boosted by the availability of public implementations
Hidding:2020ytt ; Liu:2022chg ; Armadillo:2022ugh . Nonetheless, also within
this approach to the solution, simplifying as much as possible the form of the
DEs is crucial to an efficient and stable evaluation of the solution. In
particular, it is desirable for the connection matrices to depend polynomially
on $\eps$, and for the degree in $\eps$ to be as low as possible. More
explicitly, in the generalisation of the canonical DEs we consider in this
work for topology ${\rm PB}_{B}$, the connection matrix has the form

${\mathrm{d}}A^{({\rm PB}_{B})}(\vec{x},\eps)=\sum_{k=0}^{k_{\rm
max}}\sum_{i}\eps^{k}\,c^{({\rm PB}_{B})}_{k,i}\,\omega_{i}(\vec{x})\,,$ (29)

where $\omega_{i}(\vec{x})$ are (linearly independent) _one-forms_ , $k_{\rm
max}\in\mathbb{N}$, and $c^{({\rm PB}_{B})}_{k,i}$ are matrices of rational
numbers. More explicitly, the one-forms have the form

$\displaystyle\omega_{i}(\vec{x})=\mathcal{C}_{i}(\vec{x})\sum_{j=1}^{6}f_{ij}(\vec{x})\,{\mathrm{d}}x_{j}\,,$
(30)

where $f_{ij}(\vec{x})$ are rational functions, and $\mathcal{C}_{i}(\vec{x})$
is either $1$ or a square root (possibly a product of square roots). A subset
of the one-forms $\\{\omega_{i}(\vec{x})\\}$ may be logarithmic, i.e., for
some one-form $\omega(\vec{x})$ there may exist $W(\vec{x})$ such that
$\omega(\vec{x})={\mathrm{d}}\log W(\vec{x})$. In this case, we recall that
$W(\vec{x})$ is called _letter_.

Equation (30) implies that each one-form (including the logarithmic ones)
possesses a property called _charge_ with respect to the square roots of the
problem. With respect to each square root $\mathcal{C}$, a one-form $\omega$
is either _even_ , if it stays invariant when we flip the sign of
$\mathcal{C}$, or _odd_ , if it changes sign:

$\displaystyle\omega\bigl{|}_{\mathcal{C}\to-\mathcal{C}}=\begin{cases}\omega\,,\quad&\text{if
$\omega$ is \emph{even} w.r.t.\ ${\cal C}$}\,,\\\ -\omega\,,\quad&\text{if
$\omega$ is \emph{odd} w.r.t.\ ${\cal C}$}\,.\end{cases}$ (31)

For a logarithmic one-form ${\mathrm{d}}\log W$, this definition implies that
the letter $W$ is even (odd) with respect to a square root $\mathcal{C}$ if
$W\bigl{|}_{\mathcal{C}\to-\mathcal{C}}=W$
($W\bigl{|}_{\mathcal{C}\to-\mathcal{C}}=1/W$). For example, a typical odd
letter has the form

$\displaystyle\frac{A+\mathcal{C}}{A-\mathcal{C}}\,,$ (32)

where $\mathcal{C}$ is square root and $A$ is a rational function. The fact
that all one-forms appearing in the DEs for the MIs have a well-defined
transformation under the change of the sign of the square roots follows from
the way the square roots enter in the definition of the MIs. The scalar
integrals, i.e. the integrals of the form
$I^{(F),a_{9},a_{10},a_{11}}_{a_{1},\ldots,a_{8}}$, are by definition even
with respect to all square roots. The square roots enter the game in the
construction of MIs which satisfy DEs in the canonical form. As we will see
later in this section (e.g. see eq. (38)), the square roots appear as overall
normalisation of the MIs. More explicitly, each MI $\mathcal{I}$ has the form

$\displaystyle\mathcal{I}(\vec{x},\eps)=\mathcal{C}(\vec{x})\,\sum_{i}c_{i}(\vec{x},\eps)\,I_{i}(\vec{x},\eps)\,,$
(33)

where $\mathcal{C}(\vec{x})$ is either $1$ or a product of square roots,
$c_{i}(\vec{x},\eps)$ are rational functions, and $I_{i}$ are scalar
integrals. This way, the MIs gain a charge with respect to the square roots,
and the entries of the connection matrices — and thus the one-forms — inherit
it from the MIs. With a slight abuse of notation, we say that a MI or one-form
has charge $\mathcal{C}$ if it is odd with respect to $\mathcal{C}$.

In order to construct a basis of MIs which satisfies DEs of the previous
forms, we used the approach outlined in ref. Badger:2022mrb . For
completeness, we briefly summarise it in the following. It is a bottom-up
approach, i.e. we start from the integral sectors with the fewest number of
propagators, and we bring the DEs to the desired form sector by sector. The
simplification of the DEs is done by following a procedure which exploits a
set of heuristic criteria and draws as much as possible from known results in
the literature. In particular, we perform the following steps. Let $S$ be a
sector of topology $F$.

  * •

Step 1: We choose candidate MIs $\vec{\mathcal{L}}_{F,S}(\vec{x},\epsilon)$
for sector $S$ by requiring that the homogeneous DE222Given an integral
sector, the homogeneous DE is the subset of a DE which contains the
contributions coming only from the integrals of that sector. The corresponding
connection matrix is zero everywhere except for the diagonal square block
corresponding to the MIs of the chosen sector. As a consequence, since all the
sub-topologies do not contribute to it, the maximal cut of a certain MI is
also a solution to the homogeneous DEs of the corresponding sector
Primo:2016ebd . has the $\epsilon$ structure

${\mathrm{d}}\,\vec{\mathcal{L}}_{F,S}(\vec{x},\epsilon)=\sum_{k=0}^{k_{\rm
max}}\epsilon^{k}{\mathrm{d}}\tilde{A}^{(F,S)}_{k,\operatorname{Hom}}(\vec{x})\cdot\vec{\mathcal{L}}_{F,S}(\vec{x},\epsilon)+(\text{sub-
sectors})\,,$ (34)

where $k_{\rm max}=2$ for the elliptic sector in ${\rm PB}_{B}$, otherwise
$k_{\rm max}=1$. We neglect all contributions from the sub-sectors, which
constitute the inhomogeneous terms of the DEs. As a guiding principle, we
select MI candidates following patterns observed in previously studied cases.
Since in this step we are mostly interested in the $\epsilon$ structure of DEs
for the candidate MIs under study, we exploit finite fields technique to
perform IBP reduction and to reconstruct the DEs on a univariate
$\epsilon$-slice, i.e. we set to numbers all the kinematic invariants and
reconstruct just the analytic dependence in $\epsilon$.

  * •

Step 2: We reconstruct analytically the homogeneous DEs. Except for the two
problematic sectors shown in fig. 2, which we will discuss later, we construct
a rational transformation such that the DEs take the intermediate form

${\mathrm{d}}\,\vec{\mathcal{J}}_{F,S}(\vec{x},\epsilon)=\left({\mathrm{d}}\hat{A}^{(F,S)}_{0,\operatorname{Hom}}(\vec{x})+\epsilon\,{\mathrm{d}}\hat{A}^{(F,S)}_{1,\operatorname{Hom}}(\vec{x})\right)\cdot\vec{\mathcal{J}}_{F,S}(\vec{x},\epsilon)+(\text{sub-
sectors})\,,$ (35)

where $\hat{A}^{(F,S)}_{0,\operatorname{Hom}}$, and
$\hat{A}^{(F,S)}_{1,\operatorname{Hom}}$ are matrices of rational functions,
and $\hat{A}^{(F,S)}_{0,\operatorname{Hom}}$ is diagonal and non-zero only in
correspondence of the MIs which require a square-root normalisation. This
enables the straightforward use of finite fields techniques in the
reconstruction of DEs Peraro:2019svx .

We perform the previous steps sector by sector, starting from the lower
sectors and going up to the top sectors. Finally, we proceed with the last
step.

  * •

Step 3: We reconstruct analytically the DEs, this time including also the sub-
sectors contributions, with respect to the basis
$\vec{\mathcal{J}}_{F}=\cup_{S}\vec{\mathcal{J}}_{F,S}$, once again keeping
the square roots out of the computation. While for most of the sectors the
previous steps are enough to ensure that also the sub-sectors part of the DEs
is in $\epsilon$-factorised form, we found that some integral sectors need
further adjustments. In this case, in order to achieve an
$\epsilon$-factorised form, it is sufficient to modify the definition of the
MIs in the sector by including appropriate linear combination of the MIs of
the sub-sectors which are not in $\epsilon$-factorised form. The specific form
of the linear combinations is fixed by demanding that the DEs are
$\epsilon$-factorised off diagonal (e.g. see Gehrmann:2014bfa ).

As a result, we obtain DEs whose connection matrices are rational,
$\eps$-factorised off diagonal, and linear in $\eps$ on the diagonal. The
$\eps$-factorised form can then be obtained via the rotation
$\vec{\mathcal{I}}_{F}=N_{F}(\vec{x})\cdot\vec{\mathcal{J}}_{F}$, as

${\mathrm{d}}\,\vec{\mathcal{I}}_{F}(\vec{x},\eps)=\eps\left(N_{F}(\vec{x})\cdot{\mathrm{d}}\widehat{A}^{(F)}_{1,\operatorname{Hom}}(\vec{x})\cdot
N_{F}^{-1}(\vec{x})\right)\cdot\vec{\mathcal{I}}_{F}(\vec{x},\eps)\,,$ (36)

where $N_{F}(\vec{x})$ is a diagonal matrix which captures all square-root
normalisations and satisfies

${\mathrm{d}}\widehat{A}^{(F)}_{0,\operatorname{Hom}}+N_{F}^{-1}\cdot{\mathrm{d}}N_{F}=0\,.$
(37)

(a) The elliptic sector which we refer to as $321B$.

(b) Pentagon-triangle sector containing a nested square root. We refer to it
as $421B$

Figure 2: ‘Problematic’ sectors of topology ${\rm PB}_{B}$.

Through the strategy described above, we built a basis of MIs for topology
${\rm PB}_{C}$ which satisfies the DEs in canonical form. Contrarily, such a
form is not possible for topology ${\rm PB}_{B}$. Two sectors, shown in fig.
2, present additional challenges. These can be identified by analysing the
factorisation properties of the Picard-Fuchs operators Muller-Stach:2012tgj ;
Adams:2017tga . The sectors in fig. 2 are in fact the only ones whose MIs have
Picard-Fuchs operators with irreducible factors of degree $2$, in contrast
with all the other MIs, whose Picard-Fuchs operators factorise into linear
factors. We devote sections 3.3 and 3.5 to a thorough analysis of these
sectors, and summarise here the main conclusions. By analysing the homogeneous
DEs for the sector in fig. 2(a), we find that their solution involves elliptic
integrals. While the last few years have seen important progress in the
construction of $\eps$-factorised DEs beyond the standard
‘$\eps\times{\mathrm{d}}\log$’ case in eqs. (27) and (28) Frellesvig:2021hkr ;
Frellesvig:2023iwr ; Gorges:2023zgv , this problem is challenging in general.
Moreover, the transformation required to achieve an $\eps$-factorised form in
this case involves transcendental functions (such as elliptic ones). This
feature adds a further level of complexity for the numerical evaluation. In
other words, even if an $\eps$-factorised form of the DEs could be obtained in
this case, the numerical evaluation of the solution would remain an open
problem. The most common approach in such cases is to resort to semi-numerical
methods such as the generalised power series expansion, which can be equally
applied even without an $\eps$-factorised form. Nonetheless, we put some
effort into choosing MIs for this sector such that the connection matrices are
polynomial in $\eps$ up to order $\eps^{2}$, and all the entries which do not
involve MIs of the two problematic sectors are $\eps$-factorised. With respect
to the generic form, this makes the expression of the DEs more compact, and
improves the evaluation time of the solution using generalised power series
expansions.

The second problematic sector, shown in fig. 2(b), can be put into
$\eps$-factorised form, but at the cost of introducing a nested square root,
similar to the one encountered in ref. FebresCordero:2023gjh . While in
principle this is not a problem, in practice the available codes implementing
the generalised power series expansion method cannot handle a nested square
root. For this reason, we prefer to omit this transformation, and adopt a
basis of $3$ MIs for this sector such that the connection matrices have a
$2\times 2$ block with non-zero $\eps^{0}$ terms.

In light of the previous considerations, we built a basis of MIs for topology
${\rm PB}_{B}$ which satisfies a system of DEs as in eq. (25) with the
connection matrix ${\mathrm{d}}A^{({\rm PB}_{B})}$ of the form given in eq.
(29) with $k_{\rm max}=2$. Only the 27 entries of the connection matrix which
couple the differential of $\mathcal{I}_{{\rm PB}_{B},37}$ to MIs other than
itself are quadratic in $\eps$. These entries depend on a subset of the
kinematic invariants (see sec. 3.5). The entries which are instead non-zero at
$\eps=0$ are 16. Of these, 4 involve the MIs $\mathcal{I}_{{\rm PB}_{B},19}$
and $\mathcal{I}_{{\rm PB}_{B},20}$, and can be eliminated as discussed in
sec. 3.3 at the cost of introducing a nested square root. The remaining
entries involve at least one of the MIs of the elliptic sector (see sec. 3.5).

Before we move on to discussing the choice of MIs for the most complicated
sectors, we remark that another important aspect in the construction of the
integral bases is to minimise both the highest numerator rank and the quantity
of dotted propagators. This requirement is crucial in order to prevent the
size and number of IBP relations needed from exploding and thus making the
calculation computationally too expensive. To further ameliorate this aspect,
we used the software NeatIBP Wu:2023upw  to generate optimised IBP relations
through the solution of syzygy equations Gluza:2010ws .

###  3.1 Pentagon-box sectors

The eight-propagator pentagon-box sectors shown in figures 1(b) and 1(c)
contain three MIs for topology ${\rm PB}_{B}$ and four MIs for topology ${\rm
PB}_{C}$.

Regarding topology ${\rm PB}_{B}$, since the number of MIs is the same as in
the easier mass configurations, we then find that a canonical basis of MIs for
this sector is Badger:2022hno

$\displaystyle\begin{aligned} {\mathcal{I}}_{{\rm
PB}_{B},1}&=\epsilon^{4}\,d_{15}\,\operatorname{tr}_{5}I_{1,1,1,1,1,1,1,1}^{({\rm
PB}_{B}),[12],0,0,0}\,,\\\ {\mathcal{I}}_{{\rm
PB}_{B},2}&=\epsilon^{4}\,d_{15}\,\operatorname{tr}_{5}I_{1,1,1,1,1,1,1,1}^{({\rm
PB}_{B}),[11],0,0,0}\,,\\\ {\mathcal{I}}_{{\rm
PB}_{B},3}&=\epsilon^{4}\,d_{15}\,d_{23}\,d_{34}\left(I_{1,1,1,1,1,1,1,1}^{({\rm
PB}_{B}),1,0,0}+m_{t}^{2}\,I_{1,1,1,1,1,1,1,1}^{({\rm
PB}_{B}),0,0,0}\right)\,.\end{aligned}$ (38)

From numerical evaluations with AMFlow we observe that ${\mathcal{I}}_{{\rm
PB}_{B},1}$ vanishes up to order $\eps^{4}$.

Topology ${\rm PB}_{C}$ instead has four MIs in the top sector. We choose the
first three similarly to topology ${\rm PB}_{B}$. The construction of the
fourth MI is more complicated. We start from the scalar integral in
$d=6-2\eps$ dimensions, expressed in terms of integrals in $d=4-2\eps$
dimensions using LiteRed’s implementation of the dimension-shifting relations
Tarasov:1996br . With this choice, the DEs are linear in $\eps$. We then
construct a transformation of the fourth MI to eliminate the $\eps^{0}$ terms.
The resulting expression is however complicated, has rank-$4$ numerators, and
involves a spurious pole. We then search for a representation of this MI which
is free of these undesirable features by fitting an ansatz made of (at most)
rank-$2$ integrals on the top sector. This resulted in the following choices:

$\displaystyle\begin{aligned} {\mathcal{I}}_{{\rm
PB}_{C},1}&=\eps^{4}\,\,d_{34}\,d_{45}\left(d_{12}+m_{t}^{2}\right)I_{1,1,1,1,1,1,1,1}^{({\rm
PB}_{C}),1,0,0}\,,\\\ {\mathcal{I}}_{{\rm
PB}_{C},2}&=\eps^{4}\,\operatorname{tr}_{5}\,(d_{12}+m_{t}^{2})\,I_{1,1,1,1,1,1,1,1}^{({\rm
PB}_{C}),[11],0,0,0}\,,\\\ {\mathcal{I}}_{{\rm
PB}_{C},3}&=\eps^{4}\,\operatorname{tr}_{5}\,(d_{12}+m_{t}^{2})\,I_{1,1,1,1,1,1,1,1}^{({\rm
PB}_{C}),[12],0,0,0}\,,\\\ {\mathcal{I}}_{{\rm
PB}_{C},4}&=\eps^{4}\,\frac{\beta_{12}}{\operatorname{tr}_{5}}\,4\,d_{34}d_{45}(d_{12}+m_{t}^{2})\bigl{[}(d_{45}-d_{12}-m_{t}^{2})I_{1,1,1,1,1,1,1,1}^{({\rm
PB}_{C}),1,0,1}-d_{34}d_{45}I_{1,1,1,1,1,1,1,1}^{({\rm PB}_{C}),1,0,0}+\\\
&\phantom{{}={}}2d_{23}(d_{12}+m_{t}^{2})I_{1,1,1,1,1,1,1,1}^{({\rm
PB}_{C}),0,0,1}+(\text{sub-sectors})\bigr{]}\,.\end{aligned}$ (39)

The sub-sector terms in ${\mathcal{I}}_{{\rm PB}_{C},4}$ are rather lengthy,
and can be found in the ancillary files ancillary . Interestingly, numerical
evaluations with AMFlow show that three of the four MIs of this sector
(${\mathcal{I}}_{{\rm PB}_{C},2}$, ${\mathcal{I}}_{{\rm PB}_{C},3}$ and
${\mathcal{I}}_{{\rm PB}_{C},4}$) vanish up to order $\eps^{4}$.

###  3.2 Double-box sectors

There are four double-box sectors, two in topology ${\rm PB}_{B}$ and two in
topology ${\rm PB}_{C}$, as shown in fig. 3. The sectors (b), (c) and (d) in
fig. 3 have already been discussed in ref. Badger:2022hno . Sector (a) of
topology ${\rm PB}_{B}$ is instead new and contains six MIs. As for the
pentagon-box sectors, we were able to construct compact expressions for some
of the canonical MIs of this sector using $\mu_{ij}$ numerator insertions. The
first three integrals (${\mathcal{I}}_{{\rm PB}_{B},4}$, ${\mathcal{I}}_{{\rm
PB}_{B},5}$ and ${\mathcal{I}}_{{\rm PB}_{B},6}$) can be chosen as in sector
(b) in fig. 3 (see ref. Badger:2022hno ). Of the remaining three, we defined
two (${\mathcal{I}}_{{\rm PB}_{B},7}$ and ${\mathcal{I}}_{{\rm PB}_{B},8}$)
using $\mu_{ij}$ numerators and dotted propagators:

$\displaystyle\begin{aligned} {\mathcal{I}}_{{\rm
PB}_{B},4}&=\epsilon^{4}\,d_{15}\,d_{23}\,(d_{12}+m_{t}^{2})I_{1,1,1,0,1,1,1,1}^{({\rm
PB}_{B}),0,0,0}\,,\\\ {\mathcal{I}}_{{\rm
PB}_{B},5}&=\epsilon^{4}\,d_{15}\,d_{23}\,\left(I_{1,1,1,0,1,1,1,1}^{({\rm
PB}_{B}),0,1,0}+m_{t}^{2}\,I_{1,1,1,0,1,1,1,1}^{({\rm
PB}_{B}),0,0,0}\right)\,,\\\ {\mathcal{I}}_{{\rm
PB}_{B},6}&=\epsilon^{4}\,\operatorname{tr}_{5}\,\left(I_{1,1,1,0,1,1,1,1}^{({\rm
PB}_{B}),[12],0,0,0}+(\text{sub-sectors})\right)\,,\\\ {\mathcal{I}}_{{\rm
PB}_{B},7}&=\epsilon^{3}\,d_{23}\,\operatorname{tr}_{5}\,\left(I_{1,1,2,0,1,1,1,1}^{({\rm
PB}_{B}),[12],0,0,0}+(\text{sub-sectors})\right)\,,\\\ {\mathcal{I}}_{{\rm
PB}_{B},8}&=\epsilon^{3}\,d_{15}\,\operatorname{tr}_{5}\,I_{1,1,1,0,1,1,2,1}^{({\rm
PB}_{B}),[12],0,0,0}\,.\\\ \end{aligned}$ (40)

From numerical evaluations with AMFlow, we observe that ${\mathcal{I}}_{{\rm
PB}_{B},6}$ is zero up to order $\eps^{4}$. For the sixth MI of this sector
(${\mathcal{I}}_{{\rm PB}_{B},9}$), we could not find a compact
representation. We started from the derivative of ${\mathcal{I}}_{{\rm
PB}_{B},4}$ with respect to $d_{23}$, which leads to a linear dependence of
the connection matrices on $\eps$. We then constructed a transformation to
eliminate the $\eps^{0}$ part of the connection matrices, both in the
homogeneous part of the DEs and in the sub-sectors. This constraint amounts to
first-order DEs for the entries of the transformation matrix, which we could
solve in terms of rational functions. The resulting expression for
${\mathcal{I}}_{{\rm PB}_{B},9}$ is rather lengthy, and can be found in the
ancillary files ancillary  together with the sub-sector terms of
${\mathcal{I}}_{{\rm PB}_{B},6}$ and ${\mathcal{I}}_{{\rm PB}_{B},7}$.

(a)

(b)

(c)

(d)

Figure 3: The four five-point double-box topologies covering the MIs
${\mathcal{I}}_{{\rm PB}_{B},4}$ – ${\mathcal{I}}_{{\rm PB}_{B},9}$ (a),
${\mathcal{I}}_{{\rm PB}_{B},10}$ – ${\mathcal{I}}_{{\rm PB}_{B},12}$ (b),
${\mathcal{I}}_{{\rm PB}_{C},5}$ – ${\mathcal{I}}_{{\rm PB}_{C},8}$ (c) and
${\mathcal{I}}_{{\rm PB}_{C},15}$ – ${\mathcal{I}}_{{\rm PB}_{C},18}$ (d)
respectively.

###  3.3 Pentagon-triangle sector: a nested square root

There is only one pentagon-triangle sector, shown in figure 2(b) and dubbed
$421B$, and it appears in topology ${\rm PB}_{B}$. This sector has not been
studied previously in the literature and contains three MIs. One MI can be
chosen to take the same form as in the analogous topology of the five-point
integrals with an off-shell leg and massless internal propagators
Abreu:2020jxa . For the remaining two integrals, we made use of numerator
structures inspired by local numerators Arkani-Hamed:2010pyv  and written in
terms of Dirac traces, as

$\displaystyle\begin{aligned} {\mathcal{I}}_{{\rm
PB}_{B},18}&=\epsilon^{4}\,\operatorname{tr}_{5}\,I_{1,1,1,1,1,1,0,1}^{({\rm
PB}_{B}),[11],0,0,0}\,,\\\ {\mathcal{I}}_{{\rm
PB}_{B},19}&=\eps^{4}\,d_{45}\,I_{1,1,1,1,1,1,0,1}^{({\rm
PB}_{B}),0,0,0}\left[{\rm
tr}\bigl{(}\gamma_{5}\not{p}_{3}(\not{k}_{1}-\not{p}_{2}-\not{p}_{3})\not{p}_{4}\not{p}_{2}\bigr{)}\right]\,,\\\
{\mathcal{I}}_{{\rm PB}_{B},20}&=\eps^{4}\,d_{45}\,I_{1,1,1,1,1,1,0,1}^{({\rm
PB}_{B}),0,0,0}\left[{\rm
tr}\bigl{(}\not{p}_{3}(\not{k}_{1}-\not{p}_{2}-\not{p}_{3})\not{p}_{4}\not{p}_{2}\bigr{)}\right]\,,\end{aligned}$
(41)

where the terms in the square brackets are meant to be taken under the
integral sign. Note that the numerators in ${\mathcal{I}}_{{\rm PB}_{B},19}$
and ${\mathcal{I}}_{{\rm PB}_{B},20}$ are the parity even and odd parts of the
spinor chain $\langle 3|k_{1}-p_{2}-p_{3}|4]$, multiplied by an arbitrary
factor of $\langle 42\rangle[23]$ to cancel the helicity little-group scaling.
From numerical evaluations with AMFlow, we observe that ${\mathcal{I}}_{{\rm
PB}_{B},19}$ and ${\mathcal{I}}_{{\rm PB}_{B},20}$ start at order $\eps^{4}$,
whereas ${\mathcal{I}}_{{\rm PB}_{B},18}$ vanishes up to order $\eps^{4}$.

The choice of Dirac trace in the numerator is made in order to cancel the
potential singularity as the propagator momentum $k_{1}-p_{2}-p_{3}$ becomes
collinear to either adjacent massless leg $p_{3}$ or $p_{4}$. The two simple
numerator structures that can achieve this are $\langle
3|k_{1}-p_{2}-p_{3}|4]$ and its spinor conjugate $\langle
4|k_{1}-p_{2}-p_{3}|3]$. Setting these objects inside a trace ensures they are
free from any spinor phases. We choose ${\rm
tr}_{\pm}{[\not{p}_{3}(\not{k}_{1}-\not{p}_{2}-\not{p}_{3})\not{p}_{4}\not{p}_{2}]}$,
which have the same loop-momentum dependence but different normalisation. It
is then useful to split into parity odd and even pieces by taking linear
combinations. One can also write an alternative version of this ‘local’
numerator as $(l-l^{*})^{2}$, where $l=k_{1}-p_{2}-p_{3}$ and $l^{*}$ solves
the quadruple cut constraints

$\displaystyle\bigl{\\{}(l^{*})^{2}=0\,,\ (l^{*}-p_{4})^{2}=0\,,\
(l^{*}+p_{3})^{2}=0\,,\ (l^{*}+p_{2}+p_{3})^{2}=m_{t}^{2}\bigr{\\}}\,.$ (42)

As discussed in sec. 2, we remove the parity degree of freedom and replace the
parity-odd trace with the parity-even square root $\operatorname{tr}_{5}$
defined in eq. (18).

With the MIs in eq. (41), the homogeneous DEs for this sector are linear in
$\eps$, but not $\eps$-factorised. They take the form

$\displaystyle{\mathrm{d}}\,{\mathcal{I}}_{421B}=\left[\begin{pmatrix}0&0&0\\\
0&X&Y\\\
0&Y&X\end{pmatrix}+\mathcal{O}(\eps)\right]\cdot{\mathcal{I}}_{421B}+(\text{sub-
sectors})\,,\qquad{\mathcal{I}}_{421B}=\begin{pmatrix}{\mathcal{I}}_{{\rm
PB}_{B},18}\\\ {\mathcal{I}}_{{\rm PB}_{B},19}\\\ {\mathcal{I}}_{{\rm
PB}_{B},20}\end{pmatrix}\,.$ (43)

The omitted terms in the square brackets are proportional to $\eps$. All sub-
sector terms are $\eps$-factorised, except for the entries coupling
${\mathcal{I}}_{{\rm PB}_{B},19}$ and ${\mathcal{I}}_{{\rm PB}_{B},20}$ to the
elliptic sector (${\mathcal{I}}_{{\rm PB}_{B},i}$ for $i=35,36,37$). As
anticipated in the introduction to this section, the factorisation of $\eps$
in the entire diagonal block requires a transformation involving a nested
square root.

The higher complexity of this sector can be detected by analysing the
corresponding Picard-Fuchs operators as proposed in ref. Adams:2017tga .
First, we reduce the multi-scale problem to a single-scale one by defining a
univariate phase-space slice,

$\displaystyle x_{i}=a_{i}+b_{i}\,\lambda\,,\qquad\forall\,i=1,\ldots,6\,,$
(44)

with $a_{i},b_{i}\in\mathbb{Q}$, and viewing the MIs as functions of
$\lambda$. The constants $a_{i}$ and $b_{i}$ in eq. (44) are chosen randomly,
but one must make sure that no denominator factors of the connection matrices
vanish on the univariate slice, so as to avoid singular points. Furthermore,
we work modulo sub-topologies and modulo $\eps$-corrections, i.e. we focus on
the $3\times 3$ block of the DEs corresponding to this sector and set
$\eps=0$. The first MI of this sector, $\mathcal{I}_{{\rm PB}_{B},18}$,
decouples, as the corresponding DE-entries are already $\eps$-factorised. In
order to decouple $\mathcal{I}_{{\rm PB}_{B},19}$ and $\mathcal{I}_{{\rm
PB}_{B},20}$ we need to differentiate one more time. In other words, the
remaining $2\times 2$ block is equivalent to a second-order ordinary
differential equation in $\lambda$ for each of the integrals separately:

$\displaystyle L_{i}\,{\rm MaxCut}\left[\mathcal{I}_{{\rm
PB}_{B},i}\right]_{\eps=0}=0\,,\qquad\quad
L_{i}=\sum_{k=0}^{2}c_{i,k}(\lambda)\frac{{\mathrm{d}}^{k}}{{\mathrm{d}}\lambda^{k}}\,,\qquad\quad\forall\,i=19,20\,.$
(45)

The differential operators $L_{i}$ which annihilate the MIs are called Picard-
Fuchs operators. The factorisation properties of the Picard-Fuchs operators
encode useful information for the factorisation of $\eps$. Whenever they
factor completely into linear factors, the strategy of ref. Adams:2017tga
allows one to construct a transformation which puts the DEs in
$\eps$-factorised form. The Picard-Fuchs operators of $\mathcal{I}_{{\rm
PB}_{B},19}$ and $\mathcal{I}_{{\rm PB}_{B},20}$ are however second-order and
irreducible.333We used the Maple command DFactor to factorise the differential
operators Hoeij1997FactorizationOD . Note that this algorithm factors into
differential operators with rational function coefficients, and may therefore
miss factorisations involving algebraic coefficients. The appearance of an
irreducible factor of order greater than one in the factorisation of a Picard-
Fuchs operator is an indication that the Feynman integral cannot be expressed
in terms of MPLs Weinzierl:2022eaz . In this case, however, we find that the
solutions to the Picard-Fuchs equations in eq. (56) — equivalently, the
solution to the homogeneous DEs for this block — do not contain elliptic
integrals, but rather a nested square root. Indeed, we can actually put the
DEs in $\eps$-factorised form with an algebraic change of basis, as we discuss
below. It may therefore be possible to factor these Picard-Fuchs operators
into linear factors by allowing for algebraic functions in the coefficients.
In section 3.5 we will see that also the Picard-Fuchs operators of the
integrals of the sector shown in fig. 2(a) contain second-order irreducible
factors. In that case we will however find elliptic integrals in the
solutions. The Picard-Fuchs operators (modulo $\eps$ corrections and sub-
sectors) for all the other MIs of both topology ${\rm PB}_{B}$ and ${\rm
PB}_{C}$ are instead first-order.

We now proceed to put the homogeneous DEs for this sector in $\eps$-factorised
form. Thanks to the particularly symmetric form of the $2\times 2$ $\eps^{0}$
block corresponding to ${\mathcal{I}}_{{\rm PB}_{B},19}$ and
${\mathcal{I}}_{{\rm PB}_{B},20}$ (see eq. (43)), the off-diagonal non-zero
entries can be removed by simply replacing these MIs by their sum and
difference, i.e.

$\displaystyle{\mathcal{I}}_{421B}\quad\longrightarrow\quad{\mathcal{I}}^{\prime}_{421B}=T_{1}\cdot{\mathcal{I}}_{421B}\,,\qquad
T_{1}=\begin{pmatrix}1&0&0\\\ 0&1&1\\\ 0&1&-1\end{pmatrix}\,.$ (46)

Note that this transformation mixes integrals of different
$\operatorname{tr}_{5}$-charge (${\mathcal{I}}_{{\rm PB}_{B},19}$ is even,
${\mathcal{I}}_{{\rm PB}_{B},20}$ is odd), and gives $\langle
3|(k_{1}-p_{2}-p_{3})42|3]$ and its parity conjugate as numerators for
${\mathcal{I}}^{\prime}_{{\rm PB}_{B},19}$ and ${\mathcal{I}}^{\prime}_{{\rm
PB}_{B},20}$. The resulting homogeneous DEs have non-zero $\eps^{0}$ terms
only on the diagonal, which can thus be removed by proper normalisation. The
required normalisation factors are reciprocal of the solutions to the Picard-
Fuchs equations in eq. (56), and involve a nested square root:

$\displaystyle N_{\pm}=\frac{\sqrt{n_{\pm}}}{d_{45}\,r_{2}}\,,$ (47)

where

$\displaystyle n_{\pm}=d_{23}^{2}\operatorname{tr}_{5}^{2}-8r_{2}r_{4}\pm
4d_{23}r_{3}\operatorname{tr}_{5}\,,$ (48)

with

$\displaystyle\begin{aligned}
&r_{2}=2d_{23}(d_{23}+d_{34}-d_{15})+d_{34}m_{t}^{2}\,,\\\
&r_{3}=d_{23}\left[d_{12}(d_{23}-d_{15})-d_{23}d_{34}+(d_{34}-d_{15})d_{45}\right]+\left[d_{23}(d_{15}-d_{23}-2d_{34})+\
d_{34}d_{45}\right]m_{t}^{2}-d_{34}m_{t}^{4}\,,\\\
&r_{4}=2d_{12}d_{15}d_{23}d_{45}+\left[2d_{12}d_{23}(d_{23}-d_{15})-d_{34}(2d_{23}^{2}-2d_{23}d_{45}+d_{45}^{2})\right]m_{t}^{2}+\\\
&\
\phantom{r_{4}=}2(d_{23}^{2}+d_{34}d_{45}-d_{15}d_{23})m_{t}^{4}\,.\end{aligned}$
(49)

We emphasise that $n_{\pm}$ involves the square root $\operatorname{tr}_{5}$,
and that $n_{+}$ and $n_{-}$ are related by swapping the sign of
$\operatorname{tr}_{5}$. While the expression of the normalisation factors in
eq. (47) is fairly intricate, it is straightforward to verify that they do not
introduce spurious singularities. In other words, the product $n_{+}n_{-}$
factorises in terms of the same factors present in the denominators of the
connection matrices prior to the transformation. Furthermore, we note that the
factorisation of expressions involving square roots is not unique, and
therefore a simpler representation of these normalisation factors may exist.
In particular, we find that this nested square root can also be written
compactly in terms of traces of gamma matrices,

$\displaystyle\begin{aligned} N_{\pm}=\sqrt{8\left(\frac{{\rm
tr}_{\mp}(\not{p}_{3}\not{p}_{2}\not{p}_{12}\not{p}_{5})}{{\rm
tr}_{\mp}(\not{p}_{4}\not{p}_{5}\not{p}_{3}\not{p}_{2})}+\frac{4m_{t}^{2}}{s_{45}}\right)\frac{{\rm
tr}_{\mp}(\not{p}_{3}\not{p}_{2}\not{p}_{12}\not{p}_{5})}{{\rm
tr}_{\mp}(\not{p}_{4}\not{p}_{5}\not{p}_{3}\not{p}_{2})}}\,,\end{aligned}$
(50)

in which it is clear that the interior square root could be rationalised when
using a momentum-twistor representation Hodges:2009hk ; Badger:2022mrb  with
rational parameterisations for the spinor products, and that the argument of
the outer square root becomes a perfect square in the massless limit. We have
checked that $n_{\pm}$ cannot be expressed as a perfect square of the form
$(a+b\operatorname{tr}_{5})^{2}$, for some rational functions $a$ and $b$,
which would allow us to remove the exterior square root. Finally, the
normalisation factors in eq. (47) do not have well-defined behaviour under
swapping the sign of $\operatorname{tr}_{5}$. It is instead desirable that all
MIs are either even or odd with respect to this operation. We therefore apply
another transformation of the form of eq. (46) to restore this property. In
conclusion, the basis of this sector which puts the DEs in $\eps$-factorised
form is

$\displaystyle{\mathcal{I}}^{\prime\prime}_{421B}=\begin{pmatrix}1&0&0\\\
0&1&1\\\ 0&1&-1\end{pmatrix}\cdot\begin{pmatrix}1&0&0\\\ 0&N_{+}&0\\\
0&0&N_{-}\end{pmatrix}\cdot\begin{pmatrix}1&0&0\\\ 0&1&1\\\
0&1&-1\end{pmatrix}\cdot{\mathcal{I}}_{421B}\,.$ (51)

Keeping in mind that we did not put in $\eps$-factorised form the elliptic
sector (see sec. 3.5), and that DiffExp cannot handle nested square roots, we
prefer to omit this transformation from our chosen basis of MIs.

###  3.4 Pentagon-bubble sectors

There are two sectors in the form of a pentagon with a bubble insertion, shown
in figure 4: one in topology ${\rm PB}_{B}$, with two MIs, and one in topology
${\rm PB}_{C}$, with three MIs.

(a)

(b)

Figure 4: The two pentagon-bubble sectors for topologies ${\rm PB}_{B}$ (sub-
figure (a)) and ${\rm PB}_{C}$ (sub-figure (b)).

Both sectors present a new mass configuration. However, since the number of
MIs for the pentagon-bubble integrals in topology ${\rm PB}_{B}$ is the same
as in previously studied cases, we can make a similar choice for the canonical
basis:

$\displaystyle\begin{aligned} {\mathcal{I}}_{{\rm
PB}_{B},45}&=\epsilon^{3}\,\operatorname{tr}_{5}\,I_{1,1,1,1,0,1,0,2}^{({\rm
PB}_{B}),[11],0,0,0}\,,\\\ {\mathcal{I}}_{{\rm
PB}_{B},46}&=\epsilon^{3}(1-2\epsilon)\,d_{23}\,d_{34}I_{1,1,1,1,0,1,0,1}^{({\rm
PB}_{B}),0,0,0}\,.\end{aligned}$ (52)

Regarding the pentagon-bubble sector in topology ${\rm PB}_{C}$, we find a
canonical basis where all three MIs involve dotted propagators:

$\displaystyle\begin{aligned} {\mathcal{I}}_{{\rm
PB}_{C},32}&=\epsilon^{3}\,\operatorname{tr}_{5}\,I_{1,1,1,1,0,2,0,1}^{({\rm
PB}_{C}),[11],0,0,0}\,,\\\ {\mathcal{I}}_{{\rm
PB}_{C},33}&=\epsilon^{3}\,\operatorname{tr}_{5}\,I_{1,1,1,1,0,1,0,2}^{({\rm
PB}_{C}),[11],0,0,0}\,,\\\ {\mathcal{I}}_{{\rm
PB}_{C},34}&=\epsilon^{3}\,d_{34}\,d_{45}\left(I_{1,1,1,1,0,2,0,1}^{({\rm
PB}_{C}),1,0,0}+m_{t}^{2}I_{1,1,1,1,0,2,0,1}^{({\rm
PB}_{C}),0,0,0}\right)\,.\end{aligned}$ (53)

###  3.5 Elliptic sector

The most complicated sector belongs to topology ${\rm PB}_{B}$ and is shown in
fig. 2(a). It is a four-point sector, and thus its integrals depends on four
variables only ($d_{12}$, $d_{34}$, $d_{15}$ and $m_{t}^{2}$). There are $3$
MIs, which we choose as

$\displaystyle\begin{aligned} {\mathcal{I}}_{{\rm
PB}_{B},35}&=\eps^{4}\,d_{15}\,(d_{12}+m_{t}^{2})\,I_{1,1,0,1,1,1,0,1}^{({\rm
PB}_{B}),0,0,0}\,,\\\ {\mathcal{I}}_{{\rm
PB}_{B},36}&=\eps^{4}\,\sqrt{(d_{15}-d_{34})^{2}-2d_{34}m_{t}^{2}}\,I_{1,1,0,1,1,1,0,1}^{({\rm
PB}_{B}),0,0,0}\left[2k_{1}\cdot p_{1}\right]\,,\\\ {\mathcal{I}}_{{\rm
PB}_{B},37}&=\eps^{4}\,(-1+2\eps)\,d_{15}\,I_{1,1,0,1,1,1,0,1}^{({\rm
PB}_{B}),0,0,0}\left[2k_{2}\cdot p_{2}\right]\,,\end{aligned}$ (54)

where we recall that the terms in the square brackets are meant to be taken
under the integral sign. From numerical evaluations with AMFlow, we observe
that the chosen MIs of this sector are non-zero only starting from order
$\eps^{4}$. While the DEs for this sector are not $\eps$-factorised, the MIs
above simplify them substantially with respect to an arbitrary choice. The
normalisation factor of ${\mathcal{I}}_{{\rm PB}_{B},36}$ is chosen so as to
factorise $\eps$ in the corresponding diagonal entry of the DEs. The factor of
$(-1+2\eps)$ in ${\mathcal{I}}_{{\rm PB}_{B},37}$ is inserted to remove all
$\eps$-dependent factors from the denominators of the connection matrices. The
remaining kinematic-dependent normalisation factors of ${\mathcal{I}}_{{\rm
PB}_{B},35}$ and ${\mathcal{I}}_{{\rm PB}_{B},37}$ ensure that all MIs in the
basis have the same dimensionality, and are chosen heuristically as they lead
to more compact connection matrices. The DEs have the following structure:

$\displaystyle{\mathrm{d}}\,\mathcal{I}_{321B}=\begin{pmatrix}*+\eps\,*&\eps\,*&*\\\
*+\eps\,*&\eps\,*&*\\\
(1-2\eps)(*+\eps\,*)\qquad&(1-2\eps)\eps\,*\qquad&*+\eps\,*\\\
\end{pmatrix}\cdot\,\mathcal{I}_{321B}+(\text{sub-sectors})\,,$ (55)

where $\mathcal{I}_{321B}=\left({\mathcal{I}}_{{\rm
PB}_{B},35},\,{\mathcal{I}}_{{\rm PB}_{B},36},\,{\mathcal{I}}_{{\rm
PB}_{B},37}\right)^{\top}$, and each asterisk denotes a distinct one-form. The
sub-sectors follow the same pattern as the diagonal block shown above: the
entries coupling ${\mathcal{I}}_{{\rm PB}_{B},35}$ and ${\mathcal{I}}_{{\rm
PB}_{B},36}$ to the sub-sectors are linear in $\eps$, while for
${\mathcal{I}}_{{\rm PB}_{B},37}$ they are quadratic.

The analysis of the Picard-Fuchs operators $L_{i}$ allows us to better
characterise the complexity of this sector. Following the procedure outlined
in sec. 3.3, we construct the differential operators $L_{i}$ such that

$\displaystyle L_{i}\,{\rm MaxCut}\left[\mathcal{I}_{{\rm
PB}_{B},i}\right]_{\eps=0}=0\,,\qquad\quad
L_{i}=\sum_{k=0}^{r_{i}}c_{i,k}(\lambda)\frac{{\mathrm{d}}^{k}}{{\mathrm{d}}\lambda^{k}}\,,\qquad\quad\forall\,i=35,36,37\,,$
(56)

with $r_{35}=r_{37}=2$ and $r_{36}=3$. We recall that we work on a random
univariate phase-space slice (see eq. (44)), modulo sub-sectors and
$\eps$-corrections. $L_{35}$ and $L_{37}$ are second-order, irreducible
operators. The solution to $L_{35}$ involves the elliptic integral of the
first kind ${\rm K}(x)$, while the solution to $L_{37}$ also contains
derivatives of the latter.444We obtained the solutions to the second-order
Picard-Fuchs operators with the Maple command hypergeometricsols. We thank
Christoph Dlapa for suggesting this. $L_{36}$ is instead third-order, and
factorises into the product of a second- and a first-order operator. The
first-order factor is simply ${\mathrm{d}}/{\mathrm{d}}\lambda$, which follows
from the fact that the reciprocal of the chosen normalisation factor of
$\mathcal{I}_{{\rm PB}_{B},36}$ is a solution to $L_{36}$. In appendix A we
determine the elliptic curve underlying the elliptic integrals appearing in
this analysis, and identify the solutions to the Picard-Fuchs operator
$L_{35}$ with the periods of said elliptic curve.

##  4 ‘${\mathrm{d}}\log$’ and ‘one-form’ representation of the differential
equations

Having obtained the connection matrices of the differential equations using
the integral bases described in the previous sections, we present compact
analytic expressions in terms of independent ‘${\mathrm{d}}\log$’ and ‘one-
form’ structures, as in eqs. (28) and (29). There are a variety of methods
proposed in the literature to determine the alphabet of ‘${\mathrm{d}}\log$’
forms where they exist. In our case, rather than constructing an ansatz of
possible alphabet letters a priori, we first established a set of linearly
independent one-forms organised according to the square root charges defined
earlier in section 3. As discussed there, each entry of the connection
matrices after the reconstruction is a rational function of the invariants,
$\vec{x}$. We then add back the square root normalisations of each master
integral to rotate the DEs into $\eps$-factorised form (or as far as possible,
in the case of topology ${\rm PB}_{B}$). After this stage each entry takes the
form

$\displaystyle\begin{aligned}
{\mathrm{d}}A^{(F)}_{ij}(\vec{x},\eps)&=\sum_{a}\eps^{a}{\mathrm{d}}A^{(F),a}_{ij}(\vec{x})\,,\\\
{\mathrm{d}}A^{(F),a}_{ij}(\vec{x})&=\mathcal{C}^{(F)}_{ij}(\vec{x})\sum_{k=1}^{6}{\mathrm{d}}x_{k}\,f^{(F),a}_{ijk}(\vec{x})\,,\end{aligned}$
(57)

where $\mathcal{C}_{ij}^{(F)}(\vec{x})$ are monomials in the possible square
roots (including $\mathcal{C}_{ij}^{(F)}(\vec{x})=1$) that we have called
charges (see sec. 3), and $f^{(F),a}_{ijk}$ are rational functions of the
invariants $\vec{x}$. The sum over the orders in $\eps$ (index $a$) runs from
$0$ to $2$ for ${\rm PB}_{B}$, while ${\rm PB}_{A}$ and ${\rm PB}_{C}$ are in
$\eps$-factorised form (i.e., $a=1$). There are 49 charges in total, with 23
appearing in the leading colour topologies.

At this stage, we determine the linear relations amongst each set of entries
${\mathrm{d}}A^{(F),a}_{ij}$ and of ${\mathrm{d}}\log W$’s that share the same
charge $\mathcal{C}^{(F)}_{ij}$, where we use the notation $W$ to indicate the
letters of our alphabet. Since the square roots are an overall factor, we can
divide them out and use finite field techniques to determine the linear
relations. We then solve the linear relations by ordering with respect to the
polynomial degree of the entries ${\mathrm{d}}A^{(F),a}_{ij}$ and preferring
${\mathrm{d}}\log$ structures. After iterating over all possible charges, all
entries ${\mathrm{d}}A^{(F),a}_{ij}$ are expressed in terms ${\mathrm{d}}\log$
structures of our alphabet if possible (as in eq. (28)), and additionally a
minimal set of simple one-forms otherwise (as in eq. (29), with some
$\omega_{i}$’s being ${\mathrm{d}}\log$’s). For concreteness, we spell out the
form of the connection matrix for ${\rm PB}_{B}$:

$\displaystyle{\mathrm{d}}A^{({\rm
PB}_{B})}(\vec{x},\eps)=\sum_{k=0}^{2}\eps^{k}\Biggl{[}\sum_{i}\,c^{({\rm
PB}_{B})}_{k,i}\,{\mathrm{d}}\log\left(W_{i}(\vec{x})\right)+\sum_{j}\,d^{({\rm
PB}_{B})}_{k,j}\,\omega_{j}(\vec{x})\Biggr{]}\,,$ (58)

where $c^{({\rm PB}_{B})}_{k,i}$ and $d^{({\rm PB}_{B})}_{k,j}$ are matrices
of rational numbers. We stress that the ${\mathrm{d}}\log$’s are one-forms
themselves, but in this context we call one-forms only those we could not
express in terms of ${\mathrm{d}}\log$’s. Furthermore, we emphasise that some
of these one-forms are not closed, which implies they are not exact
either.555A differential form $\omega$ is closed if ${\mathrm{d}}\omega=0$. A
one-form $\omega$ is exact if there exists a function $f$ such that
$\omega={\mathrm{d}}f$. An exact form is thus by definition closed. In other
words, they cannot be expressed as the differential of any function, let alone
of a logarithm. Table 2 shows a summary of each topology.

The rational letters (i.e. those with $\mathcal{C}^{(F)}_{ij}=1$) are easy to
determine: they are the (algebraically independent) factors in the
denominators of the connection matrices. Some of these letters are extremely
simple, linear combinations of invariants, while others are high degree
polynomials in the invariants $\vec{x}$ (up to degree $5$). We determine that
a large number of these high degree polynomials can be identified as Gram and
Cayley determinants of the external kinematics. Others can be written in
compact notation by using traces of gamma matrices. In total we find 58
rational letters for all one-loop and planar two-loop cases, although only 46
appear in the leading colour topologies.

The algebraic letters are all written in the manifestly odd form

$\frac{A+\mathcal{C}}{A-\mathcal{C}}\,,$ (59)

where $\mathcal{C}$ is one of the charges (apart from $1$) and $A$ is a
rational function. There are a variety of methods proposed to determine the
form of such letters, for example Zoia:2021zmb ; Heller:2019gkq ;
FebresCordero:2023gjh ; Jiang:2024eaj ; Henn:2024ngj . In our case, we were
able to find all the required expressions by comparison with ansätze for the
rational function $A$. Rather than printing the list of charges and letters
here we present all the relevant definitions in the ancillary files ancillary
described in appendix C. In the case of $\mathcal{C}=\operatorname{tr}_{5}$,
we find that compact representations can be found by using ratios of
${\operatorname{tr}_{+}}/{\operatorname{tr}_{-}}$ with arguments of either
four or six gamma matrices. Explicit forms are given in the ancillary files,
although this is only an aesthetic consideration.

In some cases, having found ${\mathrm{d}}\log$ representations, we have
subsequently imposed symmetries amongst the letters. This is particularly
relevant for letters related to charges involving the three-mass triangle Gram
determinants, $\Delta_{3,i}$. We note that the charge $\Delta_{3,4}$ appears
in topologies ${\rm PB}_{A}$ and ${\rm PB}_{C}$ with 3 independent one-forms.
We write the one-forms in terms of 4 letters in order to preserve these
symmetries. This would be convenient in case one was looking for a set of
special functions closed under permutations. In addition, for ${\rm PB}_{A}$
there are two additional letters containing $\operatorname{tr}_{5}$ in our
alphabet than independent entries of the connection matrix.

As a final remark, we notice that the mixed ${\mathrm{d}}\log$ and one-form
expression of the DEs for topology ${\rm PB}_{B}$ allows us to clearly
separate the features in common with the standard canonical cases (${\rm
PB}_{A}$ and ${\rm PB}_{C}$) from the new, more complicated ones. Moreover, it
makes the expression of the DEs roughly $50$ times more compact than if one
stored the connection matrices separately for each derivative, as in the form
generated by DiffExp for the semi-numerical solution. It would be interesting
to find a way to exploit this simplification of the form of the DEs in the
method of generalised power series expansions.

topology | charges |  ${\mathrm{d}}\log$ letters | one-forms
---|---|---|---
${\rm PB}_{A}$ | 17 | 74 | 0
${\rm PB}_{B}$ | 16 | 72 | 63
${\rm PB}_{C}$ | 21 | 80 | 0
2-loop LC | 23 | 98 | 63
${\rm P}_{A}$ | 8 | 41 | 0
${\rm P}_{B}$ | 10 | 45 | 0
${\rm P}_{C}$ | 17 | 55 | 0
${\rm P}_{D}$ | 29 | 64 | 0
1-loop | 40 | 105 | 0
Table 2: Summary of the ‘${\mathrm{d}}\log$’ and ‘one-form’ structures
appearing in the differential equations for each topology, and cumulatively
for all two-loop leading colour (LC) topologies and all 1-loop topologies.
${\rm PB}_{A}$ and ${\rm PB}_{C}$ contain three and one more letter than
linearly independent entries of the connection matrices respectively. This is
due to the imposed symmetries described in the text.

##  5 Numerical evaluation using generalised series expansions

In this section, we discuss the numerical evaluation of the master integrals.
We obtain a semi-analytic solution to the system of differential equations
associated with the master integrals by means of the method of generalised
power series expansions Francesco:2019yqt , by exploiting the Mathematica
package DiffExp Hidding:2020ytt . We obtain the required boundary values
numerically with the package AMFlow Liu:2022chg , which implements the
auxiliary mass flow method Liu:2017jxz ; Liu:2021wks ; Liu:2022tji . We aim to
evaluate the integrals in the physical scattering region relevant for
phenomenology. We begin by defining this region. We then motivate our choice
for the boundary point, and discuss a number of interesting features of the
boundary values. Next, we present a number of checks we performed to validate
our results. Finally, we comment on the performance of the numerical
evaluation.

###  5.1 Physical scattering region

We restrict our analysis to the physical phase-space region corresponding to
the $s_{45}$ scattering channel ($45\rightarrow 123$). All other $2\to 3$
channels relevant for $t\bar{t}j$ production can be obtained from the $s_{45}$
channel through suitable permutations of the momenta, and our numerical
evaluation procedure can therefore be straightforwardly generalised to them as
well. The $s_{45}$ channel is defined by the following linear constraints on
the kinematic invariants,

$\displaystyle\begin{gathered}p_{1}^{2}>0\,,\quad p_{2}^{2}>0\,,\quad
p_{4}\cdot p_{5}>0\,,\quad p_{1}\cdot p_{2}>0\,,\quad p_{1}\cdot
p_{3}>0\,,\quad p_{2}\cdot p_{3}>0\,,\\\ p_{4}\cdot p_{1}<0\,,\quad p_{4}\cdot
p_{2}<0\,,\quad p_{4}\cdot p_{3}<0\,,\quad p_{5}\cdot p_{1}<0\,,\quad
p_{5}\cdot p_{2}<0\,,\quad p_{5}\cdot p_{3}<0\,,\end{gathered}$ (62)

complemented by the following higher-order constraints coming from Gram
determinants,

$\displaystyle\mathrm{det}\,G(p_{i},p_{j})<0\,,\qquad\mathrm{det}\,G(p_{i},p_{j},p_{k})>0\,,\qquad\mathrm{det}\,G(p_{i},p_{j},p_{k},p_{l})<0\,,$
(63)

where $i,j,k,l$ take distinct values in $\\{1,\ldots,5\\}$. The Gram
determinants involving two and four momenta give only one constraint each:

$\displaystyle(d_{12}-m_{t}^{2})(d_{12}+m_{t}^{2})>0\,,\qquad\qquad\operatorname{tr}_{5}^{2}<0\,.$
(64)

Those involving three momenta give a number of polynomial constraints on the
kinematic invariants which we do not spell out. All constraints defining the
$s_{45}$ channel can be found in ancillary files ancillary .

###  5.2 Boundary values

The evaluation of the master integrals through the method of generalised power
series expansions entails the integration of the DEs along a path connecting
the target point with a starting point at which the values are known. We wish
to confine such paths to the $s_{45}$ channel in order to avoid analytic
continuation, which is non-trivial to determine for multi-variable problems
and increases the evaluation time. In this view, we choose a boundary point
$\vec{x}_{0}$ in the $s_{45}$ channel. The choice is arbitrary, but has an
impact on the performance of the evaluation. We choose

$\vec{x}_{0}=\bigl{\\{}2,1,-1,5,-2,1\bigr{\\}}\,,$ (65)

following the criteria set in ref. Chicherin:2021dyp , which we recall here.

  1. 1.

The point $\vec{x}_{0}$ is invariant under the symmetries of the $s_{45}$
channel, i.e. the exchanges of the external momenta $p_{1}\leftrightarrow
p_{2}$ and $p_{4}\leftrightarrow p_{5}$.

  2. 2.

The point $\vec{x}_{0}$ introduces a minimal number of distinct prime factors.

  3. 3.

The point $\vec{x}_{0}$ lies on the spurious singularity $d_{23}+d_{34}=0$.

The first two criteria reduce the number of independent transcendental
constants in the values of the master integrals at $\vec{x}_{0}$. This is
useful in view of a future analytic solution of the DEs. The second criterion
also makes the numerical evaluation using AMFlow faster, as it reduces the
number of prime fields required to reconstruct the DEs with respect to the
auxiliary mass. We recall in fact that we interface AMFlow to FiniteFlow in
order to solve the IBP relations over finite fields.

In order to understand the third criterion, we first need to define spurious
singularities. A _spurious singularity_ is a singularity of the connection
matrix of the DEs which is not a singularity of the solution to the DEs once
the boundary values are taken into account. The possible singularities of the
MIs correspond to the factors in the denominators of the connection matrices.
For DEs where a canonical ‘${\mathrm{d}}\log$ form’ is possible, this is
equivalent to the rational letters of the alphabet. We find that two of them,
$d_{23}+d_{34}$ and $d_{12}+d_{15}+m_{t}^{2}$, can vanish within the $s_{45}$
channel.666We can prove analytically that most of the denominator factors have
fixed sign in the $s_{45}$ channel. For a few, we have statistical evidence
based on $100$K random phase-space points generated by sampling uniformly a
parameterisation of the momenta in terms of energies and angles. By choosing
the boundary point such that $d_{23}+d_{34}=0$, we ensure that this spurious
singularity is never crossed by a straight path starting from $\vec{x}_{0}$.
This improves the speed of the integration with the generalised power series
expansion methods, and is useful in view of a future solution in terms of one-
fold integral representations along the lines of ref. Caron-Huot:2014lda . We
emphasise however that the presence of spurious singularities is not an
obstacle for the generalised power series expansion method.

We observe that, for a small fraction of phase-space points in the $s_{45}$
channel, the straight line connecting them to $\vec{x}_{0}$ leaves the
$s_{45}$ channel. A similar structure of the physical phase space is discussed
in ref. Chicherin:2021dyp . In such a case, one can either perform the
analytic continuation, or choose a different starting point.

We use AMFlow to obtain the values of all MIs at $\vec{x}_{0}$ with, at least,
$32$-digit precision. This gives an upper bound on the achievable precision in
the numerical evaluation of the MIs using the results provided in this work.
We expect this level of precision to be sufficient for the application in
phenomenology at NNLO in QCD, based on the available experience with two-loop
scattering amplitudes with kinematics of similar complexity. However, higher
precision can be easily achieved, if needed, as evaluating the MIs in a single
point with AMFlow does not represent a bottleneck for this computation. We
provide the boundary values in the supplementary material ancillary .

For the topologies ${\rm PB}_{A}$ and ${\rm PB}_{C}$, the canonical form of
the DEs implies that the values at order $\eps^{0}$ are rational, and can be
determined up to the overall normalisation by imposing ‘first-entry’
conditions Gaiotto:2011dt . In other words, only logarithms of the following
arguments can appear at order $\eps^{1}$,

$\displaystyle\bigl{\\{}m_{t}^{2}\,,\ 2(d_{12}+m_{t}^{2})\,,\ 2d_{23}\,,\
2d_{34}\,,\ 2d_{45}\,,\ 2d_{15}\bigr{\\}}\,,$ (66)

and this gives linear constraints on the $\mathcal{O}(\eps^{0})$ boundary
values. The allowed logarithm arguments can be read off from the graph
polynomial $\mathcal{F}$. We then fix the overall normalisation by
rationalising the values from AMFlow. Similarly, at order
$\mathcal{O}(\eps^{1})$ we verify using the PSLQ algorithm PSLQ  that the
values are $\mathbb{Q}$-linear combinations of logarithms of the functions in
eq. (66) evaluated at $\vec{x}_{0}$.

As for topology ${\rm PB}_{B}$, although the DEs are not in canonical form, we
still observe that the boundary values are rational at order $\eps^{0}$, and
linear combinations of the logarithms above at order $\eps^{1}$. This follows
from the fact that the ‘problematic’ MIs, namely those of the sector involving
the nested square root ($\mathcal{I}_{{\rm PB}_{B},19}$ and $\mathcal{I}_{{\rm
PB}_{B},20}$, see sec. 3.3) and those of the sector involving elliptic
integrals ($\mathcal{I}_{{\rm PB}_{B},35}$, $\mathcal{I}_{{\rm PB}_{B},36}$
and $\mathcal{I}_{{\rm PB}_{B},37}$, see sec. 3.5), are non-zero only starting
from order $\eps^{4}$.

###  5.3 Checks

The rationality of the boundary values at order $\eps^{0}$ and the first-entry
conditions discussed in the previous subsection are already non-trivial checks
of our results. In order to validate more robustly the numerical evaluation of
the MIs, we compared the numerical values obtained by integrating the DEs with
DiffExp starting from $\vec{x}_{0}$ against numerical evaluations performed
with AMFlow at a number of points in the $s_{45}$ channel. We found full
agreement within the accuracy estimated by DiffExp and AMFlow. In particular,
we perform this check at the following random point,

$\vec{x}_{1}=\left\\{\frac{4602}{57095},\frac{217}{8151},-\frac{8513}{67193},\frac{7}{22},-\frac{14291}{77626},\frac{1701}{90164}\right\\}\,,$
(67)

which lies in the $s_{45}$ channel on the other side of the spurious
singularity $d_{12}+d_{15}+m_{t}^{2}=0$ with respect to the boundary point
$\vec{x_{0}}$. This allows us to check the stability of the numerical
evaluation when integrating the DEs over a path that crosses this spurious
singularity. We provide the values of the master integrals at $\vec{x}_{1}$ as
benchmarks in ancillary files ancillary .

###  5.4 Performance analysis

We finish this section by making some considerations regarding the performance
of the numerical evaluation of the MIs. First of all, we want to clarify that
our evaluation strategy is not optimised for phenomenological applications,
and we therefore refrain from making absolute statements about the evaluation
time. The latter in fact depends strongly on the segmentation of the path
within the generalised power series expansion method Francesco:2019yqt ;
Hidding:2020ytt . The number of segments, in turn, depends on the chosen
endpoints of the path, and on the location of the nearest singularities. An
evaluation strategy aimed at a large number of points should therefore
minimise the number of segments in the evaluations by re-using iteratively the
values obtained with previous evaluations (see e.g. ref. Abreu:2020jxa ). We
leave this to future work.

Nonetheless, it is still interesting to compare the relative performance of
the evaluation of the MIs of topologies ${\rm PB}_{A}$, ${\rm PB}_{C}$ and
${\rm PB}_{B}$, as the $\epsilon$ structure of the corresponding DEs is
different. We recall that for topology ${\rm PB}_{A}$ and ${\rm PB}_{C}$ we
have canonical DEs where $\eps$ is factorised, whereas the entries of the
connection matrix of topology ${\rm PB}_{B}$ are degree-2 polynomials in
$\eps$. We can estimate the impact of the additional terms in $\epsilon$ for
topology ${\rm PB}_{B}$ by comparing the evaluation time per segment between
topologies for the same values of all parameters. We performed this analysis
on a sample of $100$ phase-space points in the $s_{45}$ channel, amounting to
$\approx 1.5$ K segments starting from the boundary point $\vec{x}_{0}$, with
a target accuracy of $10^{-16}$ in DiffExp. The results of this analysis are
shown in figure 5.777All the evaluations are performed on an Intel(R) Xeon(R)
Gold 5218 2.30 GHz CPU. The different $\epsilon$ structure for the DEs,
together with the larger number of MIs, results in a higher evaluation time of
${\rm PB}_{B}$ with respect to ${\rm PB}_{A}$ and ${\rm PB}_{C}$. Indeed, the
average evaluation time per segment is $\approx 12$ seconds for ${\rm
PB}_{A}$, $\approx 16$ seconds for ${\rm PB}_{C}$, and $\approx 58$ seconds
for ${\rm PB}_{B}$.

Figure 5: Histogram showing the distribution of the evaluation time per
segment for topologies ${\rm PB}_{A}$, ${\rm PB}_{B}$ and ${\rm PB}_{C}$. Note
that the number of segments depends on the singularity structure and is
therefore different for each topology.

We finish this section with some comments. First, the data shown in figure 5
have to be intended as a qualitative analysis and not as a serious attempt of
performing a statistical study of the evaluation time performance of the
method. Furthermore, if we take into account the earlier analysis, bearing in
mind that our evaluation strategy is not tailored for phenomenological
applications, we find the current results to be promising. Finally, it would
be interesting to compare the evaluation time of our solution with a possible
generalised power series implementation for the case where the connection
matrix contains elliptic kernels. Indeed, assuming to be able to build an
integral basis in which the DEs are factorised also in the elliptic case, the
resulting connection matrix will then depend on elliptic functions. Therefore,
it will be interesting to weigh whether is the polynomial dependence on
$\epsilon$, or the presence of elliptic functions, to have the most
significant impact on the numerical evaluation of the MIs.

##  6 Conclusion

In this work, we have presented compact differential equations for the master
integrals of all integral topologies required to describe the production of a
pair of top quarks in association with a jet at hadron colliders at NNLO in
leading colour QCD. There were two new pentagon-box topologies to consider,
${\rm PB}_{B}$ (fig. 1(b)) and ${\rm PB}_{C}$ (fig. 1(c)). The latter followed
a pattern similar to the one observed in ref. Badger:2022mrb  for ${\rm
PB}_{A}$ (fig. 1(a)) and we obtain the DEs in the canonical form, i.e., such
that the connection matrices are given by an overall factor of $\eps$ and
$\mathbb{Q}$-linear combinations of logarithmic one-forms
(${\mathrm{d}}\log$’s). We showed that the former, ${\rm PB}_{B}$, has new
features that prevent it from having the same canonical form. While the one-
loop topologies were previously studied and presented in $\eps$-factorised
form in ref. Badger:2022mrb , we present them here for completeness in terms
of ${\mathrm{d}}\log$’s.

Topology ${\rm PB}_{B}$ displays two new features. First, the sector shown
graphically in fig. 2(b) involves a nested square root. Secondly, the sector
in fig. 2(a) involves elliptic integrals, which we identify as the periods of
the elliptic curve in eq. (75). Obtaining canonical DEs in the presence of
elliptic integrals is at the forefront of current research and the very notion
of ‘canonical’ in such cases is still under debate. From a practical view
point, even if such a generalised canonical form is achieved, the numerical
evaluation of the solution remains challenging, and so an alternative route to
efficient and stable numerical evaluation was taken. We obtain a compact
representation of the differential equation in this case by making
${\mathrm{d}}\log$ choices for all master integrals except those in the
complicated sectors. For the problematic sectors, we find choices in which the
differential equation is at most quadratic in $\eps$. This form is compatible
with rational reconstruction over finite fields using optimised IBP relations
in the same fashion as the other topologies and with easily manageable
computation times. By using a set of linearly independent ${\mathrm{d}}\log$
and non-logarithmic one-forms we express the DEs in a compact form where the
number of non-logarithmic structures is minimised.

We evaluate the master integrals numerically by solving the corresponding DEs
by means of DiffExp Hidding:2020ytt , a Mathematica implementation of the
method of generalised power series expansions Francesco:2019yqt . We focus our
analysis on the physical phase-space region relevant for $pp\to t\bar{t}$+jet,
the $s_{45}$ channel. We obtain numerical boundary values in this region with
AMFlow Liu:2022chg , which implements the auxiliary mass flow method
Liu:2017jxz ; Liu:2021wks ; Liu:2022tji . Interestingly, we observe that the
master integrals related to the problematic features (that is, the nested
square root and the elliptic curve) are non-zero only starting from order
$\eps^{4}$. We perform a number of checks that the numerical evaluation in the
$s_{45}$ channel is reliable and computationally feasible. We find that the
quadratic $\epsilon$-dependence of the DEs for topology ${\rm PB}_{B}$ results
in an increase of the evaluation time of the method with respect to the
$\epsilon$-factorised DEs for topologies ${\rm PB}_{A}$ and ${\rm PB}_{C}$.
Nevertheless, given the level of optimisation of our solution, we find the
performance analysis promising for future phenomenology applications.

There are still a number of issues to be addressed before the method can be
applied in the context of amplitude computations. Proceeding without a
representation of the integrals in terms of a basis of special functions order
by order in $\eps$ means that the poles could not be removed analytically. The
inability to perform an expansion in four dimensions may also mean that
certain simplifications in the amplitude are not observed. Nevertheless, a
strategy of performing the integration-by-parts reduction to master integrals
numerically over a rationalised phase space with modular arithmetic appears to
be achievable without major new technological developments. While combining
this with an optimised strategy for the evaluation of the master integrals
over the full phase space is left for future work, we remark that the
application of the generalised series expansion method for phenomenological
studies has been successful for other processes Becchetti:2020wof ;
Armadillo:2022bgm ; Bonciani:2021zzf ; Becchetti:2023yat . It would also be
interesting to continue the search for $\eps$-factorised DEs for ${\rm
PB}_{B}$ following the latest line of research in this area, which could
potentially lead to even more efficient numerical evaluations. This would
however require some substantial new developments in the available theoretical
tools.

Our work paves the way to the analytic computation of the two-loop amplitudes
for $pp\to t\bar{t}$+jet in the leading colour approximation, the main
bottleneck towards obtaining predictions for this important process at NNLO in
QCD.

###### Acknowledgements.

We are grateful to Heribertus Bayu Hartanto, Colomba Brancaccio, Ekta Chaubey
and Christoph Dlapa for many enlightening discussions. We also thank Colomba
Brancaccio and Xuhang Jiang for helpful comments on this paper. This project
has received funding from the European Union’s Horizon Europe research and
innovation programme under the Marie Skłodowska-Curie grant agreement No.
101105486, and ERC Starting Grant No. 101040760 _FFHiggsTop_. This work has
received funding from the Italian Ministry of Universities and Research
through FARE grant R207777C4R. This research was supported in part by the
Swiss National Science Foundation (SNF) under contract 200021_212729. SB has
been partially supported by the Italian Ministry of Universities and Research
(MUR) through grant PRIN 2022BCXSW9.

##  Appendix A Maximal-cut analysis of the elliptic sector

In this section, we analyse the maximal cut of the following scalar integral
in the sector $321B$ of topology ${\rm PB}_{B}$ (see sec. 3.5 and fig. 2(a)):

$\displaystyle J=I^{({\rm PB}_{B}),0,0,0}_{1,1,0,1,1,1,0,1}\,.$ (68)

This is related to the basis integral $\mathcal{I}_{{\rm PB}_{B},35}$ through
eq. (54). The maximal cut is solution to the homogeneous DEs Primo:2016ebd ,
and therefore contains precious information about the analytic structure of
the full solution. We adopt the Baikov parametrisation Baikov:1996iu ;
Frellesvig:2017aai , and show that the maximal cut features the square root of
a degree-4 polynomial. The latter defines an elliptic curve, whose periods
(suitably normalised) are solutions to the Picard-Fuchs operator of
$\mathcal{I}_{{\rm PB}_{B},35}$ discussed in sec. 3.5.

Figure 6: Four-point double-box super-sector of the elliptic sector in
topology ${\rm PB}_{B}$.

Instead of starting from the Baikov parametrisation of the pentagon-box top
sector, we view this sector as a sub-sector of the four-point double-box shown
in fig. 6. This allows us to get rid of two Baikov integration variables,
corresponding to irreducible scalar products of the pentagon-box. Furthermore,
we have one fewer integration variable by adopting a loop-by-loop Baikov
parametrisation Frellesvig:2017aai , starting from the $k_{2}$
loop.888Starting from the $k_{1}$ loop would lead to an equivalent result.
Neglecting the overall kinematic-independent factors and the integration
domain, our parametrisation is given by

$\displaystyle J\propto\left[{\rm
det}\,G(p_{1},p_{2},p_{5})\right]^{\eps}\int$
$\displaystyle\frac{{\mathrm{d}}z_{1}\,{\mathrm{d}}z_{2}\,{\mathrm{d}}z_{4}\,{\mathrm{d}}z_{5}\,{\mathrm{d}}z_{6}\,{\mathrm{d}}z_{8}}{z_{1}\,z_{2}\,z_{4}\,z_{5}\,z_{6}\,z_{8}}\,{\mathrm{d}}z_{7}\,{\mathrm{d}}z_{9}\times$
$\displaystyle\left[{\rm
det}\,G(k_{1},p_{1},p_{2},p_{5})\right]^{-\frac{1}{2}-\eps}\left[{\rm
det}\,G(k_{1},p_{1},p_{5})\right]^{\eps}\left[{\rm
det}\,G(k_{2},k_{1},p_{1},p_{5})\right]^{-\frac{1}{2}-\eps}\,,$ (69)

where we recall that $G$ is the Gram matrix, i.e.
$G_{ij}(v_{1},\ldots,v_{n})=v_{i}\cdot v_{j}$ for $i,j=1,\ldots,n$, and the
scalar products $k_{i}\cdot k_{j}$ and $k_{i}\cdot p_{j}$ are understood as
written in terms of inverse propagators $z_{i}$. The maximal cut, namely the
residue of the integrand at $z_{i}=0$ for $i=1,2,4,5,6,8$ in $d=4$ dimensions
is given by

$\displaystyle{\rm
MaxCut}\left[J\right]\bigl{|}_{\eps=0}\propto\int\frac{{\mathrm{d}}z_{7}\,{\mathrm{d}}z_{9}}{\sqrt{f_{1}(z_{9})\,f_{2}(z_{9})\,f_{3}(z_{7},z_{9})}}\,,$
(70)

where the $f_{i}$’s are irreducible polynomials. In particular,
$f_{1}(z_{9})=z_{9}+m_{t}^{2}$, $f_{2}(z_{9})$ is a degree-2 polynomial in
$z_{9}$, and $f_{3}(z_{7},z_{9})$ is a polynomial of degree 1 in $z_{9}$ and 2
in $z_{7}$. The dependence on the kinematic invariants is omitted.

We now wish to rewrite the integrand on the maximal cut as a linear
combination of ${\mathrm{d}}\log$ forms. The coefficients of the
${\mathrm{d}}\log$ forms are called leading singularities. If this is
possible, the integral — normalised so that its leading singularities are
constant — is expected to satisfy canonical DEs Arkani-Hamed:2010pyv . For
example, in the massless case (i.e., we set $m_{t}^{2}=0$ under the integral
sign), we can write999We omit the exterior product between the differential
forms to simplify the discussion.

$\displaystyle{\rm MaxCut}\left[J_{\rm
massless}\right]\bigl{|}_{\eps=0}\propto\frac{1}{d_{12}d_{15}}\int{\mathrm{d}}\log\left(z_{7}-2d_{15}\right)\,{\mathrm{d}}\log\left(1-\frac{2d_{12}d_{15}}{z_{9}(d_{15}-d_{34})}\right)\,,$
(71)

which indicates that $d_{12}d_{15}\,J_{\rm massless}$ — modulo sub-sector
corrections — is a good candidate for a canonical basis. With a non-zero top-
quark mass, however, only one differential form can be expressed as a
${\mathrm{d}}\log$,

$\displaystyle{\rm
MaxCut}\left[J\right]\bigl{|}_{\eps=0}\propto\int{\mathrm{d}}\log\left[\alpha(z_{7},z_{9})\right]\,\frac{{\mathrm{d}}z_{9}}{\sqrt{\mathcal{P}(z_{9})}}\,,$
(72)

where $\alpha(z_{7},z_{9})$ is an irrelevant algebraic function, and
$\mathcal{P}(z_{9})$ is a degree-4 polynomial,

$\displaystyle\mathcal{P}(z_{9})=\,$
$\displaystyle(z_{9}+m_{t}^{2})(z_{9}-3m_{t}^{2})\bigl{(}\mathcal{P}_{0}+\mathcal{P}_{1}\,z_{9}+\mathcal{P}_{2}\,z_{9}^{2}\bigr{)}\,,$
(73)

with

$\displaystyle\begin{aligned}
\mathcal{P}_{0}&=4d_{12}^{2}d_{15}^{2}+4d_{12}d_{15}\left(d_{15}+d_{34}\right)m_{t}^{2}+\left(d_{15}^{2}+4d_{12}d_{34}+2d_{15}d_{34}-3d_{34}^{2}\right)m_{t}^{4}+2d_{34}m_{t}^{6}\,,\\\
\mathcal{P}_{1}&=4d_{12}d_{15}\left(d_{34}-d_{15}\right)-2\left(d_{15}^{2}-2d_{12}d_{34}+d_{34}^{2}\right)m_{t}^{2}\,,\\\
\mathcal{P}_{2}&=\left(d_{15}-d_{34}\right)^{2}-2d_{34}m_{t}^{2}\,.\end{aligned}$
(74)

For generic values of the kinematic invariants, this degree-4 polynomial has
four distinct roots, and thus defines an elliptic curve:

$\displaystyle y^{2}=\mathcal{P}(x)\,.$ (75)

The differential form which cannot be expressed as a ${\mathrm{d}}\log$ in eq.
(72) is the holomorphic differential form of the first kind on this elliptic
curve, which is one of the integration kernels defining the elliptic multiple
polylogarithms Broedel:2018qkq  (see also ref. Gorges:2023zgv  for a similar
example). In the massless limit ($m_{t}^{2}=0$), instead, two of the roots
degenerate, so that $\mathcal{P}(z_{9})$ becomes a perfect square and the
elliptic curve reduces to a genus-zero surface. This is consistent with the
fact that the massless pentagon-box integrals are of polylogarithmic type
Gehrmann:2015bfy ; Papadopoulos:2015jft .

The elliptic curve in eq. (75) contains useful information for putting the DEs
in $\eps$-factorised form, and is tightly connected to the Picard-Fuchs
operator of the integral $\mathcal{I}_{{\rm PB}_{B},35}$ discussed in sec.
3.5. In particular, we expect that the Feynman integral $J$ normalised by a
period of the elliptic curve in eq. (75) is a good candidate for a canonical
integral Gorges:2023zgv . While we leave this study for further work, we
discuss here the connection with the Picard-Fuchs operator.

For this purpose, we must first spell out the periods of the elliptic curve.
The four roots are given by101010In the $s_{45}$ channel we have that
$e_{1}<0$, $e_{4}>0$, $\delta<0$, while the real parts of $e_{2}$ and $e_{3}$
do not have fixed sign.

$\displaystyle\begin{aligned} &e_{1}=-m_{t}^{2}\,,\\\
&e_{2}=\frac{2d_{12}d_{15}(d_{34}-d_{15})-(d_{15}^{2}+d_{34}^{2}-2d_{12}d_{34})m_{t}^{2}+d_{34}\sqrt{\delta}}{{\rm
det}\,G(p_{2},p_{1}+p_{5})}\,,\\\
&e_{3}=e_{2}\bigl{|}_{\sqrt{\delta}\to-\sqrt{\delta}}\,,\\\
&e_{4}=3m_{t}^{2}\,,\end{aligned}$ (76)

where

$\displaystyle\delta=-4\,m_{t}^{2}\,{\rm det}\,G(p_{1},p_{2},p_{5})\,.$ (77)

The periods can then be chosen as Weinzierl:2022eaz

$\displaystyle\psi_{1}=\frac{4\,{\rm
K}(\kappa^{2})}{\sqrt{(e_{3}-e_{1})(e_{4}-e_{2})}}\,,\qquad\quad\psi_{2}=\frac{4\,{\mathrm{i}}\,{\rm
K}(1-\kappa^{2})}{\sqrt{(e_{3}-e_{1})(e_{4}-e_{2})}}\,,$ (78)

where ${\rm K}(x)$ is the elliptic integral of the first kind, and $\kappa$ is
the modulus of the elliptic curve,

$\displaystyle\kappa^{2}=\frac{(e_{3}-e_{2})(e_{4}-e_{1})}{(e_{3}-e_{1})(e_{4}-e_{2})}\,.$
(79)

Note that $e_{2}$ and $e_{3}$ involve a square root, and thus the periods
feature another nested square root, in addition to the one discussed in sec.
3.3. We find that — with a suitable algebraic normalisation — the periods
$\psi_{i}$ are solutions to the Picard-Fuchs operator $L_{35}$ of the Feynman
integral $\mathcal{I}_{{\rm PB}_{B},35}$ (modulo sub-sectors and
$\eps$-corrections). Explicitly, we have that

$\displaystyle
L_{35}\,\left[\frac{d_{15}(d_{12}+m_{t}^{2})\,\psi_{i}}{\sqrt{(d_{15}-d_{34})^{2}-2d_{34}m_{t}^{2}}}\right]=0\,,$
(80)

for $i=1,2$, where the terms in the square brackets are understood as
evaluated on the univariate phase-space slice used to derive $L_{35}$. The
rational factors multiplying the periods follow from the chosen normalisation
of $\mathcal{I}_{{\rm PB}_{B},35}$ (see eq. (54)), while the square root is
the coefficient of the degree-$4$ monomial in the polynomial $\mathcal{P}$
defining the elliptic curve. With this, we have established that the elliptic
integrals found in the solutions to the Picard-Fuchs operator $L_{35}$ in sec.
3.5 and thus in the homogeneous solution to the DEs are associated with the
elliptic curve given in eq. (75).

##  Appendix B All one-loop pentagon integrals

We update the representation of the differential equations for all one-loop
pentagon integrals contributing to $pp\to t{\bar{t}}+$jet previously
considered in reference Badger:2022mrb . We present them in
‘${\mathrm{d}}\log$’ form using an alphabet which shares as much in common
with the two-loop planar integrals as possible. For completeness we include
here the definitions of the integral families which are shown graphically in
figure 7.

(a) Topology ${\rm P}_{A}$.

(b) Topology ${\rm P}_{B}$.

(c) Topology ${\rm P}_{C}$.

(d) Topology ${\rm P}_{D}$.

Figure 7: Graphical representation of the four one-loop pentagon topologies.

The integrals of the pentagon topology $F\in\\{{\rm P}_{A},{\rm P}_{B},{\rm
P}_{C},{\rm P}_{D}\\}$ have the form

$I_{a_{1},a_{2},a_{3},a_{4},a_{5}}^{\left(F\right)}=\int\mathcal{D}^{d}k_{1}\,\frac{1}{D_{F,1}^{a_{1}}\cdots
D_{F,5}^{a_{5}}}\,,$ (81)

where we recall that the integration measure $\mathcal{D}^{d}k_{1}$ is defined
in eq. (2). The inverse propagators $D_{F,i}$ are defined in table 3.

| ${\rm P}_{A}$ | ${\rm P}_{B}$ | ${\rm P}_{C}$ | ${\rm P}_{D}$
---|---|---|---|---
$D_{F,1}$ | $k_{1}^{2}$ | $k_{1}^{2}$ | $k_{1}^{2}-m_{t}^{2}$ | $k_{1}^{2}-m_{t}^{2}$
$D_{F,2}$ | $(k_{1}-p_{1})^{2}-m_{t}^{2}$ | $(k_{1}-p_{1})^{2}-m_{t}^{2}$ | $(k_{1}-p_{1})^{2}$ | $(k_{1}-p_{1})^{2}$
$D_{F,3}$ | $(k_{1}-p_{1}-p_{2})^{2}$ | $(k_{1}-p_{1}-p_{3})^{2}-m_{t}^{2}$ | $(k_{1}-p_{1}-p_{3})^{2}$ | $(k_{1}-p_{1}-p_{2})^{2}-m_{t}^{2}$
$D_{F,4}$ | $(k_{1}+p_{4}+p_{5})^{2}$ | $(k_{1}+p_{4}+p_{5})^{2}$ | $(k_{1}+p_{4}+p_{5})^{2}-m_{t}^{2}$ | $(k_{1}+p_{4}+p_{5})^{2}-m_{t}^{2}$
$D_{F,5}$ | $(k_{1}+p_{5})^{2}$ | $(k_{1}+p_{5})^{2}$ | $(k_{1}+p_{5})^{2}-m_{t}^{2}$ | $(k_{1}+p_{5})^{2}-m_{t}^{2}$
Table 3: Inverse propagators $D_{F,i}$ of the pentagon topologies shown in
figure 7.

##  Appendix C Description of the ancillary files

The ancillary files can be downloaded from ref. ancillary . The symbols in the
files are given in the notation of the article as follows:

eps $\displaystyle=\eps\,,$ W[i] $\displaystyle=W_{i}\,,$ dij
$\displaystyle=d_{ij}\,,$ j[F,a1,…,a5]
$\displaystyle=I^{(F)}_{a_{1},\ldots,a_{5}}\,,$ mt2
$\displaystyle=m_{t}^{2}\,,$ j[F,a1,…,a11]
$\displaystyle=I^{(F),-a_{9},-a_{10},-a_{11}}_{a_{1},\ldots,a_{8}}\,,$ s[i,j]
$\displaystyle=s_{ij}\,,$ dlog[W[i]] $\displaystyle={\rm
d}\log\left(W_{i}\right)\,,$ s[i] $\displaystyle=p_{i}^{2}\,,$ of[i]
$\displaystyle=\omega_{i}\,,$ pi $\displaystyle=p_{i}\,,$ GramDet[P1,…,Pn]
$\displaystyle={\rm det}\,G(P_{1},\ldots,P_{n})\,,$ ki
$\displaystyle=k_{i}\,,$ trp[i1,…,in] $\displaystyle={\rm
tr}_{+}\left[\not{p}_{i_{1}}\ldots\not{p}_{i_{n}}\right]\,,$ p[i,j]
$\displaystyle=p_{ij}\,,$ trm[i1,…,in] $\displaystyle={\rm
tr}_{-}\left[\not{p}_{i_{1}}\ldots\not{p}_{i_{n}}\right]\,,$ spAB[i,…,j]
$\displaystyle=\langle i|\ldots|j]\,,$ sqrtratio[a,b]
$\displaystyle=\frac{a+b}{a-b}\,.$

The following files give global definitions of the square roots, the alphabet
letters and the physical scattering region under consideration.

  * •

square_roots.m – definition of the square roots in the format

$\texttt{\\{ \ldots, S -> Sqrt[expr] , \ldots\\}}\,,$

where S is the square-root label,

tr5 $\displaystyle={\rm tr}_{5},\quad\texttt{beta12}$
$\displaystyle=\beta_{12},\quad\texttt{beta34}$
$\displaystyle=\beta_{34},\quad\texttt{beta45}$
$\displaystyle=\beta_{45},\quad\texttt{Delta31}$ $\displaystyle=\Delta_{3,1},$
Delta32 $\displaystyle=\Delta_{3,2},\quad\texttt{Delta33}$
$\displaystyle=\Delta_{3,3},\quad\texttt{Delta34}$
$\displaystyle=\Delta_{3,4},\quad\texttt{Lambda1}$
$\displaystyle=\Lambda_{1},\quad\texttt{Lambda2}$ $\displaystyle=\Lambda_{2},$
Lambda3 $\displaystyle=\Lambda_{3},\quad\texttt{Lambda4}$
$\displaystyle=\Lambda_{4},\quad\texttt{Lambda5}$
$\displaystyle=\Lambda_{5},\quad\texttt{Lambda6}$ $\displaystyle=\Lambda_{6},$

and expr is a rational function of the invariants $\vec{x}$. See sec. 2 for a
representation in terms of Gram and Cayley determinants.

  * •

alphabet.m – definition of the alphabet letters $W_{i}$ in the format

$\texttt{\\{\ldots, W[i] -> expr, \ldots\\}}\,,$

where expr is given in terms of invariants $\vec{x}$, $s_{ij}$ variables, Gram
determinants, spinor-helicity chains, and gamma-matrix traces.

  * •

alphabet_dij.m – definition of the alphabet letters W[i] in the same format as
in alphabet.m, but with expr given in terms of invariants $\vec{x}$ and square
roots.

  * •

alphabet_charges.m – charges of the letters in the format

$\texttt{\\{\ldots, C -> \\{W[i1],W[i2],\ldots\\},\ldots\\}}\,,$

meaning that the letters $\\{W_{i_{1}},W_{i_{2}},\ldots\\}$ are odd with
respect to the square-root monomial C. The letters with C$\ =1$ are rational,
and thus even with respect to all square roots. See the introduction of sec. 3
for the notion of charge.

  * •

s45-channel.m – the inequalities in the $d_{ij}$ variables which define the
$s_{45}$ channel (see sec. 5.1).

For each integral topology <fam> studied in this work there is a folder with
the same name containing the definition of the propagators and of the integral
basis, the differential equations, and the boundary values. The integral
topologies are denoted by

$\displaystyle\texttt{PttjA}={\rm P}_{A}\,,\qquad\texttt{PttjB}={\rm
P}_{B}\,,\qquad\texttt{PttjC}={\rm P}_{C}\,,\qquad\texttt{PttjD}={\rm
P}_{D}\,,$ $\displaystyle\texttt{PBttjA}={\rm
PB}_{A}\,,\qquad\texttt{PBttjB}={\rm PB}_{B}\,,\qquad\texttt{PBttjC}={\rm
PB}_{C}\,.$

The folder of the family <fam> contains the following files.

  * •

<fam>_propagators.m – the inverse propagators $D_{F,i}$.

  * •

<fam>_basis_definitions.m – the definition of the (square-root free) integral
basis in the format

$\displaystyle\texttt{\\{\ldots, mi[<fam>,i] -> expr, \ldots\\}}\,,$

where mi[<fam>,i] denotes the $i$-th master integral of the topology <fam>,
and expr is a linear combination of scalar integrals with coefficients given
by rational functions of the invariants $\vec{x}$ and $\eps$.

  * •

<fam>_basis_norm.m – the square-root normalisations that should be applied in
addition to the basis integrals. More explicitly, the $i$-th master integral
in <fam>_basis_definitions.m has to be multiplied by the $i$-th entry of
<fam>_basis_norm.m. The separation between rational expressions and square
roots is motivated in sec. 3.

  * •

<fam>_de_dlogs.m (or PBttjB_de_one-forms.m where the
${\mathrm{d}}\log$-representation is not possible) – the connection matrices
of the differential equations for the master integrals expressed in terms of
${\mathrm{d}}\log$’s and, for ${\rm PB}_{B}$, one-forms $\omega_{i}$.

  * •

<fam>_boundary_values_s45.m – phase-space point and numerical values of the
master integrals in the $s_{45}$ scattering region with at least 32-digit
precision. The format is {X,values}. The first entry, X, gives the boundary
point $\vec{x}_{0}$ defined in eq. (65), in the format

$\displaystyle\texttt{<|d12->2, d23->1, d34->-1, d45->5, d15 ->-2,
mt2->1|>}\,.$

The second entry, values, is a $n_{\rm MI}\times 5$ array, where $n_{\rm MI}$
is the number of master integrals of the family, and $5$ is the number of
orders in $\eps$ (from $\eps^{0}$ to $\eps^{4}$). The entry $(i,j)$ of values
gives the value of the $\mathcal{O}\left(\eps^{j-1}\right)$ coefficient of the
$i$-th master integral.

The folder for topology ${\rm PB}_{B}$ contains two additional files, due to
the fact that the differential equations cannot be expressed in terms of
${\mathrm{d}}\log$’s only.

  * •

PBttjB_one-forms_definitions.m – definition of the one-forms $\omega_{i}$ in
the format

$\displaystyle\texttt{\\{\ldots, of[i] -> \\{r1,r2,\ldots,r6\\},
\ldots\\}}\,,$

meaning that

$\displaystyle\omega_{i}(\vec{x})=\sum_{j=1}^{6}r_{j}(\vec{x})\,{\mathrm{d}}x_{j}\,,$

where $r_{j}=\ $rj are functions of the invariants $\vec{x}$ and of the square
roots.

  * •

PBttjB_one-forms_charges.m – charges of the one-forms in the format

$\displaystyle\texttt{\\{\ldots,C -> \\{of[i1],of[i2],\ldots\\},\ldots\\}}\,,$

meaning that the one-forms $\\{\omega_{i_{1}}$, $\omega_{i_{2}},\ldots\\}$ are
odd with respect to the square-root monomial C. The one-forms with C$\ =1$ are
rational, and thus even with respect to all square roots.

The folder benchmarks/ contains benchmark values of all basis integrals.

  * •

point_x1_s45.m – phase-space point $\vec{x}_{1}$ in the $s_{45}$ channel
defined in eq. (67), in the format

$\displaystyle\texttt{\\{d12->\\#val, d23->\\#val, d34->\\#val, d45->\\#val,
d15 ->\\#val, mt2->\\#val\\}}\,,$ (82)

where #val denotes a numerical value.

  * •

<fam>_values_x1.m – values of the basis integrals of family <fam> at the
phase-space point $\vec{x}_{1}$ with at least 32-digit precision. The format
is the same as in <fam>_boundary_values_s45.m.

A Mathematica script, DiffExp_run.wl, is provided to evaluate numerically the
basis integrals through DiffExp Hidding:2020ytt . The path to DiffExp.m has to
be specified through the variable PathToDiffExp. The script may be used via
the command line as

math -script DiffExp_run.wl [-family <fam>] [-target <file>] [-storepiecewise]

The file indicated to be the target phase-space point is expected in the
format shown in eq. (82). If the option -storepiecewise is given, the script
also stores the analytic expression of the generalised power series solution
to the differential equations. Otherwise, only the numerical values of the
basis integrals at the target point are saved, in the benchmarks/ folder. The
integration of the differential equations occurs along a straight path
connecting the target point to the boundary point $\vec{x}_{0}$ in eq. (65).
If the path leaves the $s_{45}$ channel, the evaluation is aborted. This may
occur even if the target point is in the $s_{45}$ channel. In such a case, one
should either choose a different starting point for the integration path, or
study the analytic continuation from the $s_{45}$ channel to the region of
interest. The latter is instead mandatory if the target point is outside of
the $s_{45}$ channel. The accuracy goal is hard-coded to 16 digits. If the
target point is $\vec{x}_{1}$, the evaluations are compared against the
provided benchmark values.

## References

  * (1) CMS collaboration, V. Khachatryan et al., _Measurement of differential cross sections for top quark pair production using the lepton+jets final state in proton-proton collisions at 13 TeV_ , _Phys. Rev. D_ 95 (2017) 092001, [1610.04191].
  * (2) ATLAS collaboration, M. Aaboud et al., _Measurements of differential cross sections of top quark pair production in association with jets in ${pp}$ collisions at $\sqrt{s}=13$ TeV using the ATLAS detector_, _JHEP_ 10 (2018) 159, [1802.06572].
  * (3) CMS collaboration, A. M. Sirunyan et al., _Measurement of the cross section for $\text{t}\bar{\text{t}}$ production with additional jets and b jets in pp collisions at $\sqrt{s}=$ 13 TeV_, _JHEP_ 07 (2020) 125, [2003.06467].
  * (4) CMS collaboration, _Differential cross section measurements for the production of top quark pairs and of additional jets using dilepton events from pp collisions at $\sqrt{s}$ = 13 TeV_, 2402.08486.
  * (5) S. Alioli, P. Fernandez, J. Fuster, A. Irles, S.-O. Moch, P. Uwer et al., _A new observable to measure the top-quark mass at hadron colliders_ , _Eur. Phys. J. C_ 73 (2013) 2438, [1303.6415].
  * (6) G. Bevilacqua, H. B. Hartanto, M. Kraus, M. Schulze and M. Worek, _Top quark mass studies with $t\overline{t}j$ at the LHC_, _JHEP_ 03 (2018) 169, [1710.07515].
  * (7) S. Alioli, J. Fuster, M. V. Garzelli, A. Gavardi, A. Irles, D. Melini et al., _Phenomenology of $t\overline{t}j$ \+ X production at the LHC_, _JHEP_ 05 (2022) 146, [2202.07975].
  * (8) ATLAS collaboration, G. Aad et al., _Measurement of the top-quark mass in $t\bar{t}+1$-jet events collected with the ATLAS detector in $pp$ collisions at $\sqrt{s}=8$ TeV_, _JHEP_ 11 (2019) 150, [1905.02302].
  * (9) J. L. Bourjaily et al., _Functions Beyond Multiple Polylogarithms for Precision Collider Physics_ , in _Snowmass 2021_ , 3, 2022.  2203.07088.
  * (10) S. Dittmaier, P. Uwer and S. Weinzierl, _NLO QCD corrections to t anti-t \+ jet production at hadron colliders_ , _Phys. Rev. Lett._ 98 (2007) 262002, [hep-ph/0703120].
  * (11) S. Dittmaier, P. Uwer and S. Weinzierl, _Hadronic top-quark pair production in association with a hard jet at next-to-leading order QCD: Phenomenological studies for the Tevatron and the LHC_ , _Eur. Phys. J. C_ 59 (2009) 625–646, [0810.0452].
  * (12) K. Melnikov and M. Schulze, _NLO QCD corrections to top quark pair production in association with one hard jet at hadron colliders_ , _Nucl. Phys. B_ 840 (2010) 129–159, [1004.3284].
  * (13) S. Alioli, S.-O. Moch and P. Uwer, _Hadronic top-quark pair-production with one jet and parton showering_ , _JHEP_ 01 (2012) 137, [1110.5251].
  * (14) M. Czakon, H. B. Hartanto, M. Kraus and M. Worek, _Matching the Nagy-Soper parton shower at next-to-leading order_ , _JHEP_ 06 (2015) 033, [1502.00925].
  * (15) G. Bevilacqua, H. B. Hartanto, M. Kraus and M. Worek, _Top Quark Pair Production in Association with a Jet with Next-to-Leading-Order QCD Off-Shell Effects at the Large Hadron Collider_ , _Phys. Rev. Lett._ 116 (2016) 052003, [1509.09242].
  * (16) G. Bevilacqua, H. B. Hartanto, M. Kraus and M. Worek, _Off-shell Top Quarks with One Jet at the LHC: A comprehensive analysis at NLO QCD_ , _JHEP_ 11 (2016) 098, [1609.01659].
  * (17) C. Gütschow, J. M. Lindert and M. Schönherr, _Multi-jet merged top-pair production including electroweak corrections_ , _Eur. Phys. J. C_ 78 (2018) 317, [1803.00950].
  * (18) J. M. Campbell et al., _Event Generators for High-Energy Physics Experiments_ , in _Snowmass 2021_ , 3, 2022.  2203.11110.
  * (19) S. Badger, M. Becchetti, E. Chaubey, R. Marzucca and F. Sarandrea, _One-loop QCD helicity amplitudes for pp → $t\overline{t}j$ to O($\varepsilon$2)_, _JHEP_ 06 (2022) 066, [2201.12188].
  * (20) S. Badger, M. Becchetti, E. Chaubey and R. Marzucca, _Two-loop master integrals for a planar topology contributing to pp → $t\overline{t}j$_, _JHEP_ 01 (2023) 156, [2210.17477].
  * (21) A. V. Kotikov, _Differential equations method: New technique for massive Feynman diagrams calculation_ , _Phys. Lett. B_ 254 (1991) 158–164.
  * (22) E. Remiddi, _Differential equations for Feynman graph amplitudes_ , _Nuovo Cim. A_ 110 (1997) 1435–1452, [hep-th/9711188].
  * (23) J. M. Henn, _Multiloop integrals in dimensional regularization made simple_ , _Phys. Rev. Lett._ 110 (2013) 251601, [1304.1806].
  * (24) A. von Manteuffel and R. M. Schabinger, _A novel approach to integration by parts reduction_ , _Phys. Lett. B_ 744 (2015) 101–104, [1406.4513].
  * (25) T. Peraro, _Scattering amplitudes over finite fields and multivariate functional reconstruction_ , _JHEP_ 12 (2016) 030, [1608.01902].
  * (26) T. Peraro, _FiniteFlow: multivariate functional reconstruction using finite fields and dataflow graphs_ , _JHEP_ 07 (2019) 031, [1905.08019].
  * (27) T. Gehrmann, J. Henn and N. Lo Presti, _Analytic form of the two-loop planar five-gluon all-plus-helicity amplitude in QCD_ , _Phys. Rev. Lett._ 116 (2016) 062001, [1511.05409].
  * (28) T. Gehrmann, J. M. Henn and N. A. Lo Presti, _Pentagon functions for massless planar scattering amplitudes_ , _JHEP_ 10 (2018) 103, [1807.09812].
  * (29) D. Chicherin and V. Sotnikov, _Pentagon Functions for Scattering of Five Massless Particles_ , _JHEP_ 20 (2020) 167, [2009.07803].
  * (30) D. Chicherin, V. Sotnikov and S. Zoia, _Pentagon functions for one-mass planar scattering amplitudes_ , _JHEP_ 01 (2022) 096, [2110.10111].
  * (31) S. Abreu, D. Chicherin, H. Ita, B. Page, V. Sotnikov, W. Tschernow et al., _All Two-Loop Feynman Integrals for Five-Point One-Mass Scattering_ , 2306.15431.
  * (32) C. G. Papadopoulos, D. Tommasini and C. Wever, _The Pentabox Master Integrals with the Simplified Differential Equations approach_ , _JHEP_ 04 (2016) 078, [1511.09404].
  * (33) S. Abreu, L. J. Dixon, E. Herrmann, B. Page and M. Zeng, _The two-loop five-point amplitude in $\mathcal{N}=4$ super-Yang-Mills theory_, _Phys. Rev. Lett._ 122 (2019) 121603, [1812.08941].
  * (34) D. Chicherin, T. Gehrmann, J. M. Henn, P. Wasser, Y. Zhang and S. Zoia, _All Master Integrals for Three-Jet Production at Next-to-Next-to-Leading Order_ , _Phys. Rev. Lett._ 123 (2019) 041603, [1812.11160].
  * (35) S. Abreu, H. Ita, F. Moriello, B. Page, W. Tschernow and M. Zeng, _Two-Loop Integrals for Planar Five-Point One-Mass Processes_ , _JHEP_ 11 (2020) 117, [2005.04195].
  * (36) D. D. Canko, C. G. Papadopoulos and N. Syrrakos, _Analytic representation of all planar two-loop five-point Master Integrals with one off-shell leg_ , _JHEP_ 01 (2021) 199, [2009.13917].
  * (37) S. Abreu, H. Ita, B. Page and W. Tschernow, _Two-loop hexa-box integrals for non-planar five-point one-mass processes_ , _JHEP_ 03 (2022) 182, [2107.14180].
  * (38) S. Badger, C. Brønnum-Hansen, H. B. Hartanto and T. Peraro, _Analytic helicity amplitudes for two-loop five-gluon scattering: the single-minus case_ , _JHEP_ 01 (2019) 186, [1811.11699].
  * (39) S. Abreu, F. Febres Cordero, H. Ita, B. Page and V. Sotnikov, _Planar Two-Loop Five-Parton Amplitudes from Numerical Unitarity_ , _JHEP_ 11 (2018) 116, [1809.09067].
  * (40) D. Chicherin, T. Gehrmann, J. Henn, P. Wasser, Y. Zhang and S. Zoia, _Analytic result for a two-loop five-particle amplitude_ , _Phys. Rev. Lett._ 122 (2019) 121602, [1812.11057].
  * (41) S. Abreu, J. Dormans, F. Febres Cordero, H. Ita and B. Page, _Analytic Form of Planar Two-Loop Five-Gluon Scattering Amplitudes in QCD_ , _Phys. Rev. Lett._ 122 (2019) 082002, [1812.04586].
  * (42) S. Badger, D. Chicherin, T. Gehrmann, G. Heinrich, J. Henn, T. Peraro et al., _Analytic form of the full two-loop five-gluon all-plus helicity amplitude_ , _Phys. Rev. Lett._ 123 (2019) 071601, [1905.03733].
  * (43) S. Abreu, J. Dormans, F. Febres Cordero, H. Ita, B. Page and V. Sotnikov, _Analytic Form of the Planar Two-Loop Five-Parton Scattering Amplitudes in QCD_ , _JHEP_ 05 (2019) 084, [1904.00945].
  * (44) H. B. Hartanto, S. Badger, C. Brønnum-Hansen and T. Peraro, _A numerical evaluation of planar two-loop helicity amplitudes for a W-boson plus four partons_ , _JHEP_ 09 (2019) 119, [1906.11862].
  * (45) D. Chicherin, T. Gehrmann, J. M. Henn, P. Wasser, Y. Zhang and S. Zoia, _The two-loop five-particle amplitude in $\mathcal{N}$ = 8 supergravity_, _JHEP_ 03 (2019) 115, [1901.05932].
  * (46) S. Abreu, B. Page, E. Pascual and V. Sotnikov, _Leading-Color Two-Loop QCD Corrections for Three-Photon Production at Hadron Colliders_ , _JHEP_ 01 (2021) 078, [2010.15834].
  * (47) H. A. Chawdhry, M. Czakon, A. Mitov and R. Poncelet, _Two-loop leading-color helicity amplitudes for three-photon production at the LHC_ , _JHEP_ 06 (2021) 150, [2012.13553].
  * (48) S. Caron-Huot, D. Chicherin, J. Henn, Y. Zhang and S. Zoia, _Multi-Regge Limit of the Two-Loop Five-Point Amplitudes in $\mathcal{N}=4$ Super Yang-Mills and $\mathcal{N}=8$ Supergravity_, _JHEP_ 10 (2020) 188, [2003.03120].
  * (49) S. Abreu, F. F. Cordero, H. Ita, M. Klinkert, B. Page and V. Sotnikov, _Leading-Color Two-Loop Amplitudes for Four Partons and a W Boson in QCD_ , 2110.07541.
  * (50) S. Badger, H. B. Hartanto, J. Kryś and S. Zoia, _Two-loop leading-colour QCD helicity amplitudes for Higgs boson production in association with a bottom-quark pair at the LHC_ , _JHEP_ 11 (2021) 012, [2107.14733].
  * (51) B. Agarwal, F. Buccioni, A. von Manteuffel and L. Tancredi, _Two-loop leading colour QCD corrections to $q\bar{q}\to\gamma\gamma g$ and $qg\to\gamma\gamma q$_, _JHEP_ 04 (2021) 201, [2102.01820].
  * (52) S. Badger, C. Brønnum-Hansen, D. Chicherin, T. Gehrmann, H. B. Hartanto, J. Henn et al., _Virtual QCD corrections to gluon-initiated diphoton plus jet production at hadron colliders_ , _JHEP_ 11 (2021) 083, [2106.08664].
  * (53) H. A. Chawdhry, M. Czakon, A. Mitov and R. Poncelet, _Two-loop leading-colour QCD helicity amplitudes for two-photon plus jet production at the LHC_ , _JHEP_ 07 (2021) 164, [2103.04319].
  * (54) S. Badger, H. B. Hartanto and S. Zoia, _Two-Loop QCD Corrections to $Wb\bar{b}$ Production at Hadron Colliders_, _Phys. Rev. Lett._ 127 (2021) 012001, [2102.02516].
  * (55) S. Abreu, F. F. Cordero, H. Ita, B. Page and V. Sotnikov, _Leading-color two-loop QCD corrections for three-jet production at hadron colliders_ , _JHEP_ 07 (2021) 095, [2102.13609].
  * (56) B. Agarwal, F. Buccioni, A. von Manteuffel and L. Tancredi, _Two-Loop Helicity Amplitudes for Diphoton Plus Jet Production in Full Color_ , _Phys. Rev. Lett._ 127 (2021) 262001, [2105.04585].
  * (57) S. Badger, H. B. Hartanto, J. Kryś and S. Zoia, _Two-loop leading colour helicity amplitudes for W ±$\gamma$ \+ j production at the LHC_, _JHEP_ 05 (2022) 035, [2201.04075].
  * (58) H. B. Hartanto, R. Poncelet, A. Popescu and S. Zoia, _NNLO QCD corrections to $Wb\bar{b}$ production at the LHC_, 2205.01687.
  * (59) S. Abreu, G. De Laurentis, H. Ita, M. Klinkert, B. Page and V. Sotnikov, _Two-loop QCD corrections for three-photon production at hadron colliders_ , _SciPost Phys._ 15 (2023) 157, [2305.17056].
  * (60) G. De Laurentis, H. Ita and V. Sotnikov, _Double-Virtual NNLO QCD Corrections for Five-Parton Scattering: The Quark Channels_ , 2311.18752.
  * (61) S. Badger, M. Czakon, H. B. Hartanto, R. Moodie, T. Peraro, R. Poncelet et al., _Isolated photon production in association with a jet pair through next-to-next-to-leading order in QCD_ , _JHEP_ 10 (2023) 071, [2304.06682].
  * (62) G. De Laurentis, H. Ita, M. Klinkert and V. Sotnikov, _Double-Virtual NNLO QCD Corrections for Five-Parton Scattering: The Gluon Channel_ , 2311.10086.
  * (63) B. Agarwal, F. Buccioni, F. Devoto, G. Gambuti, A. von Manteuffel and L. Tancredi, _Five-Parton Scattering in QCD at Two Loops_ , 2311.09870.
  * (64) H. A. Chawdhry, M. L. Czakon, A. Mitov and R. Poncelet, _NNLO QCD corrections to three-photon production at the LHC_ , _JHEP_ 02 (2020) 057, [1911.00479].
  * (65) S. Kallweit, V. Sotnikov and M. Wiesemann, _Triphoton production at hadron colliders in NNLO QCD_ , _Phys. Lett. B_ 812 (2021) 136013, [2010.04681].
  * (66) H. A. Chawdhry, M. Czakon, A. Mitov and R. Poncelet, _NNLO QCD corrections to diphoton production with an additional jet at the LHC_ , _JHEP_ 09 (2021) 093, [2105.06940].
  * (67) S. Badger, T. Gehrmann, M. Marcoli and R. Moodie, _Next-to-leading order QCD corrections to diphoton-plus-jet production through gluon fusion at the LHC_ , _Phys. Lett. B_ 824 (2022) 136802, [2109.12003].
  * (68) X. Chen, T. Gehrmann, E. W. N. Glover, A. Huss and M. Marcoli, _Automation of antenna subtraction in colour space: gluonic processes_ , _JHEP_ 10 (2022) 099, [2203.13531].
  * (69) L. Buonocore, S. Devoto, S. Kallweit, J. Mazzitelli, L. Rottoli and C. Savoini, _Associated production of a W boson and massive bottom quarks at next-to-next-to-leading order in QCD_ , _Phys. Rev. D_ 107 (2023) 074032, [2212.04954].
  * (70) H. B. Hartanto, R. Poncelet, A. Popescu and S. Zoia, _Flavour anti- $k_{\text{T}}$ algorithm applied to $Wb\bar{b}$ production at the LHC_, 2209.03280.
  * (71) J. Mazzitelli, V. Sotnikov and M. Wiesemann, _Next-to-next-to-leading order event generation for Z-boson production in association with a bottom-quark pair_ , 2404.08598.
  * (72) J. Henn, T. Peraro, Y. Xu and Y. Zhang, _A first look at the function space for planar two-loop six-particle Feynman integrals_ , _JHEP_ 03 (2022) 056, [2112.10605].
  * (73) J. M. Henn, A. Matijašić, J. Miczajka, T. Peraro, Y. Xu and Y. Zhang, _A computation of two-loop six-point Feynman integrals in dimensional regularization_ , 2403.19742.
  * (74) F. Febres Cordero, G. Figueiredo, M. Kraus, B. Page and L. Reina, _Two-Loop Master Integrals for Leading-Color $pp\to t\bar{t}H$ Amplitudes with a Light-Quark Loop_, 2312.08131.
  * (75) F. Buccioni, P. A. Kreer, X. Liu and L. Tancredi, _One loop QCD corrections to $gg\to t\overline{t}H$ at $\mathcal{O}(\epsilon^{2})$_, 2312.10015.
  * (76) B. Agarwal, G. Heinrich, S. P. Jones, M. Kerner, S. Y. Klein, J. Lang et al., _Two-loop amplitudes for ttH production: the quark-initiated Nf-part_ , 2402.03301.
  * (77) S. Catani, S. Devoto, M. Grazzini, S. Kallweit, J. Mazzitelli and C. Savoini, _Higgs Boson Production in Association with a Top-Antitop Quark Pair in Next-to-Next-to-Leading Order QCD_ , _Phys. Rev. Lett._ 130 (2023) 111902, [2210.07846].
  * (78) G. Wang, T. Xia, L. L. Yang and X. Ye, _Two-loop QCD amplitudes for $t\bar{t}H$ production from boosted limit_, 2402.00431.
  * (79) 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].
  * (80) J. Broedel, C. Duhr, F. Dulat, B. Penante and L. Tancredi, _Elliptic Feynman integrals and pure functions_ , _JHEP_ 01 (2019) 023, [1809.10698].
  * (81) 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].
  * (82) J. Broedel, C. Duhr, F. Dulat, B. Penante and L. Tancredi, _Elliptic polylogarithms and Feynman parameter integrals_ , _JHEP_ 05 (2019) 120, [1902.09971].
  * (83) H. Frellesvig and S. Weinzierl, _On $\varepsilon$-factorised bases and pure Feynman integrals_, 2301.02264.
  * (84) 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].
  * (85) F. Moriello, _Generalised power series expansions for the elliptic planar families of Higgs + jet production at two loops_ , _JHEP_ 01 (2020) 150, [1907.13234].
  * (86) F. V. Tkachov, _A Theorem on Analytical Calculability of Four Loop Renormalization Group Functions_ , _Phys. Lett._ 100B (1981) 65–68.
  * (87) 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–204.
  * (88) R. N. Lee, _Presenting LiteRed: a tool for the Loop InTEgrals REDuction_ , 1212.2685.
  * (89) R. N. Lee, _LiteRed 1.4: a powerful tool for reduction of multiloop integrals_ , _J. Phys. Conf. Ser._ 523 (2014) 012059, [1310.1145].
  * (90) Z. Wu, J. Boehm, R. Ma, H. Xu and Y. Zhang, _NeatIBP 1.0, A package generating small-size integration-by-parts relations for Feynman integrals_ , 2305.08783.
  * (91) J. Gluza, K. Kajda and D. A. Kosower, _Towards a Basis for Planar Two-Loop Integrals_ , _Phys. Rev. D_ 83 (2011) 045012, [1009.0472].
  * (92) S. Laporta, _High precision calculation of multiloop Feynman integrals by difference equations_ , _Int. J. Mod. Phys._ A15 (2000) 5087–5159, [hep-ph/0102033].
  * (93) T. Peraro, _Analytic multi-loop results using finite fields and dataflow graphs with FiniteFlow_ , in _14th International Symposium on Radiative Corrections: Application of Quantum Field Theory to Phenomenology_ , 12, 2019.  1912.03142.  DOI.
  * (94) S. Pozzorini and E. Remiddi, _Precise numerical evaluation of the two loop sunrise graph master integrals in the equal mass case_ , _Comput. Phys. Commun._ 175 (2006) 381–387, [hep-ph/0505041].
  * (95) 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, [0705.2616].
  * (96) 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, [1709.07525].
  * (97) R. N. Lee, A. V. Smirnov and V. A. Smirnov, _Evaluating ‘elliptic’ master integrals at special kinematic values: using differential equations and their solutions via expansions near singular points_ , _JHEP_ 07 (2018) 102, [1805.00227].
  * (98) R. Bonciani, G. Degrassi, P. P. Giardino and R. Gröber, _A Numerical Routine for the Crossed Vertex Diagram with a Massive-Particle Loop_ , _Comput. Phys. Commun._ 241 (2019) 122–131, [1812.02698].
  * (99) M. Fael, F. Lange, K. Schönwald and M. Steinhauser, _A semi-analytic method to compute Feynman integrals applied to four-loop corrections to the $\overline{\mathrm{MS}}$-pole quark mass relation_, _JHEP_ 09 (2021) 152, [2106.05296].
  * (100) M. Fael, F. Lange, K. Schönwald and M. Steinhauser, _Massive Vector Form Factors to Three Loops_ , _Phys. Rev. Lett._ 128 (2022) 172003, [2202.05276].
  * (101) M. Hidding, _DiffExp, a Mathematica package for computing Feynman integrals in terms of one-dimensional series expansions_ , 2006.05510.
  * (102) X. Liu and Y.-Q. Ma, _AMFlow: a Mathematica Package for Feynman integrals computation via Auxiliary Mass Flow_ , 2201.11669.
  * (103) T. Armadillo, R. Bonciani, S. Devoto, N. Rana and A. Vicini, _Evaluation of Feynman integrals with arbitrary complex masses via series expansions_ , _Comput. Phys. Commun._ 282 (2023) 108545, [2205.03345].
  * (104) X. Liu, Y.-Q. Ma and C.-Y. Wang, _A Systematic and Efficient Method to Compute Multi-loop Master Integrals_ , _Phys. Lett. B_ 779 (2018) 353–357, [1711.09572].
  * (105) X. Liu and Y.-Q. Ma, _Multiloop corrections for collider processes using auxiliary mass flow_ , 2107.01864.
  * (106) Z.-F. Liu and Y.-Q. Ma, _Automatic computation of Feynman integrals containing linear propagators via auxiliary mass flow_ , 2201.11636.
  * (107) K. G. Chetyrkin, A. L. Kataev and F. V. Tkachov, _Higher Order Corrections to Sigma-t (e+ e- — $>$ Hadrons) in Quantum Chromodynamics_, _Phys. Lett. B_ 85 (1979) 277–279.
  * (108) S. Laporta, _High-precision calculation of multiloop Feynman integrals by difference equations_ , _Int. J. Mod. Phys. A_ 15 (2000) 5087–5159, [hep-ph/0102033].
  * (109) N. Arkani-Hamed, J. L. Bourjaily, F. Cachazo, S. Caron-Huot and J. Trnka, _The All-Loop Integrand For Scattering Amplitudes in Planar N=4 SYM_ , _JHEP_ 01 (2011) 041, [1008.2958].
  * (110) N. Arkani-Hamed, J. L. Bourjaily, F. Cachazo and J. Trnka, _Local Integrals for Planar Scattering Amplitudes_ , _JHEP_ 06 (2012) 125, [1012.6032].
  * (111) S. Badger, G. Mogull and T. Peraro, _Local integrands for two-loop all-plus Yang-Mills amplitudes_ , _JHEP_ 08 (2016) 063, [1606.02244].
  * (112) G. Barucchi and G. Ponzano, _Differential equations for one-loop generalized Feynman integrals_ , _J. Math. Phys._ 14 (1973) 396–401.
  * (113) A. V. Kotikov, _Differential equations method: The Calculation of vertex type Feynman diagrams_ , _Phys. Lett. B_ 259 (1991) 314–322.
  * (114) T. Gehrmann and E. Remiddi, _Differential equations for two loop four point functions_ , _Nucl. Phys. B_ 580 (2000) 485–518, [hep-ph/9912329].
  * (115) Z. Bern, L. J. Dixon and D. A. Kosower, _Dimensionally regulated pentagon integrals_ , _Nucl. Phys. B_ 412 (1994) 751–816, [hep-ph/9306240].
  * (116) K.-T. Chen, _Iterated path integrals_ , _Bull. Am. Math. Soc._ 83 (1977) 831–879.
  * (117) A. B. Goncharov, M. Spradlin, C. Vergu and A. Volovich, _Classical Polylogarithms for Amplitudes and Wilson Loops_ , _Phys. Rev. Lett._ 105 (2010) 151605, [1006.5703].
  * (118) S. Badger, J. Kryś, R. Moodie and S. Zoia, _Lepton-pair scattering with an off-shell and an on-shell photon at two loops in massless QED_ , _JHEP_ 11 (2023) 041, [2307.03098].