††institutetext: aLaboratoire de Physique Théorique et Hautes Energies
(LPTHE), UMR 7589, Sorbonne Université et CNRS, 4 place Jussieu, 75252 Paris
Cedex 05, France
bBethe Center for Theoretical Physics, Universität Bonn, 53115 Bonn, Germany

# Two-loop massive QCD and QED helicity amplitudes for light-by-light
scattering

Ajjath A. H.a, Ekta Chaubeyb, Hua-Sheng Shaoa<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We present the analytic and compact two-loop helicity amplitudes for QCD and
QED corrections to the light-by-light scattering process with massive internal
fermions. We express the master integrals either in terms of multiple
polylogarithms or in terms of iterated integrals with $\mathrm{dlog}$ one-
forms. We also elaborate on optimizing the analytic results for each phase-
space region. This makes the numerical evaluation of the scattering amplitudes
fast, stable and suitable for phenomenological applications.

##  1 Introduction

Light-by-light (LbL) scattering is a fundamental quantum electrodynamics (QED)
process resulting from the quantum effects predicted in the 1930s
Heisenberg:1934pza ; Euler:1935zz ; Euler:1935qgl ; Heisenberg:1936nmg . It is
a loop-induced process via one-loop Feynman diagrams with charged internal
particles in the loop, as illustrated in the left panel of figure 1. The
lowest order of cross section is $\mathcal{O}(\alpha^{4})$, where $\alpha$ is
the electromagnetic fine structure constant. In the low energy limit, the
process can be effectively described by the Euler-Heisenberg Lagrangian
Heisenberg:1936nmg . The first complete calculation in QED was carried out by
Karplus and Neuman Karplus:1950zz  in 1951. The interest in studying LbL
revived recently, in particular from the point-of-view of searching for the
elusive beyond the Standard Model (BSM) signals. It has been acknowledged that
LbL can be used to probe the quartic anomalous gauge couplings
dEnterria:2013zqi , the axion-like particles Knapen:2016moh , the graviton-
like particles dEnterria:2023npy ; Atag:2010bh , the nonlinear Born-Infeld
extension of QED as well as in the context of the Standard Model (SM)
Ellis:2017edi , the photon-photon self-interaction from the noncommutative QED
Horvat:2020ycy , large extra dimensions Cheung:1999ja ; Davoudiasl:1999di ,
and supersymmetric particles Greiner:1992fz . LbL is also a background for
looking for new particles in the SM, such as ditauonium dEnterria:2022ysg ;
dEnterria:2023yao  and glueballs Greiner:1992fz .

In spite of the early theoretical proposals, the direct experimental
observations 111LbL has also been indirectly tested by the precise
measurements of the anomalous magnetic moments of the electron and the muon
and by the elastic scattering of a photon in Coulomb fields of nuclei via the
so-called Delbrück scattering. of the LbL scattering process
$\gamma\gamma\to\gamma\gamma$ were only achieved recently in the
ultraperipheral heavy-ion collisions (UPCs) at the LHC ATLAS:2017fur ;
CMS:2018erd ; ATLAS:2019azn ; ATLAS:2020hii . Its experimental feasibility,
first suggested in ref. dEnterria:2013zqi , can be mostly attributed to the
large coherent photon flux carried by the ultra-relativistic nucleus thanks to
its large charge number $Z$. For instance, the lead (Pb) beam used by the LHC
has $Z=82$. With UPCs, we have the largest magnitude of electromagnetic fields
produced in the laboratory in an experimentally controlled way. It enhances
the LbL cross section by $Z^{4}\approx 4.5\cdot 10^{7}$ in Pb-Pb collisions
with respect to the counterparts in proton-proton ($pp$) and electron-positron
($e^{-}e^{+}$) collisions. In a UPC exclusive process, the two initial beam
hadrons do not break, which can be identified with detectors at very forward
angles, Roman Pots and Zero Degree Calorimeters for proton-proton and nucleus-
nucleus collisions respectively. This usually provides us with distinct
signatures compared to the central and peripheral collisions. The UPC
exclusive processes have the virtue of low backgrounds due to their much
smaller particle multiplicities than in the inclusive scatterings, and thus
the experimental triggers can be loosened. This allows us to probe a different
kinematic regime with respect to the inclusive reactions that breakup the beam
hadrons.

The LHC measurements have been compared with the theoretical calculations by
the two Monte Carlo event generators: SuperChic Harland-Lang:2020veo  and
gamma-UPC Shao:2022cly  222Strictly speaking, gamma-UPC is a library for
calculating the coherent photon-photon flux in UPC. It has been integrated
into the two event generators HELAC-Onia Shao:2012iz ; Shao:2015vga  and
MadGraph5_aMC@NLO Alwall:2014hca  to enable the event generations of exclusive
UPC processes.. Both of them are based on the leading order (i.e., one-loop)
matrix element and yield the consistent LbL cross sections. However, the
theoretical cross sections are about 2 standard deviations below the ATLAS
measured value. It is therefore interesting to investigate the origin of this
discrepancy. On the theory side, it would be a natural step 333Another
possibly important missing effect is the higher order Coulomb
$\mathcal{O}(Z\alpha)$ corrections, which is however beyond the scope of our
paper. to check the size of the next-to-leading order (NLO) quantum
corrections that involve two-loop helicity amplitudes. In the low-energy
limit, the two-loop corrections to the Euler-Heisenberg Lagrangian have been
computed with the string-inspired approach Reuter:1996zm ; Fliegner:1997ra ;
Schubert:2001he ; Martin:2003gb . On the other hand, in the high energy limit,
two-loop amplitudes with massless internal fermions in QCD and QED Bern:2001dg
and in supersymmetric QED Binoth:2002xg  have also been available for two
decades. The aim of this paper is to present the analytic two-loop helicity
amplitudes with general massive internal fermions in the theories of QCD and
QED, where some representative graphs are shown in figure 1. This will allow
us to lift the above two approximations, which will also be relevant for the
theory-experiment comparison. We invest efforts in simplifying the initial
lengthy expressions into more compact and concise forms. This is not only for
the purpose of aesthetics, but also to improve the numerical stability of
helicity amplitudes at two loops. This optimisation becomes particularly
important when dealing with the substantial cancellations between individual
Feynman graphs as a consequence of the gauge invariance and infrared (IR)
finiteness. The corresponding two-loop Feynman integrals with massive internal
lines have also been studied previously in refs. Caron-Huot:2014lda ;
Xu:2018eos ; Maltoni:2018zvp ; Wang:2020nnr , though the results are only
partially available in the literature.

$\gamma$$\gamma$$\gamma$$\gamma$$f$$\gamma$$\gamma$$\gamma$$\gamma$$f$$\gamma$$\gamma$$\gamma$$\gamma$$f$$\gamma$$\gamma$$\gamma$$\gamma$$f$
Figure 1: Representative Feynman diagrams of $\gamma\gamma\to\gamma\gamma$ at
one loop (left) and two loops (rest). The internal wavy line can be either a
gluon or a photon.

The rest of the paper is organised as follows. We derive the general structure
of scattering amplitude to any loop order in section 2. The two-loop master
integrals are discussed in section 3. The analytic expressions of the one- and
two-loop helicity amplitudes with full fermion mass dependence are presented
in section 4. We draw our conclusion in section 5. The appendices collect the
additional information. We define the kinematic regions of our rationalised
variables as well as the analytic continuation in appendix A. The analytic
boundary conditions for the non-rationalised two-loop master integrals are
given in appendix B. We define the one-loop master integrals in appendix C.
Finally, we report the one-loop helicity amplitudes of the $W^{\pm}$ boson
contribution in appendix D, and collect the one- and two-loop helicity
amplitudes in the low-energy limit in appendix E.

##  2 General structure of scattering amplitude

In this section, we derive the independent form factors and helicity
amplitudes. A similar derivation exists in the literature Binoth:2002xg . In
addition, an alternative approach of obtaining helicity amplitudes is also
described in ref. Peraro:2020sfm . However, our derivation is neither
(completely) identical to ref. Binoth:2002xg  nor following ref.
Peraro:2020sfm . Therefore, for completeness, we present our derivation in
detail in this section.

Let us consider the process

$\gamma(p_{1},\lambda_{1})+\gamma(p_{2},\lambda_{2})+\gamma(p_{3},\lambda_{3})+\gamma(p_{4},\lambda_{4})\rightarrow
0,$ (1)

where all the momenta $p_{i}$ of the photons are incoming and their helicities
are denoted as $\lambda_{i}$. The amplitude is expressed as

$\mathcal{M}_{\lambda_{1}\lambda_{2}\lambda_{3}\lambda_{4}}=\left(\prod_{i=1}^{4}{\varepsilon_{\lambda_{i},\mu_{i}}(p_{i})}\right)\mathcal{M}^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}(p_{1},p_{2},p_{3},p_{4})$
(2)

with $\varepsilon_{\lambda_{i},\mu_{i}}(p_{i})$ being the polarisation vectors
of external photons. The Lorentz decomposition of the scattering tensor
$\mathcal{M}^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ is given as: 444The tensorial
structure needs to be extended with the Levi-Civita tensor when considering
the chiral-dependent electroweak corrections.

$\displaystyle\mathcal{M}^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}=$ $\displaystyle
A_{1}g^{\mu_{1}\mu_{2}}g^{\mu_{3}\mu_{4}}+A_{2}g^{\mu_{1}\mu_{3}}g^{\mu_{2}\mu_{4}}+A_{3}g^{\mu_{1}\mu_{4}}g^{\mu_{2}\mu_{3}}$
$\displaystyle+\sum_{j_{1},j_{2}=1}^{3}(B_{j_{1}j_{2}}^{1}g^{\mu_{1}\mu_{2}}p_{j_{1}}^{\mu_{3}}p_{j_{2}}^{\mu_{4}}+B^{2}_{j_{1}j_{2}}g^{\mu_{1}\mu_{3}}p_{j_{1}}^{\mu_{2}}p_{j_{2}}^{\mu_{4}}+B^{3}_{j_{1}j_{2}}g^{\mu_{1}\mu_{4}}p_{j_{1}}^{\mu_{2}}p_{j_{2}}^{\mu_{3}}$
$\displaystyle+B^{4}_{j_{1}j_{2}}g^{\mu_{2}\mu_{3}}p_{j_{1}}^{\mu_{1}}p_{j_{2}}^{\mu_{4}}+B^{5}_{j_{1}j_{2}}g^{\mu_{2}\mu_{4}}p_{j_{1}}^{\mu_{1}}p_{j_{2}}^{\mu_{3}}+B^{6}_{j_{1}j_{2}}g^{\mu_{3}\mu_{4}}p_{j_{1}}^{\mu_{1}}p_{j_{2}}^{\mu_{2}})$
$\displaystyle+\sum_{j_{1},j_{2},j_{3},j_{4}=1}^{3}C_{j_{1}j_{2}j_{3}j_{4}}p_{j_{1}}^{\mu_{1}}p_{j_{2}}^{\mu_{2}}p_{j_{3}}^{\mu_{3}}p_{j_{4}}^{\mu_{4}}.$
(3)

The coefficients $A_{i}$, $B_{jk}^{i}$ and $C_{ijkl}$ are functions of the
Mandelstam variables
$s=(p_{1}+p_{2})^{2},\;t=(p_{2}+p_{3})^{2},\;u=(p_{1}+p_{3})^{2}$, as well as
the masses of the particles in loops. Therefore, they should be understood as
$A_{i}(s,t,u),B_{jk}^{i}(s,t,u),C_{ijkl}(s,t,u)$. The on-shell conditions and
the conservation of momenta guarantee that the sum $s+t+u=0$. Thanks to the
parity invariance in QCD and QED, the term proportional to the Levi-Civita
tensor is forbidden. We have $3+6\times 3\times 3+3^{4}=138$ independent
functions in (3). However, many of the coefficients are related by the
following three symmetries:

  1. 1.

Tranversity: $\varepsilon_{\lambda_{i}}\cdot p_{i}=0$ for $i=1,2,3,4$. This
reduces to 57 functions to be considered.

  2. 2.

Bose symmetry: the tensor amplitude
$\mathcal{M}^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ is invariant under the exchange
of any two index pair $(p_{i},\mu_{i})$ with $i=1,2,3,4$, which gives us $7$
independent functions.

  3. 3.

Gauge symmetry: the gauge symmetry guarantees the following four Ward
identities

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\mathcal{M}^{p_{1}\varepsilon_{\lambda_{2}}(p_{2})\varepsilon_{\lambda_{3}}(p_{3})\varepsilon_{\lambda_{4}}(p_{4})}=\mathcal{M}^{\varepsilon_{\lambda_{1}}(p_{1})p_{2}\varepsilon_{\lambda_{3}}(p_{3})\varepsilon_{\lambda_{4}}(p_{4})}$
(4) $\displaystyle=$
$\displaystyle\mathcal{M}^{\varepsilon_{\lambda_{1}}(p_{1})\varepsilon_{\lambda_{2}}(p_{2})p_{3}\varepsilon_{\lambda_{4}}(p_{4})}=\mathcal{M}^{\varepsilon_{\lambda_{1}}(p_{1})\varepsilon_{\lambda_{2}}(p_{2})\varepsilon_{\lambda_{3}}(p_{3})p_{4}}.$

They lead to a few more identities to reduce to the 3 independent functions,
which are

$A_{1},\quad\Delta B_{11}^{1}=B_{11}^{1}-B_{12}^{1},\quad\Delta
C_{2111}=C_{2111}-C_{2112}.$ (5)

Our derivation and all the relations have been cross checked with explicit
calculations. The number of independent form factors (not functions) in
$\mathcal{M}^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ is $7$, instead of $15$ in
section 3.1 of ref. Binoth:2002xg . They are
$A_{1}(s,t,u),\;A_{1}(t,s,u),\;A_{1}(u,s,t)$, $\Delta
B_{11}^{1}(s,t,u),\;\Delta B_{11}^{1}(t,s,u),$ $\Delta B_{11}^{1}(u,s,t),$ and
$\Delta C_{2111}(s,t,u)$.

The reduced tensor structure becomes

$\displaystyle\mathcal{M}^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ $\displaystyle=$
$\displaystyle
A^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}+B^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}+C^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}},$
(6)

where

$\displaystyle A^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ $\displaystyle=$
$\displaystyle
A_{1}(s,t,u)\left(g^{\mu_{1}\mu_{2}}-\frac{p_{1}^{\mu_{2}}p_{2}^{\mu_{1}}}{p_{1}\cdot
p_{2}}\right)\left(g^{\mu_{3}\mu_{4}}+\frac{p_{1}^{\mu_{3}}+p_{2}^{\mu_{3}}}{p_{1}\cdot
p_{2}}p_{3}^{\mu_{4}}\right)$
$\displaystyle+[(\mu_{1},p_{1})\leftrightarrow(\mu_{3},p_{3}),s\leftrightarrow
t]$
$\displaystyle+[(\mu_{1},p_{1})\to(\mu_{3},p_{3}),(\mu_{2},p_{2})\to(\mu_{1},p_{1}),(\mu_{3},p_{3})\to(\mu_{2},p_{2}),s\to
u,t\to s,u\to t],$ $\displaystyle B^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$
$\displaystyle=$ $\displaystyle\Delta
B_{11}^{1}(s,t,u)\left[-\left(g^{\mu_{1}\mu_{2}}-\frac{p_{1}^{\mu_{2}}p_{2}^{\mu_{1}}}{p_{1}\cdot
p_{2}}\right)\left(p_{1}^{\mu_{3}}-\frac{u}{t}p_{2}^{\mu_{3}}\right)\left(p_{2}^{\mu_{4}}-\frac{u}{s}p_{3}^{\mu_{4}}\right)\right.$
$\displaystyle\left.+\left(p_{3}^{\mu_{1}}-\frac{u}{s}p_{2}^{\mu_{1}}\right)\left(p_{1}^{\mu_{2}}-\frac{s}{t}p_{3}^{\mu_{2}}\right)\left(g^{\mu_{3}\mu_{4}}+\frac{2p_{2}^{\mu_{3}}p_{1}^{\mu_{4}}}{t}+\frac{2p_{1}^{\mu_{3}}p_{2}^{\mu_{4}}}{u}\right)\right]$
$\displaystyle+[(\mu_{1},p_{1})\leftrightarrow(\mu_{3},p_{3}),s\leftrightarrow
t]$
$\displaystyle+[(\mu_{1},p_{1})\to(\mu_{3},p_{3}),(\mu_{2},p_{2})\to(\mu_{1},p_{1}),(\mu_{3},p_{3})\to(\mu_{2},p_{2}),s\to
u,t\to s,u\to t],$ $\displaystyle C^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$
$\displaystyle=$ $\displaystyle\Delta
C_{2111}(s,t,u)\frac{1}{st^{2}u}\left(sp_{3}^{\mu_{1}}-up_{2}^{\mu_{1}}\right)\left(sp_{3}^{\mu_{2}}-tp_{1}^{\mu_{2}}\right)\left(up_{2}^{\mu_{3}}-tp_{1}^{\mu_{3}}\right)\left(sp_{2}^{\mu_{4}}-up_{3}^{\mu_{4}}\right).$
(9)

In the above eqs.(2),(2), we have used the short-hands for the replacement
rules “$\to$” and the exchange rules “$\leftrightarrow$” with respect to the
first terms on the right-hand sides. The tensor
$C^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ is invariant under the exchange
$[(\mu_{1},p_{1})\leftrightarrow(\mu_{3},p_{3}),s\leftrightarrow t]$ and the
replacement
$[(\mu_{1},p_{1})\to(\mu_{3},p_{3}),(\mu_{2},p_{2})\to(\mu_{1},p_{1}),(\mu_{3},p_{3})\to(\mu_{2},p_{2}),s\to
u,t\to s,u\to t]$.

In the next step, we need to define the projection operators to obtain the
coefficients $A_{1}$, $\Delta B_{11}^{1}$ and $\Delta C_{2111}$. We can use
the following tensor

$\displaystyle\mathcal{P}^{\mu\nu}$
$\displaystyle=g^{\mu\nu}-\sum_{j=1}^{3}{\sum_{k=1}^{3}{p_{j}^{\mu}\left[\mathbb{G}^{-1}\right]_{jk}p_{k}^{\nu}}},\quad\mathcal{R}_{j}^{\nu}=\sum_{k=1}^{3}{\left[\mathbb{G}^{-1}\right]_{jk}p_{k}^{\nu}},$
(10)

where
$\mathbb{G}=\mathbb{G}(p_{1},p_{2},p_{3})=\left(\begin{array}[]{ccc}p_{1}^{2}&p_{1}\cdot
p_{2}&p_{1}\cdot p_{3}\\\ p_{1}\cdot p_{2}&p_{2}^{2}&p_{2}\cdot p_{3}\\\
p_{1}\cdot p_{3}&p_{2}\cdot p_{3}&p_{3}^{2}\\\ \end{array}\right)$ is the Gram
matrix. $\mathcal{P}$ is the projector that projects onto the
$(d-3)$-dimensional subspace perpendicular to the space spanned by the
external momenta (i.e.,
$p_{i,\mu}\mathcal{P}^{\mu\nu}=\mathcal{P}^{\mu\nu}p_{i,\nu}=0$).
$\mathcal{R}_{j,\nu}$ is the dual vector to $p_{j}^{\nu}$ relative to this
space (i.e., $\mathcal{R}_{j,\nu}p_{i}^{\nu}=\delta_{ji}$). In addition, we
have the following relations

$\displaystyle\mathcal{P}^{\mu}_{\mu}$ $\displaystyle=$ $\displaystyle d-3,$
$\displaystyle\mathcal{P}^{\mu\rho}\mathcal{P}_{\rho}^{\nu}$ $\displaystyle=$
$\displaystyle
g^{\mu\nu}-2\sum_{j=1}^{3}{\sum_{k=1}^{3}{p_{j}^{\mu}\left[\mathbb{G}^{-1}\right]_{jk}p_{k}^{\nu}}}+\sum_{j,k,l,n=1}^{3}{p_{j}^{\mu}\left[\mathbb{G}^{-1}\right]_{jk}p_{k}^{\rho}p_{l,\rho}\left[\mathbb{G}^{-1}\right]_{ln}p_{n}^{\nu}}$
$\displaystyle=$ $\displaystyle
g^{\mu\nu}-2\sum_{j=1}^{3}{\sum_{k=1}^{3}{p_{j}^{\mu}\left[\mathbb{G}^{-1}\right]_{jk}p_{k}^{\nu}}}+\sum_{j,k,l,n=1}^{3}{p_{j}^{\mu}\left[\mathbb{G}^{-1}\right]_{jk}\left[\mathbb{G}\right]_{kl}\left[\mathbb{G}^{-1}\right]_{ln}p_{n}^{\nu}}$
$\displaystyle=$ $\displaystyle\mathcal{P}^{\mu\nu},$
$\displaystyle\mathcal{R}_{i}^{\mu}\mathcal{R}_{j,\mu}$ $\displaystyle=$
$\displaystyle\sum_{k,l=1}^{3}{\left[\mathbb{G}^{-1}\right]_{ik}p_{k}^{\mu}\left[\mathbb{G}^{-1}\right]_{jl}p_{l,\mu}}=\left[\mathbb{G}^{-1}\right]_{ji}=\left[\mathbb{G}^{-1}\right]_{ij},$
$\displaystyle\mathcal{P}^{\mu\nu}\mathcal{R}_{j,\nu}$ $\displaystyle=$
$\displaystyle\mathcal{R}_{j}^{\mu}-\sum_{k,l,n=1}^{3}{p_{k}^{\mu}\left[\mathbb{G}^{-1}\right]_{kl}p_{l}^{\nu}\left[\mathbb{G}^{-1}\right]_{jn}p_{n,\nu}}=0=\mathcal{R}_{j,\mu}\mathcal{P}^{\mu\nu},$
(11)

where we have denoted the space-time dimension by $d=4-2\epsilon$ with
$\epsilon$ being the dimensional regularisation parameter that ought to be
taken infinitesimal at the end.

We can obtain the $A_{i}$ tensor coefficients using the following operations

$\displaystyle\tilde{A}_{1}$
$\displaystyle=\frac{1}{(d-1)(d-3)}\mathcal{P}_{\mu_{1}\mu_{2}}\mathcal{P}_{\mu_{3}\mu_{4}}\mathcal{M}^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}=\frac{1}{(d-1)}\left[(d-3)A_{1}+A_{2}+A_{3}\right],$
(12) $\displaystyle\tilde{A}_{2}$
$\displaystyle=\frac{1}{(d-1)(d-3)}\mathcal{P}_{\mu_{1}\mu_{3}}\mathcal{P}_{\mu_{2}\mu_{4}}\mathcal{M}^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}=\frac{1}{(d-1)}\left[A_{1}+(d-3)A_{2}+A_{3}\right],$
(13) $\displaystyle\tilde{A}_{3}$
$\displaystyle=\frac{1}{(d-1)(d-3)}\mathcal{P}_{\mu_{1}\mu_{4}}\mathcal{P}_{\mu_{2}\mu_{3}}\mathcal{M}^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}=\frac{1}{(d-1)}\left[A_{1}+A_{2}+(d-3)A_{3}\right].$
(14)

In other words, we have

$\displaystyle\begin{pmatrix}\tilde{A}_{1}\\\ \tilde{A}_{2}\\\
\tilde{A}_{3}\end{pmatrix}=\frac{1}{(d-1)}\begin{pmatrix}d-3&1&1\\\ 1&d-3&1\\\
1&1&d-3\end{pmatrix}\begin{pmatrix}A_{1}\\\ A_{2}\\\ A_{3}\end{pmatrix}\,,$
(15)

