aainstitutetext:  Dipartimento di Fisica and Arnold-Regge Center, Università
di Torino, and INFN, Sezione di Torino, Via P. Giuria 1, I-10125 Torino, Italy
bbinstitutetext:  Institute for Particle Physics Phenomenology, Department of
Physics, Durham University, Durham DH1 3LE, United Kingdom

# Lepton-pair scattering with an off-shell and an on-shell photon at two loops
in massless QED

Simon Badger a,b Jakub Kryś a Ryan Moodie a Simone Zoia
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We compute the two-loop QED helicity amplitudes for the scattering of a lepton
pair with an off-shell and an on-shell photon,
$0\to\ell\bar{\ell}\gamma\gamma^{*}$, using the approximation of massless
leptons. We express all master integrals relevant for the scattering of four
massless particles with a single external off-shell leg up to two loops in a
basis of algebraically independent multiple polylogarithms, which guarantees
an efficient numerical evaluation and compact analytic representations of the
amplitudes. Analytic forms of the amplitudes are reconstructed from numerical
evaluations over finite fields. Our results complete the amplitude-level
ingredients contributing to the N3LO predictions of electron-muon scattering
$e\mu\to e\mu$, which are required to meet the precision goal of the future
MUonE experiment.

DE

    differential equation
HVP

    hadronic vacuum polarisation
IBP

    integration-by-parts
IR

    infrared
ISP

    irreducible scalar product
ISR

    initial state radiation
MI

    master integral
MPL

    multiple polylogarithm
N3LO

    next-to-next-to-next-to-leading order
NLO

    next-to-leading order
NNLO

    next-to-next-to-leading order
QCD

    quantum chromodynamics
QED

    quantum electrodynamics
RVV

    real-double-virtual
SM

    Standard Model
UV

    ultraviolet
VV

    double-virtual

##  1 Introduction

The MUonE experiment Abbiendi:2016xup ; Abbiendi:2022oks ; Spedicato:2022qtw ;
Pilato:2022wvg  will measure the hadronic running of the electromagnetic
coupling $\alpha$ using low-energy elastic electron-muon scattering, $e\mu\to
e\mu$. This will enable a new and precise determination of the hadronic vacuum
polarisation (HVP) contribution $a_{\mu}^{\text{HVP}}$ CarloniCalame:2015obs ;
Balzani:2021del  to the muon anomalous magnetic moment $a_{\mu}$. This is
required in light of the recent tensions between experimental Muong-2:2021ojo
, Standard Model (SM) data-driven Aoyama:2020ynm , and lattice quantum
chromodynamics (QCD) Borsanyi:2020mff  results for $a_{\mu}$. Increasing the
precision of the theoretical predictions for $e\mu\to e\mu$ scattering is a
high priority for the planned MUonE experiment Banerjee:2020tdt ;
Budassi:2022dco  and has seen good progress in the last few years
Alacevich:2018vez ; Budassi:2021twh ; Fael:2019nsf ; Fael:2018dmz ;
CarloniCalame:2020yoz . The recent completion of full next-to-next-to-leading
order (NNLO) quantum electrodynamics (QED) corrections Broggio:2022htr
indicates that next-to-next-to-next-to-leading order (N3LO) corrections in
differential distributions are required to meet MUonE’s precision goal of 10
parts per million. Electron-line corrections, meaning corrections to the
subprocess with the muon line stripped off ($e\to e\gamma^{*}$), are the
dominant corrections Broggio:2022htr , and a collaborative project was started
to perform their fixed-order calculation at N3LO Durham:n3lo . With the
triple-virtual corrections now available Fael:2022rgm ; Fael:2022miw ;
Fael:2023zqr , the main missing ingredient is the real-double-virtual (RVV)
matrix element ($e\to e\gamma\gamma^{*}$) at two loops. While these
contributions could be extracted from amplitudes in the literature
Garland:2001tf ; Garland:2002ak ; Gehrmann:2011ab , our direct computation
provides the massless RVV contribution in a complete and compact form.

Another application of the $0\to\ell\bar{\ell}\gamma\gamma^{*}$ amplitudes is
in electron-positron annihilation experiments precisionsm . They are required
for initial-state corrections in predictions of the ratio of hadron-to-muon
production in $e^{+}e^{-}$ collisions, which is an important input for
existing SM predictions of $a^{\text{HVP}}_{\mu}$ Abbiendi:2022liz . The two-
loop amplitudes contribute to RVV corrections to $e^{+}e^{+}\to\gamma^{*}$ in
direct scan measurements, while radiative return measurements concern
corrections to $e^{+}e^{-}\to\gamma\gamma^{*}$ Aoyama:2020ynm . In the latter
configuration, the $e^{+}e^{-}$ beam has a fixed centre-of-mass energy of a
few GeV and the on-shell photon originates from initial state radiation (ISR).
The energy lost to the ISR photon is used to effectively scan over the
energies of the decay of the off-shell photon. A differential cross section
of, for example, $\gamma^{*}\to\text{hadrons}$ with respect to the centre-of-
mass energy of the decay, $\text{d}\sigma/\text{d}s$, can be extracted from
measurements of the differential cross section with respect to the energy of
the ISR photon, $\text{d}\sigma/\text{d}E_{\gamma}$. State-of-the-art
predictions for these measurements are currently at next-to-leading order
(NLO) Abbiendi:2022liz . We provide the two-loop
$e^{+}e^{-}\to\gamma\gamma^{*}$ amplitudes required for the double-virtual
(VV) corrections at NNLO, although the bottleneck remains in the hadronic
decay.

Our amplitudes are calculated in the approximation of massless leptons. In the
NNLO massive $e\mu\to e\mu$ cross section calculation Broggio:2022htr , the
authors obtain photonic corrections (those with no closed fermion loops),
using a small-mass expansion Penin:2005eh ; Becher:2007cu ; Engel:2018fsb
applied to the two-loop amplitudes with massless electrons for the VV
corrections. This approximation relies on the electron mass being much smaller
than any other scale, which is valid in the bulk of phase space. Further
splitting the photonic corrections, they take the subset of electron-line
corrections and find that the relative difference to the true massive NNLO
differential cross section is generally around $10^{-3}\alpha^{2}$, where
$\alpha$ is the fine-structure constant, which is negligible compared to the
$10^{-5}$ precision goal. The approximation breaks down in soft and collinear
regions, where they treat the amplitudes using infrared (IR) factorisation
Banerjee:2021mty ; Engel:2021ccn ; Engel:2023ifn , and is not used for
contributions including closed fermion loops Engel:2018fsb ; Engel:2019nfw .
Our amplitudes can be used analogously for the RVV corrections at N3LO.

Our computation uses the modern technology developed for QCD amplitudes with
many scales. The high-multiplicity amplitude frontier in massless QCD lies
with two-loop five-particle processes, with leading-colour Abreu:2019odu ;
Abreu:2020cwb ; Chawdhry:2020for ; Agarwal:2021grm ; Abreu:2021oya ;
Chawdhry:2021mkw  and full-colour Agarwal:2021vdh ; Badger:2021imn ;
Badger:2023mgf ; Abreu:2023bdp  results in a form ready for phenomenological
application becoming available over the past few years. Recently, the first
single-external-mass calculations are also appearing Hartanto:2019uvl ;
Badger:2021nhg ; Badger:2021ega ; Badger:2022ncb . These computations have
made extensive use of finite-field arithmetic to sidestep large intermediate
expressions. This technology has had a considerable impact for solutions of
systems of integration-by-parts (IBP) identities vonManteuffel:2014ixa ;
Klappert:2020nbg ; Smirnov:2019qkx  but also applies more widely to scattering
amplitude computations Peraro:2016wsq ; Peraro:2019svx . Motivated by the
improved algorithms, we choose to implement a complete finite-field based
reduction for the $2\to 2$ processes with an off-shell leg. Since the
kinematics are relatively simple in comparison with other high-multiplicity
configurations, this technology is not essential. It does, however, provide an
opportunity to review the new techniques for readers who are not familiar with
them.

A key ingredient for computing the scattering amplitudes are analytic
expressions for the required Feynman integrals. Complete analytic results up
to two loops are already available in the literature gehrmann:2000zt ;
gehrmann:2001ck ; Gehrmann:2002zr ; Duhr:2012fh ; Gehrmann:2023etk .
Expansions of these integrals up to higher orders in the dimensional
regularisation parameter $\epsilon$ have also been reconsidered recently
Gehrmann:2023etk , in view of their usage for N3LO corrections to $2\to 2$
processes in QCD Gehrmann:2022vuk ; Gehrmann:2023zpz . The state of the art
for integrals with this kinematic configuration has reached three loops
DiVita:2014pza ; Canko:2020gqp ; Canko:2021xmn ; Henn:2023vbd . We revisit the
computation of the one- and two-loop integrals following the approach of refs.
Gehrmann:2018yef ; Chicherin:2020oor ; Badger:2021nhg ; Chicherin:2021dyp ;
Abreu:2023rco  based on the construction of a _basis_ of independent special
functions, which gives a unique and uniform representation of all the required
Feynman integrals up to transcendental weight four. This enables a more
efficient computation of the amplitudes using the modern workflow based on
finite-field arithmetic, and leads to more compact expressions. We give
explicit expressions for the basis functions in terms of multiple
polylogarithms which can be evaluated in an efficient and stable way
throughout the physical phase space. We compute all crossings of all massless
one- and two-loop four-particle Feynman integrals with an external off-shell
leg, so that our results for the integrals may be of use for any scattering
process with these kinematics.

Our paper is organised as follows. In section 2, we describe our decomposition
of the helicity amplitudes and detail how we express the off-shell currents.
In section 3, we discuss our computation of analytic amplitudes by numerical
evaluations over finite fields. In section 4, we present the computation of
the Feynman integrals in terms of a basis of special functions. We draw our
conclusions in section 5. We provide useful technical details in appendices.
We define the relevant families of Feynman integrals in appendix A. In
appendix B, we discuss in detail how we handle permutations of the integral
families in the IBP reduction. In appendix C, we describe our rational
parametrisation of the kinematics. Appendix D is devoted to the ultraviolet
(UV) renormalisation and IR factorisation which determine the pole structure
of the amplitudes. In appendix E we discuss the analytic continuation of the
special functions to the physical kinematic regions.

##  2 Structure of the amplitude

We calculate the one- and two-loop QED corrections to the process

$\displaystyle
0\to\ell(p_{1},h_{1})+\bar{\ell}(p_{2},h_{2})+\gamma(p_{3},h_{3})+\gamma^{*}(p_{4})\,,$
(1)

which we call $0\to\ell\bar{\ell}\gamma\gamma^{*}$ for short. Here, $\ell$
denotes an on-shell massless lepton and $\gamma$ ($\gamma^{*}$) an on-shell
(off-shell) photon, while $h_{i}$ and $p_{i}$ are the helicity and momentum of
the $i\textsuperscript{th}$ particle. We take the external momenta $p_{i}$ to
be all outgoing. They satisfy the following momentum-conservation and on-shell
conditions:

$\displaystyle\sum_{i=1}^{4}p_{i}^{\mu}=0\,,\qquad\qquad\qquad
p_{i}^{2}=0\quad\forall\,i=1,2,3\,.$ (2)

The single-off-shell four-particle phase space is described by three
independent scalar invariants, which we choose as

$\displaystyle\vec{s}\coloneqq\\{s_{12},s_{23},s_{4}\\}\,,$ (3)

where $s_{i\ldots j}\coloneqq(p_{i}+\ldots+p_{j})^{2}$. We use dimensional
regularisation in the ’t Hooft-Veltman scheme Gnendiger:2017pys , with
$D=4-2\epsilon$ spacetime dimensions (where $\epsilon$ is the dimensional
regulator) and four-dimensional external momenta.

Because of the off-shell photon in the process, the helicity amplitudes
$\mathcal{A}^{\mu}(1_{\ell},2_{\bar{\ell}},3_{\gamma},4_{\gamma^{*}})$ are
actually off-shell currents carrying a free Lorentz index. We consider the
perturbative QED expansion of the helicity amplitudes,

$\displaystyle\mathcal{A}^{\mu}(1_{\ell},2_{\bar{\ell}},3_{\gamma},4_{\gamma^{*}})=g_{e}^{2}\sum_{L\geq
0}\left(n_{\epsilon}\frac{\alpha}{4\pi}\right)^{L}{\mathcal{A}^{(L)}}^{\mu}(1_{\ell},2_{\bar{\ell}},3_{\gamma},4_{\gamma^{*}})\,,$
(4)

with prefactor
$n_{\epsilon}=\mathrm{i}(4\pi)^{\epsilon}\mathrm{e}^{-\epsilon\gamma_{E}}$,
electromagnetic coupling $g_{e}$, and $\alpha=g_{e}^{2}/(4\pi)$. We truncate
the expansion at $L=2$ loops. We set the renormalisation scale $\mu_{R}$ to
one throughout the computation and restore the dependence on it in the final
analytic result by dimensional analysis. For the bare amplitudes we have that

$\displaystyle{\mathcal{A}^{(L)}}^{\mu}(\mu_{R})$
$\displaystyle=\mu_{R}^{2\epsilon L}{\mathcal{A}^{(L)}}^{\mu}(\mu_{R}=1)\,.$
(5)

There are two independent helicity configurations $(h_{1},h_{2},h_{3})$, which
we take as

$\displaystyle\\{-+-\,,\ \ -++\\}\,.$ (6)

We derive the analytic expressions for these helicity amplitudes. We obtain
the remaining helicity configurations, $\\{+-+\,,\ \ +--\\}$, through parity
transformation (see appendix C of ref. Badger:2023mgf ).

(a) ${\mathcal{A}^{(1)}_{0,0}}^{\mu}$

(b) ${\mathcal{A}^{(1)}_{1,1}}^{\mu}$

(c) ${\mathcal{A}^{(2)}_{0,0}}^{\mu}$

(d) ${\mathcal{A}^{(2)}_{1,0}}^{\mu}$

(e) ${\mathcal{A}^{(2)}_{1,1}}^{\mu}$

(f) ${\mathcal{A}^{(2)}_{1,2}}^{\mu}$

(g) ${\mathcal{A}^{(2)}_{2,1}}^{\mu}$

Figure 1:  Representative Feynman diagrams for the subamplitudes defined in
eq. 7. The off-shell external leg is indicated by a bold line.

We decompose the loop-level helicity amplitudes ${\mathcal{A}^{(L)}}^{\mu}$
into gauge-invariant subamplitudes ${\mathcal{A}^{(L)}_{i,j}}^{\mu}$, where
the subscript $i$ counts the number of closed massless fermion loops and $j$
the number of external photons attached to closed fermion loops. The non-zero
contributions are

$\displaystyle{\mathcal{A}^{(1)}}^{\mu}$
$\displaystyle={\mathcal{A}^{(1)}_{0,0}}^{\mu}+n_{l}\,{\mathcal{A}^{(1)}_{1,1}}^{\mu}\,,$
(7a) $\displaystyle{\mathcal{A}^{(2)}}^{\mu}$
$\displaystyle={\mathcal{A}^{(2)}_{0,0}}^{\mu}+n_{l}\left({\mathcal{A}^{(2)}_{1,0}}^{\mu}+{\mathcal{A}^{(2)}_{1,1}}^{\mu}+{\mathcal{A}^{(2)}_{1,2}}^{\mu}\right)+n_{l}^{2}{\mathcal{A}^{(2)}_{2,1}}^{\mu}\,,$
(7b)

where $n_{l}$ denotes the number of charged lepton flavours running in the
loops. Representative Feynman diagrams contributing to these subamplitudes are
illustrated in figure 1. Amplitudes with a closed fermion loop attached to an
odd number of photons vanish by Furry’s theorem.

We decompose the amplitude and subamplitude currents as

$\displaystyle{\mathcal{A}^{(L)}}^{\mu}=\sum_{k=1}^{4}a_{k}^{(L)}\,q_{k}^{\mu}\,,\qquad\qquad{\mathcal{A}^{(L)}_{i,j}}^{\mu}=\sum_{k=1}^{4}a_{i,j;k}^{(L)}\,q_{k}^{\mu}\,,$
(8)

using the following basis written with the spinor-helicity formalism:

$\displaystyle q_{k}^{\mu}$
$\displaystyle=p_{k}^{\mu}\quad\forall\,k=1,2,3\,,$ $\displaystyle
q_{4}^{\mu}$ $\displaystyle=\frac{\langle 2|p_{3}p_{1}\sigma^{\mu}|2]-\langle
1|p_{3}p_{2}\sigma^{\mu}|1]}{2s_{12}}\,.$ (9)

Readers not familiar with the spinor-helicity formalism may like to consult
one of the many good reviews on the subject Mangano:1990by ; Dixon:1996wi ;
Badger:2023eqz . Note that $q_{4}$ is orthogonal to the momenta $p_{i}$ by
construction; one can in fact show that
$q_{4}^{\mu}\propto\varepsilon^{\mu\nu\rho\sigma}{q_{1}}_{\nu}{q_{2}}_{\rho}{q_{3}}_{\sigma}$.
The subamplitude coefficients $a_{i,j;k}^{(L)}$ can be related to the
amplitude ones $a_{k}^{(L)}$ through eq. 7.

The scattering amplitudes $\mathcal{M}^{(L)}$ for fully on-shell processes
(for instance, for $0\to e^{-}e^{+}\gamma\mu^{-}\mu^{+}$) are obtained by
contracting the amplitude currents ${\mathcal{A}^{(L)}}^{\mu}$ (for $0\to
e^{-}e^{+}\gamma\gamma^{*}$) with a suitable decay current $\mathcal{V}_{\mu}$
(in this example, $\gamma^{*}\to\mu^{-}\mu^{+}$), as

$\displaystyle\mathcal{M}^{(L)}$
$\displaystyle\coloneqq\mathcal{A}^{(L)}\cdot\mathcal{V}=\sum_{k=1}^{4}a_{k}^{(L)}\,\left(q_{k}\cdot\mathcal{V}\right)\,.$
(10)

In this manner, the on-shell amplitudes $\mathcal{M}^{(L)}$ are given by the
scalar product between the vector of coefficients
$(a_{1}^{(L)},\ldots,a_{4}^{(L)})$, and that of decay-vector contractions
$(q_{1}\cdot\mathcal{V},\ldots,q_{4}\cdot\mathcal{V})$. The coefficients
$a^{(L)}_{k}$ depend on the helicities of the three on-shell particles in eq.
1, while the decay vector $\mathcal{V}_{\mu}$ depends on the helicities of the
particles the off-shell photon decays to. The helicity-summed interference
between the $L_{1}$-loop and the $L_{2}$-loop matrix elements is then given by

$\displaystyle\mathcal{M}^{(L_{1},L_{2})}$
$\displaystyle=\frac{1}{4}\sum_{\vec{h}}{\mathcal{M}^{(L_{1})}_{\vec{h}}}^{*}\mathcal{M}^{(L_{2})}_{\vec{h}}\,,$
(11)

where the subscripts $\vec{h}$ indicates the helicities of all on-shell
particles — that is, including the decay products of the off-shell photon —
and the overall constant factor averages over the helicities of the incoming
particles.

The output of the computation described in section 3 is the set of four
projections $\mathcal{A}^{(L)}_{i,j}\cdot q_{k}$ for each helicity
configuration listed in eq. 6. From these, we determine the subamplitude
coefficients $a_{i,j;k}^{(L)}$ by inverting eq. 8, as

$\displaystyle a_{i,j;k}^{(L)}$
$\displaystyle=\sum_{m=1}^{4}\left(\mathrm{G}^{-1}\right)_{km}\left(\mathcal{A}^{(L)}_{i,j}\cdot
q_{m}\right)\,,$ (12)

where $\mathrm{G}$ is the Gram matrix of the vectors $q_{i}$, that is, the
matrix of entries $\mathrm{G}_{ij}\coloneqq q_{i}\cdot q_{j}$ for
$i,j=1,\ldots,4$. At loop level, we write the subamplitude coefficients as

$\displaystyle a_{i,j;k}^{(L)}$
$\displaystyle=\sum_{w=-2L}^{4-2L}\sum_{r}\,c_{r,w}\,\text{mon}_{r}(F)\,\epsilon^{w},$
(13)

where $\text{mon}_{r}(F)$ are monomials of special functions $F$ (see section
4), and the coefficients $c_{r,w}$ are rational functions of the kinematics.
We drop the dependence on $i$, $j$, $k$, and $L$ on the right-hand side of eq.
13 for compactness. We truncate the Laurent expansion around $\epsilon=0$ to
the orders required for computing NNLO predictions. We express the
coefficients $c_{r,w}$ as $\mathbb{Q}$-linear combinations of a smaller set of
linearly-independent coefficients (see section 3). The analytic expressions of
the latter are given explicitly in terms of momentum twistor variables (see
appendix C). We simplify these expressions through a multivariate partial
fraction decomposition using MultivariateApart Heller:2021qkz , and by
collecting the common factors.

In the ancillary files zenodo , the directory amplitudes/ contains Mathematica
files describing the bare helicity subamplitude currents
${\mathcal{A}^{(L)}_{i,j}}^{\mu}$ by their coefficients $a_{i,j;k}^{(L)}$ in
the form of eq. 13. The Mathematica script current.m is a reference
implementation of the numerical evaluation of the bare amplitude coefficients
$a_{k}^{(L)}$ in eq. 13, including summation of subamplitudes in eq. 7,
treatment of dependent helicities, and renormalisation scale restoration in
eq. 5. The Mathematica script evaluation.wl demonstrates the construction of
the five-particle on-shell amplitudes in eq. 10 for the process $0\to
e^{-}e^{+}\gamma\mu^{-}\mu^{+}$, and their helicity-summation to obtain the
squared matrix elements in eq. 11. The results of the script are checked
against a reference point included in reference_point.json.

We perform the following checks of our amplitudes.

Ward identity

We verify the gauge invariance of the subamplitudes
${\mathcal{A}^{(L)}_{i,j}}^{\mu}$ by checking that they vanish on replacing
the on-shell photon’s polarisation vector with its momentum.

One-loop crosscheck

We successfully crosscheck our one-loop $n_{l}=0$ helicity-summed matrix
element contracted with the decay $\gamma^{*}\to\mu^{-}\mu^{+}$ against the
QED NLO electron-line corrections for $e\mu\to e\mu\gamma$ obtained with
McMule Banerjee:2020rww ; ulrich_yannick_2022_6046769 .

Finite remainder

We verify that the $\epsilon$-poles of the bare amplitudes have the structure
predicted by UV renormalisation and IR factorisation Catani:1998bh ;
Gardi:2009qi ; Gardi:2009zv ; Becher:2009cu ; Becher:2009qa . We then subtract
the expected poles and define finite remainders at one and two loops as

$\displaystyle{\mathcal{F}^{(1)}}^{\mu}$
$\displaystyle=\left[{\mathcal{A}^{(1)}}^{\mu}-\frac{3}{2}\frac{\beta_{0}}{\epsilon}{\mathcal{A}^{(0)}}^{\mu}\right]-Z^{(1)}{\mathcal{A}^{(0)}}^{\mu}\,,$
(14a) $\displaystyle{\mathcal{F}^{(2)}}^{\mu}$
$\displaystyle=\left[{\mathcal{A}^{(2)}}^{\mu}-\frac{5}{2}\frac{\beta_{0}}{\epsilon}{\mathcal{A}^{(1)}}^{\mu}-\left(-\frac{15}{8}\frac{\beta_{0}^{2}}{\epsilon^{2}}+\frac{3}{4}\frac{\beta_{1}}{\epsilon}\right){\mathcal{A}^{(0)}}^{\mu}\right]-Z^{(2)}{\mathcal{A}^{(0)}}^{\mu}-Z^{(1)}{\mathcal{F}^{(1)}}^{\mu}\,,$
(14b)

where the square brackets separate renormalisation of UV poles from
subtraction of IR poles. We present the derivation of these formulae in
appendix D.

##  3 Setup of the calculation

In this section, we outline the workflow we used to calculate our amplitudes.
Firstly, we generate all Feynman diagrams contributing to eq. 1 using QGRAF
Nogueira:1991ex . Each diagram is then replaced with the corresponding Feynman
rules for vertices, propagators, and external states, leading to a collection
of $D$-dimensional Feynman integrals. Next, we filter the integrals according
to eqs. 4 and 7 using a collection of Mathematica and FORM scripts
Kuipers:2012rf ; Ruijl:2017dtg . Within each subamplitude
${\mathcal{A}^{(L)}_{i,j}}^{\mu}$, we then collect the integrals according to
their topology, by which we mean a unique set of denominators. For example,
the diagrams in figures 1(d) and 1(e) belong to different topologies, but
those in figures 1(c) and 1(f) belong to the same topology (under the
assumption of massless lepton propagators). At this point, the subamplitudes
are sums of Feynman integrals over distinct integral topologies, with the
numerators given by linear combinations of monomials that depend on the loop
as well as the external momenta. To work with the projected helicity
subamplitudes $\mathcal{A}^{(L)}_{i,j}\cdot q_{k}$, we specify the
polarisations of external particles according to eq. 6, as well as the
projector $q_{k}^{\mu}$ of the off-shell photon from eq. 9. It is natural to
express helicity-dependent objects using the spinor-helicity formalism. Then,
the monomials of loop momenta contain the following scalar products and spinor
strings:

$\left\\{k_{i}\cdot k_{j},\,k_{i}\cdot p_{j},\,\langle
ij\rangle,\,[ij],\,\langle i|k_{i}|j],\,\langle i|p_{4}|j],\,\langle
i|k_{i}p_{4}|j\rangle,\,[i|k_{i}p_{4}|j]\right\\}\,.$ (15)