and

$\displaystyle\begin{pmatrix}A_{1}\\\ A_{2}\\\
A_{3}\end{pmatrix}=\frac{1}{(d-4)}\begin{pmatrix}d-2&-1&-1\\\ -1&d-2&-1\\\
-1&-1&d-2\end{pmatrix}\begin{pmatrix}\tilde{A}_{1}\\\ \tilde{A}_{2}\\\
\tilde{A}_{3}\end{pmatrix}\,.$ (16)

In the above, the right hand side appears to introduce the spurious pole
$1/(d-4)=-1/(2\epsilon)$. On the other hand, the sum

$\displaystyle A_{S}=A_{1}+A_{2}+A_{3}$ $\displaystyle=$
$\displaystyle\tilde{A}_{1}+\tilde{A}_{2}+\tilde{A}_{3}$ (17)

is free of the spurious pole. As we discuss later, the amplitude only depends
on $A_{S}$. Hence, there are no spurious poles.

It is straightforward to define the projectors for the coefficients
$B^{1}_{kl}$ via

$\displaystyle\tilde{B}^{1}_{kl}$ $\displaystyle=$
$\displaystyle\underbrace{\frac{1}{(d-3)}\mathcal{P}_{\mu_{1}\mu_{2}}\mathcal{R}_{k,\mu_{3}}\mathcal{R}_{l,\mu_{4}}\mathcal{M}^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}}_{=\hat{B}^{1}_{kl}}-\left[\mathbb{G}^{-1}\right]_{kl}A_{1}$
(18) $\displaystyle=$
$\displaystyle\frac{1}{(d-3)}\sum_{j_{1},j_{2}=1}^{3}{\left[(d-3)B^{1}_{j_{1}j_{2}}\delta_{kj_{1}}\delta_{lj_{2}}\right]}=B^{1}_{kl}.$

Similarly, we have

$\displaystyle B^{2}_{kl}$ $\displaystyle=$
$\displaystyle\frac{1}{(d-3)}\mathcal{P}_{\mu_{1}\mu_{3}}\mathcal{R}_{k,\mu_{2}}\mathcal{R}_{l,\mu_{4}}\mathcal{M}^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}-\left[\mathbb{G}^{-1}\right]_{kl}A_{2},$
$\displaystyle B^{3}_{kl}$ $\displaystyle=$
$\displaystyle\frac{1}{(d-3)}\mathcal{P}_{\mu_{1}\mu_{4}}\mathcal{R}_{k,\mu_{2}}\mathcal{R}_{l,\mu_{3}}\mathcal{M}^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}-\left[\mathbb{G}^{-1}\right]_{kl}A_{3},$
$\displaystyle B^{4}_{kl}$ $\displaystyle=$
$\displaystyle\frac{1}{(d-3)}\mathcal{P}_{\mu_{2}\mu_{3}}\mathcal{R}_{k,\mu_{1}}\mathcal{R}_{l,\mu_{4}}\mathcal{M}^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}-\left[\mathbb{G}^{-1}\right]_{kl}A_{3},$
$\displaystyle B^{5}_{kl}$ $\displaystyle=$
$\displaystyle\frac{1}{(d-3)}\mathcal{P}_{\mu_{2}\mu_{4}}\mathcal{R}_{k,\mu_{1}}\mathcal{R}_{l,\mu_{3}}\mathcal{M}^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}-\left[\mathbb{G}^{-1}\right]_{kl}A_{2},$
$\displaystyle B^{6}_{kl}$ $\displaystyle=$
$\displaystyle\frac{1}{(d-3)}\mathcal{P}_{\mu_{3}\mu_{4}}\mathcal{R}_{k,\mu_{1}}\mathcal{R}_{l,\mu_{2}}\mathcal{M}^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}-\left[\mathbb{G}^{-1}\right]_{kl}A_{1}\,.$
(19)

Using Bose symmetry and Ward identities, we can derive the following relations
between the form factors

$\displaystyle\Delta B^{1}_{11}(s,t,u)$ $\displaystyle=$ $\displaystyle
B^{1}_{11}(s,t,u)-B^{1}_{12}(s,t,u)=\underbrace{\left(\hat{B}^{1}_{11}(s,t,u)-\hat{B}^{1}_{12}(s,t,u)\right)}_{=\Delta\hat{B}^{1}_{11}(s,t,u)}-\frac{A_{1}(s,t,u)}{u}.$

With the same spirit, we can define the projectors for the coefficients
$C_{ijkl}$ via

$\displaystyle\tilde{C}_{ijkl}$ $\displaystyle=$
$\displaystyle\underbrace{\mathcal{R}_{i,\mu_{1}}\mathcal{R}_{j,\mu_{2}}\mathcal{R}_{k,\mu_{3}}\mathcal{R}_{l,\mu_{4}}\mathcal{M}^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}}_{=\hat{C}_{ijkl}}-\left\\{\left[\mathbb{G}^{-1}\right]_{ij}B^{1}_{kl}+\left[\mathbb{G}^{-1}\right]_{ik}B^{2}_{jl}+\left[\mathbb{G}^{-1}\right]_{il}B^{3}_{jk}\right.$
(21)
$\displaystyle\left.+\left[\mathbb{G}^{-1}\right]_{jk}B^{4}_{il}+\left[\mathbb{G}^{-1}\right]_{jl}B^{5}_{ik}+\left[\mathbb{G}^{-1}\right]_{kl}B^{6}_{ij}\right\\}$
$\displaystyle-\left\\{\left[\mathbb{G}^{-1}\right]_{ij}\left[\mathbb{G}^{-1}\right]_{kl}A_{1}+\left[\mathbb{G}^{-1}\right]_{ik}\left[\mathbb{G}^{-1}\right]_{jl}A_{2}+\left[\mathbb{G}^{-1}\right]_{il}\left[\mathbb{G}^{-1}\right]_{jk}A_{3}\right\\}$
$\displaystyle=$ $\displaystyle C_{ijkl}\,.$

Furthermore, by using the relations between form factors from Bose symmetry
and Ward identities, we obtain

$\displaystyle\Delta C_{2111}(s,t,u)$
$\displaystyle=C_{2111}(s,t,u)-C_{2112}(s,t,u)$
$\displaystyle=\underbrace{\left(\hat{C}_{2111}(s,t,u)-\hat{C}_{2112}(s,t,u)\right)}_{=\Delta\hat{C}_{2111}(s,t,u)}+\;\frac{A_{S}(s,t,u)}{su}.$

With the above preparations, the helicity-summed amplitude square takes the
form

$\displaystyle\sum_{\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{4}=\pm}{\left|\mathcal{M}_{\lambda_{1}\lambda_{2}\lambda_{3}\lambda_{4}}\right|^{2}}=2\left|A_{S}(s,t,u)\right|^{2}+2\sum_{i=1}^{3}{\left|A_{i}(s,t,u)\right|^{2}}$
$\displaystyle\qquad-2su\Re{\left(\Delta
C_{2111}(s,t,u)A_{S}^{*}(s,t,u)\right)}+s^{2}u^{2}\left|\Delta
C_{2111}(s,t,u)\right|^{2}$ $\displaystyle\qquad+4\left[u\Re{\left(\Delta
B^{1}_{11}(s,t,u)A_{1}^{*}(s,t,u)\right)}+s\Re{\left(\Delta
B^{1}_{11}(t,u,s)A_{1}^{*}(t,u,s)\right)}+t\Re{\left(\Delta
B^{1}_{11}(u,s,t)A_{1}^{*}(u,s,t)\right)}\right]$
$\displaystyle\qquad+2u^{2}\left|\Delta
B^{1}_{11}(s,t,u)\right|^{2}+2s^{2}\left|\Delta
B^{1}_{11}(t,u,s)\right|^{2}+2t^{2}\left|\Delta B^{1}_{11}(u,s,t)\right|^{2}$
(23)
$\displaystyle\quad=\left|A_{S}(s,t,u)\right|^{2}+2u^{2}\left|\Delta\hat{B}^{1}_{11}(s,t,u)\right|^{2}+2s^{2}\left|\Delta\hat{B}^{1}_{11}(t,u,s)\right|^{2}+2t^{2}\left|\Delta\hat{B}^{1}_{11}(u,s,t)\right|^{2}$
(24)
$\displaystyle\qquad+s^{2}u^{2}\left|\Delta\hat{C}_{2111}(s,t,u)\right|^{2},$

where we have used the relations

$\displaystyle s\Delta B^{1}_{11}(t,u,s)$ $\displaystyle=$ $\displaystyle
u\Delta B^{1}_{11}(t,s,u),\quad
s\Delta\hat{B}^{1}_{11}(t,u,s)=u\Delta\hat{B}^{1}_{11}(t,s,u),\quad
A_{1}(t,u,s)=A_{1}(t,s,u),$ $\displaystyle A_{2}(s,t,u)$ $\displaystyle=$
$\displaystyle A_{1}(u,s,t),\quad A_{3}(s,t,u)=A_{1}(t,s,u).$ (25)

As we see from eq.(24), at the level of helicity-summed amplitude square, we
only need to calculate $5$ independent form factors $A_{S}(s,t,u)$,
$\Delta\hat{B}^{1}_{11}(s,t,u),\;\Delta\hat{B}^{1}_{11}(u,s,t),\;\Delta\hat{B}^{1}_{11}(t,u,s)$,
and $\Delta\hat{C}_{2111}(s,t,u)$. Among these, there are no interference
terms. The form factors $A_{S}(s,t,u),su\Delta C_{2111}(s,t,u),$ and
$su\Delta\hat{C}_{2111}(s,t,u)$ are fully symmetric in the Mandelstam
variables $s,t,u$.

Furthermore, the helicity amplitudes can be expressed as a linear combination
of the form factors. Like the case of form factors, there are $5$ independent
helicity amplitudes.555Note that our convention of helicities is different
with respect to ref. Bern:2001dg , because we take all momenta as incoming.
The helicities of the final particles in ref. Bern:2001dg  should be flipped
in order to match to our case. Up to a global phase, we can write them as:

$\displaystyle\mathcal{M}_{++++}$ $\displaystyle=$ $\displaystyle
A_{S}(s,t,u)+\frac{u}{2}\Delta B^{1}_{11}(s,t,u)+\frac{s}{2}\Delta
B^{1}_{11}(t,u,s)+\frac{t}{2}\Delta B^{1}_{11}(u,s,t)-\frac{su}{4}\Delta
C_{2111}(s,t,u)$ $\displaystyle=$
$\displaystyle\frac{1}{4}A_{S}(s,t,u)+\frac{u}{2}\Delta\hat{B}^{1}_{11}(s,t,u)+\frac{s}{2}\Delta\hat{B}^{1}_{11}(t,u,s)+\frac{t}{2}\Delta\hat{B}^{1}_{11}(u,s,t)-\frac{su}{4}\Delta\hat{C}_{2111}(s,t,u),$
$\displaystyle\mathcal{M}_{-+++}$ $\displaystyle=$
$\displaystyle\frac{su}{4}\Delta
C_{2111}(s,t,u)=\frac{1}{4}A_{S}(s,t,u)+\frac{su}{4}\Delta\hat{C}_{2111}(s,t,u),$
(27) $\displaystyle\mathcal{M}_{--++}$ $\displaystyle=$ $\displaystyle
A_{1}(s,t,u)+\frac{u}{2}\Delta B^{1}_{11}(s,t,u)-\frac{s}{2}\Delta
B^{1}_{11}(t,u,s)-\frac{t}{2}\Delta B^{1}_{11}(u,s,t)-\frac{su}{4}\Delta
C_{2111}(s,t,u)$ $\displaystyle=$
$\displaystyle\frac{1}{4}A_{S}(s,t,u)+\frac{u}{2}\Delta\hat{B}^{1}_{11}(s,t,u)-\frac{s}{2}\Delta\hat{B}^{1}_{11}(t,u,s)-\frac{t}{2}\Delta\hat{B}^{1}_{11}(u,s,t)-\frac{su}{4}\Delta\hat{C}_{2111}(s,t,u),$
$\displaystyle\mathcal{M}_{+-+-}$ $\displaystyle=$ $\displaystyle
A_{1}(u,s,t)-\frac{u}{2}\Delta B^{1}_{11}(s,t,u)-\frac{s}{2}\Delta
B^{1}_{11}(t,u,s)+\frac{t}{2}\Delta B^{1}_{11}(u,s,t)-\frac{su}{4}\Delta
C_{2111}(s,t,u)$ $\displaystyle=$
$\displaystyle\frac{1}{4}A_{S}(s,t,u)-\frac{u}{2}\Delta\hat{B}^{1}_{11}(s,t,u)-\frac{s}{2}\Delta\hat{B}^{1}_{11}(t,u,s)+\frac{t}{2}\Delta\hat{B}^{1}_{11}(u,s,t)-\frac{su}{4}\Delta\hat{C}_{2111}(s,t,u),$
$\displaystyle\mathcal{M}_{+--+}$ $\displaystyle=$ $\displaystyle
A_{1}(t,s,u)-\frac{u}{2}\Delta B^{1}_{11}(s,t,u)+\frac{s}{2}\Delta
B^{1}_{11}(t,u,s)-\frac{t}{2}\Delta B^{1}_{11}(u,s,t)-\frac{su}{4}\Delta
C_{2111}(s,t,u)$ $\displaystyle=$
$\displaystyle\frac{1}{4}A_{S}(s,t,u)-\frac{u}{2}\Delta\hat{B}^{1}_{11}(s,t,u)+\frac{s}{2}\Delta\hat{B}^{1}_{11}(t,u,s)-\frac{t}{2}\Delta\hat{B}^{1}_{11}(u,s,t)-\frac{su}{4}\Delta\hat{C}_{2111}(s,t,u).$

In addition, if we take into account permutations, the following relations
hold between the $3$ independent helicity amplitudes:

$\displaystyle\mathcal{M}_{+-+-}$ $\displaystyle=$
$\displaystyle\left.\mathcal{M}_{--++}\right|_{s\leftrightarrow
u}=\left.\mathcal{M}_{+--+}\right|_{t\leftrightarrow u}\,.$ (31)

While the helicity amplitudes $\mathcal{M}_{++++}$ and $\mathcal{M}_{-+++}$
are fully symmetric in the arguments $(s,t,u)$, the remaining ones
$\\{\mathcal{M}_{--++},~{}\mathcal{M}_{+-+-},~{}\mathcal{M}_{+--+}\\}$ are
symmetric in $\\{(t,u),(s,t),(s,u)\\}$ respectively. The amplitudes in other
helicity configurations can be related to the above mentioned $5$ as follows:

$\displaystyle\mathcal{M}_{----}$ $\displaystyle=$
$\displaystyle\mathcal{M}_{++++},\quad\mathcal{M}_{++--}=\mathcal{M}_{--++},\quad\mathcal{M}_{-+-+}=\mathcal{M}_{+-+-},\quad\mathcal{M}_{-++-}=\mathcal{M}_{+--+},$
$\displaystyle\mathcal{M}_{+-++}$ $\displaystyle=$
$\displaystyle\mathcal{M}_{++-+}=\mathcal{M}_{+++-}=\mathcal{M}_{+---}=\mathcal{M}_{-+--}=\mathcal{M}_{--+-}=\mathcal{M}_{---+}=\mathcal{M}_{-+++}.$

Summing these helicity amplitudes, it is straightforward to verify that
eqs.(23) and (24) hold.

##  3 Two-loop master integrals

It is well known that a form factor or an amplitude at $L$-loop can be written
as a linear combination of a small set of $L$-loop Feynman integrals, known as
master integrals. In this section, we first present the details for obtaining
analytic expressions of two-loop master integrals for LbL scattering, before
discussing in detail their numerical evaluations. The one-loop master
integrals are provided in appendix C. In dimensional regularisation, the two-
loop QED amplitude for a given fermion species with mass $m_{f}$, generated by
Qgraf Nogueira:1991ex  and FeynArts Hahn:2000kx , can be reduced into 103
$d$-dimensional master integrals across all the permutations by applying
integration-by-parts (IBP) identities with the help of Kira Maierhofer:2017gsa
; Klappert:2020nbg  or FiniteFlow Peraro:2019svx +LiteRed Lee:2012cn . Modulo
the permutations in the Mandelstam variables $s,t,u$, there are $30$ master
integrals. One among them is a product of a one-loop box integral and a one-
loop tadpole integral, which can be solved separately. The remaining $29$
integrals form a closed system in the method of differential equations
Kotikov:1990kg ; Remiddi:1997ny ; Gehrmann:1999as ; Argeri:2007up  in $d$
dimensions. We now explain how to analytically solve these 29 master
integrals.

The two-loop master integrals for LbL in the top-sector have been investigated
in ref. Caron-Huot:2014lda . In this context, the new ingredients of this
article is the presentation of the analytic results in a form that is most
suitable for phenomenological applications, facilitating the fast numerical
evaluation of the scattering amplitudes across all the physical phase-space
regions. In addition, to ensure the comprehensiveness of this article, we
elaborate on all the essential steps needed to replicate our results, even at
the risk of providing details that may already be well understood. We express
the analytic results of the master integrals in terms of Chen iterated
integrals Chen:1977oja  with logarithmic one-forms, and cast the results in a
form suitable for phenomenological applications. The family of master
integrals of interest is expressed as

$\displaystyle I_{a_{1},\cdots,a_{9}}^{(2)}(s_{ij},s_{jk},s_{ik})$
$\displaystyle=\left(\frac{e^{\epsilon\gamma_{E}}m_{f}^{2\epsilon}}{i\pi^{\frac{d}{2}}}\right)^{2}\int{\mathrm{d}}^{d}\ell_{1}{\mathrm{d}}^{d}\ell_{2}\frac{1}{D_{1}^{a_{1}}D_{2}^{a_{2}}D_{3}^{a_{3}}D_{4}^{a_{4}}D_{5}^{a_{5}}D_{6}^{a_{6}}D_{7}^{a_{7}}D_{8}^{a_{8}}D_{9}^{a_{9}}}\,,$
(33)

where $a_{i}\in\mathbb{Z}$. We represent the loop momenta through $\ell_{i}$
and Euler-Mascheroni constant by $\gamma_{E}$. The inverse propagators are

$\displaystyle D_{1}$ $\displaystyle=$
$\displaystyle\ell_{1}^{2}-m_{f}^{2}+i0^{+},\quad
D_{2}=(\ell_{1}+p_{i})^{2}-m_{f}^{2}+i0^{+},\quad
D_{3}=(\ell_{1}+p_{i}+p_{j})^{2}-m_{f}^{2}+i0^{+},$ $\displaystyle D_{4}$
$\displaystyle=$
$\displaystyle(\ell_{1}+p_{i}+p_{j}+p_{k})^{2}-m_{f}^{2}+i0^{+},\qquad
D_{5}=\ell_{2}^{2}-m_{f}^{2}+i0^{+},$ $\displaystyle D_{6}$ $\displaystyle=$
$\displaystyle(\ell_{2}+p_{i})^{2}-m_{f}^{2}+i0^{+},\qquad
D_{7}=(\ell_{2}+p_{i}+p_{j})^{2}-m_{f}^{2}+i0^{+},$ $\displaystyle D_{8}$
$\displaystyle=$
$\displaystyle(\ell_{2}+p_{i}+p_{j}+p_{k})^{2}-m_{f}^{2}+i0^{+},\qquad
D_{9}=(\ell_{1}-\ell_{2})^{2}+i0^{+}.$

The introduction of $i0^{+}$ ensures the prescription of the Feynman
propagators. We keep the general external momentum dependence with $i,j,k\in\
\left\\{1,2,3,4\right\\}$ and $i\neq j,i\neq k,j\neq k$. The arguments of
$I_{a_{1},\cdots,a_{9}}^{(2)}$ are defined as
$s_{ij}=(p_{i}+p_{j})^{2},s_{jk}=(p_{j}+p_{k})^{2},s_{ik}=(p_{i}+p_{k})^{2}$,
which are simply a permutation of $s,t,u$.

As mentioned above, to obtain the analytic results of the master integrals we
use the method of differential equations. Within this approach, it is
convenient to choose a canonical basis of master integrals that allow us to
set up the differential equations into the so-called $\epsilon$-form
Henn:2013pwa , which enables to solve the system afterwards. Our choice of the
canonical basis for the two-loop master integrals is guided by the one
outlined in ref. Caron-Huot:2014lda , although not completely identical. For
completeness of this article, we present our chosen basis
$\vec{f}^{(2)}=(f^{(2)}_{1},\cdots,f^{(2)}_{29})$ as follows:

$\displaystyle f^{(2)}_{1}$ $\displaystyle=$
$\displaystyle\epsilon^{2}I^{(2)}_{0,0,2,0,0,0,0,2,0},\quad
f^{(2)}_{2}=\epsilon^{2}\sqrt{(s_{ij}-4m_{f}^{2})s_{ij}}I^{(2)}_{2,0,1,0,2,0,0,0,0},$
$\displaystyle f^{(2)}_{3}$ $\displaystyle=$
$\displaystyle-\epsilon^{2}(s_{ij}-4m_{f}^{2})s_{ij}I^{(2)}_{2,0,1,0,2,0,1,0,0},\quad
f^{(2)}_{4}=-\epsilon^{3}s_{ij}I^{(2)}_{2,0,0,0,1,0,1,1,0},$ $\displaystyle
f^{(2)}_{5}$ $\displaystyle=$
$\displaystyle\epsilon^{3}s_{ij}\sqrt{(s_{ij}-4m_{f}^{2})s_{ij}}I^{(2)}_{2,0,1,0,1,0,1,1,0},\quad
f^{(2)}_{6}=\epsilon^{2}s_{ij}I^{(2)}_{2,0,0,0,0,0,2,0,1},$ $\displaystyle
f_{7}^{(2)}$ $\displaystyle=$
$\displaystyle\epsilon^{2}\sqrt{(s_{ij}-4m_{f}^{2})s_{ij}}\left(2I^{(2)}_{2,0,0,0,0,0,1,0,2}+I^{(2)}_{2,0,0,0,0,0,2,0,1}\right),\quad
f^{(2)}_{8}=\epsilon^{2}m_{f}^{2}s_{ij}I^{(2)}_{3,0,0,0,0,0,1,1,1},$
$\displaystyle f^{(2)}_{9}$ $\displaystyle=$
$\displaystyle\epsilon^{2}m_{f}^{2}s_{ij}I^{(2)}_{1,0,1,0,0,0,0,3,1},\quad
f^{(2)}_{10}=-\epsilon^{3}s_{ij}I^{(2)}_{1,0,1,0,0,0,0,2,1},$ $\displaystyle
f^{(2)}_{11}$ $\displaystyle=$
$\displaystyle\frac{3}{2}\epsilon^{3}\sqrt{(s_{ij}-4m_{f}^{2})s_{ij}}I^{(2)}_{1,0,1,0,0,0,0,2,1}+\epsilon^{2}m_{f}^{2}\sqrt{(s_{ij}-4m_{f}^{2})s_{ij}}(I^{(2)}_{1,0,1,0,0,0,0,3,1}+I^{(2)}_{2,0,1,0,0,0,0,2,1}),$
$\displaystyle f^{(2)}_{12}$ $\displaystyle=$
$\displaystyle\epsilon^{4}s_{ij}I^{(2)}_{1,0,1,0,0,0,1,1,1},\quad
f^{(2)}_{13}=-\epsilon^{4}s_{ij}\sqrt{(s_{ij}-4m_{f}^{2})s_{ij}}I^{(2)}_{1,0,1,0,1,0,1,1,1},$
$\displaystyle f^{(2)}_{14}$ $\displaystyle=$
$\displaystyle\epsilon^{2}s_{jk}I^{(2)}_{0,2,0,0,0,0,0,2,1},\quad
f^{(2)}_{15}=\epsilon^{2}\sqrt{(s_{jk}-4m_{f}^{2})s_{jk}}(2I^{(2)}_{0,2,0,0,0,0,0,1,2}+I^{(2)}_{0,2,0,0,0,0,0,2,1}),$
$\displaystyle f^{(2)}_{16}$ $\displaystyle=$
$\displaystyle\epsilon^{2}m_{f}^{2}s_{jk}I^{(2)}_{1,1,0,0,0,0,0,3,1},\quad
f^{(2)}_{17}=\epsilon^{4}s_{jk}I^{(2)}_{1,1,0,0,1,0,0,1,1},$ $\displaystyle
f^{(2)}_{18}$ $\displaystyle=$
$\displaystyle-\epsilon^{3}\sqrt{s_{ij}(s_{ij}(s_{jk}-m_{f}^{2})^{2}-4m_{f}^{2}s_{jk}^{2})}I^{(2)}_{1,1,1,0,0,0,0,2,1},$
$\displaystyle f^{(2)}_{19}$ $\displaystyle=$
$\displaystyle\epsilon^{2}\sqrt{s_{ij}s_{jk}(s_{ij}s_{jk}-4m_{f}^{2}(s_{ij}+s_{jk}))}(\epsilon
I^{(2)}_{1,1,1,0,0,0,0,2,1}+m_{f}^{2}I^{(2)}_{1,1,1,0,0,0,0,3,1}),$
$\displaystyle f^{(2)}_{20}$ $\displaystyle=$
$\displaystyle-\epsilon^{3}s_{ij}(I^{(2)}_{1,1,1,-1,0,0,0,2,1}+m_{f}^{2}I^{(2)}_{1,1,1,0,0,0,0,2,1}),\quad
f^{(2)}_{21}=-\epsilon^{4}s_{ik}I^{(2)}_{0,1,1,0,1,0,0,1,1},$ $\displaystyle
f^{(2)}_{22}$ $\displaystyle=$
$\displaystyle\epsilon^{3}\sqrt{s_{ij}s_{jk}(s_{ij}s_{jk}-4m_{f}^{2}(s_{ij}+s_{jk}))}I^{(2)}_{0,1,1,0,1,0,0,1,2},$
$\displaystyle f^{(2)}_{23}$ $\displaystyle=$
$\displaystyle-\epsilon^{3}m_{f}^{2}s_{ik}(I^{(2)}_{0,1,2,0,1,0,0,1,1}+I^{(2)}_{0,2,1,0,1,0,0,1,1}),\quad
f^{(2)}_{24}=\epsilon^{4}s_{ij}^{2}I^{(2)}_{1,1,1,0,1,0,1,1,0},$
$\displaystyle f^{(2)}_{25}$ $\displaystyle=$
$\displaystyle\epsilon^{4}\sqrt{s_{ij}s_{jk}(s_{ij}s_{jk}-4m_{f}^{2}(s_{ij}+s_{jk}))}I^{(2)}_{1,1,0,0,1,0,1,1,1},$
$\displaystyle f^{(2)}_{26}$ $\displaystyle=$
$\displaystyle\epsilon^{4}s_{jk}I^{(2)}_{0,1,1,0,1,0,0,1,1}+\epsilon^{4}s_{ij}I^{(2)}_{1,1,0,0,1,-1,1,1,1},$
$\displaystyle f^{(2)}_{27}$ $\displaystyle=$
$\displaystyle\epsilon^{4}\sqrt{(s_{ij}-4m_{f}^{2})s_{ij}}\sqrt{s_{ij}s_{jk}(s_{ij}s_{jk}-4m_{f}^{2}(s_{ij}+s_{jk}))}I^{(2)}_{1,1,1,0,1,0,1,1,1},$
$\displaystyle f^{(2)}_{28}$ $\displaystyle=$
$\displaystyle\epsilon^{4}\sqrt{(s_{ij}-4m_{f}^{2})s_{ij}}(2s_{jk}I^{(2)}_{1,1,0,0,1,0,1,1,1}+s_{ij}I^{(2)}_{1,1,1,0,1,-1,1,1,1}),$
$\displaystyle f^{(2)}_{29}$ $\displaystyle=$
$\displaystyle-\epsilon^{4}s_{ij}\left[I^{(2)}_{1,0,1,0,1,0,1,0,1}+2I^{(2)}_{1,1,0,0,1,-1,1,1,1}-I^{(2)}_{1,1,1,-1,1,-1,1,1,1}\right.$
(35)
$\displaystyle\left.+s_{ij}\left(I^{(2)}_{1,0,1,0,1,0,1,1,1}-I^{(2)}_{1,1,1,0,1,-1,1,1,1}\right)+s_{jk}\left(I^{(2)}_{1,1,1,0,1,0,1,1,0}-2I^{(2)}_{1,1,0,0,1,0,1,1,1}\right)\right],$

where we have suppressed the arguments $s_{ij},s_{jk},s_{ik}$ in
$I^{(2)}_{a_{1},\cdots,a_{9}}$ and the arguments of $f_{n}^{(2)}$ that are
$x_{ij}=s_{ij}/m_{f}^{2},x_{jk}=s_{jk}/m_{f}^{2},x_{ik}=s_{ik}/m_{f}^{2}$.
Each integral in the canonical basis (35) is dimensionless and has uniform
transcendental (UT) weight properties. Moreover, the system of differential
equations of $\vec{f}^{(2)}$ attains the form

${\mathrm{d}}\vec{f}^{(2)}=\epsilon\;{\mathrm{d}}A^{(2)}\vec{f}^{(2)},$ (36)

where the differential equation matrix $A^{(2)}$ does not depend on
$\epsilon$. Due to the inclusion of massive loop propagators, the presence of
square roots is unavoidable. The square roots present in our case are
$\sqrt{s_{ij}(s_{ij}-4m_{f}^{2})}$,$\sqrt{s_{jk}(s_{jk}-4m_{f}^{2})}$,
$\sqrt{s_{ij}s_{jk}(s_{ij}s_{jk}-4m_{f}^{2}(s_{ij}+s_{jk}))}$ and
$\sqrt{s_{ij}(m_{f}^{4}s_{ij}-2m_{f}^{2}s_{jk}(s_{ij}+2s_{jk})+s_{ij}s_{jk}^{2})}$.
We can rationalise the first three square roots simultaneously by choosing the
same variables as suggested in ref. Caron-Huot:2014lda :

$\displaystyle\frac{s_{ij}}{m_{f}^{2}}$ $\displaystyle=$
$\displaystyle-\frac{4(w-z)^{2}}{(1-w^{2})(1-z^{2})},\quad\frac{s_{jk}}{m_{f}^{2}}=-\frac{(w-z)^{2}}{wz},\quad
s_{ik}=-s_{ij}-s_{jk}.$ (37)

However, we still have one square root
$\rho=\sqrt{-2wz+z^{2}+w^{4}z^{2}-2w^{3}z^{3}+w^{2}(1+z^{2}+z^{4})}$ left,
which we are unable to rationalise further. There are in general multiple
(complex) solutions of eq.(37) for $w,z$ in terms of $s_{ij},s_{jk}$. Our
choice of the solutions, depending on the real values of $s_{ij}$ and
$s_{jk}$, is given in appendix A.

The next step is to integrate the system of differential equations to obtain
the analytic results. Our goal is to express the analytic results in terms of
iterated integrals with polylogarithmic kernels, with explicit representation
in terms of Goncharov polylogarithms (GPLs) goncharov2011multiple  [i.e.,
multiple polylogarithms (MPLs)] wherever possible. We start with rewriting the
one-forms as logarithmic one-forms and casting the differential equations in a
more compact form. Feynman integrals expressed in terms of iterated integrals
can be efficiently described via their symbols Goncharov:2010jf . In the
variables defined in eq. (37), the alphabet $\mathbb{A}$ consists of the
following $16$ letters:

$\displaystyle\mathbb{A}=$
$\displaystyle\left\\{l_{1}(w,z),\cdots,l_{16}(w,z)\right\\}$ $\displaystyle=$
$\displaystyle\left\\{1-w,1+w,1-wz,w-z,w,w+z,1+w-z+wz,1-w+z+wz,\right.$
$\displaystyle 1+wz,1-z,1+z,z,\frac{-\rho+w-z-wz+w^{2}z-wz^{2}}{\rho+w-z-
wz+w^{2}z-wz^{2}},\frac{-\rho+w^{2}-3wz+z^{2}}{\rho+w^{2}-3wz+z^{2}},$
$\displaystyle\left.\frac{-1-\rho+w^{2}z-wz^{2}}{-1+\rho+w^{2}z-wz^{2}},\frac{1-\rho+w^{2}z-wz^{2}}{1+\rho+w^{2}z-wz^{2}}\right\\},$
(38)

where the last four depend on the square root $\rho$. In other words, the
matrix $A^{(2)}$ in the differential equation (36) can be expressed in terms
of these $16$ one-forms:

${\mathrm{d}}A_{i,j}^{(2)}=\sum_{k=1}^{16}{c_{ijk}\mathrm{dlog}(l_{k}(w,z))},\qquad
c_{ijk}\in\mathbb{Q}.$ (39)

Furthermore, this matrix can be rewritten in terms of the following $15$
independent linear combinations of $\mathrm{dlog}$ one-forms:

$\displaystyle{\mathrm{d}}g_{1}$ $\displaystyle=$
$\displaystyle-\mathrm{dlog}(l_{1})+\mathrm{dlog}(l_{2})+\mathrm{dlog}(l_{10})-\mathrm{dlog}(l_{11}),$
$\displaystyle{\mathrm{d}}g_{2}$ $\displaystyle=$
$\displaystyle\mathrm{dlog}(l_{1})+\mathrm{dlog}(l_{2})-2\mathrm{dlog}(l_{3})+\mathrm{dlog}(l_{10})+\mathrm{dlog}(l_{11}),$
$\displaystyle{\mathrm{d}}g_{3}$ $\displaystyle=$
$\displaystyle-\mathrm{dlog}(l_{1})-\mathrm{dlog}(l_{2})+2\mathrm{dlog}(l_{4})-\mathrm{dlog}(l_{10})-\mathrm{dlog}(l_{11}),$
$\displaystyle{\mathrm{d}}g_{4}$ $\displaystyle=$ $\displaystyle
2\mathrm{dlog}(l_{4})-\mathrm{dlog}(l_{5})-\mathrm{dlog}(l_{12}),$
$\displaystyle{\mathrm{d}}g_{5}$ $\displaystyle=$
$\displaystyle-\mathrm{dlog}(l_{5})+\mathrm{dlog}(l_{12}),$
$\displaystyle{\mathrm{d}}g_{6}$ $\displaystyle=$
$\displaystyle-2\mathrm{dlog}(l_{4})+4\mathrm{dlog}(l_{5})-6\mathrm{dlog}(l_{6})+4\mathrm{dlog}(l_{12}),$
$\displaystyle{\mathrm{d}}g_{7}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\mathrm{dlog}(l_{14}),$
$\displaystyle{\mathrm{d}}g_{8}$ $\displaystyle=$
$\displaystyle-\mathrm{dlog}(l_{15})+\mathrm{dlog}(l_{16}),$
$\displaystyle{\mathrm{d}}g_{9}$ $\displaystyle=$
$\displaystyle-2\mathrm{dlog}(l_{13})+\mathrm{dlog}(l_{14})+\mathrm{dlog}(l_{15})+\mathrm{dlog}(l_{16}),$
$\displaystyle{\mathrm{d}}g_{10}$ $\displaystyle=$
$\displaystyle\mathrm{dlog}(l_{4})-\mathrm{dlog}(l_{5})+\frac{1}{2}\mathrm{dlog}(l_{7})+\frac{1}{2}\mathrm{dlog}(l_{8})-\mathrm{dlog}(l_{12}),$
$\displaystyle{\mathrm{d}}g_{11}$ $\displaystyle=$ $\displaystyle
2\mathrm{dlog}(l_{13})-\mathrm{dlog}(l_{14}),$
$\displaystyle{\mathrm{d}}g_{12}$ $\displaystyle=$
$\displaystyle-\mathrm{dlog}(l_{5})-\mathrm{dlog}(l_{12}),$
$\displaystyle{\mathrm{d}}g_{13}$ $\displaystyle=$
$\displaystyle-\mathrm{dlog}(l_{7})+\mathrm{dlog}(l_{8}),$
$\displaystyle{\mathrm{d}}g_{14}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\mathrm{dlog}(l_{1})-\frac{1}{2}\mathrm{dlog}(l_{2})+\frac{1}{2}\mathrm{dlog}(l_{10})-\frac{1}{2}\mathrm{dlog}(l_{11}),$
$\displaystyle{\mathrm{d}}g_{15}$ $\displaystyle=$
$\displaystyle-2\mathrm{dlog}(l_{4})+2\mathrm{dlog}(l_{5})+2\mathrm{dlog}(l_{7})+2\mathrm{dlog}(l_{8})-6\mathrm{dlog}(l_{9})+2\mathrm{dlog}(l_{12}).$
(40)

Here we have dropped the arguments $(w,z)$ in the one-forms. The non-zero
matrix elements of $A^{(2)}$ are then compactly expressed as the following:

$\displaystyle{\mathrm{d}}A^{(2)}_{2,1}$ $\displaystyle=$
$\displaystyle{\mathrm{d}}A^{(2)}_{9,11}={\mathrm{d}}A^{(2)}_{20,11}={\mathrm{d}}A^{(2)}_{28,3}={\mathrm{d}}A^{(2)}_{28,4}={\mathrm{d}}A^{(2)}_{29,5}={\mathrm{d}}A^{(2)}_{29,28}={\mathrm{d}}g_{1},$
$\displaystyle{\mathrm{d}}A^{(2)}_{2,2}$ $\displaystyle=$
$\displaystyle{\mathrm{d}}A^{(2)}_{5,5}={\mathrm{d}}A^{(2)}_{13,13}={\mathrm{d}}g_{2},\quad{\mathrm{d}}A^{(2)}_{3,2}={\mathrm{d}}A^{(2)}_{7,1}={\mathrm{d}}A^{(2)}_{13,9}={\mathrm{d}}A^{(2)}_{27,22}={\mathrm{d}}A^{(2)}_{28,10}=-2{\mathrm{d}}g_{1},$
$\displaystyle{\mathrm{d}}A^{(2)}_{3,3}$ $\displaystyle=$ $\displaystyle
2{\mathrm{d}}g_{2},\quad{\mathrm{d}}A^{(2)}_{4,2}={\mathrm{d}}A^{(2)}_{5,3}={\mathrm{d}}A^{(2)}_{5,4}={\mathrm{d}}A^{(2)}_{6,7}={\mathrm{d}}A^{(2)}_{11,9}={\mathrm{d}}A^{(2)}_{13,3}={\mathrm{d}}A^{(2)}_{13,4}={\mathrm{d}}A^{(2)}_{28,20}=-{\mathrm{d}}g_{1},$
$\displaystyle{\mathrm{d}}A^{(2)}_{6,6}$ $\displaystyle=$
$\displaystyle{\mathrm{d}}A^{(2)}_{10,10}={\mathrm{d}}g_{3},\quad{\mathrm{d}}A^{(2)}_{7,6}=6{\mathrm{d}}g_{1},\quad{\mathrm{d}}A^{(2)}_{7,7}=3{\mathrm{d}}g_{2}-{\mathrm{d}}g_{3},\quad{\mathrm{d}}A^{(2)}_{8,6}=-\frac{3}{4}{\mathrm{d}}g_{3},$
$\displaystyle{\mathrm{d}}A^{(2)}_{8,7}$ $\displaystyle=$
$\displaystyle\frac{1}{4}{\mathrm{d}}g_{1},\quad{\mathrm{d}}A^{(2)}_{9,2}={\mathrm{d}}A^{(2)}_{12,5}={\mathrm{d}}A^{(2)}_{26,5}=-\frac{1}{2}{\mathrm{d}}g_{1},\quad{\mathrm{d}}A^{(2)}_{9,9}={\mathrm{d}}A^{(2)}_{12,8}=-{\mathrm{d}}g_{3},\quad{\mathrm{d}}A^{(2)}_{9,10}=\frac{3}{2}{\mathrm{d}}g_{3},$
$\displaystyle{\mathrm{d}}A^{(2)}_{10,11}$ $\displaystyle=$
$\displaystyle{\mathrm{d}}A^{(2)}_{13,8}={\mathrm{d}}A^{(2)}_{24,5}={\mathrm{d}}A^{(2)}_{28,8}={\mathrm{d}}A^{(2)}_{28,9}={\mathrm{d}}A^{(2)}_{28,29}=2{\mathrm{d}}g_{1},\quad{\mathrm{d}}A^{(2)}_{11,2}={\mathrm{d}}A^{(2)}_{12,10}=-\frac{1}{2}{\mathrm{d}}g_{3},$
$\displaystyle{\mathrm{d}}A^{(2)}_{11,11}$ $\displaystyle=$
$\displaystyle{\mathrm{d}}g_{2}-{\mathrm{d}}g_{3},\quad{\mathrm{d}}A^{(2)}_{12,6}=-\frac{1}{4}{\mathrm{d}}g_{3},\quad{\mathrm{d}}A^{(2)}_{12,13}={\mathrm{d}}A^{(2)}_{26,13}=\frac{1}{2}{\mathrm{d}}g_{1},\quad{\mathrm{d}}A^{(2)}_{13,6}=-\frac{3}{2}{\mathrm{d}}g_{1},$
$\displaystyle{\mathrm{d}}A^{(2)}_{13,10}$ $\displaystyle=$ $\displaystyle
3{\mathrm{d}}g_{1},\quad{\mathrm{d}}A^{(2)}_{14,14}={\mathrm{d}}g_{4},\quad{\mathrm{d}}A^{(2)}_{14,15}={\mathrm{d}}g_{5},\quad{\mathrm{d}}A^{(2)}_{15,1}=2{\mathrm{d}}g_{5},\quad{\mathrm{d}}A^{(2)}_{15,14}=-6{\mathrm{d}}g_{5},$
$\displaystyle{\mathrm{d}}A^{(2)}_{15,15}$ $\displaystyle=$
$\displaystyle{\mathrm{d}}g_{6},\quad{\mathrm{d}}A^{(2)}_{16,14}=-\frac{3}{4}{\mathrm{d}}g_{4},\quad{\mathrm{d}}A^{(2)}_{16,15}=-\frac{1}{4}{\mathrm{d}}g_{5},\quad{\mathrm{d}}A^{(2)}_{17,14}=-\frac{1}{2}{\mathrm{d}}g_{4},\quad{\mathrm{d}}A^{(2)}_{17,16}=-2{\mathrm{d}}g_{4},$
$\displaystyle{\mathrm{d}}A^{(2)}_{18,4}$ $\displaystyle=$
$\displaystyle{\mathrm{d}}A^{(2)}_{26,18}={\mathrm{d}}g_{7},\quad{\mathrm{d}}A^{(2)}_{18,9}={\mathrm{d}}A^{(2)}_{29,18}={\mathrm{d}}g_{8},\quad{\mathrm{d}}A^{(2)}_{18,10}=-\frac{3}{2}{\mathrm{d}}g_{8},\quad{\mathrm{d}}A^{(2)}_{18,11}={\mathrm{d}}g_{9},$
$\displaystyle{\mathrm{d}}A^{(2)}_{18,14}$ $\displaystyle=$
$\displaystyle-6{\mathrm{d}}g_{7}+\frac{3}{2}{\mathrm{d}}g_{8},\quad{\mathrm{d}}A^{(2)}_{18,16}=-8{\mathrm{d}}g_{7},\quad{\mathrm{d}}A^{(2)}_{18,18}={\mathrm{d}}A^{(2)}_{20,4}={\mathrm{d}}g_{10},\quad{\mathrm{d}}A^{(2)}_{18,19}={\mathrm{d}}g_{11},$
$\displaystyle{\mathrm{d}}A^{(2)}_{18,20}$ $\displaystyle=$
$\displaystyle-3{\mathrm{d}}g_{7},\quad{\mathrm{d}}A^{(2)}_{19,2}={\mathrm{d}}A^{(2)}_{25,5}={\mathrm{d}}A^{(2)}_{29,27}={\mathrm{d}}g_{12},$
$\displaystyle{\mathrm{d}}A^{(2)}_{19,4}$ $\displaystyle=$
$\displaystyle{\mathrm{d}}A^{(2)}_{20,19}={\mathrm{d}}A^{(2)}_{25,12}={\mathrm{d}}A^{(2)}_{25,17}={\mathrm{d}}A^{(2)}_{25,26}={\mathrm{d}}g_{13},$
$\displaystyle{\mathrm{d}}A^{(2)}_{19,9}$ $\displaystyle=$
$\displaystyle{\mathrm{d}}A^{(2)}_{25,16}={\mathrm{d}}A^{(2)}_{27,28}=2{\mathrm{d}}g_{13},\quad{\mathrm{d}}A^{(2)}_{19,10}={\mathrm{d}}A^{(2)}_{19,14}={\mathrm{d}}A^{(2)}_{19,20}=-\frac{3}{2}{\mathrm{d}}g_{13},$
$\displaystyle{\mathrm{d}}A^{(2)}_{19,11}$ $\displaystyle=$
$\displaystyle{\mathrm{d}}A^{(2)}_{27,24}={\mathrm{d}}A^{(2)}_{27,14}=2{\mathrm{d}}g_{12},\quad{\mathrm{d}}A^{(2)}_{19,15}={\mathrm{d}}g_{14},\quad{\mathrm{d}}A^{(2)}_{19,16}={\mathrm{d}}A^{(2)}_{27,13}=-2{\mathrm{d}}g_{13},$
$\displaystyle{\mathrm{d}}A^{(2)}_{19,18}$ $\displaystyle=$
$\displaystyle{\mathrm{d}}A^{(2)}_{25,18}=-\frac{1}{2}{\mathrm{d}}g_{11},\quad{\mathrm{d}}A^{(2)}_{19,19}={\mathrm{d}}g_{15},\quad{\mathrm{d}}A^{(2)}_{20,9}={\mathrm{d}}g_{3}-2{\mathrm{d}}g_{4}+2{\mathrm{d}}g_{10},$
$\displaystyle{\mathrm{d}}A^{(2)}_{20,10}$ $\displaystyle=$
$\displaystyle{\mathrm{d}}A^{(2)}_{20,14}={\mathrm{d}}A^{(2)}_{26,10}=-\frac{3}{2}{\mathrm{d}}g_{3}+{\mathrm{d}}g_{4}-{\mathrm{d}}g_{10},\quad{\mathrm{d}}A^{(2)}_{20,18}=-{\mathrm{d}}g_{7},$
$\displaystyle{\mathrm{d}}A^{(2)}_{20,20}$ $\displaystyle=$
$\displaystyle{\mathrm{d}}A^{(2)}_{26,20}=-{\mathrm{d}}g_{10},\quad{\mathrm{d}}A^{(2)}_{21,6}=-\frac{1}{2}{\mathrm{d}}g_{3}+\frac{1}{2}{\mathrm{d}}g_{4},\quad{\mathrm{d}}A^{(2)}_{21,8}=-2{\mathrm{d}}g_{3}+2{\mathrm{d}}g_{4},$
$\displaystyle{\mathrm{d}}A^{(2)}_{21,14}$ $\displaystyle=$
$\displaystyle\frac{1}{2}{\mathrm{d}}g_{3}-\frac{1}{2}{\mathrm{d}}g_{4},\quad{\mathrm{d}}A^{(2)}_{21,16}=2{\mathrm{d}}g_{3}-2{\mathrm{d}}g_{4},$
$\displaystyle{\mathrm{d}}A^{(2)}_{21,21}$ $\displaystyle=$
$\displaystyle{\mathrm{d}}A^{(2)}_{21,23}={\mathrm{d}}A^{(2)}_{26,21}={\mathrm{d}}A^{(2)}_{26,23}=-4{\mathrm{d}}g_{4}+4{\mathrm{d}}g_{10},\quad{\mathrm{d}}A^{(2)}_{22,6}={\mathrm{d}}A^{(2)}_{22,14}=-6{\mathrm{d}}g_{13},$
$\displaystyle{\mathrm{d}}A^{(2)}_{22,7}$ $\displaystyle=$
$\displaystyle{\mathrm{d}}A^{(2)}_{25,13}={\mathrm{d}}A^{(2)}_{27,6}=-{\mathrm{d}}g_{12},\quad{\mathrm{d}}A^{(2)}_{22,8}={\mathrm{d}}A^{(2)}_{22,16}=-8{\mathrm{d}}g_{13},\quad{\mathrm{d}}A^{(2)}_{22,15}=2{\mathrm{d}}g_{14},$
$\displaystyle{\mathrm{d}}A^{(2)}_{22,21}$ $\displaystyle=$
$\displaystyle{\mathrm{d}}A^{(2)}_{22,23}=12{\mathrm{d}}g_{13},\quad{\mathrm{d}}A^{(2)}_{22,22}=-2{\mathrm{d}}g_{4}+4{\mathrm{d}}g_{10}+{\mathrm{d}}g_{15},\quad{\mathrm{d}}A^{(2)}_{23,6}={\mathrm{d}}A^{(2)}_{23,14}=-\frac{3}{2}{\mathrm{d}}g_{4}+\frac{3}{2}{\mathrm{d}}g_{10},$
$\displaystyle{\mathrm{d}}A^{(2)}_{23,8}$ $\displaystyle=$ $\displaystyle
2{\mathrm{d}}g_{3}-4{\mathrm{d}}g_{4}+2{\mathrm{d}}g_{10},\quad{\mathrm{d}}A^{(2)}_{23,16}=-2{\mathrm{d}}g_{3}+{\mathrm{d}}g_{10},\quad{\mathrm{d}}A^{(2)}_{23,21}={\mathrm{d}}A^{(2)}_{23,23}=6{\mathrm{d}}g_{4}-6{\mathrm{d}}g_{10},$
$\displaystyle{\mathrm{d}}A^{(2)}_{23,22}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}{\mathrm{d}}g_{13},\quad{\mathrm{d}}A^{(2)}_{25,10}={\mathrm{d}}A^{(2)}_{25,14}={\mathrm{d}}A^{(2)}_{25,20}=\frac{1}{2}{\mathrm{d}}g_{13},\quad{\mathrm{d}}A^{(2)}_{25,21}={\mathrm{d}}A^{(2)}_{25,23}=-4{\mathrm{d}}g_{13},$
$\displaystyle{\mathrm{d}}A^{(2)}_{25,25}$ $\displaystyle=$
$\displaystyle\frac{2}{3}{\mathrm{d}}g_{10}+\frac{1}{3}{\mathrm{d}}g_{15},\quad{\mathrm{d}}A^{(2)}_{26,6}=-\frac{1}{4}{\mathrm{d}}g_{3}+\frac{1}{2}{\mathrm{d}}g_{4},\quad{\mathrm{d}}A^{(2)}_{26,8}=-{\mathrm{d}}g_{3}+2{\mathrm{d}}g_{4},$
$\displaystyle{\mathrm{d}}A^{(2)}_{26,14}$ $\displaystyle=$
$\displaystyle\frac{1}{2}{\mathrm{d}}g_{4}-{\mathrm{d}}g_{10},\quad{\mathrm{d}}A^{(2)}_{26,16}={\mathrm{d}}A^{(2)}_{29,17}=2{\mathrm{d}}g_{4}-4{\mathrm{d}}g_{10},\quad{\mathrm{d}}A^{(2)}_{27,3}=-2{\mathrm{d}}g_{12},$
$\displaystyle{\mathrm{d}}A^{(2)}_{27,8}$ $\displaystyle=$
$\displaystyle{\mathrm{d}}A^{(2)}_{27,10}={\mathrm{d}}A^{(2)}_{27,17}=4{\mathrm{d}}g_{12},\quad{\mathrm{d}}A^{(2)}_{27,16}=8{\mathrm{d}}g_{12},\quad{\mathrm{d}}A^{(2)}_{27,19}={\mathrm{d}}A^{(2)}_{28,26}=4{\mathrm{d}}g_{1},$
$\displaystyle{\mathrm{d}}A^{(2)}_{27,21}$ $\displaystyle=$
$\displaystyle-12{\mathrm{d}}g_{12},\quad{\mathrm{d}}A^{(2)}_{27,23}=-8{\mathrm{d}}g_{12},\quad{\mathrm{d}}A^{(2)}_{27,27}={\mathrm{d}}g_{2}+{\mathrm{d}}g_{3}+{\mathrm{d}}g_{15},\quad{\mathrm{d}}A^{(2)}_{27,29}=-4{\mathrm{d}}g_{12},$
$\displaystyle{\mathrm{d}}A^{(2)}_{28,6}$ $\displaystyle=$
$\displaystyle\frac{5}{2}{\mathrm{d}}g_{1},\quad{\mathrm{d}}A^{(2)}_{28,13}={\mathrm{d}}A^{(2)}_{29,29}=-{\mathrm{d}}g_{3}+2{\mathrm{d}}g_{10},\quad{\mathrm{d}}A^{(2)}_{28,18}=-{\mathrm{d}}g_{9},\quad{\mathrm{d}}A^{(2)}_{28,21}={\mathrm{d}}A^{(2)}_{28,23}=-4{\mathrm{d}}g_{1},$
$\displaystyle{\mathrm{d}}A^{(2)}_{28,27}$ $\displaystyle=$
$\displaystyle-{\mathrm{d}}g_{13},\quad{\mathrm{d}}A^{(2)}_{28,28}={\mathrm{d}}g_{2}-{\mathrm{d}}g_{3}+2{\mathrm{d}}g_{4}-2{\mathrm{d}}g_{10},\quad{\mathrm{d}}A^{(2)}_{29,6}=\frac{3}{2}{\mathrm{d}}g_{3}-{\mathrm{d}}g_{4},$
$\displaystyle{\mathrm{d}}A^{(2)}_{29,8}$ $\displaystyle=$ $\displaystyle
6{\mathrm{d}}g_{3}-4{\mathrm{d}}g_{4},\quad{\mathrm{d}}A^{(2)}_{29,10}=2{\mathrm{d}}g_{3}-2{\mathrm{d}}g_{10},\quad{\mathrm{d}}A^{(2)}_{29,14}={\mathrm{d}}A^{(2)}_{29,24}={\mathrm{d}}g_{4}-2{\mathrm{d}}g_{10},$
$\displaystyle{\mathrm{d}}A^{(2)}_{29,16}$ $\displaystyle=$ $\displaystyle
4{\mathrm{d}}g_{4}-8{\mathrm{d}}g_{10},\quad{\mathrm{d}}A^{(2)}_{29,20}=-{\mathrm{d}}g_{3}+2{\mathrm{d}}g_{4}-2{\mathrm{d}}g_{10},\quad{\mathrm{d}}A^{(2)}_{29,21}=-4{\mathrm{d}}g_{3}+2{\mathrm{d}}g_{4}+4{\mathrm{d}}g_{10},$
$\displaystyle{\mathrm{d}}A^{(2)}_{29,23}$ $\displaystyle=$
$\displaystyle-4{\mathrm{d}}g_{3}+4{\mathrm{d}}g_{4},\quad{\mathrm{d}}A^{(2)}_{29,26}=-2{\mathrm{d}}g_{3}+4{\mathrm{d}}g_{4}-4{\mathrm{d}}g_{10}.$
(41)