Their coefficients, on the other hand, are composed of the same type of
objects, but do not contain any dependence on loop momenta $k_{i}$. We express
these coefficients using the rational parametrisation of the kinematics
discussed in appendix C.

This marks the start of our finite-field sampling procedure Peraro:2016wsq .
The goal of this approach is to sidestep the algebraic complexity which
typically plagues the intermediate stages of symbolic computations by
evaluating numerically all rational coefficients. Using integers modulo some
large prime number — which constitute a finite field — for the numerical
evaluation allows us to avoid the loss of accuracy inherent to floating-point
numbers, as well as the expensive arbitrary-precision arithmetic required by
rational numbers. Manipulations needed to further process the rational
coefficients are a completely separate problem from the calculation of the
integrals or special functions that these coefficients multiply. In fact, they
can be implemented as a series of rational transformations over finite fields.
We stress that this is the methodology we follow at each step of the
computation described below. In particular, we use the package FiniteFlow
Peraro:2019svx , which is conveniently interfaced to Mathematica. The analytic
form of the coefficients is not known at any intermediate step. It is
reconstructed from the finite-field samples only at the very end of the
workflow.

(a) Penta-triangles

(b) Double-box

(c) Crossed double-boxes

Figure 2:  The two-loop amplitudes include six ordered integral families:
three penta-triangles, a double-box, and two crossed double-boxes. All are
planar except for the crossed double-boxes. The off-shell external leg is
indicated by a bold line. External legs have outgoing momenta.

Firstly, we note that not all integral topologies are independent: some of
them can be written as subtopologies of others. For this reason, we define the
set of maximal topologies, i.e. topologies with the maximum number of
propagators allowed for $L$-loop, $n$-particle diagrams. In figure 2, we
present the maximal topologies for the process under consideration in an
arbitrary ordering of the external momenta (we give their explicit definitions
in appendix A). Several orderings of the external momenta are relevant for the
amplitudes, and we treat them as distinct families. Next, we map all
topologies present so far onto one of these maximal topologies. The loop
momenta dependent objects of eq. 15 are then expressed through the nine
inverse propagators and irreducible scalar products associated with the chosen
maximal topology. In this way, each subamplitude is now a sum of integrals
compatible with IBP reduction Tkachov:1981wb ; Chetyrkin:1981qh , while their
coefficients depend on the external kinematics and $\epsilon$. We generate the
required IBP relations using LiteRed Lee:2012cn . The resulting IBP system is
then solved using the Laporta algorithm Laporta:2001dd  with FiniteFlow’s
linear solver to yield the reduction of all the integrals present within our
maximal topologies onto a much smaller subset of master integrals. It is
beneficial to choose the MIs such that they satisfy differential equations in
canonical form Henn:2013pwa  (see section 4.1). This property allows us to
make the most of the optimisation measures we discuss below. We stress that
the IBP reduction is also done numerically over finite fields, since the
coefficients of the IBP relations are rational functions of external
kinematics and the dimensional regulator $\epsilon$. This is an important
simplification, since analytic IBP reduction often proves to be the bottleneck
of amplitude computations. For many amplitude applications, multiple
permutations of the ordered topologies can appear. We outline a strategy to
optimise the reduction in such situations in appendix B.

At this point, each projected helicity subamplitude
$\mathcal{A}^{(L)}_{i,j}\cdot q_{k}$ is written as a linear combination of MIs
multiplied by rational coefficients of $\epsilon$ and the kinematic variables.
We now write the MIs in terms of a basis of special functions up to the
required order in $\epsilon$ (see section 4). Finally, we Laurent expand the
amplitude around $\epsilon=0$, the deepest pole being $1/\epsilon^{2L}$ at $L$
loops. The only task left is to reconstruct the rational coefficients of the
special-function monomials from their samples over finite fields. In general,
this might be a daunting challenge and its complexity stems from two separate
factors. The workflow described so far is a series of rational operations
chained together within a so-called dataflow graph Peraro:2019svx . As such,
we essentially have a black-box algorithm which takes numerical values of the
kinematic variables as input, and returns the corresponding numerical values
of the rational coefficients of the special-function monomials. The first
factor is that several sample points are necessary to infer the analytic
expression of these coefficients from their values in the finite fields. The
required number is correlated with the polynomial degrees of the rational
functions viewed as ratios of polynomials: the higher the degree, the more
sample points are required. The second factor affecting the reconstruction
complexity is the time it takes to obtain the values of the coefficients at
each sample point. The more complicated the dataflow graph is, i.e. the more
operations are chained together and the more difficult each operation is, the
longer it will take to run the black-box algorithm. The most expensive
operation in this regard is the evaluation of the solution to the IBP system.
The total reconstruction time can thus be estimated as:

$\text{reconstruction time}\approx(\text{number of sample
points})\times(\text{evaluation time per point})\,.$ (16)

We emphasise that the sample evaluations can be run in parallel. For a
detailed discussion of various strategies to improve the reconstruction time,
see section 4 of ref. Badger:2021imn  and section 3.5 of ref. Badger:2022ncb .
Here, we give a brief overview of the tools that proved sufficient for this
work.

First, we look for $\mathbb{Q}$-linear relations among the rational
coefficients of each helicity subamplitude. This typically requires few sample
points with respect to the full reconstruction. We then solve these linear
relations to express all coefficients in terms of a minimal subset of
independent ones. Only the latter need to be reconstructed. Choosing them so
that they have the lowest degrees often leads to a decrease in the complexity
of the reconstruction.

The second strategy we employ is to match the rational coefficients with
factorised ansätze informed by the singularity structure of Feynman integrals.
The singularities of Feynman integrals can in fact be read off from the DEs
they satisfy. For each coefficient we then write an ansatz made of the
following factors:

$\left\\{\langle ij\rangle,\,[ij],\,\langle
i|p_{4}|j],\,s_{ij}-s_{k4},\,s_{i4}-s_{4},\,s_{4}\right\\}\,,$ (17)

for all $i,j,k=1,2,3$ such that $i\neq j\neq k$. This list includes
denominator factors of the DEs satisfied by the MIs (listed by eq. 19), as
well as spinor structures aimed at capturing the phase information of helicity
amplitudes. We then determine the exponents of the ansätze by comparing them
to the coefficients reconstructed on a univariate slice of the kinematic
variables Abreu:2019odu , which are very cheap to obtain. We find that with
this ansatz it is possible to determine all denominator factors — which indeed
are linked to the singularity structure of the amplitude — and sometimes also
some numerator factors. As a result, the undetermined functions yet to be
reconstructed are of lower degree and require fewer sample points. We
reconstruct the analytic form of the remaining rational functions using
FiniteFlow’s built-in multivariate functional reconstruction algorithm.

Finally, we note that, for more computationally demanding processes, further
ansatz-based techniques — for instance based on partial fraction
decompositions or informed by the singularity structure of the amplitudes —
may be used to optimise the functional reconstruction; see, for example, refs.
Badger:2021imn ; Badger:2021ega ; Badger:2022ncb ; DeLaurentis:2022otd ;
Badger:2023mgf ; Abreu:2023bdp ; Liu:2023cgs .

##  4 Computation of the master integrals

The MIs for the relevant integral families were first computed analytically in
refs. gehrmann:2000zt ; gehrmann:2001ck . The basis of special functions for
four-particle kinematics with an off-shell leg has been subsequently studied
in great detail with attention to the compactness of representations in terms
of polylogarithms and numerical evaluation across a complete phase-space
Duhr:2012fh ; Gehrmann:2013vga ; Gehrmann:2023etk  (see also ref.
Gehrmann:2002zr  for a thorough discussion of the analytic continuation). We
revisit this computation to obtain expressions for the MIs which are well
suited for the amplitude-computation workflow discussed in section 3. To this
end, we compute the MIs for _all permutations_ of the external legs in terms
of a _basis_ of special functions, following the approach of refs.
Gehrmann:2018yef ; Chicherin:2020oor ; Badger:2021nhg ; Chicherin:2021dyp ;
Abreu:2023rco . In other words, we express all the Feynman integrals
contributing to the amplitudes in terms of a set of special functions which
are (algebraically) independent. Having such a unified and unique
representation for all permutations of the integral families allows for
simplifications and cancellations among different permutations of the Feynman
integrals. This leads to a simpler expression of the amplitudes and to a more
efficient functional reconstruction in the finite-field setup presented in
section 3. We emphasise that our results cover all MIs required for computing
_any_ two-loop four-particle amplitude with a single external off-shell leg,
and not just the ones required for the amplitudes presented in this work.

We discuss the construction of the basis in section 4.1, and how we express it
in terms of multiple polylogarithms (MPLs) in section 4.2. Finally, we give
some details about the numerical evaluation and the checks we performed in
section 4.3.

###  4.1 Construction of the special function basis

We follow the strategy presented in ref. Abreu:2023rco . The starting point
are the DEs satisfied by the MIs for each family Barucchi:1973zm ;
KOTIKOV1991158 ; KOTIKOV1991123 ; Gehrmann:1999as ; Bern:1993kr . Let $\tau$
label an integral family, e.g. the double-box in figure 2(b) for an arbitrary
permutation of the external massless momenta. We choose a basis of _pure_ MIs
$\vec{g}_{\tau}$, that is, a basis which satisfies DEs in the canonical form
Henn:2013pwa

$\displaystyle\text{d}\,\vec{g}_{\tau}(\vec{s};\epsilon)=\epsilon\,\left(\sum_{i=1}^{7}\mathrm{A}^{(\tau)}_{i}\,\text{d}\log
W_{i}(\vec{s})\right)\cdot\vec{g}_{\tau}(\vec{s};\epsilon)\,.$ (18)

Here, d is the total differential,
$\text{d}f\coloneqq\text{d}s_{12}\,\partial_{s_{12}}f+\text{d}s_{23}\,\partial_{s_{23}}f+\text{d}s_{4}\,\partial_{s_{4}}f$,
$\mathrm{A}^{(\tau)}_{i}$ are constant $n_{\tau}\times n_{\tau}$ matrices,
with $n_{\tau}$ the number of MIs of the family $\tau$, and

$\displaystyle\begin{aligned}
&W_{1}=s_{12}\,,&&W_{2}=s_{23}\,,&&W_{3}=s_{12}+s_{23}\,,\qquad&&W_{4}=s_{12}-s_{4}\,,\\\
&W_{5}=s_{23}-s_{4}\,,\qquad&&W_{6}=s_{12}+s_{23}-s_{4}\,,\qquad&&W_{7}=s_{4}\,&&\end{aligned}$
(19)

are called _letters_. We emphasise that this alphabet covers all permutations
of the relevant integrals. Specific integral families may contain only subsets
of it. Canonical bases for the relevant integral families are already
available in the literature Canko:2021xmn ; Gehrmann:2023etk ; Henn:2023vbd .
Given that, by today’s standards, finding canonical bases for these integral
families is simple, we re-derived them using a mixture of methods: the package
DlogBasis Henn:2020lye , the analysis of results in the literature for related
integral families (massless two-loop five-point planar integrals
Gehrmann:2015bfy  and two-loop four-point integrals with two massive external
legs Henn:2014lfa ; Caola:2014lpa ), and a set of heuristic rules (see e.g.
ref. Dlapa:2022nct ). We normalise the MIs such that their expansion around
$\epsilon=0$ starts from order $\epsilon^{0}$,

$\displaystyle\vec{g}_{\tau}(\vec{s};\epsilon)=\sum_{w\geq
0}\epsilon^{w}\,\vec{g}^{(w)}_{\tau}(\vec{s})\,.$ (20)

For the purpose of computing two-loop scattering amplitudes up to their finite
part (i.e., up to order $\epsilon^{0}$), it suffices to restrict our attention
to $w\leq 4$. Since the MIs satisfy canonical DEs 18, the $\epsilon$-order of
the MI coefficients $\vec{g}^{(w)}_{\tau}(\vec{s})$ equals their
_transcendental weight_ Henn:2013pwa . We compute the derivatives of the MIs
using FiniteFlow Peraro:2019svx  and LiteRed Lee:2012cn . We do so only for
the integral families with the ordering of the external momenta shown in
figure 2, and obtain those for all other orderings of the external massless
legs by permutation. We provide the definition of the pure MIs and the
corresponding DEs for all one- and two-loop four-point one-mass families in
figure 2 in the folder pure_mi_bases/ of our ancillary files zenodo .

In order to solve the DEs 18 we need boundary values, i.e., values of all MIs
up to order $\epsilon^{4}$ at a phase-space point. Due to the simplicity — by
today’s standards — of the integrals under consideration, an arbitrary (non-
singular) phase-space point would do. Nonetheless, we make a more refined
choice following some of the criteria of refs. Chicherin:2020oor ;
Chicherin:2021dyp . We choose the following point in the $s_{12}$ channel (see
appendix E),