In the next step, we take the series expansion over the dimensional regulator
$\epsilon$ in the master integrals
$\vec{f}^{(2)}=\sum_{\omega=0}{\epsilon^{\omega}\vec{f}^{(2,\omega)}}$ and the
differential equation system eq.(36) becomes

$\displaystyle{\mathrm{d}}\vec{f}^{(2,\omega+1)}$ $\displaystyle=$
$\displaystyle{\mathrm{d}}A^{(2)}\vec{f}^{(2,\omega)},\quad\omega\in\mathbb{N}.$
(42)

Since our choice of $\vec{f}^{(2)}$ is a UT basis, all terms contributing to
$\vec{f}^{(2,\omega)}$ should have the same transcendental weight $\omega$.
With a proper boundary condition, we can solve eq.(42) iteratively from lower
to higher weights. The solution can be written in terms of iterated integrals
with logarithmic one-forms. Up to weight $\omega=4$ that we need, most of the
master integrals $f_{n}^{(2,\omega)}(x_{ij},x_{jk},x_{ik})$ can be expressed
through Goncharov polylogarithm $G$’s, while only the iterated integrals for
the following master integrals depend on the kernels with the square root
$\rho$:

$f_{18}^{(2,3)},f_{18}^{(2,4)},f_{19}^{(2,4)},f_{20}^{(2,4)},f_{25}^{(2,4)},f_{26}^{(2,4)},f_{28}^{(2,4)},f_{29}^{(2,4)}.$
(43)

In our two-loop amplitudes, we actually only need 5 integrals
$f_{18}^{(2,3)},f_{25}^{(2,4)},f_{26}^{(2,4)},f_{28}^{(2,4)}$, and
$f_{29}^{(2,4)}$ which cannot be written in $G$’s. For the integrals with
$G$’s, 666Their concrete expressions are available in the ancillary files.
their numerical evaluations are rather straightforward with the algorithm
outlined in ref. Vollinga:2004sn . In this paper, we use FastGPL Wang:2021imw
and take HandyG Naterop:2019xaf  as a cross check.

Now, as a new ingredient from this section, we outline the steps taken to
prepare all the master integrals present in the scattering amplitudes for
phenomenological applications. We explain how to numerically evaluate iterated
integrals with algebraic kernels using an idea that has been recently
successfully utilised for other processes (see, e.g., ref. Chicherin:2020oor
). We also find different analytic expressions for all the master integrals
valid in different regions of phase space in order to avoid performing
analytic continuations. Moreover, we identify relations among the Feynman
integrals in the scattering amplitudes, which reduces the number of iterated
integrals to be computed numerically at each point of the phase space. For the
numerical evaluations of our master integrals with square roots, we choose a
technique elucidated in refs. Caron-Huot:2014lda ; Chicherin:2020oor . The
kernels in the first two-fold integration can be rationalised and therefore
can be expressed through logarithms and classical polylogarithms
$\mathrm{Li}_{n}$ in terms of $x_{ij},x_{jk},x_{ik}$. 777We choose the
classical polylogarithms rather than the GPLs here because the numerical
evaluations of the former are much faster than the latter. This is possible
thanks to the methods outlined in ref. Duhr:2011zq  by matching symbols in
both the $(w,z)$ as well as $(s_{ij},s_{jk})$ coordinate systems. However, the
Riemann sheets for the multi-valued functions should be carefully picked up.
For the weight-$3$ integral $f_{18}^{(2,3)}$, we are only left with a single
integration to be carried out numerically, while there are two-fold
integrations left for the weight-$4$ integrals. For the latter, the basic idea
is to convert the last two-fold integration into a one-fold integration, and
to perform the last one-fold integration numerically. For numerical
integration in all phase-space regions in appendix A, we also use an approach
that makes the numerical evaluation of all the master integrals faster as we
avoid crossing any physical singularity. Instead of integrating our results
within a single region and then relying on analytic continuation to cross the
boundaries into other regions, we obtain distinct analytic results that are
specifically valid within each of these regions. Hence, our numerical
evaluations of the integrals are valid within specific phase-space regions. We
also obtain analytic boundary constants $f^{(2,\omega)}_{n}(x_{ij},0,-x_{ij})$
and $f^{(2,\omega)}_{n}(0,x_{jk},-x_{jk})$ valid in all these regions. For
these analytic boundary constants for the $5$ iterated integrals with roots,
refer to appendix B. For the last one-fold integration, the convergence of the
numerical integration strongly depends on the path. We opt for the following
path parameterisation from the starting point $(s_{ij,0},s_{jk,0})$ to the end
point $(s_{ij,1},s_{jk,1})$

$\displaystyle s_{ij}(r)$ $\displaystyle=$ $\displaystyle
s_{ij,0}(1-r)^{5}\left(1+5r+15r^{2}+35r^{3}+70r^{4}\right)$
$\displaystyle+s_{ij,1}r^{5}\left(126-420r+540r^{2}-315r^{3}+70r^{4}\right),$
$\displaystyle s_{jk}(r)$ $\displaystyle=$ $\displaystyle
s_{jk,0}(1-r)^{5}\left(1+5r+15r^{2}+35r^{3}+70r^{4}\right)$ (44)
$\displaystyle+s_{jk,1}r^{5}\left(126-420r+540r^{2}-315r^{3}+70r^{4}\right),$

where $r\in[0,1]$. It is easy to verify that
$(s_{ij}(0),s_{jk}(0))=(s_{ij,0},s_{jk,0})$ and
$(s_{ij}(1),s_{jk}(1))=(s_{ij,1},s_{jk,1})$. In this article, we choose the
following boundary conditions:

$\displaystyle\left(s_{ij,0},s_{jk,0}\right)$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}(0,0),&\quad{\rm
if}~{}s_{ij},s_{jk}<0\\\ (s_{ij},0),&\quad{\rm if}~{}s_{ij}>0,s_{jk}<0\\\
(0,s_{jk}),&\quad{\rm if}~{}s_{ij}<0,s_{jk}>0.\end{array}\right.\,$ (48)

The integrals have been cross checked numerically for a few phase space points
in each region against AMFlow Liu:2022chg .

As the final step in making the integrals pheno-ready and in order to simplify
the final helicity amplitudes, it is necessary to find all possible relations
among the master integrals needed for our amplitudes. Otherwise, not only can
we not check pole cancellations at the analytic level, as a similar finding in
ref. Badger:2021owl  for another process, but large numerical cancellations
will also happen among different pieces, which prevents us from taking any
interesting kinematic limit. In this work, we use the symbol techniques
Duhr:2019tlz  and the shuffle algebra that the iterated integrals satisfy to
find the following relations, by including the one-loop master integrals
$f_{n}^{(1,\omega)}$ defined in appendix C. The relations can be categorised
into three parts:

  * •

The following integrals are identical over permutations of external momenta:

$\displaystyle f_{n}^{(1,\omega)}(x_{ij},x_{jk},x_{ik})=$
$\displaystyle~{}f_{n}^{(1,\omega)}(x_{ik},x_{jk},x_{ij})\quad\text{for}~{}n\in\\{3,5\\}\,,$
$\displaystyle f_{n}^{(1,\omega)}(x_{ij},x_{jk},x_{ik})=$
$\displaystyle~{}f_{n}^{(1,\omega)}(x_{ij},x_{ik},x_{jk})\quad\text{for}~{}n\in\\{2,4\\}\,,$
$\displaystyle f_{6}^{(1,\omega)}(x_{ij},x_{jk},x_{ik})=$
$\displaystyle~{}f_{6}^{(1,\omega)}(x_{jk},x_{ij},x_{ik})\,,$ $\displaystyle
f_{n}^{(2,\omega)}(x_{ij},x_{jk},x_{ik})=$
$\displaystyle~{}f_{n}^{(2,\omega)}(x_{ij},x_{ik},x_{jk})\quad\text{for}~{}n\in\\{2,3,4,5,6,7,8,9,10,11,12,13,24\\}\,,$
$\displaystyle f_{n}^{(2,\omega)}(x_{ij},x_{jk},x_{ik})=$
$\displaystyle~{}f_{n}^{(2,\omega)}(x_{ik},x_{jk},x_{ij})\quad\text{for}~{}n\in\\{14,15,16,17\\}\,,$
$\displaystyle f_{19}^{(2,2)}(x_{ij},x_{jk},x_{ik})=$
$\displaystyle~{}f_{19}^{(2,2)}(x_{jk},x_{ij},x_{ik})\,,$ $\displaystyle
f_{n}^{(2,\omega)}(x_{ij},x_{jk},x_{ik})=$
$\displaystyle~{}f_{n}^{(2,\omega)}(x_{jk},x_{ij},x_{ik})\quad\text{for}~{}n\in\\{21,22,23\\}\,.$
(49)

  * •

Some integrals with same transcendental weight can be related:

$\displaystyle f_{2}^{(1,\omega)}(x_{ij},x_{jk},x_{ik})=$
$\displaystyle~{}f_{3}^{(1,\omega)}(x_{jk},x_{ij},x_{ik}),\quad
f_{4}^{(1,\omega)}(x_{ij},x_{jk},x_{ik})=~{}f_{5}^{(1,\omega)}(x_{jk},x_{ij},x_{ik})\,,$
$\displaystyle f_{6}^{(2,\omega)}(x_{ij},x_{jk},x_{ik})=$
$\displaystyle~{}f_{14}^{(2,\omega)}(x_{jk},x_{ij},x_{ik}),\quad
f_{7}^{(2,\omega)}(x_{ij},x_{jk},x_{ik})=~{}f_{15}^{(2,\omega)}(x_{jk},x_{ij},x_{ik})\,,$
$\displaystyle f_{8}^{(2,\omega)}(x_{ij},x_{jk},x_{ik})=$
$\displaystyle~{}f_{16}^{(2,\omega)}(x_{jk},x_{ij},x_{ik}),\quad
f_{9}^{(2,2)}(x_{ij},x_{jk},x_{ik})=f_{16}^{(2,2)}(x_{jk},x_{ij},x_{ik})\,,$
$\displaystyle f_{15}^{(2,1)}(x_{ij},x_{jk},x_{ik})=$
$\displaystyle-2f_{2}^{(2,1)}(x_{jk},x_{ij},x_{ik}),\quad
f_{3}^{(2,2)}(x_{ij},x_{jk},x_{ik})=4f_{16}^{(2,2)}(x_{jk},x_{ij},x_{ik})\,,$
$\displaystyle f_{4}^{(2,2)}(x_{ij},x_{jk},x_{ik})=$
$\displaystyle~{}2f_{16}^{(2,2)}(x_{jk},x_{ij},x_{ik}),\quad
f_{14}^{(2,2)}(x_{ij},x_{jk},x_{ik})=-4f_{16}^{(2,2)}(x_{ij},x_{jk},x_{ik})\,,$
$\displaystyle f_{22}^{(2,2)}(x_{ij},x_{jk},x_{ik})=$
$\displaystyle~{}2f_{19}^{(2,2)}(x_{ij},x_{jk},x_{ik}),\quad
f_{4}^{(2,3)}(x_{ij},x_{jk},x_{ik})=-f_{5}^{(1,3)}(x_{jk},x_{ij},x_{ik})\,,$
$\displaystyle f_{2}^{(2,\omega)}(x_{ij},x_{jk},x_{ik})=$ $\displaystyle-
f_{2}^{(1,\omega)}(x_{ij},x_{jk},x_{ik})\qquad\text{for}~{}\omega\in\\{1,2\\}\,,$
$\displaystyle f_{4}^{(2,\omega)}(x_{ij},x_{jk},x_{ik})=$ $\displaystyle-
f_{4}^{(1,\omega)}(x_{ij},x_{jk},x_{ik})\qquad\text{for}~{}\omega\in\\{2,3\\}\,,$
$\displaystyle f_{11}^{(2,2)}(x_{ij},x_{jk},x_{ik})=$
$\displaystyle~{}\frac{3}{2}f_{3}^{(1,2)}(x_{jk},x_{ij},x_{ik})-\frac{1}{4}f_{15}^{(2,2)}(x_{jk},x_{ij},x_{ik})\,,$
$\displaystyle f_{8}^{(2,3)}(x_{ij},x_{jk},x_{ik})=$
$\displaystyle-\frac{3}{4}f_{6}^{(2,3)}(x_{ij},x_{jk},x_{ik})+f_{10}^{(2,3)}(x_{ij},x_{jk},x_{ik})-3f_{4}^{(1,3)}(x_{ij},x_{jk},x_{ik})\,,$
$\displaystyle f_{9}^{(2,3)}(x_{ij},x_{jk},x_{ik})=$
$\displaystyle~{}\frac{1}{4}f_{14}^{(2,3)}(x_{jk},x_{ij},x_{ik})-2f_{5}^{(1,3)}(x_{jk},x_{ij},x_{ik})\,.$
(50)

  * •

Some higher transcendental weight two-loop master integrals can be expressed
into lower weight integrals:

$\displaystyle f_{5}^{(1,2)}(x_{ij},x_{jk},x_{ik})=$
$\displaystyle~{}\frac{1}{2}\left(f_{3}^{(1,1)}(x_{ij},x_{jk},x_{ik})\right)^{2},$
$\displaystyle f_{24}^{(2,4)}(x_{ij},x_{jk},x_{ik})=$
$\displaystyle~{}\frac{1}{4}\left(f_{3}^{(1,1)}(x_{jk},x_{ij},x_{ik})\right)^{4}\,,$
$\displaystyle f_{16}^{(2,2)}(x_{ij},x_{jk},x_{ik})=$
$\displaystyle-\frac{1}{4}\left(f_{3}^{(1,1)}(x_{ij},x_{jk},x_{ik})\right)^{2}\,,$
$\displaystyle f_{5}^{(2,3)}(x_{ij},x_{jk},x_{ik})=$
$\displaystyle-\frac{1}{2}\left(f_{3}^{(1,1)}(x_{jk},x_{ij},x_{ik})\right)^{3},$
$\displaystyle f_{3}^{(2,3)}(x_{ij},x_{jk},x_{ik})=$
$\displaystyle-2f_{3}^{(1,1)}(x_{jk},x_{ij},x_{ik})f_{3}^{(1,2)}(x_{jk},x_{ij},x_{ik})\,,$
$\displaystyle f_{5}^{(2,4)}(x_{ij},x_{jk},x_{ik})=$ $\displaystyle-
f_{2}^{(1,2)}(x_{ij},x_{jk},x_{ik})f_{4}^{(1,2)}(x_{ij},x_{jk},x_{ik})\;-$
$\displaystyle\hskip
128.0374ptf_{2}^{(1,1)}(x_{ij},x_{jk},x_{ik})f_{4}^{(1,3)}(x_{ij},x_{jk},x_{ik})\,,$
$\displaystyle f_{10}^{(2,3)}(x_{ij},x_{jk},x_{ik})=$
$\displaystyle~{}\frac{1}{3}f_{14}^{(2,3)}(x_{jk},x_{ij},x_{ik})-2f_{5}^{(1,3)}(x_{jk},x_{ij},x_{ik})\;+$
$\displaystyle\frac{1}{6}f_{3}^{(1,1)}(x_{jk},x_{ij},x_{ik})\left(-6f_{3}^{(1,2)}(x_{jk},x_{ij},x_{ik})+f_{15}^{(2,2)}(x_{jk},x_{ij},x_{ik})\right)\,,$
$\displaystyle f_{16}^{(2,3)}(x_{ij},x_{jk},x_{ik})=$
$\displaystyle~{}f_{5}^{(1,3)}(x_{ij},x_{jk},x_{ik})-\frac{5}{12}f_{14}^{(2,3)}(x_{ij},x_{jk},x_{ik})\;+$
$\displaystyle\frac{1}{6}f_{3}^{(1,1)}(x_{ij},x_{jk},x_{ik})\left(-6f_{3}^{(1,2)}(x_{ij},x_{jk},x_{ik})+f_{15}^{(2,2)}(x_{ij},x_{jk},x_{ik})\right)\,,$
$\displaystyle f_{23}^{(2,3)}(x_{ij},x_{jk},x_{ik})=$
$\displaystyle-\frac{1}{4}f_{10}^{(2,3)}(x_{ij},x_{jk},x_{ik})-\frac{1}{2}f_{10}^{(2,3)}(x_{jk},x_{ij},x_{ik})-f_{20}^{(2,3)}(x_{ij},x_{jk},x_{ik})$
$\displaystyle+\frac{3}{2}f_{11}^{(2,2)}(x_{ij},x_{jk},x_{ik})f_{2}^{(2,1)}(x_{ij},x_{jk},x_{ik})$
$\displaystyle+f_{11}^{(2,2)}(x_{jk},x_{ij},x_{ik})f_{2}^{(2,1)}(x_{jk},x_{ij},x_{ik})\,.$
(51)

We remind the readers that
$x_{ij}=s_{ij}/m_{f}^{2},\;x_{jk}=s_{jk}/m_{f}^{2},$ and
$x_{ik}=s_{ik}/m_{f}^{2}$. Using the above relations, we manage to cancel the
poles analytically. Furthermore, the relations also considerably reduce the
number of independent integrals at a transcendental weight from around 300 to
84 in the two-loop helicity amplitudes.

##  4 Analytic expressions of helicity amplitudes

In this section, we present our main results, i.e., the analytic expressions
of helicity amplitudes for LbL. To obtain the helicity amplitudes we use Qgraf
and FeynArts to generate the Feynman diagrams, and further manipulations such
as performing $d$-dimensional Lorentz and Dirac algebra are done in FORM
Vermaseren:2000nd ; Ruijl:2017dtg  or Mathematica.

We first present the one-loop form factors and helicity amplitudes as a warm-
up exercise. In this section, we will only consider a given fermion loop with
the fermion mass $m_{f}$, while the one-loop expressions for the $W^{\pm}$
loop can be found in appendix D. Throughout this section, let us define the
following dimensionless variables

$x_{s}=\frac{s}{m_{f}^{2}},\quad x_{t}=\frac{t}{m_{f}^{2}},\quad
x_{u}=\frac{u}{m_{f}^{2}}$ (52)

with $x_{s}+x_{t}+x_{u}=0$. Let us assume the fermion has charge $Q_{f}$ and
colour charge $N_{c,f}$. 888In the SM, the quarks have $N_{c,f}=3$ and the
leptons have $N_{c,f}=1$. The up-, down-type quarks have $Q_{f}=2/3$ and
$Q_{f}=-1/3$, and the charged leptons have $Q_{f}=-1$. In addition, we also
denote the LO (one-loop) helicity amplitudes with the fermion $f$ as
$\mathcal{M}_{\lambda_{1}\lambda_{2}\lambda_{3}\lambda_{4}}^{(0,0,f)}$, while
the two-loop QCD and QED amplitudes with the fermion are
$\mathcal{M}_{\lambda_{1}\lambda_{2}\lambda_{3}\lambda_{4}}^{(1,0,f)}$ and
$\mathcal{M}_{\lambda_{1}\lambda_{2}\lambda_{3}\lambda_{4}}^{(0,1,f)}$
respectively.

###  4.1 One-loop results

The one-loop form factors are

$\displaystyle-
iA_{S}(s,t,u)=4N_{c,f}Q_{f}^{4}\alpha^{2}\left\\{2+\mathop{\sum_{(i,j,k)=(s,t,u),}}_{(t,u,s),(u,s,t)}{\left[r^{(1)}_{1}(x_{i},x_{j},x_{k})\frac{f_{2}^{(1,1)}(x_{i},x_{j},x_{k})}{\sqrt{x_{i}(x_{i}-4)}}\right.}\right.$
$\displaystyle\qquad\qquad\qquad\left.\left.-r^{(1)}_{2}(x_{i},x_{j},x_{k})f_{4}^{(1,2)}(x_{i},x_{j},x_{k})-r^{(1)}_{3}(x_{i},x_{j},x_{k})\frac{f_{6}^{(1,2)}(x_{i},x_{j},x_{k})}{\sqrt{x_{i}x_{j}(x_{i}x_{j}+4x_{k})}}\right]\right\\},$
$\displaystyle-i\Delta\hat{B}^{1}_{11}(s,t,u)=\frac{4N_{c,f}Q_{f}^{4}\alpha^{2}}{u}\left\\{2-r^{(1)}_{1}(x_{s},x_{t},x_{u})\frac{f_{2}^{(1,1)}(x_{s},x_{t},x_{u})}{\sqrt{x_{s}(x_{s}-4)}}+r^{(1)}_{4}(x_{t},x_{u},x_{s})\frac{f_{2}^{(1,1)}(x_{t},x_{u},x_{s})}{\sqrt{x_{t}(x_{t}-4)}}\right.$
$\displaystyle\qquad\qquad\qquad-r^{(1)}_{4}(x_{u},x_{s},x_{t})\frac{f_{2}^{(1,1)}(x_{u},x_{s},x_{t})}{\sqrt{x_{u}(x_{u}-4)}}+\left(r^{(1)}_{5}(x_{s},x_{t},x_{u})+r^{(1)}_{6}(x_{s},x_{t},x_{u})\right)f_{4}^{(1,2)}(x_{s},x_{t},x_{u})$
$\displaystyle\qquad\qquad\qquad+r^{(1)}_{5}(x_{t},x_{u},x_{s})f_{4}^{(1,2)}(x_{t},x_{u},x_{s})-r^{(1)}_{5}(x_{u},x_{s},x_{t})f_{4}^{(1,2)}(x_{u},x_{s},x_{t})$
$\displaystyle\qquad\qquad\qquad-r^{(1)}_{7}(x_{s},x_{t},x_{u})\frac{f_{6}^{(1,2)}(x_{s},x_{t},x_{u})}{\sqrt{x_{s}x_{t}(x_{s}x_{t}+4x_{u})}}+\left(r^{(1)}_{7}(x_{t},x_{u},x_{s})-4x_{u}\right)\frac{f_{6}^{(1,2)}(x_{t},x_{u},x_{s})}{\sqrt{x_{t}x_{u}(x_{t}x_{u}+4x_{s})}}$
$\displaystyle\qquad\qquad\qquad\left.-\left(r^{(1)}_{7}(x_{u},x_{s},x_{t})+4(x_{u}-x_{s})\right)\frac{f_{6}^{(1,2)}(x_{u},x_{s},x_{t})}{\sqrt{x_{u}x_{s}(x_{u}x_{s}+4x_{t})}}\right\\},$
$\displaystyle-i\Delta
C_{2111}(s,t,u)=\frac{32N_{c,f}Q_{f}^{4}\alpha^{2}}{su}\left\\{1+2\left(x_{s}^{-1}+x_{t}^{-1}+x_{u}^{-1}\right)\mathop{\sum_{(i,j,k)=(s,t,u),}}_{(t,u,s),(u,s,t)}{f_{4}^{(1,2)}(x_{i},x_{j},x_{k})}\right.$
$\displaystyle\qquad\qquad\qquad\left.-\mathop{\sum_{(i,j,k)=(s,t,u),}}_{(t,u,s),(u,s,t)}{r^{(1)}_{8}(x_{i},x_{j},x_{k})\frac{f_{6}^{(1,2)}(x_{i},x_{j},x_{k})}{\sqrt{x_{i}x_{j}(x_{i}x_{j}+4x_{k})}}}\right\\},$
(53)

where the basis of rational coefficients is as follows:

$\displaystyle r^{(1)}_{1}(x_{i},x_{j},x_{k})$ $\displaystyle=$ $\displaystyle
2(x_{i}-4)\left(\frac{x_{i}}{x_{k}}-\frac{x_{k}}{x_{j}}\right)$ $\displaystyle
r^{(1)}_{2}(x_{i},x_{j},x_{k})$ $\displaystyle=$
$\displaystyle-2\left(\frac{4-x_{i}}{x_{i}}+\frac{x_{i}(x_{j}^{3}-6x_{j}x_{k}+x_{k}^{3})}{x_{j}^{2}x_{k}^{2}}\right),$
$\displaystyle r^{(1)}_{3}(x_{i},x_{j},x_{k})$ $\displaystyle=$
$\displaystyle\frac{32x_{k}^{2}-x_{i}x_{j}(x_{i}^{2}+x_{j}^{2}-16x_{k})}{x_{k}^{2}},$
$\displaystyle r^{(1)}_{4}(x_{i},x_{j},x_{k})$ $\displaystyle=$
$\displaystyle\frac{2x_{i}(x_{i}-4)(x_{j}-x_{k})}{x_{j}x_{k}},$ $\displaystyle
r^{(1)}_{5}(x_{i},x_{j},x_{k})$ $\displaystyle=$
$\displaystyle\frac{2(2x_{j}x_{k}(x_{j}-x_{k})+x_{i}(x_{j}^{3}-x_{k}^{3}))}{x_{j}^{2}x_{k}^{2}},$
$\displaystyle r^{(1)}_{6}(x_{i},x_{j},x_{k})$ $\displaystyle=$
$\displaystyle\frac{2(x_{i}^{2}+x_{j}^{2}-4x_{k})}{x_{k}^{2}},$ $\displaystyle
r^{(1)}_{7}(x_{i},x_{j},x_{k})$ $\displaystyle=$
$\displaystyle\frac{x_{i}(4(x_{j}^{2}-x_{i}^{2})+x_{j}(x_{i}^{2}+x_{j}^{2})}{x_{k}^{2}},$
$\displaystyle r^{(1)}_{8}(x_{i},x_{j},x_{k})$ $\displaystyle=$
$\displaystyle\frac{2(x_{i}x_{j}+2x_{k})}{x_{k}}.$ (54)

The one-loop master integrals, defined in appendix C, can be expressed through
GPLs:

$\displaystyle f_{2}^{(1,1)}(x_{s},x_{t},x_{u})$ $\displaystyle=$
$\displaystyle\left(-1\right)^{\theta(x_{s}-4)}\left[-G(-1;w)+G(-1;z)+G(1;w)-G(1;z)\right],$
$\displaystyle f_{4}^{(1,2)}(x_{s},x_{t},x_{u})$ $\displaystyle=$
$\displaystyle\left(G(1;w)-G(-1;w)\right)\left(G(-1;z)-G(1;z)\right)$
$\displaystyle+\sum_{i,j=0}^{1}{(-1)^{i+j}\left[G((-1)^{i},(-1)^{j};w)+G((-1)^{i},(-1)^{j};z)\right]},$
$\displaystyle f_{6}^{(1,2)}(x_{s},x_{t},x_{u})$ $\displaystyle=$
$\displaystyle\left(-1\right)^{\theta(x_{s})+\theta(x_{t})}\left[G(0;w)\left[G(1;z)-G(-1;z)\right]+G(0;z)\left[G(-1;w)-G(1;w)\right]\right.$
(55) $\displaystyle+~{}G(1,0;w)-G(-1,0;w)+G(0,-1;w)-G(0,1;w)$
$\displaystyle\left.-~{}G(1,0;z)+G(-1,0;z)-G(0,-1;z)+G(0,1;z)\right].$

The prefactors $\left(-1\right)^{\theta(x_{s}-4)}$ and
$\left(-1\right)^{\theta(x_{s})+\theta(x_{t})}$ ensure the correct analytic
continuation when using the convention in section A, so that the multi-valued
square-root function in eq.(53) takes the values in the first Riemann sheet.

Alternatively, the one-loop helicity amplitudes with a given massive fermion
$f$ loop
$\mathcal{M}_{\lambda_{1}\lambda_{2}\lambda_{3}\lambda_{4}}^{(0,0,f)}$ are

$\displaystyle-i\mathcal{M}^{(0,0,f)}_{++++}$ $\displaystyle=$ $\displaystyle
8N_{c,f}Q_{f}^{4}\alpha^{2}\left[1-4\mathop{\sum_{(i,j,k)=(s,t,u),}}_{(t,u,s),(u,s,t)}{\frac{f_{6}^{(1,2)}(x_{i},x_{j},x_{k})}{\sqrt{x_{i}x_{j}(x_{i}x_{j}+4x_{k})}}}\right],$
$\displaystyle-i\mathcal{M}^{(0,0,f)}_{-+++}$ $\displaystyle=$ $\displaystyle
8N_{c,f}Q_{f}^{4}\alpha^{2}\left[1+2\left(x_{s}^{-1}+x_{t}^{-1}+x_{u}^{-1}\right)\mathop{\sum_{(i,j,k)=(s,t,u),}}_{(t,u,s),(u,s,t)}{f_{4}^{(1,2)}(x_{i},x_{j},x_{k})}\right.$
$\displaystyle\qquad\qquad\qquad\left.-\mathop{\sum_{(i,j,k)=(s,t,u),}}_{(t,u,s),(u,s,t)}{r^{(1)}_{8}(x_{i},x_{j},x_{k})\frac{f_{6}^{(1,2)}(x_{i},x_{j},x_{k})}{\sqrt{x_{i}x_{j}(x_{i}x_{j}+4x_{k})}}}\right],$
$\displaystyle-i\mathcal{M}^{(0,0,f)}_{--++}$ $\displaystyle=$ $\displaystyle
8N_{c,f}Q_{f}^{4}\alpha^{2}\left[-1+r_{9}^{(1)}(x_{t},x_{u},x_{s})\frac{f_{2}^{(1,1)}(x_{t},x_{u},x_{s})}{\sqrt{x_{t}(x_{t}-4)}}+r_{9}^{(1)}(x_{u},x_{t},x_{s})\frac{f_{2}^{(1,1)}(x_{u},x_{t},x_{s})}{\sqrt{x_{u}(x_{u}-4)}}\right.$
(56) $\displaystyle-
r_{10}^{(1)}(x_{s},x_{t},x_{u})\left(f_{4}^{(1,2)}(x_{t},x_{u},x_{s})+f_{4}^{(1,2)}(x_{u},x_{s},x_{t})\right)+2(x_{s}-2)\frac{f_{6}^{(1,2)}(x_{s},x_{t},x_{u})}{\sqrt{x_{s}x_{t}(x_{s}x_{t}+4x_{u})}}$
$\displaystyle\left.+2(x_{s}-2)\frac{f_{6}^{(1,2)}(x_{u},x_{s},x_{t})}{\sqrt{x_{s}x_{u}(x_{s}x_{u}+4x_{t})}}-r_{11}^{(1)}(x_{s},x_{t},x_{u})\frac{f_{6}^{(1,2)}(x_{t},x_{u},x_{s})}{\sqrt{x_{t}x_{u}(x_{t}x_{u}+4x_{s})}}\right],$

where

$\displaystyle r_{9}^{(1)}(x_{i},x_{j},x_{k})$ $\displaystyle=$
$\displaystyle\frac{(x_{i}-4)(x_{i}-x_{j})}{x_{k}},$ $\displaystyle
r_{10}^{(1)}(x_{i},x_{j},x_{k})$ $\displaystyle=$ $\displaystyle
1-\frac{4}{x_{i}}-\frac{2x_{j}x_{k}}{x_{i}^{2}},$ $\displaystyle
r_{11}^{(1)}(x_{i},x_{j},x_{k})$ $\displaystyle=$ $\displaystyle
4-2x_{i}-x_{j}x_{k}+\frac{2x_{j}x_{k}(x_{j}x_{k}+4x_{i})}{x_{i}^{2}}.$ (57)

It is interesting to note that, while the conversion between the helicity
amplitudes and the form factors is merely a linear transformation (see
eq.(LABEL:eq:ampinff)), the expressions for the amplitudes are shorter than
those for the form factors. Such an observation motivates us to simplify at
the level of the helicity amplitude rather than at the form factor. Note that
the one-loop amplitudes have been given in ref. Bardin:2009gq  in terms of
one-loop scalar integrals, rather than as presented here in terms of a UT
basis. The UT basis has the advantage of rendering shorter expressions. Our
one-loop results have been numerically checked against MadLoop Hirschi:2011pa
; Alwall:2014hca .

###  4.2 Two-loop helicity amplitudes

We are now ready to present the UV-renormalised two-loop helicity amplitudes
999At this order, we only need to renormalise the internal fermion mass
$m_{f}$. In this article, we choose the on-shell renormalisation condition for
the mass. with the full fermion mass dependence. Following the IBP reduction,
we obtain large rational coefficients in front of the master integrals. We put
dedicated efforts to simplify these coefficients into more compact and well-
organised forms as given below. This is achieved with the help of FiniteFlow
Peraro:2019svx  and MultivariateApart Heller:2021qkz . Also, the relations
between the master integrals at different weights as given in eqs.(49),(50)
and (51) were one of the crucial ingredients for simplifying the final
expressions. The size of the amplitude file is reduced from $\sim 300$ MB just
after IBP to around $180$ KB counting with Mathematica.

After simplification, the two-loop QCD helicity amplitudes are given as
follows, with $C_{F,f}=(N_{c,f}^{2}-1)/(2N_{c,f})$. The QED counterparts are
obtained with a simple rescaling
$\mathcal{M}^{(0,1,f)}_{\lambda_{1}\lambda_{2}\lambda_{3}\lambda_{4}}/\mathcal{M}^{(1,0,f)}_{\lambda_{1}\lambda_{2}\lambda_{3}\lambda_{4}}=\alpha
Q_{f}^{2}/(\alpha_{s}C_{F,f})$. It is implicitly understood that the two-loop
QCD amplitudes only apply to quarks, while the QED amplitudes are non-zero for
all charged fermions, i.e., quarks and charged leptons in the SM.

The all-plus and one-minus amplitudes are

$\displaystyle i\mathcal{M}^{(1,0,f)}_{++++}$
$\displaystyle=4N_{c,f}Q_{f}^{4}\alpha^{2}\frac{\alpha_{s}}{\pi}C_{F,f}\left\\{3-4\mathop{\sum_{(i,j,k)=(s,t,u),}}_{(t,u,s),(u,s,t)}\left[\left(\frac{2}{x_{i}}+{2\over
x_{k}}-{1\over x_{j}}\right)f_{5}^{(1,2)}(x_{i},x_{j},x_{k})\right.\right.$
$\displaystyle-$ $\displaystyle{1\over
x_{i}^{2}}f_{24}^{(2,4)}(x_{i},x_{j},x_{k})+\frac{2}{\sqrt{x_{i}x_{j}(x_{i}x_{j}+4x_{k})}}\left(\frac{x_{i}x_{j}+4x_{k}}{x_{k}}f_{19}^{(2,2)}(x_{i},x_{j},x_{k})+f_{22}^{(2,3)}(x_{i},x_{j},x_{k})\right)$
$\displaystyle+$
$\displaystyle\frac{1}{2x_{i}x_{j}}\left(6f_{3}^{(2,3)}(x_{i},x_{j},x_{k})-f_{14}^{(2,3)}(x_{k},x_{i},x_{j})-{4}f_{16}^{(2,3)}(x_{k},x_{i},x_{j})-{8}f_{20}^{(2,3)}(x_{i},x_{j},x_{k})\right.$
$\displaystyle\left.\left.+2f_{15}^{(2,2)}(x_{k},x_{i},x_{j})f_{3}^{(1,1)}(x_{k},x_{i},x_{j})\right)\right]$
$\displaystyle+\left.\mathop{\sum_{(i,j,k)=(s,t,u),(s,u,t),(t,s,u),}}_{(t,u,s),(u,s,t),(u,t,s)}\left[\frac{8~{}f_{19}^{(2,3)}(x_{i},x_{j},x_{k})}{\sqrt{x_{i}x_{j}(x_{i}x_{j}+4x_{k})}}-\frac{2(x_{i}-2)f_{27}^{(2,4)}(x_{i},x_{j},x_{k})}{\sqrt{x_{i}(x_{i}-4)}\sqrt{x_{i}x_{j}(x_{i}x_{j}+4x_{k})}}\right]\right\\}\,,$
(58)

and

$\displaystyle i\mathcal{M}^{(1,0,f)}_{-+++}$
$\displaystyle=N_{c,f}Q_{f}^{4}\alpha^{2}\frac{\alpha_{s}}{\pi}C_{F,f}\left\\{\mathop{\sum_{(i,j,k)=(s,t,u),}}_{(t,u,s),(u,s,t)}\left[\frac{1}{\sqrt{x_{j}(x_{j}-4)}}\left(2r_{1}^{(1)}(x_{j},x_{k},x_{i})f_{3}^{(1,1)}(x_{i},x_{j},x_{k})\right.\right.\right.$
$\displaystyle\left.+r_{1}^{(2)}(x_{i},x_{j},x_{k})\Big{(}f_{3}^{(1,2)}(x_{i},x_{j},x_{k})-\frac{1}{6}f_{15}^{(2,2)}(x_{i},x_{j},x_{k})\Big{)}\right)-\frac{2x_{i}x_{j}f_{6}^{(1,2)}(x_{i},x_{j},x_{k})}{\sqrt{x_{i}x_{j}(x_{i}x_{j}+4x_{k})}}$
$\displaystyle+r_{2}^{(2)}(x_{i},x_{j},x_{k})\left(f_{3}^{(2,3)}(x_{i},x_{j},x_{k})+{1\over
3}f_{3}^{(1,1)}(x_{k},x_{i},x_{j})f_{15}^{(2,2)}(x_{k},x_{i},x_{j})\right)$
$\displaystyle+\frac{r_{3}^{(2)}(x_{i},x_{j},x_{k})}{\sqrt{x_{i}(x_{i}-4)}}f_{5}^{(2,3)}(x_{i},x_{j},x_{k})+\frac{r_{4}^{(2)}(x_{i},x_{j},x_{k})}{\sqrt{x_{i}(x_{i}-4)}}f_{13}^{(2,4)}(x_{i},x_{j},x_{k})$
$\displaystyle+r_{5}^{(2)}(x_{i},x_{j},x_{k})f_{14}^{(2,3)}(x_{i},x_{j},x_{k})+r_{6}^{(2)}(x_{i},x_{j},x_{k})f_{5}^{(1,2)}(x_{i},x_{j},x_{k})$
$\displaystyle+r_{7}^{(2)}(x_{i},x_{j},x_{k})f_{5}^{(1,3)}(x_{i},x_{j},x_{k})-\left({6\over
x_{i}}+{4\over x_{j}}+{6\over x_{k}}\right)f_{17}^{(2,4)}(x_{i},x_{j},x_{k})$
$\displaystyle+\left({4\over x_{i}}+{4\over x_{j}}+{8\over
x_{k}}\right)\left(f_{21}^{(2,4)}(x_{i},x_{j},x_{k})-{1\over
2}f_{24}^{(2,4)}(x_{k},x_{i},x_{j})\right)+r_{8}^{(2)}(x_{i},x_{j},x_{k})f_{20}^{(2,3)}(x_{i},x_{j},x_{k})$
$\displaystyle\left.+r_{10}^{(2)}(x_{i},x_{j},x_{k})\frac{f_{19}^{(2,2)}(x_{i},x_{j},x_{k})}{\sqrt{x_{i}x_{j}(x_{i}x_{j}+4x_{k})}}-\left(32+\frac{16x_{i}x_{j}}{x_{k}}\right)\frac{f_{22}^{(2,3)}(x_{i},x_{j},x_{k})}{\sqrt{x_{i}x_{j}(x_{i}x_{j}+4x_{k})}}\right]$
$\displaystyle-\mathop{\sum_{(i,j,k)=(s,t,u),(s,u,t)}}_{(t,u,s),(t,s,u),(u,s,t),(u,t,s)}\left[\left({4\over
x_{i}}+{6\over x_{j}}+{6\over
x_{k}}\right)f_{26}^{(2,4)}(x_{i},x_{j},x_{k})+\left({2\over x_{j}}-{2\over
x_{k}}\right)f_{29}^{(2,4)}(x_{i},x_{j},x_{k})\right.$
$\displaystyle-\frac{r_{9}^{(2)}(x_{i},x_{j},x_{k})f_{18}^{(2,3)}(x_{i},x_{j},x_{k})}{\sqrt{x_{i}\left(x_{i}(x_{j}-1)^{2}-4x_{j}^{2}\right)}}-\left(32+\frac{16x_{i}x_{j}}{x_{k}}\right)\frac{f_{19}^{(2,3)}(x_{i},x_{j},x_{k})}{\sqrt{x_{i}x_{j}(x_{i}x_{j}+4x_{k})}}$
$\displaystyle\left.-r_{11}^{(2)}(x_{i},x_{j},x_{k})\frac{f_{25}^{(2,4)}(x_{i},x_{j},x_{k})}{\sqrt{x_{i}x_{j}(x_{i}x_{j}+4x_{k})}}+\left(2+\frac{x_{i}x_{j}}{x_{k}}\right)\frac{4(x_{i}-2)}{\sqrt{x_{i}(x_{i}-4)}}\frac{f_{27}^{(2,4)}(x_{i},x_{j},x_{k})}{\sqrt{x_{i}x_{j}(x_{i}x_{j}+4x_{k})}}\right.$
$\displaystyle\left.\left.-\left(\frac{4(x_{i}-2)}{x_{k}}-\frac{2x_{k}}{x_{j}}\right)\frac{f_{28}^{(2,4)}(x_{i},x_{j},x_{k})}{\sqrt{x_{i}(x_{i}-4)}}\right]\right\\}\,,$
(59)

where the rational coefficients are defined as

$\displaystyle r_{1}^{(2)}(x_{i},x_{j},x_{k})$
$\displaystyle=-3(x_{j}-4)\left(\frac{x_{i}-x_{k}-x_{i}x_{j}}{x_{j}(x_{i}-1)^{2}-4x_{i}^{2}}+\frac{x_{k}-x_{i}-x_{k}x_{j}}{x_{j}(x_{k}-1)^{2}-4x_{k}^{2}}\right)\,,$
$\displaystyle r_{2}^{(2)}(x_{i},x_{j},x_{k})$ $\displaystyle={2\over
x_{i}}+{9\over x_{j}}+{30\over
x_{k}}+\frac{36}{x_{i}x_{j}}-\frac{8}{x_{i}x_{k}}+\frac{20+x_{i}-2x_{j}}{x_{i}x_{k}+4x_{j}}+\frac{112+15x_{i}+13x_{j}}{2(x_{i}x_{j}+4x_{k})}\,,$
$\displaystyle r_{3}^{(2)}(x_{i},x_{j},x_{k})$
$\displaystyle=-\frac{2(x_{i}-4)^{2}(x_{k}-x_{j})^{2}}{(x_{i}x_{j}+4x_{k})(x_{i}x_{k}+4x_{j})}\,,$
$\displaystyle r_{4}^{(2)}(x_{i},x_{j},x_{k})$
$\displaystyle=2(x_{i}-4)\frac{x_{j}^{2}+x_{k}^{2}}{x_{i}x_{j}x_{k}}\,,$
$\displaystyle r_{5}^{(2)}(x_{i},x_{j},x_{k})$ $\displaystyle=-{1\over
12}\left({93\over x_{i}}+{68\over x_{j}}+{135\over
x_{k}}+\frac{88}{x_{i}x_{j}}+\frac{168}{x_{i}x_{k}}+\frac{116+8x_{i}+21x_{j}}{(x_{i}x_{j}+4x_{k})}+\frac{44+x_{i}+13x_{j}}{x_{j}x_{k}+4x_{i}}\right)\,,$
$\displaystyle r_{6}^{(2)}(x_{i},x_{j},x_{k})$ $\displaystyle=-4+{4\over
x_{j}}+\frac{4x_{j}(x_{i}^{3}-2x_{i}x_{k}+x_{k}^{3})}{x_{i}^{2}x_{k}^{2}}+\frac{2+x_{i}+x_{k}(x_{j}-3)}{x_{j}(x_{k}-1)^{2}-4x_{k}^{2}}+\frac{2+x_{k}+x_{i}(x_{j}-3)}{x_{j}(x_{i}-1)^{2}-4x_{i}^{2}}$
$\displaystyle+\frac{16-6x_{i}-6x_{j}(x_{k}-3)}{x_{k}(x_{j}-1)^{2}-4x_{j}^{2}}+\frac{16-6x_{k}-6x_{j}(x_{i}-3)}{x_{i}(x_{j}-1)^{2}-4x_{j}^{2}}\,,$
$\displaystyle r_{7}^{(2)}(x_{i},x_{j},x_{k})$
$\displaystyle=-6r_{5}^{(2)}(x_{i},x_{j},x_{k})-32\left(x_{i}^{-1}+x_{j}^{-1}+x_{k}^{-1}\right)\,,$
$\displaystyle r_{8}^{(2)}(x_{i},x_{j},x_{k})$ $\displaystyle={8\over
x_{i}}+{8\over x_{j}}-{34\over x_{k}}-{88\over
x_{i}x_{j}}-\frac{72-9x_{k}}{x_{i}x_{j}+4x_{k}}\,,$ $\displaystyle
r_{9}^{(2)}(x_{i},x_{j},x_{k})$
$\displaystyle=37-8x_{i}+\frac{(17-9x_{i})x_{i}-8}{x_{k}}-\frac{2(4+x_{k})}{x_{j}}-\frac{x_{k}(28+9x_{i})-4x_{i}}{x_{i}x_{j}+4x_{k}}-\frac{2x_{i}(1+x_{j})}{x_{i}(x_{j}-1)^{2}-4x_{j}^{2}}\,,$
$\displaystyle r_{10}^{(2)}(x_{i},x_{j},x_{k})$
$\displaystyle=56+\frac{4x_{i}x_{j}(x_{i}x_{j}-(x_{k}-4)x_{k})}{x_{k}^{2}}+\frac{4x_{j}(1+2x_{i})}{x_{j}(x_{i}-1)^{2}-4x_{i}^{2}}+\frac{4x_{i}(1+2x_{j})}{x_{i}(x_{j}-1)^{2}-4x_{j}^{2}}\,,$
$\displaystyle r_{11}^{(2)}(x_{i},x_{j},x_{k})$
$\displaystyle=4(1-x_{j})-\frac{2x_{i}^{2}}{x_{k}}-\frac{4(x_{i}x_{k}+4x_{k})}{x_{i}x_{j}+4x_{k}}.$
(60)

They are the shortest and second-shortest amplitudes, and are written in a
form that manifests the $s,t,u$ symmetry. It is also interesting to note that
in the high-energy limit, i.e., $x_{s},-x_{t},-x_{u}\gg 1$, the all-plus
amplitude should approach to a constant as shown in ref. Bern:2001dg . We can
indeed verify that eq.(58) is only left with the first constant term that
matches with ref. Bern:2001dg .

The most-lengthy amplitudes are the two-minus-two-plus helicity amplitudes.
Thanks to the relation in eq.(31), we only need to present
$\mathcal{M}^{(1,0,f)}_{--++}$, which is:

$\displaystyle i\mathcal{M}^{(1,0,f)}_{--++}$
$\displaystyle=N_{c,f}Q_{f}^{4}\alpha^{2}\frac{\alpha_{s}}{\pi}C_{F,f}\left\\{4+\left({8\over
x_{s}}+{4x_{s}\over
x_{t}x_{u}}\right)\left[\mathop{\sum_{(i,j,k)=(s,t,u),}}_{(t,u,s),(u,s,t)}\frac{x_{j}(x_{j}-4)f_{3}^{(1,1)}(x_{i},x_{j},x_{k})}{\sqrt{x_{j}(x_{j}-4)}}\right]\right.$
$\displaystyle+r_{17}^{(2)}(x_{u},x_{s},x_{t})f_{5}^{(1,2)}(x_{u},x_{s},x_{t})-\frac{8x_{s}^{2}(x_{s}^{2}-16)}{(x_{s}x_{t}+4x_{u})(x_{s}x_{u}+4x_{t})}\frac{f_{5}^{(2,3)}(x_{s},x_{t},x_{u})}{\sqrt{x_{s}(x_{s}-4)}}$
$\displaystyle+r_{29}^{(2)}(x_{s},x_{t},x_{u})\left(f_{3}^{(2,3)}(x_{s},x_{t},x_{u})+{1\over
3}f_{15}^{(2,2)}(x_{u},x_{s},x_{t})f_{3}^{(1,1)}(x_{u},x_{s},x_{t})\right)$
$\displaystyle-8(2x_{s}+x_{t}x_{u})\frac{f_{6}^{(1,2)}(x_{t},x_{u},x_{s})}{\sqrt{x_{t}x_{u}(x_{t}x_{u}+4x_{s})}}+r_{24}^{(2)}(x_{t},x_{u},x_{s})\frac{f_{19}^{(2,2)}(x_{t},x_{u},x_{s})}{\sqrt{x_{t}x_{u}(x_{t}x_{u}+4x_{s})}}$
$\displaystyle+\frac{8}{x_{s}^{2}}(8+x_{s}-x_{t}^{2}-x_{u}^{2})f_{21}^{(2,4)}(x_{t},x_{u},x_{s})-r_{25}^{(2)}(x_{t},x_{u},x_{s})\frac{f_{22}^{(2,3)}(x_{t},x_{u},x_{s})}{\sqrt{x_{t}x_{u}(x_{t}x_{u}+4x_{s})}}$
$\displaystyle+{16\over
x_{s}^{2}}f_{24}^{(2,4)}(x_{s},x_{t},x_{u})+r_{26}^{(2)}(x_{u},x_{s},x_{t})\left(f_{14}^{(2,3)}(x_{u},x_{s},x_{t})-6f_{5}^{(1,3)}(x_{u},x_{s},x_{t})\right)$
$\displaystyle+\mathop{\sum_{(i,j,k)=(s,t,u),}}_{(s,u,t)}\left[-16x_{i}\frac{f_{6}^{(1,2)}(x_{i},x_{j},x_{k})}{\sqrt{x_{i}x_{j}(x_{i}x_{j}+4x_{k})}}-16(x_{i}-2)\frac{\left(f_{19}^{(2,3)}(x_{i},x_{j},x_{k})-f_{22}^{(2,3)}(x_{i},x_{j},x_{k})\right)}{\sqrt{x_{i}x_{j}(x_{i}x_{j}+4x_{k})}}\right.$
$\displaystyle+\frac{4(x_{i}-2)^{2}}{\sqrt{x_{i}(x_{i}-4)}}\frac{f_{27}^{(2,4)}(x_{i},x_{j},x_{k})}{\sqrt{x_{i}x_{j}(x_{i}x_{j}+4x_{k})}}+r_{12}^{(2)}(x_{i},x_{j},x_{k})\left(f_{17}^{(2,4)}(x_{i},x_{j},x_{k})+f_{26}^{(2,4)}(x_{k},x_{j},x_{i})\right)$
$\displaystyle+r_{13}^{(2)}(x_{i},x_{j},x_{k})f_{20}^{(2,3)}(x_{i},x_{j},x_{k})+r_{16}^{(2)}(x_{i},x_{j},x_{k})f_{5}^{(1,2)}(x_{i},x_{j},x_{k})$
$\displaystyle+r_{18}^{(2)}(x_{i},x_{j},x_{k})\frac{f_{3}^{(1,2)}(x_{i},x_{j},x_{k})}{\sqrt{x_{j}(x_{j}-4)}}+r_{20}^{(2)}(x_{i},x_{j},x_{k})\frac{f_{15}^{(2,2)}(x_{i},x_{j},x_{k})}{\sqrt{x_{j}(x_{j}-4)}}$
$\displaystyle+\frac{8x_{i}(2x_{i}+(x_{i}-4)x_{j}^{2})}{x_{j}(x_{i}x_{j}+4x_{k})}\frac{f_{18}^{(2,3)}(x_{i},x_{j},x_{k})}{\sqrt{x_{i}(x_{i}(x_{j}-1)^{2}-4x_{j}^{2})}}$
$\displaystyle+r_{23}^{(2)}(x_{i},x_{j},x_{k})\frac{f_{19}^{(2,2)}(x_{i},x_{j},x_{k})}{\sqrt{x_{i}x_{j}(x_{i}x_{j}+4x_{k})}}+\frac{16x_{i}^{2}}{x_{i}x_{j}+4x_{k}}\frac{f_{25}^{(2,4)}(x_{i},x_{j},x_{k})}{\sqrt{x_{i}x_{j}(x_{i}x_{j}+4x_{k})}}$
$\displaystyle+r_{27}^{(2)}(x_{i},x_{j},x_{k})f_{14}^{(2,3)}(x_{i},x_{j},x_{k})+r_{28}^{(2)}(x_{i},x_{j},x_{k})f_{5}^{(1,3)}(x_{i},x_{j},x_{k})$
$\displaystyle\left.+r_{31}^{(2)}(x_{i},x_{j},x_{k})f_{15}^{(2,2)}(x_{i},x_{j},x_{k})f_{3}^{(1,1)}(x_{i},x_{j},x_{k})\right]$
$\displaystyle+\mathop{\sum_{(i,j,k)=(u,s,t),}}_{(t,s,u)}\left[-16(x_{j}-2)\frac{f_{19}^{(2,3)}(x_{i},x_{j},x_{k})}{\sqrt{x_{i}x_{j}(x_{i}x_{j}+4x_{k})}}+r_{21}^{(2)}(x_{i},x_{j},x_{k})\frac{f_{18}^{(2,3)}(x_{i},x_{j},x_{k})}{\sqrt{x_{i}(x_{i}(x_{j}-1)^{2}-4x_{j}^{2})}}\right.$
$\displaystyle-\frac{8}{x_{j}^{2}}(4+x_{j}+3x_{k}-x_{i}^{2})\left(f_{26}^{(2,4)}(x_{i},x_{j},x_{k})+\frac{1}{2}f_{29}^{(2,4)}(x_{i},x_{j},x_{k})\right)$
$\displaystyle+\frac{4(x_{j}-2)(x_{i}-2)}{\sqrt{x_{i}(x_{i}-4)}}\frac{f_{27}^{(2,4)}(x_{i},x_{j},x_{k})}{\sqrt{x_{i}x_{j}(x_{i}x_{j}+4x_{k})}}+4(4x_{j}+12x_{i}-x_{i}^{3}-x_{k}^{2})\frac{f_{28}^{(2,4)}(x_{i},x_{j},x_{k})}{x_{j}^{2}\sqrt{x_{i}(x_{i}-4)}}$
$\displaystyle\left.+r_{30}^{(2)}(x_{i},x_{j},x_{k})f_{3}^{(2,3)}(x_{i},x_{j},x_{k})\right]$
$\displaystyle+\mathop{\sum_{(i,j,k)=(t,u,s),}}_{(u,t,s)}\left[r_{14}^{(2)}(x_{i},x_{j},x_{k})f_{20}^{(2,3)}(x_{i},x_{j},x_{k})+\frac{r_{25}^{(2)}(x_{i},x_{j},x_{k})}{\sqrt{x_{i}x_{j}(x_{i}x_{j}+4x_{k})}}\left(f_{19}^{(2,3)}(x_{i},x_{j},x_{k})\right.\right.$
$\displaystyle\left.-\frac{(x_{i}-2)f_{27}^{(2,4)}(x_{i},x_{j},x_{k})}{4\sqrt{x_{i}(x_{i}-4)}}\right)-2(x_{i}^{2}+x_{j}^{2}-4x_{k})\frac{(x_{i}-2)f_{28}^{(2,4)}(x_{i},x_{j},x_{k})}{x_{k}^{2}\sqrt{x_{i}(x_{i}-4)}}$
$\displaystyle+r_{15}^{(2)}(x_{i},x_{j},x_{k})\left(f_{24}^{(2,4)}(x_{i},x_{j},x_{k})-f_{29}^{(2,4)}(x_{i},x_{j},x_{k})\right)+r_{19}^{(2)}(x_{i},x_{j},x_{k})\frac{f_{13}^{(2,4)}(x_{i},x_{j},x_{k})}{\sqrt{x_{i}(x_{i}-4)}}$
$\displaystyle\left.\left.+r_{22}^{(2)}(x_{i},x_{j},x_{k})\frac{f_{18}^{(2,3)}(x_{i},x_{j},x_{k})}{\sqrt{x_{i}(x_{i}(x_{j}-1)^{2}-4x_{j}^{2})}}\right]\right\\}\,.$
(61)

Here the additional rational coefficients are:

$\displaystyle r_{12}^{(2)}(x_{i},x_{j},x_{k})$
$\displaystyle=4+\frac{8(1+x_{j})}{x_{i}}+\frac{8(x_{j}^{2}+3x_{j}-4)}{x_{i}^{2}}\,,$
$\displaystyle r_{13}^{(2)}(x_{i},x_{j},x_{k})$
$\displaystyle=-12-\frac{108+24x_{j}}{x_{i}^{2}}+{72+24x_{j}\over
x_{i}}+{48\over x_{j}}-{16\over
x_{k}}-\frac{128}{x_{i}x_{j}}+\frac{144-36x_{i}}{x_{i}x_{j}+4x_{k}}\,,$
$\displaystyle r_{14}^{(2)}(x_{i},x_{j},x_{k})$ $\displaystyle=24-{102\over
x_{k}}+\frac{34x_{k}-64}{x_{i}x_{j}}+\frac{36-48x_{i}x_{j}}{x_{k}^{2}}\,,$
$\displaystyle r_{15}^{(2)}(x_{i},x_{j},x_{k})$
$\displaystyle=2+\frac{4+4x_{j}}{x_{k}}+\frac{4(x_{j}^{2}+3x_{j}-4)}{x_{k}^{2}}\,,$
$\displaystyle r^{(2)}_{16}(x_{i},x_{j},x_{k})$ $\displaystyle=38-{48\over
x_{j}}-\frac{152-48x_{j}}{x_{i}}-\frac{48+8x_{j}}{x_{k}}+{4x_{j}^{2}\over
x_{k}^{2}}-{8x_{j}^{2}\over
x_{i}^{2}}-\frac{4-12x_{i}}{x_{j}(x_{i}-1)^{2}-4x_{i}^{2}}$
$\displaystyle-\frac{4+2x_{i}-2(x_{j}+1)x_{k}}{x_{j}(x_{k}-1)^{2}-4x_{k}^{2}}-\frac{8\left(4+5x_{j}-5x_{i}-2x_{j}x_{k}\right)}{x_{k}(x_{j}-1)^{2}-4x_{j}^{2}}\,,$
$\displaystyle r^{(2)}_{17}(x_{i},x_{j},x_{k})$
$\displaystyle=4\left(4+{8\over x_{j}}+{x_{j}^{2}\over
x_{k}^{2}}+{x_{j}^{2}\over
x_{i}^{2}}-\frac{4(x_{j}-3)x_{j}}{x_{k}x_{i}}-\frac{4(x_{j}+2)}{x_{k}(x_{j}-1)^{2}-4x_{j}^{2}}-\frac{4(x_{j}+2)}{x_{i}(x_{j}-1)^{2}-4x_{j}^{2}}\right)\,,$
$\displaystyle r^{(2)}_{18}(x_{i},x_{j},x_{k})$
$\displaystyle=x_{j}(x_{j}-4)\left({40\over x_{i}}+{20\over
x_{j}}-\frac{3(1+x_{k})^{2}-3x_{i}(3+x_{k}^{2})}{x_{i}(x_{j}(x_{k}-1)^{2}-4x_{k}^{2})}+\frac{3(x_{i}-1)^{2}(x_{i}+1)}{x_{i}(x_{j}(x_{i}-1)^{2}-4x_{i}^{2})}\right)\,,$
$\displaystyle r_{19}^{(2)}(x_{i},x_{j},x_{k})$
$\displaystyle=(x_{k}-4)\left(2-\frac{12x_{k}}{x_{j}^{2}}+\frac{4(x_{k}-2)}{x_{j}}\right)\,,$
$\displaystyle r_{20}^{(2)}(x_{i},x_{j},x_{k})$
$\displaystyle=-\frac{1}{6}r_{18}^{(2)}(x_{i},x_{j},x_{k})-\frac{8}{3x_{i}}(x_{j}-x_{k})(x_{j}-4)\,,$
$\displaystyle r_{21}^{(2)}(x_{i},x_{j},x_{k})$
$\displaystyle=4\left(33-4x_{j}+3x_{k}+x_{i}x_{j}-2x_{i}^{2}-\frac{x_{i}(2x_{i}+9)}{x_{j}^{2}}+\frac{x_{i}(4x_{i}+19)}{x_{j}}+\frac{2x_{i}-2x_{j}(x_{j}-x_{k})}{x_{i}(x_{j}-1)^{2}-4x_{j}^{2})}\right)\,,$
$\displaystyle r_{22}^{(2)}(x_{i},x_{j},x_{k})$
$\displaystyle=8-2x_{k}-6x_{j}+2x_{k}x_{j}+6x_{j}^{2}+\frac{4x_{j}(x_{j}+1)^{2}(9+4x_{j})}{x_{k}^{2}}+\frac{4(x_{k}-x_{j})(x_{k}-x_{j}+x_{j}x_{i})}{x_{i}(x_{j}-1)^{2}-4x_{j}^{2}}$
$\displaystyle+\frac{4(9+x_{j}(5x_{j}^{2}+5x_{j}-11))}{x_{k}}\,,$
$\displaystyle r_{23}^{(2)}(x_{i},x_{j},x_{k})$
$\displaystyle=-264-4x_{i}(x_{j}-2)+8x_{j}(x_{j}+3)+\frac{8(8x_{i}^{2}-3x_{j}(1-x_{i}))}{x_{j}(x_{i}-1)^{2}-4x_{i}^{2}}-\frac{4x_{j}^{4}}{x_{k}^{2}}+\frac{4x_{j}^{2}(x_{j}+14)}{x_{k}}\,,$
$\displaystyle r_{24}^{(2)}(x_{i},x_{j},x_{k})$
$\displaystyle=4x_{i}x_{j}(4x_{k}+x_{i}x_{j})\left({2\over
x_{k}^{2}}+\frac{3(x_{k}-6)}{x_{i}x_{j}x_{k}}-\frac{x_{k}-2+x_{i}(x_{j}-1)}{x_{j}(x_{j}(x_{i}-1)^{2}-4x_{i}^{2})}-\frac{x_{k}-2+x_{j}(x_{i}-1)}{x_{i}(x_{i}(x_{j}-1)^{2}-4x_{j}^{2})}\right)\,,$
$\displaystyle r_{25}^{(2)}(x_{i},x_{j},x_{k})$
$\displaystyle=32-16x_{k}-8x_{i}x_{j}+\frac{64x_{i}x_{j}}{x_{k}}+\frac{16x_{i}^{2}x_{j}^{2}}{x_{k}^{2}}\,,$
$\displaystyle r_{26}^{(2)}(x_{i},x_{j},x_{k})$
$\displaystyle=\frac{1}{3}\left(28-{90\over x_{j}^{2}}-{74\over
x_{j}}+{19x_{j}-32\over x_{i}x_{k}}+{12(4-x_{j})\over
x_{i}x_{j}+4x_{k}}+{12(4-x_{j})\over x_{j}x_{k}+4x_{i}}\right)\,,$
$\displaystyle r_{27}^{(2)}(x_{i},x_{j},x_{k})$
$\displaystyle=\frac{1}{6}\left(37+{45\over x_{j}}-{32\over
x_{i}x_{j}}+{73(2x_{j}-1)\over x_{i}}+{2(40x_{j}^{2}-175x_{j}-63)\over
x_{i}^{2}}+{x_{i}-32\over x_{i}x_{k}}-{16(x_{i}-4)\over
x_{i}x_{j}+4x_{k}}\right)\,,$ $\displaystyle r_{28}^{(2)}(x_{i},x_{j},x_{k})$
$\displaystyle=-45+{105-162x_{j}\over x_{i}}+{3\over x_{j}}+{32\over
x_{i}x_{j}}+{126+350x_{j}-48x_{j}^{2}\over x_{i}^{2}}-{1\over x_{k}}+{32\over
x_{i}x_{k}}+{16(x_{i}-4)\over x_{i}x_{j}+4x_{k}}\,,$ $\displaystyle
r_{29}^{(2)}(x_{i},x_{j},x_{k})$ $\displaystyle=8+{252\over
x_{i}^{2}}-{52\over x_{i}}+{8(x_{i}-8)\over x_{j}x_{k}}+{30(x_{i}-4)\over
x_{i}x_{j}+4x_{k}}+{30(x_{i}-4)\over x_{i}x_{k}+4x_{j}}\,,$ $\displaystyle
r_{30}^{(2)}(x_{i},x_{j},x_{k})$ $\displaystyle=-10+{29\over x_{i}}-{32\over
x_{i}x_{j}}+{23-20x_{k}\over x_{j}}+{25\over x_{k}}-{32\over
x_{j}x_{k}}+{8(x_{j}-4)\over x_{i}x_{j}+4x_{k}}+{90-52x_{k}-40x_{k}^{2}\over
x_{j}^{2}}\,,$ $\displaystyle r_{31}^{(2)}(x_{i},x_{j},x_{k})$
$\displaystyle=\frac{1}{3}\left(-26+{59-52x_{j}\over x_{i}}+{21\over
x_{j}}-{32\over x_{i}x_{j}}+{90+52x_{j}-40x_{j}^{2}\over x_{i}^{2}}+{25\over
x_{k}}-{32\over x_{i}x_{k}}+{8(x_{i}-4)\over x_{i}x_{j}+4x_{k}}\right)\,.$
(62)