$\displaystyle\vec{s}_{0}=\left(2,\,-\frac{1}{2},\,1\right)\,,$ (21)

motivated by two principles: that it is symmetric under the permutations which
preserve the $s_{12}$ channel (i.e., swapping $p_{1}\leftrightarrow p_{2}$),
and that it contains few distinct prime factors. The first condition reduces
the number of permuted integral families we need to evaluate in order to
obtain the boundary values. The second condition reduces the number of
independent transcendental constants appearing in the boundary values, which
simplifies the construction of the basis of special functions. The
order-$\epsilon^{0}$ boundary values $\vec{g}_{\tau}^{(0)}$ are rational
constants. We obtain them up to their overall normalisation by solving the
‘first-entry conditions’ Gaiotto:2011dt , i.e., by requiring the absence of
unphysical branch cuts in the solutions. We fix the overall normalisation and
the higher-order boundary values $\vec{g}_{\tau}^{(w)}(\vec{s}_{0})$ (for
$1\leq w\leq 4$) by evaluating all MIs with AMFlow Liu:2022chg  (interfaced to
FiniteFlow Peraro:2019svx  and LiteRed Lee:2012cn ) at $\vec{s}_{0}$ with at
least $60$-digit precision. We anticipate from section 4.2 that, although we
use floating-point boundary values, our results in terms of MPLs are fully
analytic.

The canonical DEs 18 and the boundary values for all integral families are the
input for the algorithm of ref. Abreu:2023rco  for constructing a basis of
special functions. We refer to the original work for a thorough discussion.
Out of all MI coefficients up to transcendental weight $4$, the algorithm
selects a subset, denoted $F\coloneqq\\{F^{(w)}_{i}(\vec{s})\\}$, which
satisfy two constraints. First, they are _algebraically independent_ , that
is, there are no polynomial functional relations among them. Second, the MI
coefficients of all families (including all permutations of the external
massless legs) up to transcendental weight $4$ are expressed as polynomials in
the $\\{F^{(w)}_{i}(\vec{s})\\}$ and the zeta values $\zeta_{2}=\pi^{2}/6$ and
$\zeta_{3}$. For example, an arbitrary weight-$2$ MI coefficient
$g^{(2)}(\vec{s})$ has the general form

$\displaystyle
g^{(2)}(\vec{s})=\sum_{i=1}^{3}c_{i}\,F_{i}^{(2)}(\vec{s})+\sum_{i\leq
j=1}^{4}d_{ij}\,F_{i}^{(1)}(\vec{s})\,F_{j}^{(1)}(\vec{s})+e\,\zeta_{2}\,,$
(22)

with $c_{i},d_{ij},e\in\mathbb{Q}$. This special subset of MI coefficients,
$\\{F^{(w)}_{i}(\vec{s})\\}$, constitutes our special function basis. We give
the number of functions in the basis in table 1. Note that there is freedom in
the choice of which MI coefficients make up the basis. We make use of this
freedom to choose as many basis-elements as possible from the one-loop family,
then complement them with coefficients from the planar two-loop families, and
finally complete them with coefficients from the non-planar two-loop families.
In this way no two-loop MI coefficients appear in the one-loop amplitudes, and
no non-planar two-loop MI coefficients appear in those amplitudes where only
planar diagrams contribute (as is often the case in the leading colour
approximation of QCD).

Weight | Number of basis functions
---|---
$1$ | 4
$2$ | 3
$3$ | 20
$4$ | 67
Table 1: Number of functions $\\{F^{(w)}_{i}\\}$ in the basis weight by
weight.

The folder mi2func/ of our ancillary files zenodo  contains the expression of
all MI coefficients (for all one- and two-loop integral families in all
permutations of the external massless legs) up to weight $4$ in terms of our
special function basis. This result enables the efficient amplitude-
computation strategy based on finite-field arithmetic discussed in section 3.
However, at this stage the basis functions $\\{F^{(w)}_{i}\\}$ are expressed
in terms of Chen iterated integrals Chen:1977oja  and numerical boundary
values $\vec{g}^{(w)}(\vec{s}_{0})$. This representation is excellent for
investigating the analytic properties of Feynman integrals and amplitudes, but
it is not readily suitable for an efficient numerical evaluation. In the next
section we discuss how we construct a representation of the function basis in
terms of MPLs and zeta values, which is well suited for an efficient and
stable numerical evaluation.

###  4.2 Expression in terms of multiple polylogarithms

In this section we construct a representation of our function basis in terms
of MPLs. Note that, while the basis functions $\\{F^{(w)}_{i}\\}$ are by
construction algebraically independent, the MPLs appearing in their
representation constructed in this section may not be. Nonetheless, we take a
number of simple measures to reduce their number and optimise the
representation. The weight-$n$ MPL of indices $\\{a_{1},\ldots,a_{n}\\}$ and
argument $x$ is defined recursively as

$\displaystyle
G(a_{1},a_{2},\ldots,a_{n};x)\coloneqq\int_{0}^{x}\frac{\text{d}t}{t-a_{1}}\,G(a_{2},\ldots,a_{n};t)\,,\qquad
a_{n}\neq 0\,,$ (23)

starting with $G(;x)=1$. Trailing zeros, i.e., zeros in the right-most
indices, are allowed through the definition

$\displaystyle
G(\underbrace{0,\ldots,0}_{k};x)\coloneqq\frac{1}{k!}\log^{k}(x)\,.$ (24)

We refer to ref. Vollinga:2004sn  for a thorough discussion.

Since the letters in eq. 19 are rational and linear in all variables, we can
solve the canonical DEs in eq. 18 algorithmically in terms of MPLs. Order by
order in $\epsilon$, the solution is given by

$\displaystyle\vec{g}^{(w)}_{\tau}(\vec{s})=\sum_{i=1}^{7}\mathrm{A}_{i}^{(\tau)}\cdot\int_{\gamma}\text{d}\log\bigl{(}W_{i}(\vec{s}=\gamma)\bigr{)}\,\vec{g}^{(w-1)}_{\tau}\bigl{(}\vec{s}=\gamma\bigr{)}+\vec{b}^{(w)}_{\tau}\,,$
(25)

starting from the constant weight-$0$ boundary values $\vec{g}^{(0)}_{\tau}$
determined in the previous subsection. Here, $\gamma$ is a path connecting an
arbitrary base-point $\vec{s}_{\mathrm{base}}$ to the end-point $\vec{s}$. The
weight-$w$ constants $\vec{b}^{(w)}_{\tau}$ are given by the values of the
integrals at the base-point,
$\vec{b}^{(w)}_{\tau}=\vec{g}^{(w)}_{\tau}(\vec{s}_{\mathrm{base}})$. For
$\vec{s}_{\mathrm{base}}$ we may use the boundary point $\vec{s}_{0}$ in eq.
21, so that the constants $\vec{b}^{(w)}_{\tau}$ coincide with the boundary
values determined numerically in the previous section. We follow a different
approach, which allows us to trade all numerical constants in the expressions
for zeta values.

We find it convenient to change variables from $(s_{12},s_{23},s_{4})$ to
$(z_{1},z_{2},s_{4})$, with

$\displaystyle z_{1}=\frac{s_{12}}{s_{4}}\,,\qquad\qquad
z_{2}=\frac{s_{23}}{s_{4}}\,.$ (26)

This way, there is only one dimensionful variable, $s_{4}$, the dependence on
which is fixed as an overall factor by dimensional analysis. We then integrate
the canonical DEs as in eq. 25 along the following piece-wise path in the
$(z_{1},z_{2},s_{4})$ space:

$\displaystyle(0,0,0)\overset{\gamma_{1}}{\longrightarrow}(z_{1},0,0)\overset{\gamma_{2}}{\longrightarrow}(z_{1},z_{2},0)\overset{\gamma_{3}}{\longrightarrow}(z_{1},z_{2},s_{4})\,.$
(27)

Since the Feynman integrals are divergent at the chosen base-point, the latter
is understood in a _regularised_ sense (we refer to section 4 of ref.
Abreu:2022mfk  for a thorough discussion). Choosing $(0,0,0)$ as base-point
has the important benefit of removing spurious transcendental numbers that
would pollute the solution were we to choose a base-point where the integrals
are finite. As we will see below, only zeta values appear. Roughly speaking,
we define regularised, finite values
$\vec{b}^{(w)}_{\tau}\coloneqq\mathrm{Reg}\,\vec{g}^{(w)}_{\tau}(\vec{s}_{\mathrm{base}})$
by introducing a regulator and formally setting to $0$ the (divergent)
logarithms of the regulator. Since the integrals are finite at a generic end-
point $\vec{s}$, the divergences at the base-point must cancel out with
divergences arising in the integration. We can thus drop all these
divergences. Provided that we do it consistently between the integration and
the base-point values $\vec{b}^{(w)}_{\tau}$, this leads to a finite and
unique result. In practice, we fix the finite base-point values
$\vec{b}^{(w)}_{\tau}$ by matching the solution
$\vec{g}^{(w)}_{\tau}(\vec{s})$ evaluated at the boundary point $\vec{s}_{0}$
against the boundary values discussed in the previous subsection.

We therefore keep the $\vec{b}^{(w)}_{\tau}$ as symbols and integrate the
canonical DEs as in eq. 25 along the path in eq. 27 up to weight $4$. We
parameterise each piece of the path in eq. 27 linearly. For instance,
$\gamma_{2}(t)=(z_{1},t,0)$, with $t\in[0,z_{2}]$.

  * •

The $\gamma_{1}$ integration leads to MPLs with indices in $\\{0,1\\}$ and
argument $z_{1}$.

  * •

The $\gamma_{2}$ integration leads to MPLs with indices in
$\\{0,1,1-z_{1},-z_{1}\\}$ and argument $z_{2}$.

  * •

The $\gamma_{3}$ integration leads to powers of $\log(-s_{4})$, fixed by
dimensional analysis.

Once we have obtained expressions for all MIs in terms of MPLs and symbolic
constants $\vec{b}^{(w)}_{\tau}$, we equate them to the numerical boundary
values at $\vec{s}_{0}$, and solve for the $\vec{b}^{(w)}_{\tau}$. We use
GiNaC Bauer:2000cp ; Vollinga:2004sn  to evaluate the MPLs numerically.
Finally, we use the PSLQ algorithm PSLQ  to express the ensuing values of
$\vec{b}^{(w)}_{\tau}$ in terms of $\zeta_{2}$ and $\zeta_{3}$. As a result,
we obtain a fully analytic representation of all MIs — and thus of our special
function basis $\\{F^{(w)}_{i}\\}$ — in terms of MPLs and zeta values, up to
weight $4$.

Contrary to the functions in the basis $\\{F^{(w)}_{i}\\}$, the MPLs in their
representation satisfy functional relations. We make use of this freedom to
optimise our expressions in view of their numerical evaluation by reducing the
number of distinct MPLs that need to be evaluated. First, we use the shuffle
algebra of MPLs to push all trailing zeros into logarithms through eq. 24
Vollinga:2004sn . Next, we employ the scaling relation

$\displaystyle
G(a_{1},\ldots,a_{n};x)=G\left(\frac{a_{1}}{x},\ldots,\frac{a_{n}}{x};1\right)\,,$
(28)

which holds for $x,a_{n}\neq 0$. As a result, all MPLs have argument $1$ and
indices

$\displaystyle l_{0}=0\,,\qquad l_{1}=\frac{s_{4}}{s_{12}}\,,\qquad
l_{2}=\frac{s_{4}}{s_{23}}\,,\qquad l_{3}=\frac{s_{4}-s_{12}}{s_{23}}\,,\qquad
l_{4}=-\frac{s_{12}}{s_{23}}\,.$ (29)

Finally, we decompose the MPLs to _Lyndon words_ Radford1979ANR  using
PolyLogTools Duhr:2019tlz ; we refer to the latter work for a thorough
explanation, and give here only a simple example. This procedure requires that
we choose a symbolic ordering of the MPL indices. We choose $l_{0}\prec
l_{1}\prec l_{2}\prec l_{3}\prec l_{4}$, meaning that $l_{1}$ is greater than
$l_{0}$, and so on. Consider the MPL $G(l_{1},l_{0};1)$, whose indices are not
sorted according to the ordering above, since $l_{1}\succ l_{0}$. We can use
the shuffle algebra of MPLs to rewrite it in terms of MPLs whose indices are
sorted according to the chosen ordering, as

$\displaystyle G(l_{1},l_{0};1)=G(l_{0};1)\,G(l_{1};1)-G(l_{0},l_{1};1)\,.$
(30)

Doing this consistently throughout all expressions reduces the number of
higher-weight MPLs in favour of products of lower-weight ones, which are
cheaper to evaluate numerically. To maximise the impact in this sense, we
tested all possible orderings of the indices and selected the one — given
above — which minimises the number of weight-4 MPLs. The resulting
representation of the function basis contains $4$ weight-1, 6 weight-2, 19
weight-3, and 25 weight-4 MPLs, as well as $3$ logarithms:

$\displaystyle\log(s_{12}/s_{4})\,,\qquad\quad\log(s_{23}/s_{4})\,,\qquad\quad\log(-s_{4})\,.$
(31)

We write the latter in terms of logarithms rather than MPLs as they play an
important role in the factorisation of the IR divergences in the scattering
amplitudes (see appendix D for the IR structure of the amplitudes we compute
here). We stress that $\log(-s_{4})$ is the only function of a dimensionful
argument in our representation of the function basis.

We provide in the folder mi2func/ of our ancillary files zenodo  the
expression of the basis functions $\\{F^{(w)}_{i}\\}$ in terms of MPLs,
logarithms, $\zeta_{2}$ and $\zeta_{3}$.

It is important to stress that the MPLs are multi-valued functions. For unit
argument, there is a pole on the integration contour whenever one of the
indices lies between $0$ and $1$. In this case the contour must be deformed in
the complex plane, either above or below the pole, leading to different
branches. Our MPLs are thus well-defined only in the kinematic region where
all MPL indices in eq. 29 are either less than 0 or greater than 1, and we
need $s_{4}<0$ for the argument of all logarithms in eq. 31 to be positive. We
analytically continue the MPLs to the kinematic regions of interest by adding
infinitesimal imaginary parts to the indices in the numerical evaluation, as
we discuss in appendix E. The analytic continuation of these functions has
been analysed in great detail in ref. Gehrmann:2002zr , where it is achieved
by suitable changes of variables so that the imaginary parts are extracted
explicitly, and all MPLs are real and single valued. This approach could lead
to faster evaluations although it was unnecessary for the applications
considered here.

###  4.3 Performance and validation

We validated our results for the MIs of all families by crosschecking them
against values obtained with AMFlow Liu:2022chg  at several random points in
all the physical kinematic regions discussed in appendix E. Furthermore, we
find agreement with the results of ref. Gehrmann:2023etk . We employ GiNaC
Bauer:2000cp ; Vollinga:2004sn  to evaluate the MPLs.

Our results allow for an efficient and stable evaluation of the MIs, and are
thus ready for immediate deployment in phenomenology. Indeed, the amplitudes
we computed in this work have already been implemented in McMule
Banerjee:2020rww ; ulrich_yannick_2022_6046769  to provide the real-double-
virtual electron-line corrections to $e\mu\to e\mu$ scattering. The evaluation
is efficient, running at $\approx 130$ events per second in the bulk of the
phase space ulrich-radcor  using handyG Naterop:2019xaf  for the evaluation of
the MPLs.

##  5 Conclusions

In this article, we calculated analytically the two-loop QED helicity
amplitudes for the process $0\to\ell\bar{\ell}\gamma\gamma^{*}$ in terms of a
basis of multiple polylogarithms that are suitable for fast and stable
numerical evaluation. We employed modern finite-field evaluation techniques to
reconstruct the amplitudes directly in terms of the special function basis,
sidestepping the symbolic computation in all intermediate stages. As a by-
product we have recomputed all two-loop master integrals for four-point
functions with an off-shell leg up to transcendental weight four, and provide
all the necessary ingredients needed to use them in amplitude computations
with the same kinematics.

We hope these new results will now open the path to N3LO predictions that can
be used for the future MUonE experiment.

Note added: Another computation of these amplitudes Fadin:2023phc  became
available on the arXiv shortly after ours.

###### Acknowledgements.

We thank Yannick Ulrich for providing the one-loop crosschecks, correspondence
on the McMule implementation of these amplitudes, and other useful
discussions. We further thank Heribertus Bayu Hartanto and Tiziano Peraro for
collaboration on the codebase. SZ wishes to thank Dmitry Chicherin and Vasily
Sotnikov for many useful discussions. This project received funding from the
European Union’s Horizon 2020 research and innovation programme High precision
multi-jet dynamics at the LHC (grant agreement number 772099).

##  Appendix A Definition of the Feynman integral families

For each two-loop integral family $\tau$ corresponding to one of the maximal
topologies shown in figure 2, the Feynman integrals have the form

$j^{\tau}(a_{1},\ldots,a_{9})=\mathrm{e}^{2\epsilon\gamma_{E}}\int\frac{\mathrm{d}^{4-2\epsilon}k_{1}}{\mathrm{i}\pi^{2-\epsilon}}\frac{\mathrm{d}^{4-2\epsilon}k_{2}}{\mathrm{i}\pi^{2-\epsilon}}\frac{1}{D_{\tau,1}^{a_{1}}\ldots
D_{\tau,9}^{a_{9}}}\,.$ (32)

The sets $\\{D_{\tau,1},\ldots,D_{\tau,9}\\}$ contain seven (inverse)
propagators and two ISPs ($a_{8},a_{9}\leq 0$). For the maximal topologies
under consideration, they are given by111We use a naming convention analogous
to that of ref. Abreu:2020jxa .:

  * •

penta-triangle, mzz configuration:

$\displaystyle\begin{aligned}
\big{\\{}k_{1}^{2},(k_{1}+p_{1}+p_{2}+p_{3})^{2},(k_{1}+p_{2}+p_{3})^{2},(k_{1}+p_{3})^{2},k_{2}^{2},(k_{2}-p_{3})^{2},\\\
(k_{1}+k_{2})^{2},(k_{2}-p_{1}-p_{2}-p_{3})^{2},(k_{2}-p_{2}-p_{3})^{2}\big{\\}}\,,\end{aligned}$
(33)

  * •

penta-triangle, zmz configuration:

$\displaystyle\begin{aligned}
\big{\\{}k_{1}^{2},(k_{1}-p_{1})^{2},(k_{1}+p_{2}+p_{3})^{2},(k_{1}+p_{3})^{2},k_{2}^{2},(k_{2}-p_{3})^{2},(k_{1}+k_{2})^{2},\\\
(k_{2}+p_{1})^{2},(k_{2}-p_{2}-p_{3})^{2}\big{\\}}\,,\end{aligned}$ (34)

  * •

penta-triangle, zzz configuration:

$\displaystyle\begin{aligned}
\big{\\{}k_{1}^{2},(k_{1}-p_{1})^{2},(k_{1}-p_{1}-p_{2})^{2},(k_{1}-p_{1}-p_{2}-p_{3})^{2},k_{2}^{2},(k_{2}+p_{1}+p_{2}+p_{3})^{2},\\\
(k_{1}+k_{2})^{2},(k_{2}+p_{1})^{2},(k_{2}+p_{1}+p_{2})^{2}\big{\\}}\,,\end{aligned}$
(35)

  * •

planar double-box:

$\displaystyle\begin{aligned}
\big{\\{}k_{1}^{2},(k_{1}-p_{1})^{2},(k_{1}-p_{1}-p_{2})^{2},k_{2}^{2},(k_{2}+p_{1}+p_{2}+p_{3})^{2},(k_{2}+p_{1}+p_{2})^{2},\\\
(k_{1}+k_{2})^{2},(k_{1}-p_{1}-p_{2}-p_{3})^{2},(k_{2}+p_{1})^{2}\big{\\}}\,,\end{aligned}$
(36)

  * •

crossed double-box, mz configuration:

$\displaystyle\begin{aligned}
\big{\\{}k_{1}^{2},(k_{1}+p_{1}+p_{2}+p_{3})^{2},(k_{1}+p_{2}+p_{3})^{2},k_{2}^{2},(k_{2}-p_{2})^{2},(k_{1}+k_{2})^{2},\\\
(k_{1}+k_{2}+p_{3})^{2},(k_{1}+p_{3})^{2},(k_{2}-p_{1}-p_{2}-p_{3})^{2}\big{\\}}\,,\end{aligned}$
(37)

  * •

crossed double-box, zz configuration:

$\displaystyle\begin{aligned}
\big{\\{}k_{1}^{2},(k_{1}-p_{1})^{2},(k_{1}-p_{1}-p_{2})^{2},k_{2}^{2},(k_{2}-p_{3})^{2},(k_{1}+k_{2})^{2},\\\
(k_{1}+k_{2}-p_{1}-p_{2}-p_{3})^{2},(k_{1}-p_{1}-p_{2}-p_{3})^{2},(k_{2}+p_{1})^{2}\big{\\}}\,.\end{aligned}$
(38)

We also use the one-loop (one-mass) box family, made of the following
integrals:

$j^{\rm
box}(a_{1},a_{2},a_{3},a_{4})=\mathrm{e}^{\epsilon\gamma_{E}}\int\frac{\mathrm{d}^{4-2\epsilon}k}{\mathrm{i}\pi^{2-\epsilon}}\frac{1}{D_{\mathrm{box},1}^{a_{1}}D_{\mathrm{box},2}^{a_{2}}D_{\mathrm{box},3}^{a_{3}}D_{\mathrm{box},4}^{a_{4}}}\,,$
(39)

with the four inverse propagators $D_{\mathrm{box},i}$

$\big{\\{}k_{1}^{2},(k_{1}-p_{1})^{2},(k_{1}-p_{1}-p_{2})^{2},(k_{1}-p_{1}-p_{2}-p_{3})^{2}\big{\\}}\,.$
(40)

Feynman’s prescription for the imaginary parts of all propagators is implicit.

These family definitions (strictly with the ordering of inverse propagators
and ISPs shown above) correspond to the integrals j[family,$a_{1},$…] that
build the canonical MI bases provided in the pure_mi_bases/ directory of our
ancillary files zenodo . In this notation, each j[...] represents a Feynman
integral within a given integral family, while the numbers $a_{i}$ refer to
the powers of its propagators and ISPs.

##  Appendix B Optimised IBP reduction procedure for amplitudes with many
permuted integral families

An amplitude will in general have contributions from permutations of the
ordered integral families shown in figure 2. To reduce the tensor integrals in
the amplitude, IBP identities must be generated for all the permutations of
these ordered families. This can lead to a very large IBP system. The
performance of the reduction setup is extremely sensitive to the number of IBP
identities required so, to minimise the memory consumption, we choose to
generate IBP identities only for the ordered families. Next, we obtain the
reduction for any permutation of these families by permuting the ‘ordered’
reduction numerically over finite fields. The result is then given in terms of
MIs of each family permutation, but it is missing the symmetry relations that
can be found between subsectors of different families. To express the final
result in terms of a minimal set of MIs, we find such relations from a
separate computation. One may account for integral symmetries using automated
tools such as LiteRed Lee:2012cn . Since we use a pure basis of MIs, the
symmetry relations amongst them will have rational numbers as coefficients.
This is because the presence of any kinematic invariant would spoil the purity
of the canonical DEs (see section 4), and would mean that such a symmetry
relation in fact involves non-pure integrals. Therefore, the computation of
the missing symmetry relations can be performed with all kinematic invariants
set to numeric values, which significantly lowers the complexity of this task.
Finally, we note that even if symmetries amongst the MIs were missed, a
representation of the integrals in terms of a basis of special functions — as
we construct in section 4 — would automatically incorporate the extra
simplifications and so the same final result would be obtained. Nonetheless,
in practice we do find it useful to include these symmetry relations, as they
reduce the number of independent coefficients that have to be processed
further.

The procedure can be summarised as follows:

  1. 1.

Generate (analytic) IBPs for the six ordered families.

  2. 2.

Compute the mappings between permutations of the MIs of the system above.

  3. 3.

Take the tensor integrals in the amplitudes for each permutation of these
families and solve the linear system over finite fields.

  4. 4.

Apply the symmetry mappings between the MIs of each family permutation to find
the minimal set for the full system.

Since there are a few additional bits of terminology, we can consider a
concrete example to clarify everything. At one-loop, a four-point process with
a single off-shell leg can be described by a single independent integral
family which is simply the box topology (see appendix A for its explicit
definition). Following the Laporta reduction algorithm leads to a basis of
four MIs,

${\rm MI}^{\rm box}=\\{j^{\rm box}(1,1,1,1),\,j^{\rm box}(1,0,1,0),\,j^{\rm
box}(0,1,0,1),\,j^{\rm box}(1,0,0,1)\\}\,,$ (41)

which are the scalar box and scalar bubble integrals in channels
$s_{12},s_{23}$ and $s_{4}$ respectively. An amplitude will, in general, be
written in terms of three permutations of this family. Let us denote these
permutations as $j^{\rm box,1234}$, $j^{\rm box,2314}$, and $j^{\rm
box,3124}$, where $j^{\rm box,1234}=j^{\rm box}$ as above and the additional
superscript indices refer to the order of the external legs. Following our
procedure we would load one set of IBP relations generated for $j^{\rm box}$.
These identities can then be permuted numerically, for example as FiniteFlow
graphs, to reduce tensor integrals in each of the three permuted families. The
result is now in terms of twelve MIs: three boxes and nine bubbles. While the
amplitude is already in a minimal basis of box integrals, there is clearly an
over-complete set of bubbles. The independent bubbles are in the channels
$s_{12}$, $s_{23}$, $s_{13}$, and $s_{4}$, so the five additional symmetry
mappings are