The amplitude $\mathcal{M}^{(1,0,f)}_{--++}$ is symmetric in $x_{t}$ and
$x_{u}$. The expressions of all these two-loop amplitudes are also provided in
the ancillary files. We have also implemented them into a Fortran code, that
will be published elsewhere. The timing of evaluating the 2-loop amplitude
square after summing all helicity configurations is around $0.3$ second for a
random phase space point on a single core with the 2.6 GHz Intel Core i7 CPU.

We have also carried out a few checks. First, the crossing symmetry relations
of the amplitudes have been explicitly verified by plugging in the concrete
analytic expressions. Our high-energy limit approaches the massless results in
ref. Bern:2001dg  numerically, and our low-energy limit approaches the results
present in ref. Martin:2003gb  (also see appendix E). Finally, at the cross
section level, our results agree perfectly with the pure numerical approach
based on the local unitarity construction Capatti:2020xjc ; Capatti:2022tit .

##  5 Conclusion

In this paper, we have presented a compact and analytic two-loop helicity
amplitudes of the LbL process $\gamma\gamma\to\gamma\gamma$ keeping the full
dependence of the internal fermion mass. The results remain simple enough even
at two loops to be written in a few pages thanks to the numerous symmetries of
the process. They could serve as a benchmark for developing novel multi-loop
techniques, such as the recent local unitarity approach Capatti:2020xjc ;
Capatti:2022tit . The connection between the simplicity of our amplitudes and
the string-derived Bern-Kosower formalism Bern:1990cu ; Bern:1991aq ;
Strassler:1992zr  needs to be explored in the future. The phenomenological
applications and the comparisons with the local unitarity method will be
reported in a companion paper Ajjath:2023yy . As an outlook, it is intriguing
in the future to take into account the full electroweak corrections and next-
to-next-to-leading order QCD and QED corrections. Their phenomenological
relevance may show up in UPCs at the high-luminosity LHC by tagging proton(s)
with forward proton spectrometers and future $e^{+}e^{-}$ colliders.

## Acknowledgements

We thank Claude Duhr, Mathijs Fraaije, Johannes Henn, Valentin Hirschi, Kirill
Melnikov, Lukas Simon, Vasily Sotnikov, Guoxing Wang and Simone Zoia for
useful discussions. This work is supported by the grants from the ERC (grant
101041109 ‘BOSON’, grant 101043686 ‘LoCoMotive’), the French ANR (grant
ANR-20-CE31-0015 ‘PrecisOnium’), the French LIA FCPPL, the European Union’s
Horizon 2020 research and innovation program (grant 824093 STRONG-2020, EU
Virtual Access ‘NLOAccess’).

##  Appendix A Kinematic regions and analytic continuation

In this section, we determine the kinematic regions that will be needed in
two-loop amplitude calculations. In particular, we use the following
parametrisation for the Mandelstam variables for all the polylogarithmic
integrals:

$\displaystyle\frac{s_{ij}}{m_{f}^{2}}$ $\displaystyle=$ $\displaystyle
x_{ij}=-\frac{4(w-z)^{2}}{(1-w^{2})(1-z^{2})},\quad\frac{s_{jk}}{m^{2}}=x_{jk}=-\frac{(w-z)^{2}}{wz},\quad
s_{ik}=-s_{ij}-s_{jk}.$ (63)

In the above parametrisation, there are multiple solutions of $w$ and $z$ in
terms of $x_{ij}$ and $x_{jk}$. We choose a unique solution that covers the
whole $2\to 2$ scattering phase space. First of all, since there is a symmetry
in $w$ and $z$, we assume $\Re{z}\geq\Re{w}$ if $w,z\in\mathbb{R}$. we derive
the regions of $w$ and $z$ that cover the required three kinematic regions:

  * •

Region I: $s_{ij}<0,s_{jk}<0,s_{ik}>0$ $\longrightarrow$ $0<w\leq z<1$.

  * •

Region II: $s_{ij}>0,s_{jk}<0,s_{ik}<0$ $\longrightarrow$ this region can be
further divided into 2 sub-regions.

    * –

$s_{ij}>4m_{f}^{2},s_{jk}<0,s_{ik}<0$: $\sqrt{2}-1<w<1<z<\frac{1}{w}$ or
$0<w<\sqrt{2}-1,1<z<\frac{1+w}{1-w}$. It can also be written as $0<w<1<z<{\rm
min}(\frac{1}{w},\frac{1+w}{1-w})$.

    * –

$0<s_{ij}<4m_{f}^{2},s_{jk}<0,s_{ik}<0$: $w,z\in\mathbb{C}$,
$z=\frac{w}{|w|^{2}}$.
$\sqrt{2}-1<\Re{w}<1,0<\Im{w}<\sqrt{1-\left(\Re{w}\right)^{2}}$ or
$0<\Re{w}<\sqrt{2}-1,\sqrt{1-2\Re{w}-\left(\Re{w}\right)^{2}}<\Im{w}<\sqrt{1-\left(\Re{w}\right)^{2}}$.
It can also be rewritten as $0<\Re{w}<1$ and $\sqrt{{\rm
max}(0,1-2\Re{w}-\left(\Re{w}\right)^{2})}<\Im{w}<\sqrt{1-\left(\Re{w}\right)^{2}}$.

  * •

Region III: $s_{ij}<0,s_{jk}>0,s_{ik}<0$ $\longrightarrow$ this region can
also be further divided into 2 sub-regions.

    * –

$s_{ij}<0,s_{jk}>4m_{f}^{2},s_{ik}<0$: $1-\sqrt{2}<w<0<z<-w$ or
$-1<w<1-\sqrt{2},0<z<\frac{1+w}{1-w}$. It can also be written as
$-1<w<0<z<{\rm min}(-w,\frac{1+w}{1-w})$.

    * –

$s_{ij}<0,0<s_{jk}<4m_{f}^{2},s_{ik}<0$: $w,z\in\mathbb{C}$, $z=-w^{*}$.
$1-\sqrt{2}<\Re{w}<0,0<\Im{w}<\sqrt{1+2\Re{w}-\left(\Re{w}\right)^{2}}$.

The regions of $(w,z)$ in the above cover the entire phase space only once.

For $w,z\in\mathbb{R}$, we need to know the sign of their imaginary parts in
order to choose the correct Riemann sheet in GPLs. Since $s_{ij}\to
s_{ij}+i0,s_{jk}\to s_{jk}+i0,s_{ik}\to s_{ik}+i0$ and $m_{f}^{2}\to
m_{f}^{2}-i0$, it is sufficient for our purpose to take $w\to w-i0$. For the
square roots in the prefactors of the canonical basis eq.(35), we always
evaluate them in the first Riemann sheet for simplicity. However, a minus sign
stemming from the square roots should be introduced in the canonical master
integrals $f_{n}^{(2,\omega)}(x_{ij},x_{jk},x_{ik})$ with the corresponding
square roots when the following condition is satisfied