$\displaystyle\begin{aligned} j^{\rm box,2314}(1,0,1,0)&=j^{\rm
box,1234}(0,1,0,1)\,,&j^{\rm box,3124}(1,0,1,0)&=j^{\rm
box,2314}(0,1,0,1)\,,\\\ j^{\rm box,3124}(0,1,0,1)&=j^{\rm
box,1234}(1,0,1,0)\,,&j^{\rm box,2314}(1,0,0,1)&=j^{\rm
box,1234}(1,0,0,1)\,,\\\ j^{\rm box,3124}(1,0,0,1)&=j^{\rm
box,1234}(1,0,0,1)\,.\end{aligned}$ (42)

After applying these identities we arrive at the final result with seven MIs
which cover all permutations of the integral families. This approach would not
lead to any significant performance enhancements in this simple example of
course, but it can be particularly important when considering high-
multiplicity examples where the number of permutations is high.

##  Appendix C Rational parametrisation of the kinematics

Since we are applying finite-field techniques to helicity amplitudes, we
employ a rational parametrisation of the external kinematics using Hodges’s
momentum twistor formalism Hodges:2009hk . While this is not essential to
combat the algebraic complexity for the kinematics considered here, it does
provide a convenient parametrisation of the spinor products.

The single-off-shell four-particle phase space $p$ is obtained from a massless
five-particle parametrisation $q$ (defined in appendix A of ref.
Badger:2021imn  with $\\{x_{2}\leftrightarrow x_{4},x_{3}\leftrightarrow
x_{5}\\}$) under

$\displaystyle p_{i}$ $\displaystyle=q_{i}\qquad\forall i=1,2,3,$
$\displaystyle p_{4}=q_{4}+q_{5}\,.$ (43)

The momentum twistor variables $x_{i}$ for $p$ are then related to the scalar
invariants $\vec{s}$ through

$\displaystyle s_{12}$ $\displaystyle=x_{1}\,,$ $\displaystyle s_{23}$
$\displaystyle=x_{1}x_{2}\,,$ $\displaystyle s_{4}$
$\displaystyle=x_{1}x_{3}\,.$ (44)

Momentum twistors allow us to express any spinor expression as a rational
function in the variables $x_{i}$. In this representation the helicity scaling
is however obscured, as we have fixed the spinor phases in order to achieve a
parameterisation in terms of the minimal number of variables (see e.g. ref.
Badger:2016uuq ). Therefore, we need to manually restore the phase information
at the end of the computation. This can be achieved by multiplying the
momentum twistor expression by an arbitrary factor $\Phi$ with the same
helicity scaling as the helicity amplitude under consideration, divided by
that factor written in terms of momentum twistor variables. For example, for
the helicity configurations of eq. 6, we can use the phase factors

$\displaystyle\Phi(-++)$ $\displaystyle=\frac{\langle 12\rangle}{\langle
23\rangle^{2}}\,,$ $\displaystyle\Phi(-+-)$
$\displaystyle=\frac{[12]}{[13]^{2}}\,,$ (45)

which in our momentum twistor parameterisation are given by

$\displaystyle\Phi(-++)$ $\displaystyle=x_{1}^{2}\,,$ $\displaystyle\Phi(-+-)$
$\displaystyle=-\frac{1}{x_{1}(1+x_{2}-x_{3})^{2}}\,.$ (46)

We refer to appendix C of ref. Badger:2023mgf  for a thorough discussion of
how to restore the phase information in a momentum twistor parameterisation.

##  Appendix D Renormalisation and infrared structure

We renormalise the coupling constant by trading the bare coupling
$\alpha_{\mathrm{bare}}$ for the renormalised one $\alpha_{\mathrm{R}}$
through

$\displaystyle\alpha_{\mathrm{bare}}=\alpha_{\mathrm{R}}(\mu_{R})\,Z_{\alpha}\bigl{(}\alpha_{\mathrm{R}}(\mu_{R})\bigr{)}\,\mu_{R}^{2\epsilon}\,S_{\epsilon}\,,$
(47)

with $S_{\epsilon}=(4\pi)^{-\epsilon}\mathrm{e}^{\epsilon\gamma_{E}}$. The
renormalisation factor $Z_{\alpha}$ in the $\overline{\text{MS}}$ scheme is
Barnreuther:2013qvf ; Bonciani:2021okt

$\displaystyle
Z_{\alpha}(\alpha)=1-\frac{\alpha}{4\pi}\frac{\beta_{0}}{\epsilon}-\left(\frac{\alpha}{4\pi}\right)^{2}\left(-\frac{\beta_{0}^{2}}{\epsilon^{2}}+\frac{1}{2}\frac{\beta_{1}}{\epsilon}\right)+\mathcal{O}\bigl{(}\alpha^{3}\bigr{)}\,.$
(48)

The $\beta$-function is defined from the renormalised coupling as

$\displaystyle\frac{\text{d}\alpha_{\rm
R}(\mu_{R})}{\text{d}\ln\mu_{R}}=\left[-2\,\epsilon+\beta\bigl{(}\alpha_{\rm
R}(\mu_{R})\bigr{)}\right]\alpha_{\rm R}(\mu_{R})\,,$ (49)

and expanded as

$\displaystyle\beta(\alpha)=-2\,\frac{\alpha}{4\pi}\sum_{k\geq
0}\beta_{k}\left(\frac{\alpha}{4\pi}\right)^{k}\,,$ (50)

with

$\displaystyle\beta_{0}=-\frac{4}{3}n_{l}\,,\qquad\qquad\beta_{1}=-4\,n_{l}\,.$
(51)

The photon wavefunction renormalisation factor is $Z_{A}=Z_{\alpha}$, which we
include due to the external off-shell photon. The complete renormalisation
procedure then is

$\displaystyle\mathcal{A}^{\mu}_{\text{renorm}}(\alpha_{R})=Z_{A}^{\frac{1}{2}}(\alpha_{R})\,\mathcal{A}^{\mu}_{\text{bare}}(\alpha_{\text{bare}})\,,$
(52)

where $\alpha_{\text{bare}}$ is expressed in terms of $\alpha_{R}$ through eq.
47.

The IR poles of the renormalised amplitude factorise as Catani:1998bh ;
Gardi:2009qi ; Gardi:2009zv ; Becher:2009cu ; Becher:2009qa

$\displaystyle\mathcal{A}^{\mu}_{\text{renorm}}(\alpha_{R})=Z(\alpha_{R})\,\mathcal{F}^{\mu}(\alpha_{R})\,,$
(53)

so that $Z(\alpha_{R})$ captures all IR poles and $\mathcal{F}^{\mu}$ is a
finite remainder. We obtain the explicit two-loop expression of the IR factor
$Z(\alpha_{R})$ by choosing QED parameters ($C_{A}=0$, $C_{F}=1$, and
$T_{F}=1$) in the non-abelian gauge-theory expressions of ref. Becher:2009qa .
We expand it as

$\displaystyle Z(\alpha)$ $\displaystyle=\sum_{k\geq
0}Z^{(L)}\left(\frac{\alpha}{4\pi}\right)^{L}\,.$ (54)

The coefficients $Z^{(L)}$ are expressed in terms of the anomalous dimension

$\displaystyle\Gamma$
$\displaystyle=\gamma^{\text{cusp}}\ln\left(\frac{-s_{12}}{\mu^{2}}\right)+2\gamma^{l}+\gamma^{A}\,,$
(55)

and its derivative

$\displaystyle\Gamma^{\prime}$
$\displaystyle\coloneqq\frac{\partial\Gamma}{\partial\ln\mu}=-2\gamma^{\text{cusp}}\,.$
(56)

Here, $\gamma^{\text{cusp}}$ is the cusp anomalous dimension, while
$\gamma^{l}$ and $\gamma^{A}$ are the lepton’s and the photon’s collinear
anomalous dimensions, respectively. We expand all anomalous dimensions
$y\in\\{\Gamma,\gamma^{i}\\}$ as

$\displaystyle y$ $\displaystyle=\frac{\alpha}{4\pi}\sum_{k\geq
0}y_{k}\left(\frac{\alpha}{4\pi}\right)^{k}\,,$ (57)

with coefficients

$\displaystyle\gamma^{l}_{0}$ $\displaystyle=-3\,,$
$\displaystyle\gamma^{l}_{1}$
$\displaystyle=-\frac{3}{2}+2\pi^{2}-24\,\zeta_{3}+n_{l}\left(\frac{130}{27}+\frac{2}{3}\pi^{2}\right)\,,$
(58a) $\displaystyle\gamma^{A}_{0}$ $\displaystyle=-\beta_{0}\,,$
$\displaystyle\gamma^{A}_{1}$ $\displaystyle=-\beta_{1}\,,$ (58b)
$\displaystyle\gamma^{\text{cusp}}_{0}$ $\displaystyle=4\,,$
$\displaystyle\gamma^{\text{cusp}}_{1}$ $\displaystyle=-\frac{80}{9}n_{l}\,.$
(58c)

Finally, the coefficients of the IR factor $Z$ up to two loop are given by

$\displaystyle Z^{(0)}$ $\displaystyle=1\,,$ $\displaystyle Z^{(1)}$
$\displaystyle=\frac{\Gamma_{0}^{\prime}}{4\epsilon^{2}}+\frac{\Gamma_{0}}{2\epsilon}\,,$
$\displaystyle Z^{(2)}$
$\displaystyle=\frac{{Z^{(1)}}^{2}}{2}-\frac{3\beta_{0}\Gamma_{0}^{\prime}}{16\epsilon^{3}}+\frac{\Gamma_{1}^{\prime}-4\beta_{0}\Gamma_{0}}{16\epsilon^{2}}+\frac{\Gamma_{1}}{4\epsilon}\,.$
(59)

Putting together the subtraction of UV and IR poles, and expanding the
resulting finite remainder $\mathcal{F}^{\mu}(\alpha_{R})$ in $\alpha_{R}$
leads to the definitions in eq. 14.

##  Appendix E Analytic continuation

We analytically continue the MPLs by adding an infinitesimal positive (or
negative) imaginary part to the MPL indices $l_{i}$ in eq. 29 whenever they
fall between $0$ and $1$. The imaginary part of each index prescribes how to
deform the integration contour around the pole associated with it. We do
similarly for the logarithms in eq. 31. To this end, following ref.
Gehrmann:2002zr , we change variables from $(s_{12},s_{23},s_{4})$ to
$(s_{12},s_{23},s_{13})$, with $s_{4}=s_{12}+s_{23}+s_{13}$. We then add a
small positive imaginary part to the latter variables, as

$\displaystyle s_{12}\longrightarrow s_{12}+\mathrm{i}\,c_{1}\,\delta\,,\qquad
s_{23}\longrightarrow s_{23}+\mathrm{i}\,c_{2}\,\delta\,,\qquad
s_{13}\longrightarrow s_{13}+\mathrm{i}\,c_{3}\,\delta\,,$ (60)

where $c_{1}$, $c_{2}$ and $c_{3}$ are arbitrary positive constants, and
$\delta$ is a positive infinitesimal. Finally, we check whether this
substitution gives a positive or negative imaginary part to each MPL index
$l_{i}$. This depends on the domain of the kinematic variables. We focus on
three kinematic regions which are of phenomenological interest. The analytic
continuation for any other region may be obtained similarly.

#### Electron-line corrections to $e^{-}\mu^{-}\to e^{-}\mu^{-}\gamma$.

To define the domain of the kinematic variables relevant for this application,
we embed the four-particle off-shell process of eq. 1 in the five-particle
process $e^{-}\mu^{-}\to e^{-}\mu^{-}\gamma$. We then determine the kinematic
constraints for the five-particle process (see e.g. appendix A of ref.
Chicherin:2021dyp ), and from them derive the constraints on the four-point
off-shell kinematics. The result is

$\displaystyle\mathcal{P}_{e\mu\to{e}\mu\gamma}\coloneqq\\{\vec{s}\colon
s_{12}<0\,\land\,s_{23}<0\,\land\,0<s_{13}<-s_{12}-s_{23}\\}\,.$ (61)

The MPL index $l_{4}=-s_{12}/s_{23}$ is always negative in
$\mathcal{P}_{{e}\mu\to{e}\mu\gamma}$, hence no analytic continuation is
required. The other three indices may instead fall between $0$ and $1$. Let us
study $l_{1}$. Changing variables from $s_{4}$ to $s_{13}$ and adding
imaginary parts as in eq. 60 gives

$\displaystyle
l_{1}=\frac{s_{12}+s_{13}+s_{23}}{s_{12}}+\frac{\mathrm{i}\delta}{s_{12}^{2}}\left[(c_{2}+c_{3})s_{12}-c_{1}(s_{13}+s_{23})\right]+\mathcal{O}\left(\delta^{2}\right)\,.$
(62)

The imaginary part of $l_{1}$ may be either negative or positive in
$\mathcal{P}_{{e}\mu\to{e}\mu\gamma}$. However, it is strictly negative in the
subregion of $\mathcal{P}_{{e}\mu\to{e}\mu\gamma}$ where $0<l_{1}<1$. We
therefore assign a negative imaginary part to $l_{1}$ whenever $0<l_{1}<1$ in
$\mathcal{P}_{{e}\mu\to{e}\mu\gamma}$. The analysis of the other indices
follows similarly, and is summarised in table 2. The arguments of the three
logarithms in eq. 31 are positive in $\mathcal{P}_{e\mu\to{e}\mu\gamma}$.

Index | $\mathcal{P}_{{e}\mu\to{e}\mu\gamma}$ | $\mathcal{P}_{{e}\bar{{e}}\to\gamma\gamma^{*}}$ | $\mathcal{P}_{\gamma^{*}\to e\bar{e}\gamma}$
---|---|---|---
$l_{1}$ | $-$ | $+$ | $0$
$l_{2}$ | $-$ | $0$ | $0$
$l_{3}$ | $-$ | $0$ | $0$
$l_{4}$ | $0$ | $0$ | $0$
Table 2: Imaginary parts of the MPL indices defined by eq. 29 in the three
kinematic regions discussed in appendix E. The symbol $+$ ($-$) denotes a
positive (negative) imaginary part, while $0$ means no analytic continuation
is needed.

#### Corrections to $e^{-}e^{+}\to\gamma\gamma^{*}$.

The relevant domain of the kinematic variables in this case can be derived
directly for the four-point kinematics, and is typically named the $s_{12}$
channel. It is given by

$\displaystyle\mathcal{P}_{e\bar{e}\to\gamma\gamma^{*}}\coloneqq\\{\vec{s}\colon
s_{23}<0\,\land\,s_{13}<0\,\land\,s_{12}>-s_{23}-s_{13}\\}\,.$ (63)

The MPL indices $l_{2}$, $l_{3}$ and $l_{4}$ can never fall between $0$ and
$1$ in $\mathcal{P}_{e\bar{e}\to\gamma\gamma^{*}}$, and hence require no
analytic continuation. We instead need to add a positive imaginary part to
$l_{1}$. In this region also the logarithms in eq. 31 need to be analytically
continued. The argument of $\log(s_{12}/s_{4})$ is positive in
$\mathcal{P}_{e\bar{e}\to\gamma\gamma^{*}}$. By adding imaginary parts to the
arguments of the other logarithms and studying them where the arguments are
negative in $\mathcal{P}_{e\bar{e}\to\gamma\gamma^{*}}$, we determine that the
analytic continuation is achieved through the following replacements:

$\displaystyle\log(s_{23}/s_{4})\longrightarrow\log(-s_{23}/s_{4})+\mathrm{i}\pi\,,\qquad\log(-s_{4})\longrightarrow\log(s_{4})-\mathrm{i}\pi\,.$
(64)

#### Corrections to the decay $\gamma^{*}\to e^{-}e^{+}\gamma$.

The relevant domain of the kinematic variables is

$\displaystyle\mathcal{P}_{\gamma^{*}\to
e\bar{e}\gamma}\coloneqq\\{\vec{s}\colon
s_{12}>0\,\land\,s_{23}>0\,\land\,s_{13}>0\\}\,.$ (65)

All MPL indices $l_{i}$ in eq. 29 are either $l_{i}<0$ or $l_{i}>1$, hence no
analytic continuation is required. The same holds for the first two logarithms
in eq. 31, whose arguments are positive. The only function which needs to be
analytically continued is $\log(-s_{4})$. We achieve this by replacing

$\displaystyle\log(-s_{4})\longrightarrow\log(s_{4})-\mathrm{i}\pi\,.$ (66)

The information about the imaginary parts of the MPL indices can be fed into
the publicly available libraries for evaluating these functions numerically,
such as FastGPL Wang:2021imw , GiNaC Bauer:2000cp ; Vollinga:2004sn , and
handyG Naterop:2019xaf . This typically leads to longer evaluation times with
respect to MPLs which do not need analytic continuation. We find that this is
not an issue for the planned applications of our results (see section 4.3).
Nonetheless, we note that a more performant evaluation may be achieved by
tailoring the representation to the kinematic region of interest in such a way
that no MPLs require analytic continuation. We refer to ref. gehrmann:2000zt ;
gehrmann:2001ck ; Gehrmann:2002zr ; Duhr:2012fh ; Gehrmann:2023etk  for a
detailed discussion.