$\displaystyle\left\\{\begin{array}[]{ll}\sqrt{s_{ij}(s_{ij}-4m_{f}^{2})}&~{}~{}{\rm
with}~{}s_{ij}>4m_{f}^{2}\\\ \sqrt{s_{jk}(s_{jk}-4m_{f}^{2})}&~{}~{}{\rm
with}~{}0<s_{jk}<4m_{f}^{2}\\\
\sqrt{s_{ij}s_{jk}(s_{ij}s_{jk}-4m_{f}^{2}(s_{ij}+s_{jk}))}&~{}~{}{\rm
with}~{}s_{ij}>0~{}{\rm or}~{}s_{jk}>0.\end{array}\right.\,$ (67)

For instance, for the two square roots $\sqrt{s_{ij}(s_{ij}-4m_{f}^{2})}$ and
$\sqrt{s_{ij}s_{jk}(s_{ij}s_{jk}-4m_{f}^{2}(s_{ij}+s_{jk}))}$ in
$f_{27}^{(2,4)}$, we should multiply the prefactor
$(-1)^{\theta(x_{ij}-4)}(-1)^{\theta(x_{ij})+\theta(x_{jk})}$ in front of the
expression consisting of $G$’s as what is similarly done in eq.(55). We have
taken care of these prefactors in the relations given in eqs.(50) and (51),
and the remaining square roots in the amplitudes can be evaluated on the first
Riemann sheet. The same applies to the one-loop master integrals
$f_{n}^{(1,\omega)}(x_{ij},x_{jk},x_{ik})$, which are defined in appendix C.

##  Appendix B Analytic boundary conditions

In this section, we present the analytic boundary conditions in the regions II
and III defined in appendix A. For the region I, we simply choose
$x_{s}=x_{t}=x_{u}=0$ as the boundary point, where all the master integrals
are zero except the first one $f_{1}^{(2,\omega)}$. In the regions II and III,
we use the parametrisation $x_{s}=s/m^{2}=-(y_{s}-1)^{2}/y_{s}>0$ to
rationalise the square roots in the boundary conditions. The solution is

$\displaystyle y_{s}$ $\displaystyle=$ $\displaystyle
1-x_{s}/2-\sqrt{x_{s}(x_{s}-4)}/2.$ (68)

Given the analytic continuation of $x_{s}\to x_{s}+i0$, we have $y_{s}\to
y_{s}-i0$ when $x_{s}>4$. When $0<x_{s}<4$, we simply take the square root in
the first Riemann sheet in eq.(68) with the analytic continuation
$x_{s}(x_{s}-4)\to x_{s}(x_{s}-4)+i0$. We write down the independent non-zero
boundary constants that are needed for evaluating the five master integrals
$f_{18}^{(2,3)},f_{25}^{(2,4)},f_{26}^{(2,4)},f_{28}^{(2,4)},f_{29}^{(2,4)}$
with the non-rationalised square root $\rho$. They are:

$\displaystyle f_{2}^{(2,1)}(x_{s},0,-x_{s})$ $\displaystyle=$ $\displaystyle
G(0;y_{s}),\quad
f_{2}^{(2,2)}(x_{s},0,-x_{s})=-\zeta_{2}+G(0,0;y_{s})-2G(-1,0;y_{s}),$
$\displaystyle f_{4}^{(2,3)}(x_{s},0,-x_{s})$ $\displaystyle=$ $\displaystyle
3\zeta_{3}+\zeta_{2}G(0;y_{s})-G(0,0,0;y_{s})+2G(0,-1,0;y_{s}),$
$\displaystyle f_{10}^{(2,3)}(x_{s},0,-x_{s})$ $\displaystyle=$ $\displaystyle
4\zeta_{3}+2\zeta_{2}G(0;y_{s})+G(0,0,0;y_{s})-2G(0,1,0;y_{s}),$
$\displaystyle f_{11}^{(2,2)}(x_{s},0,-x_{s})$ $\displaystyle=$
$\displaystyle\zeta_{2}+\frac{1}{2}G(0,0;y_{s})-G(1,0;y_{s}),$ $\displaystyle
f_{18}^{(2.3)}(x_{s},0,-x_{s})$ $\displaystyle=$
$\displaystyle-3\zeta_{3}-\zeta_{2}G(0;y_{s})+G(0,1,0;y_{s})-G(1,0,0;y_{s}),$
$\displaystyle f_{20}^{(2,3)}(x_{s},0,-x_{s})$ $\displaystyle=$
$\displaystyle\zeta_{3}+\zeta_{2}G(0;y_{s})+G(0,0,0;y_{s})-G(0,1,0;y_{s})-G(1,0,0;y_{s}),$
$\displaystyle f_{23}^{(2,3)}(x_{s},0,-x_{s})$ $\displaystyle=$
$\displaystyle-2\zeta_{3}+G(0,0,0;y_{s})-2G(1,0,0;y_{s}),$ $\displaystyle
f_{26}^{(2,4)}(x_{s},0,-x_{s})$ $\displaystyle=$
$\displaystyle\frac{5}{2}\zeta_{4}+\zeta_{2}\left[G(0,0;y_{s})-2G(1,0;y_{s})\right]-G(0,0,1,0;y_{s})-2G(0,1,0,0;y_{s})$
$\displaystyle-3G(1,0,0,0;y_{s})+2G(1,0,1,0;y_{s})+4G(1,1,0,0;y_{s}),$
$\displaystyle f_{28}^{(2,4)}(x_{s},0,-x_{s})$ $\displaystyle=$
$\displaystyle-9\zeta_{4}-8\zeta_{3}G(0;y_{s})-4\zeta_{2}G(0,0;y_{s})+8G(0,-1,0,0;y_{s})-16G(0,0,-1,0;y_{s})$
$\displaystyle-G(0,0,0,0;y_{s})-4G(0,0,1,0;y_{s})+14G(0,1,0,0;y_{s}),$
$\displaystyle f_{29}^{(2,4)}(x_{s},0,-x_{s})$ $\displaystyle=$
$\displaystyle-5\zeta_{4}-2\zeta_{2}\left[G(0,0;y_{s})-2G(1,0;y_{s})\right]+2G(0,0,1,0;y_{s})+4G(0,1,0,0;y_{s})$
$\displaystyle+6G(1,0,0,0;y_{s})-4G(1,0,1,0;y_{s})-8G(1,1,0,0;y_{s}),$
$\displaystyle f_{14}^{(2,3)}(0,x_{s},-x_{s})$ $\displaystyle=$
$\displaystyle-6\zeta_{3}-2\zeta_{2}G(0;y_{s})-12G(0,-1,0;y_{s})+6G(0,0,0;y_{s})-4G(0,1,0;y_{s})$
$\displaystyle+4G(1,0,0;y_{s}),$ $\displaystyle
f_{16}^{(2,2)}(0,x_{s},-x_{s})$ $\displaystyle=$
$\displaystyle-\frac{1}{2}G(0,0;y_{s}),$ $\displaystyle
f_{16}^{(2,3)}(0,x_{s},-x_{s})$ $\displaystyle=$
$\displaystyle-\frac{\zeta_{3}}{2}+\frac{\zeta_{2}}{2}G(0;y_{s})+3G(0,-1,0;y_{s})-\frac{1}{2}G(0,0,0;y_{s})+G(0,1,0;y_{s})$
$\displaystyle-3G(1,0,0;y_{s}),$ $\displaystyle
f_{23}^{(2,3)}(0,x_{s},-x_{s})$ $\displaystyle=$
$\displaystyle-2\zeta_{3}+G(0,0,0;y_{s})-2G(1,0,0;y_{s}),$ $\displaystyle
f_{26}^{(2,4)}(0,x_{s},-x_{s})$ $\displaystyle=$
$\displaystyle-6\zeta_{4}-4\zeta_{3}\left[G(0;y_{s})-2G(1;y_{s})\right]+2G(0,0,0,0;y_{s})-4G(0,1,0,0;y_{s})$
(69) $\displaystyle-4G(1,0,0,0;y_{s})+8G(1,1,0,0;y_{s}).$

The other non-zero integrals can be obtained by using the relations among the
master integrals in eqs.(50) and (51). Otherwise, they are either zero or not
needed in the evaluations of the five master integrals with square roots. In
eq.(69), $\zeta_{n}$ is the usual Riemann zeta function $\zeta(n)$.

##  Appendix C One-loop master integrals

Let us define the one-loop master integrals that are needed for the one-loop
helicity amplitudes and the UV renormalisation counterterms. The family of the
one-loop master integrals in $d=4-2\epsilon$ is

$\displaystyle I_{a_{1},a_{2},a_{3},a_{4}}^{(1)}(s_{ij},s_{jk},s_{ik})$
$\displaystyle=\frac{m_{f}^{2\epsilon}}{\Gamma(1+\epsilon)}\int\frac{{\mathrm{d}}^{d}\ell_{1}}{i\pi^{\frac{d}{2}}}\frac{1}{D_{1}^{a_{1}}D_{2}^{a_{2}}D_{3}^{a_{3}}D_{4}^{a_{4}}},\quad
a_{i}\in\mathbb{Z}\,.$ (70)

The loop momentum is $\ell_{1}$ and $\Gamma$ is the Euler gamma function. The
inverse propagators are

$\displaystyle D_{1}$ $\displaystyle=$
$\displaystyle\ell_{1}^{2}-m_{f}^{2}+i0^{+},\quad
D_{2}=(\ell_{1}+p_{i})^{2}-m_{f}^{2}+i0^{+},$ $\displaystyle D_{3}$
$\displaystyle=$
$\displaystyle(\ell_{1}+p_{i}+p_{j})^{2}-m_{f}^{2}+i0^{+},\quad
D_{4}=(\ell_{1}+p_{i}+p_{j}+p_{k})^{2}-m_{f}^{2}+i0^{+}.$ (71)

The canonical basis with uniform transcendental weight properties is Caron-
Huot:2014lda

$\displaystyle f_{1}^{(1)}$ $\displaystyle=$
$\displaystyle-2m_{f}^{2}I^{(1)}_{0,0,0,3},$ $\displaystyle f_{2}^{(1)}$
$\displaystyle=$
$\displaystyle-\epsilon\sqrt{s_{ij}(s_{ij}-4m_{f}^{2})}I^{(1)}_{1,0,2,0},\quad
f_{3}^{(1)}=-\epsilon\sqrt{s_{jk}(s_{jk}-4m_{f}^{2})}I^{(1)}_{0,1,0,2},$
$\displaystyle f_{4}^{(1)}$ $\displaystyle=$
$\displaystyle\epsilon^{2}s_{ij}I^{(1)}_{1,1,1,0},\quad
f_{5}^{(1)}=\epsilon^{2}s_{jk}I^{(1)}_{1,1,0,1},$ $\displaystyle f_{6}^{(1)}$
$\displaystyle=$
$\displaystyle\frac{\epsilon^{2}}{2}\sqrt{s_{ij}s_{jk}\left(s_{ij}s_{jk}-4m_{f}^{2}(s_{ij}+s_{jk})\right)}I^{(1)}_{1,1,1,1}.$
(72)

Again, let us define the vector
$\vec{f}^{(1)}=(f_{1}^{(1)},\cdots,f_{6}^{(1)})$, which satisfies the
$\epsilon$-form of differential equations

$\displaystyle{\mathrm{d}}\vec{f}^{(1)}$ $\displaystyle=$
$\displaystyle\epsilon{\mathrm{d}}A^{(1)}\vec{f}^{(1)}.$ (73)

The matrix $A^{(1)}$ is independent of the dimensional regulator $\epsilon$.
Unlike the two-loop, in this case there are only three square roots
$\sqrt{s_{ij}(s_{ij}-4m_{f}^{2})}$, $\sqrt{s_{jk}(s_{jk}-4m_{f}^{2})}$, and
$\sqrt{s_{ij}s_{jk}(s_{ij}s_{jk}-4m_{f}^{2}(s_{ij}+s_{jk}))}$, that can be
rationalised simultaneously with the parametrisation eq.(37). In other words,
the one-loop master integrals can be expressed through GPLs up to any order in
$\epsilon$. Only $10$ $\mathrm{dlog}$ one-forms as already defined in eq.(40)
appear in the matrix $A^{(1)}$:

$\displaystyle{\mathrm{d}}A^{(1)}$ $\displaystyle=$
$\displaystyle\left(\begin{array}[]{cccccc}0&0&0&0&0&0\\\
-{\mathrm{d}}g_{1}&{\mathrm{d}}g_{2}&0&0&0&0\\\
{\mathrm{d}}g_{5}&0&\frac{{\mathrm{d}}g_{4}}{3}+\frac{{\mathrm{d}}g_{6}}{3}&0&0&0\\\
0&-{\mathrm{d}}g_{1}&0&0&0&0\\\ 0&0&{\mathrm{d}}g_{5}&0&0&0\\\
0&{\mathrm{d}}g_{12}&-2{\mathrm{d}}g_{14}&{\mathrm{d}}g_{13}&{\mathrm{d}}g_{13}&\frac{2{\mathrm{d}}g_{10}}{3}+\frac{{\mathrm{d}}g_{15}}{3}\\\
\end{array}\right)\,.$ (80)

The series expansion over the dimensional regulator $\epsilon$ in
$\vec{f}^{(1)}=\sum_{\omega=0}{\epsilon^{\omega}\vec{f}^{(1,\omega)}}$ or
equivalently
$f_{n}^{(1)}=\sum_{\omega=0}{\epsilon^{\omega}f_{n}^{(1,\omega)}}$ enables the
differential equation system

$\displaystyle{\mathrm{d}}\vec{f}^{(1,\omega+1)}$ $\displaystyle=$
$\displaystyle{\mathrm{d}}A^{(1)}\vec{f}^{(1,\omega)},\quad\omega\in\mathbb{N}$
(81)

to be solved iteratively in GPLs with the boundary condition at
$s_{ij}=s_{jk}=s_{ik}=0$. The analytic continuation should be performed as
explained in appendix A. We provide the solutions up to weight $\omega=3$ in
the ancillary files.

##  Appendix D Helicity amplitudes of the one-loop $W^{\pm}$ contribution

For completeness, we also present the one-loop helicity amplitudes with
$W^{\pm}$ boson. The one-loop helicity amplitudes with the $W^{\pm}$ boson and
its Goldstone boson from the sum of the $162$ Feynman graphs are:

$\displaystyle-i\mathcal{M}_{++++}^{(0,0,W)}$ $\displaystyle=$
$\displaystyle-12\alpha^{2}\left[1-4\mathop{\sum_{(i,j,k)=(s,t,u),}}_{(t,u,s),(u,s,t)}{\frac{f_{6}^{(1,2)}(x_{i},x_{j},x_{k})}{\sqrt{x_{i}x_{j}(x_{i}x_{j}+4x_{k})}}}\right],$
$\displaystyle-i\mathcal{M}^{(0,0,W)}_{-+++}$ $\displaystyle=$
$\displaystyle-12\alpha^{2}\left[1+2\left(x_{s}^{-1}+x_{t}^{-1}+x_{u}^{-1}\right)\mathop{\sum_{(i,j,k)=(s,t,u),}}_{(t,u,s),(u,s,t)}{f_{4}^{(1,2)}(x_{i},x_{j},x_{k})}\right.$
$\displaystyle\qquad\qquad\qquad\left.-\mathop{\sum_{(i,j,k)=(s,t,u),}}_{(t,u,s),(u,s,t)}{r^{(1)}_{8}(x_{i},x_{j},x_{k})\frac{f_{6}^{(1,2)}(x_{i},x_{j},x_{k})}{\sqrt{x_{i}x_{j}(x_{i}x_{j}+4x_{k})}}}\right],$
$\displaystyle-i\mathcal{M}^{(0,0,W)}_{--++}$ $\displaystyle=$
$\displaystyle-12\alpha^{2}\left[-1+r_{9}^{(1)}(x_{t},x_{u},x_{s})\frac{f_{2}^{(1,1)}(x_{t},x_{u},x_{s})}{\sqrt{x_{t}(x_{t}-4)}}+r_{9}^{(1)}(x_{u},x_{t},x_{s})\frac{f_{2}^{(1,1)}(x_{u},x_{t},x_{s})}{\sqrt{x_{u}(x_{u}-4)}}\right.$
(82)
$\displaystyle-\left(r_{10}^{(1)}(x_{s},x_{t},x_{u})+\frac{5}{3}\right)\left(f_{4}^{(1,2)}(x_{t},x_{u},x_{s})+f_{4}^{(1,2)}(x_{u},x_{s},x_{t})\right)$
$\displaystyle-\frac{4}{3}(x_{s}-1)(x_{s}-3)\left(\frac{f_{6}^{(1,2)}(x_{s},x_{t},x_{u})}{\sqrt{x_{s}x_{t}(x_{s}x_{t}+4x_{u})}}+\frac{f_{6}^{(1,2)}(x_{u},x_{s},x_{t})}{\sqrt{x_{s}x_{u}(x_{s}x_{u}+4x_{t})}}\right)$
$\displaystyle\left.-r_{12}^{(1)}(x_{s},x_{t},x_{u})\frac{f_{6}^{(1,2)}(x_{t},x_{u},x_{s})}{\sqrt{x_{t}x_{u}(x_{t}x_{u}+4x_{s})}}\right],$

where

$\displaystyle r_{12}^{(1)}(x_{i},x_{j},x_{k})$ $\displaystyle=$
$\displaystyle
r_{11}^{(1)}(x_{i},x_{j},x_{k})+\frac{4x_{i}^{2}-10x_{i}-5x_{j}x_{k}}{3}.$
(83)

Here, we have defined $x_{s}=\frac{s}{m_{W}^{2}},x_{t}=\frac{t}{m_{W}^{2}},$
and $x_{u}=\frac{u}{m_{W}^{2}}$ with $m_{W}$ being the mass of the $W^{\pm}$
boson. It is interesting to note that, according to eqs.(56) and (82), we have
the simple relation
$\mathcal{M}_{++++}^{(0,0,W)}/\mathcal{M}_{++++}^{(0,0,f)}=\mathcal{M}_{-+++}^{(0,0,W)}/\mathcal{M}_{-+++}^{(0,0,f)}=-3/(2N_{c,f}Q_{f}^{4})$.
However, such a simple relation does not hold for $\mathcal{M}_{--++}$, though
the form of $\mathcal{M}^{(0,0,W)}_{--++}$ is very close to
$\mathcal{M}^{(0,0,f)}_{--++}$. The one-loop amplitudes provided in ref.
Bardin:2009gq  are in terms of one-loop scalar integrals rather than a UT
basis. Our expressions expressed in terms of a UT basis is however shorter
than ref. Bardin:2009gq . The squared amplitude has also been numerically
checked with MadLoop Hirschi:2011pa ; Alwall:2014hca  for a few random phase
space points.

##  Appendix E One- and two-loop helicity amplitudes in the low-energy limit

In this appendix, we collect the one- and two-loop helicity amplitudes in the
low energy (LE) limit $m_{f}^{2}\gg s,-t,-u$ for a given fermion $f$, which
has been derived in ref. Martin:2003gb  from the Euler-Heisenberg Lagrangian.
We express them in terms of the Mandelstam variables instead of the spinor-
helicity formalism presented in ref. Martin:2003gb , with
$x_{s}=s/m_{f}^{2},x_{t}=t/m_{f}^{2}$ and $x_{u}=u/m_{f}^{2}$. Up to a global
phase, we have the one-loop helicity amplitudes in the LE limit as:

$\displaystyle-i\mathcal{M}^{(0,0,f)}_{++++,{\rm LE}}$ $\displaystyle=$
$\displaystyle-\frac{1}{15}N_{c,f}Q_{f}^{4}\alpha^{2}\left(x_{s}^{2}+x_{t}^{2}+x_{u}^{2}\right),$
$\displaystyle-i\mathcal{M}^{(0,0,f)}_{-+++,{\rm LE}}$ $\displaystyle=$
$\displaystyle 0,$ $\displaystyle-i\mathcal{M}^{(0,0,f)}_{--++,{\rm LE}}$
$\displaystyle=$
$\displaystyle\frac{11}{45}N_{c,f}Q_{f}^{4}\alpha^{2}x_{s}^{2}.$ (84)

The two-loop QCD amplitudes in the same limit are

$\displaystyle\mathcal{M}^{(1,0,f)}_{++++,{\rm LE}}$ $\displaystyle=$
$\displaystyle\frac{25}{4}\frac{\alpha_{s}}{\pi}C_{F,f}\mathcal{M}^{(0,0,f)}_{++++,{\rm
LE}},$ $\displaystyle\mathcal{M}^{(1,0,f)}_{-+++,{\rm LE}}$ $\displaystyle=$
$\displaystyle 0,$ $\displaystyle\mathcal{M}^{(1,0,f)}_{--++,{\rm LE}}$
$\displaystyle=$
$\displaystyle\frac{1955}{396}\frac{\alpha_{s}}{\pi}C_{F,f}\mathcal{M}^{(0,0,f)}_{--++,{\rm
LE}}.$ (85)

The two-loop QED amplitudes can be obtained from the relation

$\displaystyle\mathcal{M}^{(0,1,f)}_{\lambda_{1}\lambda_{2}\lambda_{3}\lambda_{4},{\rm
LE}}=(\alpha
Q_{f}^{2})/(\alpha_{s}C_{F,f})\mathcal{M}^{(1,0,f)}_{\lambda_{1}\lambda_{2}\lambda_{3}\lambda_{4},{\rm
LE}}.$ (86)

Similarly, the one-loop helicity amplitudes for the $W^{\pm}$ boson loop in
the LE limit $m_{W}^{2}\gg s,-t,-u$ with $x_{s}=s/m_{W}^{2},x_{t}=t/m_{W}^{2}$
and $x_{u}=u/m_{W}^{2}$ are:

$\displaystyle-i\mathcal{M}_{++++,{\rm LE}}^{(0,0,W)}$ $\displaystyle=$
$\displaystyle\frac{1}{10}\alpha^{2}\left(x_{s}^{2}+x_{t}^{2}+x_{u}^{2}\right),$
$\displaystyle-i\mathcal{M}_{-+++,{\rm LE}}^{(0,0,W)}$ $\displaystyle=$
$\displaystyle 0,$ $\displaystyle-i\mathcal{M}_{--++,{\rm LE}}^{(0,0,W)}$
$\displaystyle=$ $\displaystyle\frac{14}{5}\alpha^{2}x_{s}^{2}.$ (87)

## References

  * (1) W. Heisenberg, Remarks on the Dirac theory of the positron, Z. Phys. 90 no. 3-4, (1934) 209–231. [Erratum: Z.Phys. 92, 692–692 (1934)].
  * (2) H. Euler and B. Kockel, The scattering of light by light in Dirac’s theory, Naturwiss. 23 no. 15, (1935) 246–247.
  * (3) H. Euler, On the scattering of light by light according to Dirac’s theory, Annalen Phys. 26 no. 5, (1936) 398–448.
  * (4) W. Heisenberg and H. Euler, Consequences of Dirac’s theory of positrons, Z. Phys. 98 no. 11-12, (1936) 714–732, arXiv:physics/0605038.
  * (5) R. Karplus and M. Neuman, The scattering of light by light, Phys. Rev. 83 (1951) 776–784.
  * (6) D. d’Enterria and G. G. da Silveira, Observing light-by-light scattering at the Large Hadron Collider, Phys. Rev. Lett. 111 (2013) 080405, arXiv:1305.7142 [hep-ph]. [Erratum: Phys.Rev.Lett. 116, 129901 (2016)].
  * (7) S. Knapen, T. Lin, H. K. Lou, and T. Melia, Searching for Axionlike Particles with Ultraperipheral Heavy-Ion Collisions, Phys. Rev. Lett. 118 no. 17, (2017) 171801, arXiv:1607.06083 [hep-ph].
  * (8) D. d’Enterria, M. A. Tamlihat, L. Schoeffel, H.-S. Shao, and Y. Tayalati, Collider constraints on massive gravitons coupling to photons, arXiv:2306.15558 [hep-ph].
  * (9) S. Atag, S. C. Inan, and I. Sahin, Extra dimensions in ${\gamma\gamma\rightarrow\gamma\gamma}$ process at the CERN-LHC, JHEP 09 (2010) 042, arXiv:1005.4792 [hep-ph].
  * (10) J. Ellis, N. E. Mavromatos, and T. You, Light-by-Light Scattering Constraint on Born-Infeld Theory, Phys. Rev. Lett. 118 no. 26, (2017) 261802, arXiv:1703.08450 [hep-ph].
  * (11) R. Horvat, D. Latas, J. Trampetić, and J. You, Light-by-Light Scattering and Spacetime Noncommutativity, Phys. Rev. D 101 no. 9, (2020) 095035, arXiv:2002.01829 [hep-ph].
  * (12) K.-m. Cheung, Diphoton signals for low scale gravity in extra dimensions, Phys. Rev. D 61 (2000) 015005, arXiv:hep-ph/9904266.
  * (13) H. Davoudiasl, Gamma gamma —$>$ gamma gamma as a test of weak scale quantum gravity at the NLC, Phys. Rev. D 60 (1999) 084022, arXiv:hep-ph/9904425.
  * (14) M. Greiner, M. Vidovic, and G. Soff, Electromagnetic production of Higgs bosons, SUSY particles, glueballs and mesons in ultrarelativistic heavy ion collisions, Phys. Rev. C 47 (1993) 2288–2298.
  * (15) D. d’Enterria and H.-S. Shao, Observing true tauonium via two-photon fusion at e+e- and hadron colliders, Phys. Rev. D 105 no. 9, (2022) 093008, arXiv:2202.02316 [hep-ph].
  * (16) D. d’Enterria and H.-S. Shao, Prospects for ditauonium discovery at colliders, Phys. Lett. B 842 (2023) 137960, arXiv:2302.07365 [hep-ph].
  * (17) ATLAS Collaboration, M. Aaboud et al., Evidence for light-by-light scattering in heavy-ion collisions with the ATLAS detector at the LHC, Nature Phys. 13 no. 9, (2017) 852–858, arXiv:1702.01625 [hep-ex].
  * (18) CMS Collaboration, A. M. Sirunyan et al., Evidence for light-by-light scattering and searches for axion-like particles in ultraperipheral PbPb collisions at $\sqrt{s_{\mathrm{NN}}}=$ 5.02 TeV, Phys. Lett. B 797 (2019) 134826, arXiv:1810.04602 [hep-ex].
  * (19) ATLAS Collaboration, G. Aad et al., Observation of light-by-light scattering in ultraperipheral Pb+Pb collisions with the ATLAS detector, Phys. Rev. Lett. 123 no. 5, (2019) 052001, arXiv:1904.03536 [hep-ex].
  * (20) ATLAS Collaboration, G. Aad et al., Measurement of light-by-light scattering and search for axion-like particles with 2.2 nb-1 of Pb+Pb data with the ATLAS detector, JHEP 03 (2021) 243, arXiv:2008.05355 [hep-ex]. [Erratum: JHEP 11, 050 (2021)].
  * (21) L. A. Harland-Lang, M. Tasevsky, V. A. Khoze, and M. G. Ryskin, A new approach to modelling elastic and inelastic photon-initiated production at the LHC: SuperChic 4, Eur. Phys. J. C 80 no. 10, (2020) 925, arXiv:2007.12704 [hep-ph].
  * (22) H.-S. Shao and D. d’Enterria, gamma-UPC: automated generation of exclusive photon-photon processes in ultraperipheral proton and nuclear collisions with varying form factors, JHEP 09 (2022) 248, arXiv:2207.03012 [hep-ph].
  * (23) H.-S. Shao, HELAC-Onia: An automatic matrix element generator for heavy quarkonium physics, Comput. Phys. Commun. 184 (2013) 2562–2570, arXiv:1212.5293 [hep-ph].
  * (24) H.-S. Shao, HELAC-Onia 2.0: an upgraded matrix-element and event generator for heavy quarkonium physics, Comput. Phys. Commun. 198 (2016) 238–259, arXiv:1507.03435 [hep-ph].
  * (25) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli, and M. Zaro, The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations, JHEP 07 (2014) 079, arXiv:1405.0301 [hep-ph].
  * (26) M. Reuter, M. G. Schmidt, and C. Schubert, Constant external fields in gauge theory and the spin 0, 1/2, 1 path integrals, Annals Phys. 259 (1997) 313–365, arXiv:hep-th/9610191.
  * (27) D. Fliegner, M. Reuter, M. G. Schmidt, and C. Schubert, The Two loop Euler-Heisenberg Lagrangian in dimensional renormalization, Theor. Math. Phys. 113 (1997) 1442–1451, arXiv:hep-th/9704194.
  * (28) C. Schubert, Perturbative quantum field theory in the string inspired formalism, Phys. Rept. 355 (2001) 73–234, arXiv:hep-th/0101036.
  * (29) L. C. Martin, C. Schubert, and V. M. Villanueva Sandoval, On the low-energy limit of the QED N photon amplitudes, Nucl. Phys. B 668 (2003) 335–344, arXiv:hep-th/0301022.
  * (30) Z. Bern, A. De Freitas, L. J. Dixon, A. Ghinculov, and H. L. Wong, QCD and QED corrections to light by light scattering, JHEP 11 (2001) 031, arXiv:hep-ph/0109079.
  * (31) T. Binoth, E. W. N. Glover, P. Marquard, and J. J. van der Bij, Two loop corrections to light by light scattering in supersymmetric QED, JHEP 05 (2002) 060, arXiv:hep-ph/0202266.
  * (32) S. Caron-Huot and J. M. Henn, Iterative structure of finite loop integrals, JHEP 06 (2014) 114, arXiv:1404.2922 [hep-th].
  * (33) X. Xu and L. L. Yang, Towards a new approximation for pair-production and associated-production of the Higgs boson, JHEP 01 (2019) 211, arXiv:1810.12002 [hep-ph].
  * (34) F. Maltoni, M. K. Mandal, and X. Zhao, Top-quark effects in diphoton production through gluon fusion at next-to-leading order in QCD, Phys. Rev. D 100 no. 7, (2019) 071501, arXiv:1812.08703 [hep-ph].
  * (35) G. Wang, Y. Wang, X. Xu, Y. Xu, and L. L. Yang, Efficient computation of two-loop amplitudes for Higgs boson pair production, Phys. Rev. D 104 no. 5, (2021) L051901, arXiv:2010.15649 [hep-ph].
  * (36) T. Peraro and L. Tancredi, Tensor decomposition for bosonic and fermionic scattering amplitudes, Phys. Rev. D 103 no. 5, (2021) 054042, arXiv:2012.00820 [hep-ph].
  * (37) P. Nogueira, Automatic Feynman Graph Generation, J. Comput. Phys. 105 (1993) 279–289.
  * (38) T. Hahn, Generating Feynman diagrams and amplitudes with FeynArts 3, Comput. Phys. Commun. 140 (2001) 418–431, arXiv:hep-ph/0012260.
  * (39) P. Maierhöfer, J. Usovitsch, and P. Uwer, Kira—A Feynman integral reduction program, Comput. Phys. Commun. 230 (2018) 99–112, arXiv:1705.05610 [hep-ph].
  * (40) 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, arXiv:2008.06494 [hep-ph].
  * (41) T. Peraro, FiniteFlow: multivariate functional reconstruction using finite fields and dataflow graphs, JHEP 07 (2019) 031, arXiv:1905.08019 [hep-ph].
  * (42) R. N. Lee, Presenting LiteRed: a tool for the Loop InTEgrals REDuction, arXiv:1212.2685 [hep-ph].
  * (43) A. V. Kotikov, Differential equations method: New technique for massive Feynman diagrams calculation, Phys. Lett. B 254 (1991) 158–164.
  * (44) E. Remiddi, Differential equations for Feynman graph amplitudes, Nuovo Cim. A 110 (1997) 1435–1452, arXiv:hep-th/9711188.
  * (45) T. Gehrmann and E. Remiddi, Differential equations for two loop four point functions, Nucl. Phys. B 580 (2000) 485–518, arXiv:hep-ph/9912329.
  * (46) M. Argeri and P. Mastrolia, Feynman Diagrams and Differential Equations, Int. J. Mod. Phys. A 22 (2007) 4375–4436, arXiv:0707.4037 [hep-ph].
  * (47) K.-T. Chen, Iterated path integrals, Bull. Am. Math. Soc. 83 (1977) 831–879.
  * (48) J. M. Henn, Multiloop integrals in dimensional regularization made simple, Phys. Rev. Lett. 110 (2013) 251601, arXiv:1304.1806 [hep-th].
  * (49) A. B. Goncharov, Multiple polylogarithms, cyclotomy and modular complexes, 2011.
  * (50) A. B. Goncharov, M. Spradlin, C. Vergu, and A. Volovich, Classical Polylogarithms for Amplitudes and Wilson Loops, Phys. Rev. Lett. 105 (2010) 151605, arXiv:1006.5703 [hep-th].
  * (51) J. Vollinga and S. Weinzierl, Numerical evaluation of multiple polylogarithms, Comput. Phys. Commun. 167 (2005) 177, arXiv:hep-ph/0410259.
  * (52) Y. Wang, L. L. Yang, and B. Zhou, FastGPL: a C++ library for fast evaluation of generalized polylogarithms, arXiv:2112.04122 [hep-ph].
  * (53) L. Naterop, A. Signer, and Y. Ulrich, handyG —Rapid numerical evaluation of generalised polylogarithms in Fortran, Comput. Phys. Commun. 253 (2020) 107165, arXiv:1909.01656 [hep-ph].
  * (54) D. Chicherin and V. Sotnikov, Pentagon Functions for Scattering of Five Massless Particles, JHEP 20 (2020) 167, arXiv:2009.07803 [hep-ph].
  * (55) C. Duhr, H. Gangl, and J. R. Rhodes, From polygons and symbols to polylogarithmic functions, JHEP 10 (2012) 075, arXiv:1110.0458 [math-ph].
  * (56) X. Liu and Y.-Q. Ma, AMFlow: A Mathematica package for Feynman integrals computation via auxiliary mass flow, Comput. Phys. Commun. 283 (2023) 108565, arXiv:2201.11669 [hep-ph].
  * (57) S. Badger, E. Chaubey, H. B. Hartanto, and R. Marzucca, Two-loop leading colour QCD helicity amplitudes for top quark pair production in the gluon fusion channel, JHEP 06 (2021) 163, arXiv:2102.13450 [hep-ph].
  * (58) C. Duhr and F. Dulat, PolyLogTools — polylogs for the masses, JHEP 08 (2019) 135, arXiv:1904.07279 [hep-th].
  * (59) J. A. M. Vermaseren, New features of FORM, arXiv:math-ph/0010025.
  * (60) B. Ruijl, T. Ueda, and J. Vermaseren, FORM version 4.2, arXiv:1707.06453 [hep-ph].
  * (61) D. Bardin, L. Kalinovskaya, and E. Uglov, Standard Model light-by-light scattering in SANC: analytic and numeric evaluation, Phys. Atom. Nucl. 73 (2010) 1878–1888, arXiv:0911.5634 [hep-ph].
  * (62) V. Hirschi, R. Frederix, S. Frixione, M. V. Garzelli, F. Maltoni, and R. Pittau, Automation of one-loop QCD corrections, JHEP 05 (2011) 044, arXiv:1103.0621 [hep-ph].
  * (63) M. Heller and A. von Manteuffel, MultivariateApart: Generalized partial fractions, Comput. Phys. Commun. 271 (2022) 108174, arXiv:2101.08283 [cs.SC].
  * (64) Z. Capatti, V. Hirschi, A. Pelloni, and B. Ruijl, Local Unitarity: a representation of differential cross-sections that is locally free of infrared singularities at any order, JHEP 04 (2021) 104, arXiv:2010.01068 [hep-ph].
  * (65) Z. Capatti, V. Hirschi, and B. Ruijl, Local unitarity: cutting raised propagators and localising renormalisation, JHEP 10 (2022) 120, arXiv:2203.11038 [hep-ph].
  * (66) Z. Bern and D. A. Kosower, Efficient calculation of one loop QCD amplitudes, Phys. Rev. Lett. 66 (1991) 1669–1672.
  * (67) Z. Bern and D. A. Kosower, The Computation of loop amplitudes in gauge theories, Nucl. Phys. B 379 (1992) 451–561.
  * (68) M. J. Strassler, Field theory without Feynman diagrams: One loop effective actions, Nucl. Phys. B 385 (1992) 145–184, arXiv:hep-ph/9205205.
  * (69) A. A H, E. Chaubey, M. Fraaije, V. Hirschi, and H.-S. Shao, Light-by-Light Scattering at Next-to-Leading Order QCD and QED, arXiv:2312.YYYYY [hep-ph].