## References

  * (1) G. Abbiendi et al., _Measuring the leading hadronic contribution to the muon g-2 via $\mu e$ scattering_, _Eur. Phys. J. C_ 77 (2017) 139 [1609.08987].
  * (2) G. Abbiendi, _Status of the MUonE experiment_ , _Phys. Scripta_ 97 (2022) 054007 [2201.13177].
  * (3) MUonE collaboration, _Status of the MUonE experiment_ , _PoS_ EPS-HEP2021 (2022) 642.
  * (4) R.N. Pilato, _The MUonE experiment_ , _Nuovo Cim. C_ 45 (2022) 67.
  * (5) C.M. Carloni Calame, M. Passera, L. Trentadue and G. Venanzoni, _A new approach to evaluate the leading hadronic corrections to the muon $g$-2_, _Phys. Lett. B_ 746 (2015) 325 [1504.02228].
  * (6) E. Balzani, S. Laporta and M. Passera, _Hadronic vacuum polarization contributions to the muon g-2 in the space-like region_ , _Phys. Lett. B_ 834 (2022) 137462 [2112.05704].
  * (7) Muon g-2 collaboration, _Measurement of the Positive Muon Anomalous Magnetic Moment to 0.46 ppm_ , _Phys. Rev. Lett._ 126 (2021) 141801 [2104.03281].
  * (8) T. Aoyama et al., _The anomalous magnetic moment of the muon in the Standard Model_ , _Phys. Rept._ 887 (2020) 1 [2006.04822].
  * (9) S. Borsanyi et al., _Leading hadronic contribution to the muon magnetic moment from lattice QCD_ , _Nature_ 593 (2021) 51 [2002.12347].
  * (10) P. Banerjee et al., _Theory for muon-electron scattering @ 10 ppm: A report of the MUonE theory initiative_ , _Eur. Phys. J. C_ 80 (2020) 591 [2004.13663].
  * (11) E. Budassi and C.L. Del Pio, _High precision calculations for the MUonE experiment_ , _PoS_ ICHEP2022 (2022) 796.
  * (12) M. Alacevich, C.M. Carloni Calame, M. Chiesa, G. Montagna, O. Nicrosini and F. Piccinini, _Muon-electron scattering at NLO_ , _JHEP_ 02 (2019) 155 [1811.06743].
  * (13) E. Budassi, C.M. Carloni Calame, M. Chiesa, C.L. Del Pio, S.M. Hasan, G. Montagna et al., _NNLO virtual and real leptonic corrections to muon-electron scattering_ , _JHEP_ 11 (2021) 098 [2109.14606].
  * (14) M. Fael and M. Passera, _Muon-Electron Scattering at Next-To-Next-To-Leading Order: The Hadronic Corrections_ , _Phys. Rev. Lett._ 122 (2019) 192001 [1901.03106].
  * (15) M. Fael, _Hadronic corrections to $\mu$-$e$ scattering at NNLO with space-like data_, _JHEP_ 02 (2019) 027 [1808.08233].
  * (16) C.M. Carloni Calame, M. Chiesa, S.M. Hasan, G. Montagna, O. Nicrosini and F. Piccinini, _Towards muon-electron scattering at NNLO_ , _JHEP_ 11 (2020) 028 [2007.01586].
  * (17) A. Broggio et al., _Muon-electron scattering at NNLO_ , _JHEP_ 01 (2023) 112 [2212.06481].
  * (18) Y. Ulrich, “N3LO kick-off workstop/thinkstart.” https://conference.ippp.dur.ac.uk/event/1104/, Durham, Aug, 2022.
  * (19) 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].
  * (20) M. Fael, F. Lange, K. Schönwald and M. Steinhauser, _Singlet and nonsinglet three-loop massive form factors_ , _Phys. Rev. D_ 106 (2022) 034029 [2207.00027].
  * (21) M. Fael, F. Lange, K. Schönwald and M. Steinhauser, _Massive three-loop form factors: Anomaly contribution_ , _Phys. Rev. D_ 107 (2023) 094017 [2302.00693].
  * (22) L.W. Garland, T. Gehrmann, E.W.N. Glover, A. Koukoutsakis and E. Remiddi, _The Two loop QCD matrix element for $e^{+}e^{-}\to 3$ jets_, _Nucl. Phys. B_ 627 (2002) 107 [hep-ph/0112081].
  * (23) L.W. Garland, T. Gehrmann, E.W.N. Glover, A. Koukoutsakis and E. Remiddi, _Two loop QCD helicity amplitudes for $e^{+}e^{-}\to$ three jets_, _Nucl. Phys. B_ 642 (2002) 227 [hep-ph/0206067].
  * (24) T. Gehrmann and L. Tancredi, _Two-loop QCD helicity amplitudes for $q\bar{q}\to W^{\pm}\gamma$ and $q\bar{q}\to Z^{0}\gamma$_, _JHEP_ 02 (2012) 004 [1112.1531].
  * (25) “PrecisionSM database.” https://precision-sm.github.io/.
  * (26) G. Abbiendi et al., _Mini-Proceedings of the STRONG2020 Virtual Workshop on ”Space-like and Time-like determination of the Hadronic Leading Order contribution to the Muon $g-2$”_, in _STRONG2020 Virtual Workshop “Space-like and Time-like determination of the Hadronic Leading Order contribution to the Muon g $-$ 2”_, 1, 2022 [2201.12102].
  * (27) A.A. Penin, _Two-loop photonic corrections to massive Bhabha scattering_ , _Nucl. Phys. B_ 734 (2006) 185 [hep-ph/0508127].
  * (28) T. Becher and K. Melnikov, _Two-loop QED corrections to Bhabha scattering_ , _JHEP_ 06 (2007) 084 [0704.3582].
  * (29) T. Engel, C. Gnendiger, A. Signer and Y. Ulrich, _Small-mass effects in heavy-to-light form factors_ , _JHEP_ 02 (2019) 118 [1811.06461].
  * (30) P. Banerjee, T. Engel, N. Schalch, A. Signer and Y. Ulrich, _Bhabha scattering at NNLO with next-to-soft stabilisation_ , _Phys. Lett. B_ 820 (2021) 136547 [2106.07469].
  * (31) T. Engel, A. Signer and Y. Ulrich, _Universal structure of radiative QED amplitudes at one loop_ , _JHEP_ 04 (2022) 097 [2112.07570].
  * (32) T. Engel, _The LBK theorem to all orders_ , 2304.11689.
  * (33) T. Engel, A. Signer and Y. Ulrich, _A subtraction scheme for massive QED_ , _JHEP_ 01 (2020) 085 [1909.10244].
  * (34) 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].
  * (35) 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].
  * (36) 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].
  * (37) 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].
  * (38) S. Abreu, F. Febres 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].
  * (39) 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].
  * (40) 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].
  * (41) 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].
  * (42) 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_ , 2304.06682.
  * (43) 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_ , 2305.17056.
  * (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) 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].
  * (46) 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].
  * (47) 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].
  * (48) A. von Manteuffel and R.M. Schabinger, _A novel approach to integration by parts reduction_ , _Phys. Lett. B_ 744 (2015) 101 [1406.4513].
  * (49) J. Klappert, F. Lange, P. Maierhöfer and J. Usovitsch, _Integral reduction with Kira 2.0 and finite field methods_ , _Comput. Phys. Commun._ 266 (2021) 108024 [2008.06494].
  * (50) A.V. Smirnov and F.S. Chuharev, _FIRE6: Feynman Integral REduction with Modular Arithmetic_ , _Comput. Phys. Commun._ 247 (2020) 106877 [1901.07808].
  * (51) T. Peraro, _Scattering amplitudes over finite fields and multivariate functional reconstruction_ , _JHEP_ 12 (2016) 030 [1608.01902].
  * (52) T. Peraro, _FiniteFlow: multivariate functional reconstruction using finite fields and dataflow graphs_ , _JHEP_ 07 (2019) 031 [1905.08019].
  * (53) T. Gehrmann and E. Remiddi, _Two loop master integrals for $\gamma^{*}\to 3$ jets: The Planar topologies_, _Nucl. Phys. B_ 601 (2001) 248 [hep-ph/0008287].
  * (54) T. Gehrmann and E. Remiddi, _Two loop master integrals for $\gamma^{*}\to 3$ jets: The Nonplanar topologies_, _Nucl. Phys. B_ 601 (2001) 287 [hep-ph/0101124].
  * (55) T. Gehrmann and E. Remiddi, _Analytic continuation of massless two loop four point functions_ , _Nucl. Phys. B_ 640 (2002) 379 [hep-ph/0207020].
  * (56) C. Duhr, _Hopf algebras, coproducts and symbols: an application to Higgs boson amplitudes_ , _JHEP_ 08 (2012) 043 [1203.0454].
  * (57) T. Gehrmann, P. Jakubčík, C.C. Mella, N. Syrrakos and L. Tancredi, _Two-loop helicity amplitudes for $H+$jet production to higher orders in the dimensional regulator_, _JHEP_ 04 (2023) 016 [2301.10849].
  * (58) T. Gehrmann, T. Peraro and L. Tancredi, _Two-loop QCD corrections to the V → $q\overline{q}g$ helicity amplitudes with axial-vector couplings_, _JHEP_ 02 (2023) 041 [2211.13596].
  * (59) T. Gehrmann, P. Jakubčík, C.C. Mella, N. Syrrakos and L. Tancredi, _Two-loop helicity amplitudes for $V+$jet production including axial vector couplings to higher orders in $\epsilon$_, 2306.10170.
  * (60) S. Di Vita, P. Mastrolia, U. Schubert and V. Yundin, _Three-loop master integrals for ladder-box diagrams with one massive leg_ , _JHEP_ 09 (2014) 148 [1408.3107].
  * (61) D.D. Canko and N. Syrrakos, _Resummation methods for Master Integrals_ , _JHEP_ 02 (2021) 080 [2010.06947].
  * (62) D.D. Canko and N. Syrrakos, _Planar three-loop master integrals for 2 → 2 processes with one external massive particle_ , _JHEP_ 04 (2022) 134 [2112.14275].
  * (63) J.M. Henn, J. Lim and W.J. Torres Bobadilla, _First look at the evaluation of three-loop non-planar Feynman diagrams for Higgs plus jet production_ , _JHEP_ 05 (2023) 026 [2302.12776].
  * (64) T. Gehrmann, J.M. Henn and N.A. Lo Presti, _Pentagon functions for massless planar scattering amplitudes_ , _JHEP_ 10 (2018) 103 [1807.09812].
  * (65) D. Chicherin and V. Sotnikov, _Pentagon Functions for Scattering of Five Massless Particles_ , _JHEP_ 20 (2020) 167 [2009.07803].
  * (66) D. Chicherin, V. Sotnikov and S. Zoia, _Pentagon functions for one-mass planar scattering amplitudes_ , _JHEP_ 01 (2022) 096 [2110.10111].
  * (67) 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.
  * (68) C. Gnendiger et al., _To ${d}$, or not to ${d}$: recent developments and comparisons of regularization schemes_, _Eur. Phys. J. C_ 77 (2017) 471 [1705.01827].
  * (69) M.L. Mangano and S.J. Parke, _Multiparton amplitudes in gauge theories_ , _Phys. Rept._ 200 (1991) 301 [hep-th/0509223].
  * (70) L.J. Dixon, _Calculating scattering amplitudes efficiently_ , in _Theoretical Advanced Study Institute in Elementary Particle Physics (TASI 95): QCD and Beyond_ , pp. 539–584, 1, 1996 [hep-ph/9601359].
  * (71) S. Badger, J. Henn, J. Plefka and S. Zoia, _Scattering Amplitudes in Quantum Field Theory_ , 2306.05976.
  * (72) M. Heller and A. von Manteuffel, _MultivariateApart: Generalized partial fractions_ , _Comput. Phys. Commun._ 271 (2022) 108174 [2101.08283].
  * (73) S. Badger, J. Kryś, R. Moodie and S. Zoia, _Ancillary files for “Lepton-pair scattering with an off-shell and an on-shell photon at two loops in massless QED”_ , July, 2023.  10.5281/zenodo.8120540.
  * (74) P. Banerjee, T. Engel, A. Signer and Y. Ulrich, _QED at NNLO with McMule_ , _SciPost Phys._ 9 (2020) 027 [2007.01654].
  * (75) Y. Ulrich, P. Banerjee, A. Coutinho, T. Engel, A. Gurgone, F. Hagelstein et al., _The mcmule manual_ , Feb., 2022.  10.5281/zenodo.6046769.
  * (76) S. Catani, _The Singular behavior of QCD amplitudes at two loop order_ , _Phys. Lett. B_ 427 (1998) 161 [hep-ph/9802439].
  * (77) E. Gardi and L. Magnea, _Factorization constraints for soft anomalous dimensions in QCD scattering amplitudes_ , _JHEP_ 03 (2009) 079 [0901.1091].
  * (78) E. Gardi and L. Magnea, _Infrared singularities in QCD amplitudes_ , _Nuovo Cim. C_ 32N5-6 (2009) 137 [0908.3273].
  * (79) T. Becher and M. Neubert, _Infrared singularities of scattering amplitudes in perturbative QCD_ , _Phys. Rev. Lett._ 102 (2009) 162001 [0901.0722].
  * (80) T. Becher and M. Neubert, _On the Structure of Infrared Singularities of Gauge-Theory Amplitudes_ , _JHEP_ 06 (2009) 081 [0903.1126].
  * (81) P. Nogueira, _Automatic Feynman graph generation_ , _J. Comput. Phys._ 105 (1993) 279.
  * (82) J. Kuipers, T. Ueda, J.A.M. Vermaseren and J. Vollinga, _FORM version 4.0_ , _Comput. Phys. Commun._ 184 (2013) 1453 [1203.6543].
  * (83) B. Ruijl, T. Ueda and J. Vermaseren, _FORM version 4.2_ , 1707.06453.
  * (84) F.V. Tkachov, _A Theorem on Analytical Calculability of Four Loop Renormalization Group Functions_ , _Phys. Lett. B_ 100 (1981) 65.
  * (85) K.G. Chetyrkin and F.V. Tkachov, _Integration by Parts: The Algorithm to Calculate beta Functions in 4 Loops_ , _Nucl. Phys._ B192 (1981) 159.
  * (86) R.N. Lee, _Presenting LiteRed: a tool for the Loop InTEgrals REDuction_ , 1212.2685.
  * (87) S. Laporta, _High precision calculation of multiloop Feynman integrals by difference equations_ , _Int. J. Mod. Phys._ A15 (2000) 5087 [hep-ph/0102033].
  * (88) J.M. Henn, _Multiloop integrals in dimensional regularization made simple_ , _Phys. Rev. Lett._ 110 (2013) 251601 [1304.1806].
  * (89) G. De Laurentis and B. Page, _Ansätze for scattering amplitudes from p-adic numbers and algebraic geometry_ , _JHEP_ 12 (2022) 140 [2203.04269].
  * (90) X. Liu, _Reconstruction of rational functions made simple_ , 2306.12262.
  * (91) T. Gehrmann, L. Tancredi and E. Weihs, _Two-loop QCD helicity amplitudes for $g\,g\to Z\,g$ and $g\,g\to Z\,\gamma$_, _JHEP_ 04 (2013) 101 [1302.2630].
  * (92) G. Barucchi and G. Ponzano, _Differential equations for one-loop generalized feynman integrals_ , _J. Math. Phys._ 14 (1973) 396.
  * (93) A. Kotikov, _Differential equations method. new technique for massive feynman diagram calculation_ , _Physics Letters B_ 254 (1991) 158.
  * (94) A. Kotikov, _Differential equation method. the calculation of n-point feynman diagrams_ , _Physics Letters B_ 267 (1991) 123.
  * (95) T. Gehrmann and E. Remiddi, _Differential equations for two loop four point functions_ , _Nucl. Phys. B_ 580 (2000) 485 [hep-ph/9912329].
  * (96) Z. Bern, L.J. Dixon and D.A. Kosower, _Dimensionally regulated pentagon integrals_ , _Nucl. Phys. B_ 412 (1994) 751 [hep-ph/9306240].
  * (97) J. Henn, B. Mistlberger, V.A. Smirnov and P. Wasser, _Constructing d-log integrands and computing master integrals for three-loop four-particle scattering_ , _JHEP_ 04 (2020) 167 [2002.09492].
  * (98) T. Gehrmann, J.M. Henn and N.A. 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].
  * (99) J.M. Henn, K. Melnikov and V.A. Smirnov, _Two-loop planar master integrals for the production of off-shell vector bosons in hadron collisions_ , _JHEP_ 05 (2014) 090 [1402.7078].
  * (100) F. Caola, J.M. Henn, K. Melnikov and V.A. Smirnov, _Non-planar master integrals for the production of two off-shell vector bosons in collisions of massless partons_ , _JHEP_ 09 (2014) 043 [1404.5590].
  * (101) C. Dlapa, _Algorithms and techniques for finding canonical differential equations of Feynman integrals_ , Ph.D. thesis, Munich U., 2022.  10.5282/edoc.29769.
  * (102) D. Gaiotto, J. Maldacena, A. Sever and P. Vieira, _Pulling the straps of polygons_ , _JHEP_ 12 (2011) 011 [1102.0062].
  * (103) X. Liu and Y.-Q. Ma, _AMFlow: A Mathematica package for Feynman integrals computation via auxiliary mass flow_ , _Comput. Phys. Commun._ 283 (2023) 108565 [2201.11669].
  * (104) K.-T. Chen, _Iterated path integrals_ , _Bull. Am. Math. Soc._ 83 (1977) 831.
  * (105) J. Vollinga and S. Weinzierl, _Numerical evaluation of multiple polylogarithms_ , _Comput. Phys. Commun._ 167 (2005) 177 [hep-ph/0410259].
  * (106) S. Abreu, R. Britto and C. Duhr, _The SAGEX review on scattering amplitudes Chapter 3: Mathematical structures in Feynman integrals_ , _J. Phys. A_ 55 (2022) 443004 [2203.13014].
  * (107) C.W. Bauer, A. Frink and R. Kreckel, _Introduction to the GiNaC framework for symbolic computation within the C++ programming language_ , _J. Symb. Comput._ 33 (2002) 1 [cs/0004015].
  * (108) H.R.P. Ferguson and D.H. Bailey, _A Polynomial Time, Numerically Stable Integer Relation Algorithm_ , _RNR Technical Report RNR-91-032_ (1992) .
  * (109) D.E. Radford, _A natural ring basis for the shuffle algebra and an application to group schemes_ , _Journal of Algebra_ 58 (1979) 432.
  * (110) C. Duhr and F. Dulat, _PolyLogTools — polylogs for the masses_ , _JHEP_ 08 (2019) 135 [1904.07279].
  * (111) Y. Ulrich, _QED at NNLO and beyond for precision experiments_ , _RADCOR 2023 – 16th International Symposium on Radiative Corrections: Applications of Quantum Field Theory to Phenomenology_ (2023) .
  * (112) L. Naterop, A. Signer and Y. Ulrich, _handyG —Rapid numerical evaluation of generalised polylogarithms in Fortran_ , _Comput. Phys. Commun._ 253 (2020) 107165 [1909.01656].
  * (113) V.S. Fadin and R.N. Lee, _Two-loop radiative corrections to $e^{+}e^{-}\rightarrow\gamma\gamma^{*}$ cross section_, 2308.09479.
  * (114) 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].
  * (115) A. Hodges, _Eliminating spurious poles from gauge-theoretic amplitudes_ , _JHEP_ 05 (2013) 135 [0905.1473].
  * (116) S. Badger, _Automating QCD amplitudes with on-shell methods_ , _J. Phys. Conf. Ser._ 762 (2016) 012057 [1605.02172].
  * (117) P. Bärnreuther, M. Czakon and P. Fiedler, _Virtual amplitudes and threshold behaviour of hadronic top-quark pair-production cross sections_ , _JHEP_ 02 (2014) 078 [1312.6279].
  * (118) R. Bonciani et al., _Two-Loop Four-Fermion Scattering Amplitude in QED_ , _Phys. Rev. Lett._ 128 (2022) 022002 [2106.13179].
  * (119) Y. Wang, L.L. Yang and B. Zhou, _FastGPL: a C++ library for fast evaluation of generalized polylogarithms_ , 2112.04122.

  *[N3LO]: next-to-next-to-next-to-leading order
  *[HVP]: hadronic vacuum polarisation
  *[SM]: Standard Model
  *[QCD]: quantum chromodynamics
  *[NNLO]: next-to-next-to-leading order
  *[QED]: quantum electrodynamics
  *[RVV]: real-double-virtual
  *[ISR]: initial state radiation
  *[NLO]: next-to-leading order
  *[VV]: double-virtual
  *[IR]: infrared
  *[IBP]: integration-by-parts
  *[UV]: ultraviolet
  *[MIs]: master integral
  *[DEs]: differential equation
  *[MPLs]: multiple polylogarithm
  *[MI]: master integral
  *[MPL]: multiple polylogarithm
  *[ISPs]: irreducible scalar product