CERN-TH-2023-090, PSI-PR-23-16, ZU-TH 23/23

Two-Loop QCD Corrections for Three-Photon Production at Hadron Colliders

Samuel Abreu1,2, Giuseppe De Laurentis3, Harald Ita3,4, Maximillian Klinkert5,
Ben Page1 and Vasily Sotnikov6,7

1 CERN, Theoretical Physics Department, CH-1211 Geneva 23, Switzerland

2 Higgs Centre for Theoretical Physics, School of Physics and Astronomy,

The University of Edinburgh, Edinburgh EH9 3FD, Scotland, UK

3 Paul Scherrer Institut, CH-5232 Villigen PSI, Switzerland

4 ICS, University of Zurich, Winterthurerstrasse 190, 8057 Zurich, Switzerland

5 Physikalisches Institut, Albert-Ludwigs-Universität Freiburg,

Hermann-Herder. Str. 3, D-79104 Freiburg, Germany

6 Physik-Institut, University of Zurich,

Winterthurerstrasse 190, 8057 Zurich, Switzerland

7 Department of Physics and Astronomy, Michigan State University,

567 Wilson Road, East Lansing, MI 48824, USA

## Abstract

We complete the computation of the two-loop helicity amplitudes for the
production of three photons at hadron colliders, including all contributions
beyond the leading-color approximation. We reconstruct the analytic form of
the amplitudes from numerical finite-field samples obtained with the numerical
unitarity method. This method requires as input surface terms for all relevant
five-point non-planar integral topologies, which we obtain by solving the
associated syzygy problem in embedding space. The numerical samples are used
to constrain compact spinor-helicity ansätze, which are optimized by taking
advantage of the known one-loop analytic structure. We make our analytic
results available in a public C++ library, which is suitable for immediate
phenomenological applications. We estimate that the inclusion of the
subleading-color contributions will decrease the size of the two-loop
corrections by about $30\%$ to $50\%$ compared to the results in the leading-
color approximation.

###### Contents

  1. 1 Introduction
  2. 2 Notation and Conventions
  3. 3 Calculation
    1. 3.1 Overview
    2. 3.2 Surface Terms
      1. 3.2.1 Embedding-Space Formalism
      2. 3.2.2 IBP-Generating Vectors In Embedding Space
      3. 3.2.3 Computing IBP-Generating Vectors
      4. 3.2.4 Computing Surface Terms
    3. 3.3 Analytic Reconstruction In Spinor Helicity Formalism
  4. 4 Results
    1. 4.1 Efficiency Of Analytic Reconstruction
    2. 4.2 Validation
    3. 4.3 Ancillary Files
    4. 4.4 Numerical Evaluation
  5. 5 Subleading Color Corrections To Hard Function
  6. 6 Conclusion
  7. A IR Renormalization
    1. A.1 Conventions In Catani Scheme
    2. A.2 IR-Subtraction Scheme Change
  8. B Reference Values

##  1 Introduction

With the rapid development of precision studies within the physics program at
the Large Hadron Collider (LHC), there is a growing need for precise theory
predictions for many Standard Model processes. The knowledge of higher-order
radiative corrections in the strong coupling constant is essential for the
correct interpretation of the associated experimental measurements, and
provides better control over theoretical uncertainties. While many two-to-two
reactions are now known at next-to-next-to-leading order (NNLO) in
perturbation theory, NNLO theoretical predictions for two-to-three processes
constitute the current state of the art.

The study of triphoton final states at hadron colliders such as the LHC offers
excellent opportunities for several interesting investigations. For instance,
it allows to explore the consequences of potentially anomalous gauge and Higgs
couplings [1, 2, 3, 4, 5]. Furthermore, the study of triphoton production is
an important irreducible background process for the associated production of a
photon with a beyond-the-Standard-Model particle that subsequently decays into
a pair of photons [6, 7, 8]. Similar to diphoton production [9, 10, 11],
triphoton production exhibits large next-to-leading order (NLO) and NNLO
corrections [12, 13, 14]. Accordingly, NLO predictions significantly deviate
from data [15], and NNLO is the first perturbative order providing reliable
results [13, 12] (see also the related work of refs. [16, 17]). Thus, good
control of NNLO QCD corrections to this process is of great importance.

While diphoton production has been known at NNLO QCD accuracy for more than a
decade [9, 10], first results for triphoton production at the same level of
accuracy have only recently been obtained [12, 13]. This is due to the highly
challenging nature of NNLO computations for five-particle scattering
processes. It has been demonstrated that existing frameworks for the
subtraction of infrared divergences at NNLO are, in principle, capable of
handling arbitrary production process [13, 18, 19, 20]. However, the
calculation of two-loop amplitudes for five-particle processes represents the
current state of the art, and needs to be addressed on a case-by-case basis.
While NNLO QCD corrections for most massless five-particle processes have
already been considered [18, 19, 21, 12, 13, 22, 23, 24], in most of these
studies, including the NNLO QCD corrections to triphoton production [13, 12],
only the leading-color approximation of the double-virtual corrections has
been employed. The first complete cross-section calculation that does not
employ the leading-color approximation was performed in ref. [24], and a
number of complete two-loop five-point amplitudes are also available [25, 23,
24].

Subleading-color corrections are more challenging to compute because they
include non-planar Feynman graphs, which are notoriously more challenging to
handle. In pure QCD, all non-planar contributions vanish in the limit of a
large number of colors ${N_{c}}$ [26]. This still holds if the number of
fermions flavors ${N_{f}}$ is considered to be of the order of ${N_{c}}$, i.e.
when the ratio ${N_{f}}/{N_{c}}$ is kept constant as ${N_{c}}$ approaches
infinity. The latter variant of the leading-color approximation is typically
phenomenologically justified, provided errors of about 10% in the approximated
contributions are deemed acceptable. On the other hand, terms originating from
photons coupling to closed fermion loops involve nonplanar diagrams. While
these terms are still suppressed in the formal ${N_{c}}\to\infty$ limit, one
might be concerned that their contribution is not suppressed numerically.

The goal of this work is to compute analytic expressions for all two-loop
five-point amplitudes that contribute to the NNLO corrections for triphoton
production at hadron colliders in massless QCD (ignoring top loops), and
provide an efficient numerical implementation for use in phenomenological
studies. Our calculation follows the multi-loop numerical unitarity approach
[27, 28, 29, 30] as implemented in Caravel [31]. In numerical unitarity,
amplitudes are reduced to a set of master integrals by matching numerical
evaluations of generalized unitarity cuts to a parametrization of the loop
integrands. Analytic expressions can then be reconstructed using multivariate
functional reconstruction techniques [32, 33] (see also refs. [34, 35, 36, 37,
38] for related developments). In order to perform our calculation within this
framework, we make a number of theoretical developments.

Firstly, we develop a new approach to the problem of parametrising the
integrand so that as many terms as possible vanish upon integration, due to
integration-by-parts (IBP) identities [39, 40]. These terms are commonly
referred to as _surface terms_ , and it is a highly non-trivial problem to
construct them in a form that enables efficient numerical evaluations. In this
work, we present a novel method of deriving numerically efficient
representations of surface terms. We achieve this by solving an associated
syzygy problem [41, 27, 42] which, taking inspiration from refs. [43, 44, 45],
is formulated in so-called embedding space. This allows us to identify a
simpler _homogeneous_ syzygy system, whose solutions we lift to full syzygies
via linear algebra. By performing this calculation on a numerical phase-space
point, we are able to construct “skeleton” syzygies, that allow us to
numerically determine the full set of syzygies and surface terms phase-space
point by phase-space point in an efficient manner. This significantly reduces
the expression size of surface terms, allowing us to efficiently match the
integrand parametrisation to generalized unitarity cuts.

Secondly, we tackle the important problem of the large amount of samples
required to perform analytic reconstruction. Indeed, while several
improvements to the generic black-box reconstruction [32, 33] have been
explored [46, 23, 47] over the years, cutting-edge calculations (see e.g.
[47]) indicate that further developments will be important to tackle
amplitudes with an increased number of scales. For this reason, we explore new
techniques to construct ansätze whose analytic structure better exhibit the
physical properties of the scattering amplitudes, performing the analytic
reconstruction using spinor-helicity ansatz techniques [35, 36]. In contrast
to a more traditional reconstruction using a set of independent kinematic
invariants (as e.g. in [48]), this better manifests physical properties of the
amplitudes. Combined with a further optimization of the ansatz based on the
expectation that some features of the analytic structure of one-loop
amplitudes are preserved at two loops, we find a significant reduction in the
number of numerical samples that is required. In practice, we are able to
reconstruct the most complicated helicity amplitude from only about 4000
evaluations, corresponding to an order of magnitude less than what was
originally required for the reconstruction of the planar amplitudes in ref.
[48].

Alongside this paper, we provide the analytic results for the complete two-
loop triphoton production amplitudes in a collection of supplementary
material. Furthermore, in order to facilitate the applicability of our results
in phenomenological studies of triphoton production, we have implemented them
in the efficient public C++ library FivePointAmplitudes [49]. This further
allows us to analyze the important question of the impact of the subleading-
color contributions on the two-loop corrections. Our study suggests that
including these contributions will lead to a significant decrease in the size
of the two-loop corrections, reducing them by approximately 30% to 50%
compared to the results obtained in the leading-color approximation. This
effect is larger than the corrections of about 10% expected from typical
color-suppressed contributions. This confirms the concerns that nonplanar
contributions arising from the photons coupling to closed fermion loops are
not necessarily numerically suppressed. We note nevertheless that this
substantial change in the two-loop corrections should be contrasted against
the fact that the double-virtual contributions to the NNLO corrections to this
process are observed to be small [12, 13].

The paper is organized as follows. In section 2 we classify the full set of
gauge-invariant contributions to triphoton production. In section 3 we discuss
our computational approach. We review numerical unitarity, discuss our
approach to the construction of surface terms and how we construct compact
spinor-helicity ansätze. In section 4, we discuss the structure of our
results, their validation, and the format in which they are presented in
ancillary files. We also showcase the numerical performance of our C++
implementation. In section 5, we discuss the impact of the subleading-color
contributions on the double-virtual corrections. Finally, in section 6, we
present our conclusions.

##  2 Notation and Conventions

We consider the $\order{\alpha_{s}^{2}}$ corrections to the production of
three photons at hadron colliders. The loop-induced process
$gg\to\gamma\gamma\gamma$ vanishes to all orders in the combined theory of QCD
and QED due to charge-conjugation symmetry [50]. Therefore, the only
contributing partonic process is

$q(-p_{1},-h_{1})+\bar{q}(-p_{2},-h_{2})\,\to\,\gamma(p_{3},h_{3})+\gamma(p_{4},h_{4})+\gamma(p_{5},h_{5})\,,$
(1)

where $p_{i}$ and $h_{i}$ denote the momentum and the helicity of the
$i^{\text{th}}$ particle, respectively. Throughout this paper, momenta and
helicity labels are understood in the all-outgoing convention.

The process involves five massless particles. Thus, the underlying kinematic
is specified by five Mandelstam invariants, which can be chosen to be

$\displaystyle\begin{split}s_{12}=(p_{1}+&\,p_{2})^{2}\,,\quad
s_{23}=(p_{2}+p_{3})^{2}\,,\quad s_{34}=(p_{3}+p_{4})^{2}\,,\\\
&s_{45}=(p_{4}+p_{5})^{2}\,,\quad s_{15}=(p_{1}+p_{5})^{2}\,,\end{split}$ (2)

as well as the parity-odd contraction of four momenta,

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

Strictly speaking, scattering amplitudes for processes such as that of eq. (1)
cannot be expressed in terms of just the set
$\\{s_{12},s_{23},s_{34},s_{45},s_{51},\text{tr}_{5}\\}$, as this requires
removing an arbitrary little-group-dependent factor. In fact, amplitudes
depend not only on the four-momenta, but also on the helicities of the
external states.

To better represent this dependence on the helicities, we can adopt a
different set of variables, namely the two-component spinors,
$\lambda_{i}^{\alpha}$ and $\tilde{\lambda}_{i}^{\dot{\alpha}}$, with
$i\in\\{1,\dots,5\\}$. Starting from the $2\times 2$ spinors
$p_{i}^{\dot{\alpha}\alpha}$, which are given in terms of the respective four-
momenta as
$p_{i}^{\dot{\alpha}\alpha}=p_{i,\mu}\sigma^{\mu\dot{\alpha}\alpha}$, with
$\sigma^{\mu\dot{\alpha}\alpha}=(\mathbb{1},\vec{\sigma})$ and $\vec{\sigma}$
the Pauli matrices, the two-component spinors $\lambda_{i}^{\alpha}$ and
$\tilde{\lambda}_{i}^{\dot{\alpha}}$ are defined by noting that for massless
particles

$\det(\\{p_{i}^{\dot{\alpha}\alpha}\\})=0\quad\Longrightarrow\quad
p_{i}^{\dot{\alpha}\alpha}=\tilde{\lambda}_{i}^{\dot{\alpha}}\lambda_{i}^{\alpha}\,.$
(4)

Lowering of spinor indices is performed as
$\lambda_{i,\alpha}=\epsilon_{\alpha\beta}\lambda_{i}^{\beta}$ and
$\tilde{\lambda}_{i,\dot{\alpha}}=\epsilon_{\dot{\alpha}\dot{\beta}}\tilde{\lambda}_{i}^{\dot{\beta}}$,
where we make use of the Levi-Civita symbol
$\epsilon^{\alpha\beta}=\epsilon^{\dot{\alpha}\dot{\beta}}=-\epsilon_{\alpha\beta}=-\epsilon_{\dot{\alpha}\dot{\beta}}=i\sigma_{2}$.
Invariant contractions of spinors give so-called spinor brackets, which we
define as

$\langle
ij\rangle=\lambda^{\alpha}_{i}\lambda_{j,\alpha}\quad\text{and}\quad[ij]=\tilde{\lambda}_{i,\dot{\alpha}}\tilde{\lambda}_{j}^{\dot{\alpha}}\,.$
(5)

These are related to the Mandelstam invariants in eq. (2) through
$s_{ij}=\langle ij\rangle[ji]$. We also use longer spinor contractions, in
particular

$\langle i|j+k|i]=\langle ij\rangle[ji]+\langle ik\rangle[ki]\,.$ (6)

Finally, we can express ${\rm tr}_{5}$ as a polynomial in spinor brackets as

${\rm tr}_{5}=[12]\langle 23\rangle[34]\langle 41\rangle-\langle
12\rangle[23]\langle 34\rangle[41]\,.$ (7)

##### Helicity amplitudes

We closely follow the notation and conventions of ref. [48], and denote the
(renormalised) amplitudes for this process by

$\displaystyle\mathcal{M}(1_{q}^{h_{1}},2_{\bar{q}}^{h_{2}},3_{\gamma}^{h_{3}},4_{\gamma}^{h_{4}},5_{\gamma}^{h_{5}})\coloneqq
e_{q}^{3}\delta_{i_{1}i_{2}}\mathcal{A}(1_{q}^{h_{1}},2_{\bar{q}}^{h_{2}},3_{\gamma}^{h_{3}},4_{\gamma}^{h_{4}},5_{\gamma}^{h_{5}})\,,$
(8)

where $i_{1}$ and $i_{2}$ are the color indices of the external quarks and
$e_{q}$ is their electric charge. We call
$\mathcal{A}(1_{q}^{h_{1}},2_{\bar{q}}^{h_{2}},3_{\gamma}^{h_{3}},4_{\gamma}^{h_{4}},5_{\gamma}^{h_{5}})$
the _helicity amplitudes_ for the process in eq. (1), and will often suppress
their arguments for simplicity.

Helicity amplitudes satisfy relations under permutations of the photon momenta
or under charge and parity conjugation. For the process in eq. 1 there are two
independent helicity configurations, which we choose to be

$\displaystyle\mathcal{A}_{+++}(1,2,3,4,5)$
$\displaystyle\coloneqq\mathcal{A}(1_{q}^{+},2_{\bar{q}}^{-},3_{\gamma}^{+},4_{\gamma}^{+},5_{\gamma}^{+})\,,$
(9) $\displaystyle\qquad\mathcal{A}_{-++}(1,2,3,4,5)$
$\displaystyle\coloneqq\mathcal{A}(1_{q}^{+},2_{\bar{q}}^{-},3_{\gamma}^{-},4_{\gamma}^{+},5_{\gamma}^{+})\,,$

where we indexed the independent amplitudes by the photon helicities. We work
in the ’t Hooft–Veltman scheme of dimensional regularisation, setting the
space-time dimensions to $D=4-2\epsilon$, and use the definition of
dimensionally regularised helicity amplitudes with external quarks given in
ref. [46]. We perform the UV renormalisation in the $\overline{\text{MS}}$
scheme, where the amplitudes admit an expansion in terms of the renormalised
QCD coupling constant $\alpha_{s}$ of the form

$\mathcal{A}=\mathcal{A}^{(0)}+\frac{\alpha_{s}}{2\pi}\mathcal{A}^{(1)}+\quantity(\frac{\alpha_{s}}{2\pi})^{2}\mathcal{A}^{(2)}+\ldots{}\,.$
(10)

The coupling $\alpha_{s}$ is related to the bare coupling $\alpha_{s}^{0}$
through

$\alpha_{s}^{0}\mu_{0}^{2\epsilon}S_{\epsilon}=\alpha_{s}\mu^{2\epsilon}\left(1-\frac{\beta_{0}}{\epsilon}\frac{\alpha_{s}}{2\pi}+\mathcal{O}\left(\alpha_{s}^{2}\right)\right)\,,\;\quad
S_{\epsilon}=(4\pi)^{\epsilon}e^{-\epsilon\gamma_{E}}\,,$ (11)

where $\gamma_{E}$ is the Euler–Mascheroni constant, $\mu_{0}$ and $\mu$ are
the dimensional regularization and renormalization scales (which we assume to
be equal), and $\beta_{0}$ is the first coefficient of the QCD
$\beta$-function,

$\beta_{0}=\frac{11}{6}C_{A}-\frac{2}{3}T_{F}{N_{f}}\,.$ (12)

Here, $C_{A}={N_{c}}$ is the quadratic Casimir of the adjoint representation
of the $SU({N_{c}})$ group, ${N_{f}}$ is the number of massless quarks, and
$T_{F}=1/2$ is the normalization of the generators of the fundamental
representation. Below we will also need the quadratic Casimir of the
fundamental representation, $C_{F}=\frac{{N_{c}}^{2}-1}{2{N_{c}}}$. The
coefficients of the perturbative expansion of the renormalised amplitudes
$\mathcal{A}^{(\ell)}$ are related to their bare counterparts
$\mathcal{A}^{(\ell)}_{\mathcal{B}}$, which are coefficients in a perturbative
expansion in powers of $\alpha_{s}^{0}$, by

$\mathcal{A}^{(0)}_{\mathcal{B}}=\mathcal{A}^{(0)}\,,\quad\mathcal{A}^{(1)}_{\mathcal{B}}=S_{\epsilon}\mathcal{A}^{(1)}\,,\quad\mathcal{A}^{(2)}_{\mathcal{B}}=S_{\epsilon}^{2}\left(\mathcal{A}^{(2)}+\frac{\beta_{0}}{\epsilon}\mathcal{A}^{(1)}\right)\,.$
(13)

The coefficients $\mathcal{A}^{(\ell)}$ can be decomposed into individually
gauge-invariant contributions that scale differently with the number of light
quarks ${N_{f}}$, the number of colors ${N_{c}}$ and the electric charges of
the light fermions. Up to second order, we have [51, 52]

$\displaystyle\begin{split}\mathcal{A}^{(1)}&=C_{F}A^{(1)}\,,\\\
\mathcal{A}^{(2)}&=C_{F}^{2}B^{(2,0)}+C_{F}C_{A}B^{(2,1)}+C_{F}T_{F}{N_{f}}A^{(2,{N_{f}})}+C_{F}T_{F}\left(\sum_{f=1}^{{N_{f}}}Q_{f}^{2}\right)\,A^{(2,{\tilde{N}_{f}})}\,,\end{split}$
(14)

where $Q_{f}$ denotes the ratio of the charges of the light fermions running
in a closed fermion loop to the charge of the initial state quark/anti-quark
pair. We can rearrange eq. 14 as

$\displaystyle\begin{split}\mathcal{A}^{(2)}&=\frac{N_{c}^{2}}{4}\left(A^{(2,0)}-\frac{1}{{N_{c}}^{2}}(A^{(2,0)}+A^{(2,1)})+\frac{1}{{N_{c}}^{4}}A^{(2,1)}\right)\\\
&\qquad+C_{F}T_{F}{N_{f}}A^{(2,{N_{f}})}+C_{F}T_{F}\left(\sum_{f=1}^{{N_{f}}}Q_{f}^{2}\right)\,A^{(2,{\tilde{N}_{f}})}\,,\end{split}$
(15)

where we have

$\displaystyle A^{(2,0)}\coloneqq B^{(2,0)}+2B^{(2,1)},\qquad
A^{(2,1)}\coloneqq B^{(2,0)}.$ (16)

The contributions $A^{(2,0)}$ and $A^{(2,{N_{f}})}$ involve only planar
diagrams and were previously computed in refs. [12, 48, 53]. In this work we
obtain the missing contributions $A^{(2,1)}$ and $A^{(2,{\tilde{N}_{f}})}$.
Representative diagrams for each contribution are shown in figure 1.

((a)) $A^{(2,0)}$

((b)) $A^{(2,{N_{f}})}$

((c)) $A^{(2,{\tilde{N}_{f}})}$

((d)) $A^{(2,1)}$

Figure 1: Representative diagrams for the different contributions to eq. (14).

The renormalised amplitudes still contain infrared (IR) poles. They can be
predicted and subtracted in a scheme-dependent way. We use the scheme of refs.
[54, 51, 52], which we denote as the Catani scheme. The finite remainders
$\mathcal{R}$ are obtained from the UV-renormalised amplitudes $\mathcal{A}$
through

$\mathcal{R}\coloneqq\mathbf{I}\mathcal{A}\,.$ (17)

Upon expansion of eq. 17 in $\alpha_{s}$ and $\epsilon$, we obtain

$\displaystyle\begin{split}\mathcal{R}^{(0)}&=\mathcal{A}^{(0)}\,,\\\
\mathcal{R}^{(1)}&=\mathcal{A}^{(1)}+\mathbf{I}^{(1)}\mathcal{A}^{(0)}+\mathcal{O}(\epsilon)\,,\\\
\mathcal{R}^{(2)}&=\mathcal{A}^{(2)}+\mathbf{I}^{(1)}\mathcal{A}^{(1)}+\mathbf{I}^{(2)}{\mathcal{A}}^{(0)}+\mathcal{O}(\epsilon)\,,\end{split}$
(18)

where the functions $\mathbf{I}^{(1)}$, $\mathbf{I}^{(2)}$ are the
coefficients of the expansion of $\mathbf{I}$ in $\alpha_{s}/(2\pi)$, and are
given in section A.1.111We note that our definition of $\mathbf{I}^{(1)}$,
$\mathbf{I}^{(2)}$ differs by a sign from the one of ref. [48]. Up to two-
loops we can expand the finite remainders into the gauge-invariant
contributions $R^{(1)}$, $R^{(2,0)}$, $R^{(2,1)}$, $R^{(2,{N_{f}})}$ and
$R^{(2,{\tilde{N}_{f}})}$ as in eqs. 14 and 15.

For phenomenological applications, we are mostly interested in squared finite
remainders, summed over helicity and color, which we will refer to as the
_hard function_. We define it as

$\mathcal{H}=\frac{1}{\mathcal{B}}~{}\sum_{h}\absolutevalue{\mathcal{R}_{h}}^{2}\,,\qquad\mathcal{B}\coloneqq\sum_{h}\absolutevalue{\mathcal{A}^{(0)}_{h}}^{2}\,,$
(19)

where the sum is over all different choices of helicities of the involved
particles. Expanding in $\alpha_{s}$ and in color factors, we get

$\displaystyle\begin{split}\mathcal{H}^{(1)}&=C_{F}H^{(1)},\\\
\mathcal{H}^{(2)}&=\frac{N_{c}^{2}}{4}\left(H^{(2,0)}-\frac{1}{{N_{c}}^{2}}(H^{(2,0)}+H^{(2,1)})+\frac{1}{{N_{c}}^{4}}H^{(2,1)}\right)\\\
&\qquad+C_{F}T_{F}{N_{f}}H^{(2,{N_{f}})}+C_{F}T_{F}\left(\sum_{f=1}^{{N_{f}}}Q_{f}^{2}\right)\,H^{(2,{\tilde{N}_{f}})},\end{split}$
(20)

and by definition $\mathcal{H}^{(0)}=1$. We note that the hard functions are
scheme dependent, with expressions in different schemes being related by a
finite shift. We give examples of such relations in section A.2.

##  3 Calculation

###  3.1 Overview

Our calculation is done within the framework of two-loop numerical unitarity
[27, 28, 29]. We rely on the implementation of the method in the program
Caravel [31].

This approach starts from the observation that the integrands of the gauge-
invariant contributions $A^{(2,k)}$ of eq. 15 all admit a decomposition of the
form

$A(\ell)=\sum_{\Gamma}\sum_{i\in M_{\Gamma}\cup
S_{\Gamma}}c_{\Gamma,i}\frac{m_{\Gamma,i}(\ell)}{\prod_{j}\varrho_{\Gamma,j}(\ell)}\,,$
(21)

where the outer sum is over all distinct sets $\Gamma$ of inverse propagators
$\varrho_{\Gamma,j}$ contributing to the amplitude (which we call
_topologies_). For each topology $\Gamma$, the numerators $m_{\Gamma,i}(\ell)$
are polynomials in the loop momenta $\ell$ (we use a single $\ell$ to
collectively denote all loop momenta), and they are constructed such that each
$m_{\Gamma,i}(\ell)$ either corresponds to a master integral ($i\in
M_{\Gamma}$) or can be expressed as a total derivative and therefore
integrates to zero, i.e. it is a surface term ($i\in S_{\Gamma}$). The
coefficients $c_{\Gamma,i}$ in eq. 21 are unknown rational functions of the
external particles’ momenta and the dimensional regulator $\epsilon$. In the
generalized unitarity method, the coefficients of each topology $\Gamma$ are
constrained by evaluating (generalized) _cuts_ , corresponding to residues of
eq. 21 at values of $\ell=\hat{\ell}$ such that
$\varrho_{\Gamma,j}(\hat{\ell})=0~{}\forall j$. The cuts of the LHS of eq. 21
are evaluated as products of $D$-dimensional tree-amplitudes (which we call
_cut diagrams_) and matched to the integrand parametrisation of eq. 21. In
this way, _cut equations_ are generated for various loop-momentum
configurations $\hat{\ell}_{k}$ satisfying the cut conditions
$\varrho_{\Gamma,j}(\hat{\ell}_{k})=0$. Given a large enough sample of
$\hat{\ell}_{k}$, the coefficients $c_{\Gamma,i}$ in eq. 21 can be determined
as (numerical) solutions to the cut equations. Our computation is based on the
public implementation of the two-loop numerical unitarity method in Caravel
[31], which we extended to account for the non-planar topologies appearing in
eq. 21. These were absent in previous applications of the two-loop numerical
unitarity approach. In particular, we constructed the complete set of surface
terms required for the non-planar contributions in $A^{(2,1)}$ and
$A^{(2,{\tilde{N}_{f}})}$. We discuss this in more detail in section 3.2.

We note that to be compatible with generalized cuts in $D$ dimensions, we
require that the numerators $m_{\Gamma,i}(\ell)$ in eq. 21 are polynomials in
the loop momenta (i.e. additional denominator powers are not allowed). In
addition, for each $\Gamma$ the numerators $m_{\Gamma,i}(\ell)$ must be
linearly-independent on $D$-dimensional cuts that set all $\varrho_{\Gamma,j}$
to zero. While for the surface terms these condition hold by construction, the
definitions of master integrals are frequently chosen such that at least one
of the conditions is violated. Indeed, the basis of master integrals of refs.
[55, 56] violates the second condition. We find that the basis of ref. [57] is
therefore more convenient for our approach. It can easily be written in terms
of the pentagon functions of ref. [56] using modern integration-by-parts
codes, e.g. [58, 59], which among other benefits allows for an efficient
numerical evaluation of the master integrals. The decomposition of eq. 21 then
leads to a decomposition of the finite remainders in eq. 18 in terms of
pentagon functions,

$\mathcal{R}^{(\ell)}=\sum_{i}r_{i}h_{i}\,,$ (22)

where the $r_{i}$ are rational functions of external kinematics, and the
$h_{i}$ are monomials of the pentagon functions of ref. [56]. In summary, two-
loop numerical unitarity gives us a way to numerically compute the $r_{i}$. We
can then use these numerical evaluations to reconstruct their analytic form.

Before discussing some details of the construction of surface terms in section
3.2, and of the analytic reconstruction of the finite remainders in section
3.3, we close this brief overview of the approach with some technical
comments. Cut diagrams are generated with qgraf [60], and arranged into a
hierarchy of cuts with a private code. To match the cuts evaluated through
color-ordered tree amplitudes to the amplitude definitions in section 2 we
employ the unitarity-based colour decomposition of refs. [61, 62]. We
determine the $\epsilon$-dependence of cut diagrams originating from the state
sums in the loops through the dimensional reduction method developed in refs.
[63, 64, 65]. This allows us to perform the entire calculation with six-
dimensional states only.

###  3.2 Surface Terms

The two-loop numerical unitarity framework [27, 28, 29] builds on the
parametrisation of the integrand as in eq. 21. A crucial step in this
procedure is the determination of a basis of surface terms that integrate to
zero. In this section we present the method we used to construct surface terms
for the non-planar topologies.

The construction of surface terms starts from the observation that total
derivatives of Feynman integrands integrate to zero in dimensional
regularisation [39, 40, 66]. That is, we can construct surface terms from

$\displaystyle\int\text{d}^{D}\ell_{1}\cdots\text{d}^{D}\ell_{L}\frac{\partial}{\partial\ell_{a}^{\mu}}\frac{v_{a}^{\mu}(\ell)}{\varrho_{1}\cdots\varrho_{N}}=0\,.$
(23)

For an arbitrary vector $v_{a}^{\mu}$, the equation above will generate
surface terms involving integrals that are not relevant for the amplitudes we
are computing, either because they have numerators of too high degree or
because they have propagators raised to too high powers. In order to have an
efficient construction of the surface terms required for the decomposition in
eq. 21, it is beneficial to construct a minimal set of surface terms, which
contain specific propagator powers and numerators whose polynomial degree is
limited by the interactions of the underlying process.

While the numerator power in the surface terms is easy to control (because
derivatives do not increase it), care must be taken with the power of the
propagators. Propagator powers in surface terms can be controlled by requiring
that the vector $v^{\mu}_{a}$ in eq. 23 satisfies [41]

$\displaystyle\sum_{a,\mu}v^{\mu}_{a}(\ell)\frac{\partial\varrho_{i}}{\partial\ell^{\mu}_{a}}=f_{i}(\ell)\varrho_{i}\,,\quad\forall
i\,.$ (24)

where the unknowns $f_{i}$ and $v^{\mu}_{a}$ are polynomials in loop momenta.
We call such vectors unitarity-compatible integration-by-parts (IBP)
generating vectors, or simply IBP-generating vectors. The mathematical
structure of eq. 24 is well known and defines the vectors $v^{\mu}_{a}$ and
the $f_{i}$ to be elements of a syzygy module.

Solving eq. 24 allows us to construct a minimal set of surface terms. In
practice, however, obtaining analytic expressions for the surface terms can be
challenging, both due to the difficulty of solving eq. 24 and due to the size
of the final expressions (see refs. [67, 68, 69, 70, 71, 72] for related
work). Since our goal is to reconstruct analytic expressions for the two-loop
remainders from their numerical evaluations on a sufficient number of phase-
space points, we do not actually require analytic solutions to eq. 24. Indeed,
it is sufficient to obtain IBP-vectors that are analytic in the loop-momentum
variables but numerical in external momenta. Naturally, this requires solving
eq. 24 at each phase-space point, and one must therefore have an approach that
is efficient. We now present our solution to this problem: we first discuss
its formulation in embedding space [73, 43], and then discuss how unitarity-
compatible IBP vectors and surface terms are constructed at a given phase-
space point, and how we efficiently extend this to subsequent phase-space
points.

####  3.2.1 Embedding-Space Formalism

We start by reviewing the formulation of Feynman integrals in embedding space.
We will observe in the next section that this formulation simplifies the
solution of the syzygy equations in eq. 24. Formally, our goal is to discuss
how momentum space can be mapped into a subset of a projective space, commonly
called embedding space. To this end, we map each point $z^{\mu}$ in momentum
space into a line $Z$ of projectively equivalent points. In the context of a
Feynman integral with $N$ external legs, this can be done as follows [73, 43].
We first define

$q_{i}^{\mu}=\sum_{j=1}^{i-1}k_{j}^{\mu}\,,\qquad 1\leq i\leq N\,,$ (25)

where the $k_{j}^{\mu}$ denote the external momenta and the empty sum gives
the zero vector. We then arrange the loop and external momenta into the
$(D+2)$-dimensional embedding-space vectors

$Y_{a}=c_{Y_{a}}\left(\ell_{a}^{\mu},(\ell_{a})^{2},1\right)\,,\quad
X_{i}=c_{X_{i}}\left(q_{i}^{\mu},(q_{i})^{2},1\right)\,,$ (26)

where $c_{Y_{a}}$ and $c_{X_{i}}$ parametrise the points in embedding space
that are projectively equivalent. As a projective space, embedding space is
equipped with a special point, commonly referred to as the infinity point and
denoted here by $X_{0}$,

$\displaystyle X_{0}=\left(0^{\mu},1,0\right)\,.$ (27)

We also define an inner product as

$\displaystyle(AB)=c_{A}c_{B}\left(-2a^{\mu}\cdot
b_{\mu}+a^{2}+b^{2}\right)\,.$ (28)

An $L$-loop Feynman integral is then mapped into embedding space as

$\int\left[\prod_{a=1}^{L}\text{d}^{D}\ell_{a}\right]f(\ell;q)=\int\left[\prod_{a=1}^{L}\frac{d^{D+2}Y_{a}}{{\rm
vol}[{\rm
GL}(1)]}\frac{\delta[(Y_{a}Y_{a})]}{(X_{0}Y_{a})^{D}}\right]F(Y;X)\,,$ (29)

where $F(Y;X)$ is implicitly defined as the image of $f(\ell;q)$ under the
embedding procedure and ${\rm vol}[{\rm GL}(1)]$ accounts for the projective
equivalence parametrised by the $c_{Y_{a}}$.

In this construction, we seem to have increased the number of degrees of
freedom from $D$ to $D+2$. One apparent extra degree of freedom is removed by
noting that any $z^{\mu}$ in momentum space is mapped to an embedding space
point $Z$ that satisfies $(ZZ)=0$, i.e., onto the light-cone in embedding
space. Indeed, we can easily check from eq. 26 that $(Y_{a}Y_{a})=0$ and
$(X_{i}X_{i})=0$. The integration in embedding space in eq. 29 is restricted
to the light-cone by the delta function $\delta[(Y_{a}Y_{a})]$, thereby
corresponding to the integration over all loop-momentum space. Another
apparent extra degree of freedom is parametrised by the non-vanishing
parameters $c_{Y_{a}}$ and $c_{X_{i}}$. To understand how it is removed in
practice, we first note that, for an arbitrary embedding-space vector $Z$,
$c_{Z}=(X_{0}Z)$. It therefore follows that

$(\ell_{1}-\ell_{2})^{2}=\frac{(Y_{1}Y_{2})}{(X_{0}Y_{1})(X_{0}Y_{2})}\,,\quad(\ell_{a}-q_{i})^{2}=\frac{(Y_{a}X_{i})}{(X_{0}Y_{a})(X_{0}X_{i})}\,,\quad(q_{i}-q_{j})^{2}=\frac{(X_{i}X_{j})}{(X_{0}X_{i})(X_{0}X_{j})}\,.$
(30)

The function $f(\ell;q)$ in eq. 29 is a rational function of these momentum-
space inner products, which implies that $F(Y;X)$ is homogeneous of degree 0
in both $c_{Y_{a}}$ and $c_{X_{i}}$, that is

$F(Y;X)=F(Y_{1},\ldots,\lambda Y_{a},\ldots,Y_{L};X)=F(Y;X_{1},\ldots,\lambda
X_{i},\ldots,X_{N})\,,\qquad\lambda\neq 0\,.$ (31)

Since the $X_{i}$ only appear in eq. 29 within $F(Y;X)$, we can set
$c_{X_{i}}=(X_{0}X_{i})=1$. In principle we could make the same choice for
$c_{Y_{a}}$, but this would obscure some of the properties of the Feynman
integral, such as the fact that in dimensional regularization the integrand
has a branch point at $(X_{0}Y_{a})=0$ [73, 43, 45]. Equivalently, we can keep
the explicit dependence on the $(X_{0}Y_{a})$, and the division by ${\rm
vol}[{\rm GL}(1)]$ ensures that the integral is well defined despite the
invariance of $F(Y;X)$ with respect to rescalings of each $Y_{a}$.

The (inverse) propagators in a Feynman integral can all be written in terms of
the expressions given in eq. 30, where we are free to set $(X_{0}X_{i})=1$ as
discussed above. For two-loop five-point non-planar integrals it is also
convenient to define another combination of the form

$\displaystyle(\ell_{1}-\ell_{2}+k_{i})^{2}=\frac{r_{AB}Y_{1}^{A}Y_{2}^{B}}{(X_{0}Y_{1})(X_{0}Y_{2})}\,,$
(32)

where $r_{AB}$ is a matrix that depends only on external kinematics. That is,
the components of $r_{AB}$ can be written in terms of the $(X_{i}X_{j})$ and
are independent of $Y_{1}$ and $Y_{2}$. The denominator
$(\ell_{1}-\ell_{2})^{2}$ is a particular case of eq. 32. We denote the set of
propagator denominators in a two-loop five-point non-planar integral as

${\cal
P}=\big{\\{}\left(X_{1}Y_{1}\right)\,,\left(X_{2}Y_{1}\right)\,,\ldots\,,(Y_{1}Y_{2})\,,r_{AB}Y^{A}_{1}Y^{B}_{2}\big{\\}}\,,$
(33)

where we note that the last term is not required for planar integrals.

Figure 2: Non-planar diagram corresponding to the propagators in eq. 34.

For concreteness we explicitly give an example of a simple non-planar diagram
in embedding space. The propagators of the diagram in fig. 2 take the form

$\displaystyle\begin{split}\varrho_{1}&=\ell_{1}^{2}=\frac{(X_{1}Y_{1})}{(X_{0}Y_{1})}\,,\quad\varrho_{2}=(\ell_{1}-k_{1})^{2}=\frac{(X_{2}Y_{1})}{(X_{0}Y_{1})}\,,\quad\varrho_{3}=(\ell_{2}+k_{3})^{2}=\frac{(X_{3}Y_{2})}{(X_{0}Y_{2})}\,,\\\
\varrho_{4}&=(\ell_{2}+k_{2}+k_{3})^{2}=\frac{(X_{2}Y_{2})}{(X_{0}Y_{2})}\,,\quad\varrho_{5}=(\ell_{1}-\ell_{2})^{2}=\frac{(Y_{1}Y_{2})}{(X_{0}Y_{1})(X_{0}Y_{2})}\,,\\\
\varrho_{6}&=(\ell_{1}-\ell_{2}-k_{3})^{2}=\frac{1}{(X_{0}Y_{1})(X_{0}Y_{2})}\left[(Y_{1}Y_{2})+(X_{0}Y_{1})(X_{3}Y_{2})-\right.\\\
&\left.\qquad\quad(X_{0}Y_{1})(X_{1}Y_{2})+(X_{0}Y_{2})(X_{1}Y_{1})-(X_{0}Y_{2})(X_{3}Y_{1})+(X_{0}Y_{1})(X_{0}Y_{2})(X_{1}X_{3})\right]\,,\end{split}$
(34)

where $\varrho_{6}$ has the form of eq. 32 and we have written the components
of the respective matrix $r_{AB}$ explicitly in terms of the $(X_{i}X_{j})$.

####  3.2.2 IBP-Generating Vectors In Embedding Space

Let us now discuss how IBP relations arise in embedding space, following
previous work on the subject [43, 44]. To begin, we reformulate eq. 23 in
embedding space as222 In practice, it can be more convenient to replace the
$\delta\left[\left(Y_{b}Y_{b}\right)\right]$ by $1/\left(Y_{b}Y_{b}\right)$
and modify the integration contour to encircle the pole at
$\left(Y_{b}Y_{b}\right)$. This perspective can be taken to more easily derive
some of the properties of IBP vectors in embedding space.

$\displaystyle\int\prod_{a}\left[\frac{d^{D+2}Y_{a}}{{\rm vol}[{\rm
GL}(1)]}\right]\times\frac{\partial}{\partial
Y_{c}^{C}}\left[\frac{V_{c}^{C}}{\varrho_{1}\cdots\varrho_{N}}\,\prod_{b}\frac{\delta\left[\left(Y_{b}Y_{b}\right)\right]}{\left(X_{0}Y_{b}\right)^{D}}\right]=0\,,$
(35)

where the $\rho_{i}$ are now understood as their embedding-space expressions.
The IBP-generating vectors are denoted $V_{c}^{C}$, where $C$ runs over the
$D+2$ dimensions of an embedding-space vector. A vector $V_{c}$ that generates
linear relations between Feynman integrals in embedding space must satisfy a
number of non-trivial properties. First, the components of $V_{c}^{C}$ must be
homogeneous functions in $Y_{c}$ of degree one in order to compensate for the
scaling of the partial derivatives. More precisely, they are required to be
rational functions with denominators given by products of the factors
$(X_{0}Y_{a})$, so that no further propagator poles are introduced. Second, it
can be verified that

$\displaystyle\frac{\partial}{\partial
Y_{c}^{C}}\left[\frac{F_{c}Y_{c}^{C}}{\varrho_{1}\cdots\varrho_{N}}\,\prod_{b}\frac{\delta\left[\left(Y_{b}Y_{b}\right)\right]}{\left(X_{0}Y_{b}\right)^{D}}\right]=0\,,$
(36)

for any homogeneous $F_{c}$ that is degree 0 in $Y_{c}$. As such, any IBP-
generating vector $V_{c}$ generates the same surface term as another vector
obtained from it by translations in the $Y_{c}$ direction. Finally, any valid
IBP vector must only generate relations between Feynman integrals. We
therefore impose that any terms containing derivatives of the delta function
that arise from expanding the derivative in eq. 35 should cancel. This can be
achieved by requiring that

$\displaystyle V_{c}^{C}\frac{\partial}{\partial Y_{c}^{C}}(Y_{b}Y_{b})=F^{\rm
quadric}_{b}\,(Y_{b}Y_{b})\,,\qquad\forall b\,,$ (37)

i.e. by requiring that the $V_{c}$ generate translations along the light cone.
It turns out that many of the solutions to eq. 37 can be discarded as they
lead to vanishing surface terms. Indeed, it is easy to see that a subset of
the solutions is given by $V_{c}$ such that $V_{c}=\frac{1}{2}F_{c}^{\rm
quadric}Y_{c}$. By comparison with eq. 36, we see that such solutions give a
vanishing surface term. We can use this observation to simplify IBP-vector
construction. Specifically, if we have an IBP vector $V_{c}$ for which the
associated $F_{c}^{\rm quadric}\neq 0$, then the same surface term is
generated by a vector $V^{\prime}_{c}=V_{c}-\frac{1}{2}F_{c}^{\rm
quadric}Y_{c}$. It follows that the new vector $V^{\prime}_{c}$ satisfies
equation (37) with $F^{\rm quadric}_{b}=0$. Therefore, the solutions of eq. 37
with $F^{\rm quadric}_{b}=0$ generate the full set of surface terms and it is
sufficient to require that the IBP vectors $V_{c}$ satisfy

$\displaystyle V_{c}^{C}\frac{\partial}{\partial
Y_{c}^{C}}(Y_{b}Y_{b})=0\,,\qquad\forall b\,.$ (38)

Within the two-loop generalised unitarity framework, we also require that IBP-
generating vectors are unitarity compatible, that is that they satisfy eq. 24.
In embedding space, the analogous condition is that the exponents of the
denominator factors $(X_{i}Y_{a})$ and
$\big{(}r_{AB}Y^{A}_{1}Y^{B}_{2}\big{)}$ are not increased [44]. For
propagators that depend on a single $Y_{a}$, the conditions read,

$\displaystyle V_{c}^{C}\frac{\partial}{\partial
Y_{c}^{C}}\left(X_{i}Y_{b}\right)=F_{i}\,\left(X_{i}Y_{b}\right)\,,$ (39)

and for the propagators that depend on both loop momenta

$\displaystyle V_{c}^{C}\frac{\partial}{\partial
Y_{c}^{C}}\left(r_{AB}Y_{1}^{A}Y_{2}^{B}\right)=F_{r}\left(r_{AB}Y_{1}^{A}Y_{2}^{B}\right)\,,$
(40)

which also covers the case of the propagator $(Y_{1}Y_{2})$.

In summary, we must construct IBP-generating vectors satisfying eqs. 38, 39
and 40. These constrain the vectors $V_{a}^{A}$ and the polynomials $F_{i}$
and $F_{r}$ to be elements of the syzygies of a module defined by

$\displaystyle\begin{split}\left(V_{a}Y_{a}\right)&=0\,,\quad\forall a\,,\\\
\left(V_{a}X_{i}\right)&=F_{i}\left(X_{i}Y_{a}\right)\,,\quad\text{for
}\left(X_{i}Y_{a}\right)\in{\cal P}\,,\\\
r_{AB}\left[V_{1}^{A}Y^{B}_{2}+V_{2}^{B}Y^{A}_{1}\right]&=F_{r}\left(r_{AB}Y_{1}^{A}Y_{2}^{B}\right)\,,\quad\text{for
}(Y_{1}Y_{2})\text{ and }\left(r_{AB}Y_{1}^{A}Y_{2}^{B}\right)\in{\cal
P}\,,\end{split}$ (41)

where ${\cal P}$ was defined in eq. 33. The unknowns in these equations are
the $V_{a}^{A}$, $F_{i}$ and $F_{r}$. From now on, for simplicity we set
$(X_{0}Y_{a})=1$ and require the solutions to be polynomials in the
$(X_{i}Y_{a})$ and $(Y_{1}Y_{2})$, keeping in mind that those factors can be
reinstated by requiring that quantities have the correct homogeneous degree.

####  3.2.3 Computing IBP-Generating Vectors

We now turn towards solving the syzygy equations (41). In short, our strategy
is to construct a degree-bounded generating set of solutions on a single
phase-space point and then use the shape of this set to compute a generating
set of vectors on further phase space points, via linear algebra. We start by
expressing the IBP vectors in terms of the momenta in the problem,

$\displaystyle V_{a}=\sum_{i=0}^{N}v^{i}_{a}X_{i}+v^{Y}_{a}Y_{b\neq a},$ (42)

where we took into consideration eq. 36 and the surrounding discussion to
remove the $Y_{a}$ direction. We can then insert this parametrisation into eq.
41 and re-express these equations as

$\displaystyle M_{IJ}\,w^{J}=0\,,\quad
w^{J}=\\{v^{i}_{1},v^{Y}_{1},v^{i}_{2},v^{Y}_{2},F_{i},F_{r}\\}\,,$ (43)

where the vector $w^{J}$ collects all unknowns in the problem (while it
contains a single $v_{1}^{Y}$ and $v_{2}^{Y}$, there are several $v^{i}_{1}$
and $v^{i}_{2}$, several $F_{i}$ and two $F_{r}$). The direct solution of eq.
43 is in general very challenging. Since we are only interested in a subset of
solutions that generates sufficiently many surface terms to construct the
decomposition in eq. 21, which is degree bounded by the theory, we instead
solve the problem for a fixed polynomial degree [67, 29, 72]. To make this
degree manifest, we consider the ansatz

$\displaystyle w^{J}=\sum_{|\vec{m}|\leq Q}w^{J}_{\vec{m}}y^{\vec{m}}\,,\qquad
y=\\{\left(X_{1}Y_{1}\right),\ldots,\left(X_{1}Y_{2}\right),\ldots,\left(Y_{1}Y_{2}\right)\\}\,,$
(44)

where the $w^{J}_{\vec{m}}$ are $Y$-independent rational functions of the
external kinematics, i.e. of the $(X_{i}X_{j})$, that multiply associated
monomials of the $y_{k}$ variables,

$\displaystyle y^{\vec{m}}:=\prod_{k}y_{k}^{m_{k}}\,,$ (45)

where the multi-index $\vec{m}$ specifies the exponents of the monomial. The
total degree of a monomial $y^{\vec{m}}$ is denoted by $|\vec{m}|$,

$\displaystyle\quad|\vec{m}|:=\sum_{k}m_{k}\,.$ (46)

Inserting eq. 44 into eq. 43 one obtains a linear system for the coefficients
$w^{J}_{\vec{m}}$ allowing us to solve the syzygy problem with a degree bound
$Q$ as a linear algebra problem.

We can however further simplify the problem. Indeed, we note that the matrix
$M_{IJ}$ is linear in the $y_{k}$, that is

$\displaystyle M_{IJ}=M_{IJ}^{[0]}+M_{IJ}^{[1]}\,,\qquad
M_{IJ}^{[1]}=\sum_{k}M_{IJk}^{[1]}y_{k}\,.$ (47)

It is interesting to note that, as the $X_{i}$ are linearly independent,
$M_{IJ}^{[1]}$ is independent of external kinematics, i.e., the
$M_{IJk}^{[1]}$ are matrices of rational numbers. We can also decompose the
ansatz of eq. 44 into leading and subleading $y$ contributions,

$\displaystyle w^{J}$ $\displaystyle=w^{J}_{\rm max}+w^{J}_{\rm rem}\,,\qquad
w^{J}_{\rm max}=\sum_{|\vec{m}|=Q}w_{J}^{\vec{m}}y^{\vec{m}}\,.$ (48)

It then follows that

$\displaystyle M_{IJ}^{[1]}\,w^{J}_{\rm max}$ $\displaystyle=0\,,$ (49)

which is a simpler syzygy problem than eq. 43 as it is homogeneous and
independent of kinematics. For $w^{J}$ to be a solution to eq. 43, it is
necessary that its maximal degree piece solves eq. 49.

To solve eq. 49 we use the methods implemented in Singular [74], which yield a
generating set of solutions, with each generator homogeneous of a given
degree. That is, we have the generating set,

$\displaystyle\left\\{w_{k}^{J,[q_{k}]}(\vec{y})\right\\}_{k=1,\,...\,k_{\rm
max}}\,,$ (50)

where $q_{k}$ is the degree of the homogeneous generator $w_{k}^{J,[q_{k}]}$.
A basis of solutions with degree $q$ to eq. 49 is obtained from linear
combinations of the generating set of eq. 50 with coefficients that are
homogeneous polynomials of degree $q-q_{k}$.

We now use the homogeneous syzygies of eq. 49 to construct solutions to the
full syzygy problem of eq. 43. Since we require solutions only up to a maximal
polynomial degree we again consider solutions up to a given polynomial degree
$Q$. To this end, we consider a parametrization of a syzygy with degree $q$
that consists of

  1. 1.

a maximal degree piece of degree $q$ constructed from the generators eq. 50.
This consists of combinations of the generators multiplied with polynomials of
suitable polynomial degree, and,

  2. 2.

a generic polynomial-valued vector with degree less than $q$.

More precisely, such a syzygy can be written as

$\displaystyle
v^{I,[q]}(\vec{y})=\sum_{k,\,\,|\vec{m}_{k}|+q_{k}=q}b_{k,\vec{m}_{k}}\,y^{\vec{m}_{k}}\,w_{k}^{I,[q_{k}]}(\vec{y})+\sum_{|\vec{m}|<q}b^{I}_{\vec{m}}y^{\vec{m}}\,,$
(51)

for $y$-independent coefficients $b_{k,\vec{m}_{k}}$ and $b^{I}_{\vec{m}}$.
The monomials $y^{\vec{m}}$ and $y^{\vec{m}_{k}}$ are bounded by the degree
constraints $|\vec{m}|<q$ and $|\vec{m}_{k}|<q$. Inserting this ansatz into
eq. 43 then yields linear equations for $b_{k,\vec{m}_{k}}$ and
$b^{I}_{\vec{m}}$. In practice, we observe that this approach of combining
maximal degree solutions with a parametrization of lower degree terms
significantly reduces the size of the linear system compared to starting from
the generic parametrisation of eq. 44. In this way, working degree by degree
we can construct a basis of the solutions up to degree $Q$. Once we have
constructed this basis of the solutions it is easy to find a subset of the
basis elements which generate the module up to degree $Q$. Working degree by
degree, one simply removes basis elements that are polynomial combinations of
lower degree generators using linear algebra.

At this point we have obtained a generating set of the solutions of eq. 43 up
to polynomial degree $Q$ at a given numerical phase-space point. We now use
this information to streamline the construction of IBP vectors on further
phase-space points. We start with the generating set we have obtained and
express their components as a linear combination of monomials

$\displaystyle v^{I}(\vec{y})=\sum_{\vec{m}}r^{I}_{\vec{m}}y^{\vec{m}},$ (52)

where we sum over all monomials such that the associated $r^{I}_{\vec{m}}$ are
non-zero. We then replace the numerical coefficients $r^{I}_{\vec{m}}$ by
parameters $b^{I}_{\vec{m}}(\vec{p})$ that explicitly depend on the external
kinematic point $\vec{p}$, in order to obtain skeleton vectors

$\displaystyle
v^{I}_{\text{skel}}(\vec{y})=\sum_{\vec{m}}b^{I}_{\vec{m}}(\vec{p}){\vec{y}}^{\vec{m}}.$
(53)

We finally insert these skeleton vectors into (43), and obtain a linear system
of equations that constrain the $b^{I}_{\vec{m}}(\vec{p})$. Importantly, these
linear systems are analytic in the external kinematics and so can be solved
phase-space point by phase-space point. These equations may be linearly
dependent and not uniquely determine the $b^{I}_{\vec{m}}(\vec{p})$. We
discard linearly dependent equations, and choose a set of
$b^{I}_{\vec{m}}(\vec{p})$ to set to zero so that the solution is unique. By
making use of these skeleton vectors instead of fully general vectors, we
obtain equations for a minimal set of monomial coefficients. In practice, the
largest linear systems that we encounter when solving for the
$b^{I}_{\vec{m}}$ of eq. 53 are approximately $2000\times 2000$. Finally, we
note that for simpler case we find it efficient to solve these linear systems
analytically.

####  3.2.4 Computing Surface Terms

Now that we have constructed a set of IBP generating vectors up to degree $Q$,
we use them to construct a collection of surface terms. Recalling eq. 21, our
goal is to construct a set of linearly-independent surface terms $S_{\Gamma}$
for a given topology $\Gamma$ that is sufficiently large for the power
counting of our theory. This set of surface terms is a subset of the
$m_{\Gamma,i}$ of eq. 21 and they are polynomials in the loop momentum
components. The maximum degree of these polynomials is known a priori, since
it is determined by the theory describing the scattering process, and as such
it is easy to construct a basis of $\text{span}(M_{\Gamma}\cup S_{\Gamma})$.
As a basis of master integrals is also known [55, 57, 56], we have sufficient
information to construct a linearly-independent set of surface terms.

To construct a single surface term, we start with an IBP-generating vector and
a monomial of the $y_{k}$ variables. We take their product, and use this as
the vector $V_{c}^{C}$ in eq. 23. Expanding out the derivative leads to a
surface term in embedding space that is easily mapped to one in momentum space
by identifying all $y_{k}$ variables with their momentum space counterparts.
In this way, taking all pairs of IBP-vectors and monomials of $y_{k}$
variables that satisfy our power counting bounds, we construct a set of
surface terms. Many of the surface terms obtained in this way are linearly
dependent and can be discarded. As our construction of the IBP-generating
vectors was only performed up to degree $Q$, the set of surface terms
generated in this way may not be a basis of $\text{span}(S_{\Gamma})$. In
practice, to solve this problem, we repeat the procedure described in section
3.2.3, increasing the value of $Q$ in order to construct more IBP generating
vectors until our set of surface terms is a basis of
$\text{span}(S_{\Gamma})$.

###  3.3 Analytic Reconstruction In Spinor Helicity Formalism

As noted above, the two-loop numerical unitarity approach allows us to
numerically evaluate the coefficients $r_{i}$ in the finite remainder of eq.
22. We now discuss how to obtain analytic expressions using these numerical
evaluations. As suggested in e.g. refs. [75, 36], we will consider the $r_{i}$
to be rational functions of the spinor-helicity variables $\lambda$ and
$\tilde{\lambda}$ already introduced in section 2.333More precisely, the
coefficients belong to the field of fractions of the ring of independent
Lorentz invariants; see ref. [36, Section 2.2] for more details. The
coefficients admit a least common denominator representation, which reads

$r_{i}(\lambda,\tilde{\lambda})=\frac{\mathcal{N}_{i}(\lambda,\tilde{\lambda})}{\prod_{j}\mathcal{D}_{j}^{q_{ij}}(\lambda,\tilde{\lambda})}\,.$
(54)

The exponents $q_{ij}$ are allowed to take negative values, thus denoting
numerator factors. The $\mathcal{N}_{i}$ and $\mathcal{D}_{j}$ are polynomials
in spinor brackets. Like the amplitude, the $r_{i}$ are dimensionful and
transform non-trivially under the little group. For a spinor function
$\mathcal{E}(\lambda,\tilde{\lambda})$ we define the mass dimension
$[\mathcal{E}]$ and $k^{\mathrm{th}}$ little-group weight
$\\{\mathcal{E}\\}_{k}$ through

$\displaystyle\mathcal{E}(z\lambda_{1},\ldots,z\lambda_{n},z\tilde{\lambda}_{1},\ldots,z\tilde{\lambda}_{n})$
$\displaystyle=z^{2[\mathcal{E}]}\,\mathcal{E}(\lambda_{1},\ldots,\lambda_{n},\tilde{\lambda}_{1},\ldots,\tilde{\lambda}_{n})\,,$
(55)
$\displaystyle\mathcal{E}(\ldots,z\lambda_{k},\ldots,\tilde{\lambda}_{k}/z,\ldots)$
$\displaystyle=z^{\\{\mathcal{E}\\}_{k}}\,\mathcal{E}(\ldots,\lambda_{k},\ldots,\tilde{\lambda}_{k},\ldots).$
(56)

The mass dimension of an $n$-point amplitude $\mathcal{A}_{n}$ is
$[\mathcal{A}_{n}]=4-n$. The $k^{\mathrm{th}}$ little-group weight of an
amplitude depends on the $k^{\mathrm{th}}$-particles’ helicity $h_{k}$,

$\displaystyle\\{\mathcal{A}\\}_{k}=-2h_{k}\,.$ (57)

Crucially, there is by now a large collection of evidence that the
$\mathcal{D}_{j}$ in eq. 54 are known a priori. Indeed, they can be
constructed from the symbol alphabet [30]. In our normalisation of the
amplitudes (see eq. 8), they also include factors with little-group weight,
such as the tree level amplitude if it is not zero, or the leading order of
the one-loop amplitude. Given that the pentagon functions $h_{i}$ in eq. 22
are dimensionless and little-group invariant, the mass dimensions and little-
group weights of the $r_{i}$ are the same as those of the amplitudes. Since
the ${\cal D}_{j}$ are known, the mass dimensions and little-group weights of
the numerators $\mathcal{N}_{i}$ are easily determined.

In a nutshell, the reconstruction procedure amounts to constructing and
fitting an ansatz for the $\mathcal{N}_{i}$ in terms of spinor variables. An
important feature of spinor-helicity-based reconstruction methods is that the
ansatz directly captures the little-group and mass-dimension properties of the
coefficients. This is in contrast to Mandelstam-based reconstruction methods,
where the $r_{i}$ are normalized with a (dimensionful) phase factor to be
little-group invariant, and split into parity even and odd parts. This
difference has an important practical consequence, as the common denominator
form of the $r_{i}$ often simplifies in spinor variables. For example, the
symbol alphabet may factorize in spinor variables. This can lead to
denominator factors that cancel further against the coefficient’s numerators,
lowering the mass dimensions of the latter, which thus require simpler ansätze
in the reconstruction procedure.

To begin constructing our ansatz, we specify the set of denominator factors
${\cal D}_{j}$ in eq. (54). As already alluded to, these denominator factors
can be constructed from the symbol alphabet. In Mandelstam variables they
correspond to the subset of so-called even letters, which describe the
collection of singular surfaces associated to the master integrals. When
written in terms of spinor variables, these surfaces may branch. That is,
alphabet letters that appear irreducible when written in terms of Mandelstam
variables may further factorize in spinor space. Beyond alphabet letters, the
${\cal D}_{j}$ also include rational functions appearing in tree-level or one-
loop amplitudes. We therefore take these functions and the irreducible factors
of the symbol alphabet in spinor space as the set of expected denominator
factors. This analysis was already performed in ref. [36] for the amplitudes
considered in this paper, and it was found that the expected set of
denominator factors contains $35$ elements. It can be expressed as

${\cal D}=\big{\\{}\langle ij\rangle,\,[ij]\,\,:\,\,1\leq i<j\leq
5\big{\\}}\cup\left(\bigcup_{\sigma\in Z_{5}}\sigma\circ\big{\\{}\langle
1|2+3|1],\langle 1|2+4|1],\langle 1|2+5|1]\big{\\}}\right),$ (58)

where the elements of $Z_{5}$ are the cyclic permutations acting on the
momentum indices associated to each spinor. Note that, despite its presence in
the alphabet we do not include $\mathrm{tr}_{5}$ in $\mathcal{D}$, as it is
expected to cancel in the finite remainder (see e.g. refs. [57, 76, 77]).

In the following, it will be useful to consider the zero set associated to
various denominators. We denote the algebraic variety corresponding to the
common zero set associated to a list of denominator factors $\\{{\cal
D}_{i_{1}}\,,\ldots\,,{\cal D}_{i_{n}}\\}$ by

$V({\cal D}_{i_{1}}\,,\ldots\,,{\cal D}_{i_{n}})\,.$ (59)

##### Denominator exponents

Our first task is to determine the denominator exponents $q_{ij}$ in eq. (54).
We will use the standard technology of univariate-slice reconstruction [30],
combined with the approach introduced in ref. [78], based on all-line BCFW-
shifts [79, 80, 81]. Its application to functions of Mandelstam variables was
discussed in ref. [47]. Here we review the procedure and discuss its
generalization to functions of spinor variables.

In order to determine the $q_{ij}$, we must choose a set of curves in phase
space that intersect each surface $V(\mathcal{D}_{j})$ at a generic point at
least once and are not contained in any of the $V(\mathcal{D}_{j})$. We will
achieve this with two BCFW shifts, one holomorphic and one anti-holomorphic,
following the all-line shift approach of ref. [47] which we now review. Let us
begin with a generic, numeric momentum-conserving configuration of spinors
$\\{\lambda_{1},\tilde{\lambda}_{1},...,\lambda_{5},\tilde{\lambda}_{5}\\}$.
We make an all-line holomorphic shift by adjusting every $\lambda$ spinor in a
way which is proportional to a common reference spinor $\eta$. That is we
shift

$\displaystyle\lambda_{i}\rightarrow\lambda_{i}+tc_{i}\eta\,,\quad\tilde{\lambda}_{i}\rightarrow\tilde{\lambda}_{i}\,.$
(60)

Here, we introduce unknowns $c_{i}$ that we use to ensure that the shifted
kinematics satisfy momentum conservation. Specifically, we choose the $c_{i}$
to satisfy

$\displaystyle\sum_{i=1}^{5}c_{i}\tilde{\lambda}_{i}=0\,.$ (61)

This linear equation does not have a unique solution. Nevertheless, any
solution with $c_{i}\neq 0$ will guarantee a non-trivial shift that satisfies
momentum conservation. In order to ensure that a slice is generic, we pick
both the initial set of spinor variables
$\\{\lambda_{1},\tilde{\lambda}_{1},...,\lambda_{5},\tilde{\lambda}_{5}\\}$
and an independent subset of the $c_{a}$ randomly over a finite field. The
corresponding anti-holomorphic shift can be constructed by taking the parity
conjugate.

Let us now consider how a holomorphic shift affects the set of denominator
factors $\mathcal{D}$ in eq. 58. Firstly, the holomorphic spinor products are
linear functions of $t$

$\displaystyle\langle ab\rangle\rightarrow\langle
ab\rangle+t\,\big{(}c_{a}\langle\eta b\rangle+c_{b}\langle
a\eta\rangle\big{)}\,,$ (62)

whereas the anti-holomorphic spinor products remain unchanged. It is easy to
see that the curve parametrised by $t$ intersects each of the ten distinct
codimension-one varieties $V\left(\langle ab\rangle\right)$ as well as each of
the fifteen varieties $V\left(\langle a|b+c|a]\right)$. Nevertheless, it does
not intersect any of the ten varieties $V\left([ab]\right)$. By comparing the
form of the $\mathcal{D}_{j}$ on this univariate slice with the denominator of
$r_{i}$ on the slice, one can compute the exponents $q_{ij}$ corresponding to
the $\mathcal{D}_{j}$ which are not purely anti-holomorphic (i.e. the $[ab]$).
In order to determine the exponents of the $[ab]$ denominator factors, we
repeat the procedure on the anti-holomorphic shift. We note that the anti-
holomorphic shift must yield the same $q_{ij}$ for $\langle a|b+c|a]$ as the
holomorphic one, providing a consistency check. Having determined the $q_{ij}$
exponents in eq. 54, we can now determine the mass dimension and spinor weight
of the numerators $\mathcal{N}_{i}$, allowing us to write an ansatz for them
in terms of spinor-helicity variables.

##### Improved Ansatz

It is well known that the $r_{i}$ are not all independent rational functions,
and it is sufficient to reconstruct a basis of this space of functions. We
begin by sorting the $r_{i}$ by the mass dimension of their numerator
$\mathcal{N}_{i}$, and then apply standard linear-algebra techniques on
finite-field-valued evaluations of the $r_{i}$ to determine a basis of the
function space of the $r_{i}$. We therefore express the pentagon-function
coefficients as

$r_{i}=\tilde{r}_{j}M_{ji}\,,$ (63)

where $M_{ji}$ is a rectangular matrix with $M_{ji}\in\mathbb{Q}$. This
observation allows us to reduce the number of functions to reconstruct.

A further important simplification comes from the observation that the
rational coefficients in a scattering amplitude are less naturally expressed
in common-denominator form, and should in fact be cast in some partial-
fractions representation of the form

$\tilde{r}_{i}=\sum_{k}\tilde{r}_{ik},\qquad\text{with}\qquad\tilde{r}_{ik}=\sum_{k}\frac{\mathcal{N}_{ik}}{\prod_{j}\mathcal{D}_{j}^{q_{ijk}}}\,,$
(64)

where the $\mathcal{N}_{ik}$ are polynomials in spinor brackets and the
$q_{ijk}$ are integer exponents that determine the exact partial-fractions
decomposition. For analytic reconstruction approaches, this means that in
general a common-denominator ansatz is unnecessarily large. However,
systematically constructing a compact partial-fractions ansatz for a generic
amplitude remains a challenging task (see refs. [35, 36, 82] for recent
progress in this area). A practical solution to this issue is to explore
several possible partial-fractions decompositions, determine the size of the
ansatz of the $\mathcal{N}_{ik}$ for each of them, and declare that a suitable
(if not optimal) decomposition has been found when the ansatz is small enough
for the $\mathcal{N}_{ik}$ to be reconstructed. For the amplitudes we are
concerned with in this paper, however, we were able to identify a very
convenient partial-fractions decomposition that we now discuss.

To construct compact ansätze for the rational functions, we build upon the
observation that for several one-loop five-point massless amplitudes all poles
of the form $\langle c|a+b|c]^{\alpha}$ can be separated into different
fractions (see e.g. [83]). This feature is related to the fact that the poles
$\langle c|a+b|c]$ are spurious, and can be used to greatly simplify analytic
expressions for amplitudes [84]. Given this observation, we work under the
assumption that a one-loop-like partial-fractions representation exists at two
loops. The validity of this assumption will be tested during the calculation.
To take an explicit example, consider a function whose common-denominator form
reads

$\frac{\mathcal{N}}{\langle 1|2+4|1]\langle 1|2+5|1]\langle 2|1+5|2]\langle
4|1+2|4]^{2}\langle 4|1+5|4]^{2}\langle 5|1+2|5]^{2}\langle 5|1+4|5]^{2}}\,,$

where $\mathcal{N}$ is an unknown polynomial in spinor brackets to be
determined. The partial-fractions decomposition would then take the form

$\displaystyle\begin{split}&\frac{\mathcal{N}_{1}}{\langle
1|2+4|1]}+\frac{\mathcal{N}_{2}}{\langle
1|2+5|1]}+\frac{\mathcal{N}_{3}}{\langle
2|1+5|2]}+\frac{\mathcal{N}_{4}}{\langle 4|1+2|4]^{2}}\\\
&\qquad+\frac{\mathcal{N}_{5}}{\langle
4|1+5|4]^{2}}+\frac{\mathcal{N}_{6}}{\langle
5|1+2|5]^{2}}+\frac{\mathcal{N}_{7}}{\langle 5|1+4|5]^{2}}\,,\end{split}$ (65)

where the $\mathcal{N}_{i}$ are unknown polynomials in spinor brackets. The
ansätze for each of them is substantially simpler than that for the numerator
$\mathcal{N}$ in the common-denominator form, making the reconstruction
substantially more efficient.

Another layer of simplification of the reconstruction procedure comes from the
observation that, when considering a partial-fractions representation, many of
the basis functions are partially contained within the vector space spanned by
the others. More precisely, given a basis element $\tilde{r}_{n}$ expressed in
the form of eq. (64), many of the $\tilde{r}_{nk}$ belong to
$\operatorname{span}(\tilde{r}_{i\neq n})$. If we consider the reconstruction
of a basis function $\tilde{r}_{n}$, this motivates us to include other basis
functions in its ansatz, while dropping many of the terms in the partial-
fractions decomposition. That is, if we include other basis elements
$\tilde{r}_{i\neq n}$ when constructing an ansatz for $\tilde{r}_{n}$, then it
will suffice to supplement this set of functions with a set that parametrises
only the part of $\tilde{r}_{n}$ not contained in
$\operatorname{span}(\tilde{r}_{i\neq n})$. Naturally, there is a family of
such ansätze, with the members corresponding to different choices of terms in
the partial-fractions decomposition to include. Concretely, we consider a
family of ansätze, where each member is of the form

$\tilde{r}_{n}=\frac{\tilde{\mathcal{N}}_{n}}{\prod_{j}\mathcal{D}_{j}^{\tilde{q}_{nj}}}+\sum_{i\neq
n}c_{ni}\tilde{r}_{i},$ (66)

where $\tilde{\mathcal{N}}_{n}$ is an unknown polynomial in spinor brackets
and the $c_{ni}$ are unknown rational numbers. The exponents $\tilde{q}_{nj}$
in eq. (66) are chosen so that the involved denominator factors are a subset
of those in the common denominator form of $\tilde{r}_{n}$. Given our one-
loop-like motivation, we choose the $\tilde{q}_{nj}$ so that only one pole of
the form $\langle a|b+c|a]^{\alpha}$ is involved in each member of the family.
In the example of eq. (65), there are seven members of the family of ansätze,
each corresponding to a term in eq. (65). Taking each member in turn, we
sample the $\tilde{r}_{n}$ appropriately using finite-field kinematics
randomly generated with lips[85, 86]. We then attempt to fit $\tilde{r}_{n}$
with each ansatz, until we find that one is successful, validating our working
assumption. If any of the working assumptions we made were not valid we would
not find a successful fit, but in all cases considered in this work we find
that our assumption holds. We note that an important step in this procedure is
to parametrise polynomials in spinor brackets. This is a non-trivial exercise
due to momentum conservation and Schouten identities and we make use of the
systematic solution proposed in ref. [36].

Having constrained a single $\tilde{r}_{n}$ with an ansatz of the form given
in eq. (66), we have determined the $\tilde{\mathcal{N}}_{n}$, the
$\tilde{q}_{nj}$ and the $c_{ni}$. As we do not yet know the analytic form of
the $\tilde{r}_{i\neq n}$, we have not yet determined the analytic form of
$\tilde{r}_{n}$. Nevertheless, we can sidestep this problem by choosing to
change basis and replace $\tilde{r}_{n}$ with the first term in eq. (66). In
order to fully determine a complete set of linearly independent rational
functions, we then apply this procedure iteratively. Once this procedure is
finished we will need a different matrix $M_{ji}$ that relates our basis to
the original set of rational functions $r_{i}$ in eq. (63). This matrix can be
obtained in a similar way as $M_{ji}$ was determined.

Let us finish by briefly commenting on the benefits of our ansatz procedure.
We find that the dimension of each ansatz that we make is greatly reduced in
comparison to common denominator form. Indeed, the improved ansatz in eq. (66)
involves a numerator of much lower mass dimension than that in common
denominator form, supplemented only by a few extra parameters for the other
rational functions. In practice we find that the largest ansatz we use has
only $\mathcal{O}(4000)$ unknowns, which is approximately 10 times smaller
than the largest ansatz in common-denominator form. Importantly, as we
systematically walk through all possible choices of ansätze we can recycle the
numerical evaluations. Furthermore, given the small number of unknowns in the
improved ansatz, the many linear systems that we solve are performed at
negligible computational cost in comparison to the numerical sampling.
Finally, we note that once a successful small ansatz for an $\tilde{r}_{n}$ of
the form in eq. (66) has been found, we can efficiently further simplify the
analytic form of the $\tilde{r}_{n}$ by systematically trying simpler ansätze
with fewer poles and potentially partial fractioning.

##### Mandelstam reconstruction

We also perform the reconstruction computation in Mandelstam invariants
following the approach of ref. [47]. The pentagon-function coefficients are
split into parity-even and parity-odd parts, and are normalized by the
corresponding tree, or a spinor-weight factor in case the latter vanishes. The
denominators in Mandelstam invariants are obtained with the same procedure
described above, with the difference that a single univariate slice now
suffices since the denominators are little-group invariant (given by products
of alphabet letters). We sort the coefficients by ansatz size, and select the
simplest subset which constitutes a basis of the vector space of the rational
functions. The numerator ansatz for the most complicated parity-even or
parity-odd coefficient is now a polynomial of degree 32 in 5 independent
variables, corresponding to $58905$ unknown parameters. The number of
evaluations of the remainders required to fit the free parameters of this
ansatz is twice this number, since for each kinematic point we also require an
evaluation at the parity conjugate point in order to differentiate the parity-
even from the party-odd parts of the pentagon-function coefficients. In
conclusion, to fit the ansatz in common-denominator form the number of
required evaluations is $117810$ in Mandelstam invariants, compared to $29059$
in spinor variables. We verified numerically over both $\mathbb{C}$ and
$\mathbb{F}_{p}$ that the results obtained with the two computations match.

##  4 Results

###  4.1 Efficiency Of Analytic Reconstruction

In order to show the impact of our reconstruction procedure on the efficiency
of the calculation, we gather here representative data at various stages of
the calculation. In table 1, we display a summary of the complexity of the
amplitudes in common denominator form when making use of spinor-helicity
variables. Extending the notation of eq. (56), we denote the collective
little-group weights of $\mathcal{N}_{i}$ as $\\{\mathcal{N}_{i}\\}$. For
comparison, in the same table, we also reproduce the corresponding information
for the planar finite remainders previously obtained in ref. [48]. In table 2,
we summarize the size of the improved ansatz. For each remainder, we provide
the mass-dimension and little-group weights of the most complicated rational
function. We use this information to construct the ansätze with the algorithm
presented in ref. [36]. The Gröbner-basis computation is performed with
Singular [74]. Enumeration of the spinor monomials is performed with OR-Tools
CP-SAT [87]. The resulting dimension of the ansatz can therefore easily be
counted, and we provide this information in the final column of our tables. We
note that this cannot be determined from the mass dimension alone.

Contribution | $\text{dim}(\text{span}(r_{i}))$ | $\text{max}_{j}([\text{Num}(\tilde{r}_{j})]),\,\\{\text{Num}(\tilde{r}_{j_{max}})\\}$ |  | Common Den.
---
Ansatz Size
non-planar | $R^{(2,1)}_{-++}$ | 174 | 48, {1, -3, -6, 2, 2} | 29059
$R^{(2,\tilde{{N_{f}}})}_{-++}$ | 88 | 47, {4, 4, -5, 3, 4} | 24582
$R^{(2,1)}_{+++}$ | 49 | 21, {5, 4, 3, 3, 3} | 1092
$R^{(2,\tilde{{N_{f}}})}_{+++}$ | 24 | 20, {2, 4, 6, 6, 6} | 535
planar | $R^{(2,0)}_{-++}$ | 87 | 35, {-3, 0, 6, -3, -2} | 7358
$R^{(2,{N_{f}})}_{-++}$ | 29 | 15, {-2, -2, 0, -3, -3} | 378
$R^{(2,0)}_{+++}$ | 31 | 20, {-2, -4, -2, -2, -2} | 1140
$R^{(2,{N_{f}})}_{+++}$ | 6 | 8, {1, 3, 1, 1, 2} | 44
Table 1: Summary of rational function space in common denominator form for the non-planar finite remainders (this work), and, for reference, for the planar ones (originally calculated in ref. [48]). In the first column, we label the associated finite remainder. In the second, we give the number of linearly independent rational functions that arise in the remainder. In the third, we record the mass dimension and little group weights of the most complicated rational function in the basis. In the final column, we state the number of terms in the common-denominator ansatz for this function. Contribution | $\text{max}([\tilde{\mathcal{N}}]),\,\\{\tilde{\mathcal{N}}\\}$ |  | Improved ansatz Size
---
$R^{(2,1)}_{-++}$ | 30, {1, -3, -6, 2, 2} | 4003
$R^{(2,\tilde{{N_{f}}})}_{-++}$ | 31, {4, 4, -5, 3, 4} | 3810
$R^{(2,1)}_{+++}$ | 21, {5, 4, 3, 3, 3} | 1092
$R^{(2,\tilde{{N_{f}}})}_{+++}$ | 20, {2, 4, 6, 6, 6} | 535
Table 2: Summary of size of the function space of the numerator
$\tilde{\mathcal{N}}$ in the improved ansatz of eq. (66) that was used to
reconstruct the analytic results. In the second column, we record the mass
dimension and little-group weights of the most complicated numerator function
in the basis. In the last column, we state the number of terms in the
numerator ansatz for this function. By comparison to table 1, we see that the
improved ansatz has greatly decreased the number of required samples.

In summary, we see that reconstruction of the non-planar finite remainders
with the original common-denominator ansatz requires about four times more
samples than in the planar case. Furthermore, we see that the sampling
requirement for the non-planar computation is strongly reduced when we use the
improved ansatz of eq. (66). In fact, with $\mathcal{O}(4000)$ required
samples, the complexity is below that of the common-denominator ansatz of the
planar amplitudes. We note that for the all-plus configurations, the ansatz
size is unchanged with respect to table 1 as the most complex function does
not contain a $\langle c|a+b|c]$ pole and so is unaffected by the procedure of
section 3.3. After the analytic reconstruction, we perform further clean-up
following the partial-fraction strategies of refs. [35, 36] and simplify the
rational functions to the point where none has more than about $100$ free
coefficients. The file size of the final results is then dominated by the
matrices of rational numbers, while the rational functions are about one order
of magnitude smaller. We discuss the ancillary files in more detail in section
4.3.

Let us now briefly comment on the analytic properties of the amplitudes.
Firstly, we verified that all iterated integrals with the letter
$W_{31}=\tr_{5}$ cancel out in all finite remainders. Secondly, we find that
only 162 independent combinations of irreducible weight-four functions are
required to express $\mathcal{H}^{(2)}$, a number that is significantly lower
than the 472 required to span the full space of irreducible weight-four
pentagon functions [56]. It is also interesting to compare the functions
arising in the planar triphoton amplitudes with those for leading-color three-
jet production. We find that only 48 independent combinations of irreducible
weight-four functions are required to express the planar triphoton amplitudes
and that they span a proper subspace of those needed for three-jet production.

###  4.2 Validation

In order to validate our computation, we have performed various checks at
intermediate stages, as well as on the final result which we now discuss. To
verify the surface terms, we produced numerical reduction tables on one phase-
space point using FIRE6 [88] and checked that the surface terms reduce to
zero. The evaluation of the amplitudes with the numerical unitarity method was
carried out by the well-tested code Caravel [31]. The construction of finite
remainders was also carried out within Caravel, thus confirming the expected
pole-structure on every phase-space point. In particular we confirm that the
contribution $A^{(2,{\tilde{N}_{f}})}$ is finite, which is expected as its
coupling structure is non-zero for the first time at two loops. This
furthermore implies that this contribution is invariant under scale
variations, which we also confirm. To validate the reconstruction methodology,
we confirm that the analytic results reconstructed in spinor variables agree
with the results obtained using established Mandelstam-variable reconstruction
techniques. Furthermore, we check that the analytic results match numerical
evaluations from Caravel on finite fields of different characteristic. For the
planar contributions $R^{(2,0)}$ and $R^{(2,{N_{f}})}$, we find full agreement
between our result and a previous calculation [48]. Finally, we evaluated the
remainders in a number of collinear configurations to verify the expected
singular behavior. In particular, we validate the sign of the
$R^{(2,{\tilde{N}_{f}})}$ contribution by numerically checking that its
behavior in collinear limits is consistent with the universal expectation. To
this end, we make use of the known two-loop four-point amplitudes of ref.
[52].

###  4.3 Ancillary Files

In the ancillary files associated to this paper, we present all finite
remainders through two loops. All finite remainders are presented in the form

$R=\tilde{r}_{j}M_{ji}h_{i}\,,$ (67)

where $\tilde{r}_{j}$ are rational functions, $M_{ji}$ is a matrix of rational
numbers and the $h_{i}$ are monomials of pentagon functions. The ancillary
files are provided in Mathematica format, organized by gauge-invariant
contributions, with naming conventions based on couplings and helicities.

To aid with the understanding of the ancillary files, we provide assembly
scripts that evaluate the results of this paper and reproduce the benchmark
values in appendix B, which are also found in the file targets.m. For ease of
use, we provide these assembly scripts both in Mathematica (amp_eval.m) and in
Python (amp_eval.py) format. The Python tests can be run with pytest [89]. The
assembly scripts make use of a series of raw analytic files for each
amplitude. Each finite remainder is associated to a subfolder of
anc/amplitudes which contains the following files:

  * •

BasisCoefficients.m: The rational functions $\tilde{r}_{j}$ in eq. (67). Each
element of the list is either a rational function of spinor variables or a
list of integers. A list of integers represents a permutation of the spinors
in the previous element in the list. For a list of integers {i1, i2, i3, i4,
i5} the permutation to be applied is

$(i_{1}i_{2}i_{3}i_{4}i_{5}\rightarrow 12345).$

  * •

Matrix.m: The matrix of rational numbers $M_{ji}$ in eq. (67). The matrices
are given in the sparse coordinate list (COO) format. In this format, only the
non-zero values are stated and all unspecified entries are implicitly zero.
The non-zero entries in COO format are specified as {row index, column index}
-> value.

  * •

BasisPentagons.m: The pentagon functions in eq. (67).

Cached values for the pentagon-function monomials at the kinematic point in
appendix B (see eq. (83)) are given in FValues.m. The files Benchmark.m
contain target values for the individual pentagon-function coefficients
$r_{i}=\tilde{r}_{j}M_{ji}$ at this kinematic point.

##### Explicit example: $R^{(2,{\tilde{N}_{f}})}_{+++}$

In order to elucidate the structure of the results accompanying this paper, we
discuss in detail the structure of one of the non-planar finite remainders,
$R^{(2,{\tilde{N}_{f}})}_{+++}$. The dimension of the rational function space
is 24 (see table 1). The number of contributing pentagon-function monomials is
71. The sparse rational matrix $M$ has dimension $24\times 71$, with 177
nonzero entries. The basis of rational functions $\tilde{r}$ can be written as

$\displaystyle\begin{split}\tilde{r}=\Big{\\{}&\tilde{r}_{1}\,,\;\tilde{r}_{1}\big{|}_{453\rightarrow
345}\,,\;\tilde{r}_{1}\big{|}_{534\rightarrow
345}\,,\;\tilde{r}_{4}\,,\;\tilde{r}_{4}\big{|}_{453\rightarrow
345}\,,\;\tilde{r}_{4}\big{|}_{534\rightarrow
345}\,,\;\tilde{r}_{4}\big{|}_{543\rightarrow
345}\,,\;\tilde{r}_{4}\big{|}_{435\rightarrow 345}\,,\;\\\
&\,\,\tilde{r}_{4}\big{|}_{354\rightarrow
345}\,,\;\tilde{r}_{10}\,,\;\tilde{r}_{10}\big{|}_{453\rightarrow
345}\,,\;\tilde{r}_{10}\big{|}_{534\rightarrow
345}\,,\;\tilde{r}_{13}\,,\;\tilde{r}_{13}\big{|}_{453\rightarrow
345}\,,\;\tilde{r}_{13}\big{|}_{534\rightarrow 345}\,,\;\tilde{r}_{16}\,,\;\\\
&\,\,\tilde{r}_{16}\big{|}_{453\rightarrow
345}\,,\;\tilde{r}_{18}\,,\;\tilde{r}_{18}\big{|}_{453\rightarrow
345}\,,\;\tilde{r}_{18}\big{|}_{534\rightarrow
345}\,,\;\tilde{r}_{18}\big{|}_{543\rightarrow
345}\,,\;\tilde{r}_{18}\big{|}_{435\rightarrow
345}\,,\;\tilde{r}_{23}\,,\;\tilde{r}_{24}\Big{\\}}\,,\end{split}$ (68)

where we make explicit that many of the basis functions are obtained by
permuting momenta in previous elements. It is interesting that owing to the
symmetries of the rational function space it suffices to give 8 explicit
spinor-helicity expressions. The basis elements that generate the others under
permutations are given by

$\displaystyle\begin{split}\tilde{r}_{1}&=16\frac{[13]^{2}}{[12]\langle
45\rangle^{2}}\,,\\\\[4.2679pt] \tilde{r}_{4}&=-8\frac{\langle
12\rangle[13][34](\langle 34\rangle\langle 12\rangle+2\langle 24\rangle\langle
13\rangle)}{\langle 14\rangle^{2}[14]\langle 15\rangle\langle 25\rangle\langle
34\rangle}\,,\\\\[4.2679pt] \tilde{r}_{10}^{D}&=-24\frac{\langle
12\rangle\langle 24\rangle[34]}{\langle 14\rangle\langle 15\rangle\langle
25\rangle\langle 34\rangle}\,,\quad\tilde{r}_{10}^{S}=-16\frac{\langle
12\rangle[15][34]}{\langle 15\rangle\langle 34\rangle\langle
1|2+5|1]}\,,\\\\[4.2679pt]
\tilde{r}_{10}&=\tilde{r}_{10}^{D}+\tilde{r}_{10}^{D}\big{|}_{345\rightarrow
435}+\tilde{r}_{10}^{S}\,,\\\\[4.2679pt] \tilde{r}_{13}^{D}&=24\frac{\langle
12\rangle[15]}{\langle 13\rangle\langle 34\rangle\langle
45\rangle}-48\frac{[15]\langle 25\rangle^{2}}{\langle 24\rangle\langle
35\rangle^{2}\langle 45\rangle}+24\frac{\langle 12\rangle^{2}[15]}{\langle
13\rangle\langle 14\rangle\langle 23\rangle\langle 45\rangle}\,,\\\\[4.2679pt]
\tilde{r}_{13}^{S}&=48\frac{[15]\langle 24\rangle}{\langle
34\rangle^{2}\langle 45\rangle}-24\frac{[15]\langle 25\rangle}{\langle
34\rangle\langle 35\rangle\langle 45\rangle}\,,\\\\[4.2679pt]
\tilde{r}_{13}&=\tilde{r}_{13}^{D}+\tilde{r}_{13}^{D}\big{|}_{345\rightarrow
435}+\tilde{r}_{13}^{S}\,,\\\\[4.2679pt] \tilde{r}_{16}^{D}&=-48\frac{\langle
23\rangle^{2}[23]}{\langle 15\rangle\langle 34\rangle^{2}\langle
35\rangle}+24\frac{[23]\langle 24\rangle\langle 25\rangle}{\langle
15\rangle\langle 34\rangle\langle 45\rangle^{2}}\,,\\\\[4.2679pt]
\tilde{r}_{16}^{S}&=16\frac{\langle 12\rangle[13][45]}{\langle
13\rangle\langle 45\rangle\langle 1|2+3|1]}-48\frac{\langle 12\rangle\langle
23\rangle[23]}{\langle 14\rangle\langle 15\rangle\langle 34\rangle\langle
35\rangle}\,,\\\\[4.2679pt]
\tilde{r}_{16}&=\tilde{r}_{16}^{D}+\tilde{r}^{D}_{16}\big{|}_{345\rightarrow
354}+\tilde{r}_{16}^{S}\,,\end{split}$ (69)

where the superscripts refer to doublets (D) and singlets (S) under the
permutation transformation used to express the coefficient. We omit the three
functions $\tilde{r}_{18},\tilde{r}_{23}$ and $\tilde{r}_{24}$ as they are too
complicated to print in the text.

###  4.4 Numerical Evaluation

To facilitate the use of our results in phenomenological applications, we have
implemented our analytic expressions into the C++ library FivePointAmplitudes
[49], making use of PentagonFunctions++ [56] for the numerical evaluation of
the pentagon functions. For use in such applications, it is important that our
implementation is capable of producing numerically-stable results while
maintaining reasonable evaluation times. To demonstrate the numerical
performance of our implementation, we study the (helicity) finite remainders
and the hard function $\mathcal{H}$ on a sample of $100k$ phase-space points
generated by MATRIX v2 [90]. We adopt the phase-space definition from ref.
[13], and we set the renormalization scale to $\mu=m_{\gamma\gamma\gamma}$.

((a)) Distributions of correct digits for full $\mathcal{H}^{(2)}$ with five
active flavours, i.e. ${N_{c}}=3,$ ${N_{f}}=5$. The dashed line shows the
cumulative distribution.

((b)) Relative error of the contributions defined in eq. (20) to the hard
function $\mathcal{H}^{(2)}$.

Figure 3: Distributions of correct digits $R$ (see eq. (70)) that characterize
the numerical performance of the C++ implementation of our analytic results in
FivePointAmplitudes [49]. The phase-space sample of $100k$ phase-space points
is defined in ref. [13] and generated by MATRIX v2 [90]. The renormalization
scale is set to $\mu=m_{\gamma\gamma\gamma}$.

We characterize the numerical stability of our implementation in fig. 3, where
we plot the distributions of correct digits $R$,

$R\coloneqq-\log_{10}\absolutevalue{1-\frac{X_{\text{double}}}{X_{\text{quad}}}}\,,$
(70)

for various quantities $X$. In this way, we use a quadruple-precision
evaluation ($X_{\text{quad}}$) to calculate the accuracy of the double-
precision evaluation ($X_{\text{double}}$) on each point. To catch and correct
unstable evaluations we recycle the precision rescue system developed in ref.
[91] for the amplitudes describing three-jet production at hadron colliders.
We note however that in this case it triggers only on a few points from the
whole sample, so its effect on the evaluation time is insignificant. In fig.
3(a) we show the $R$-distribution for the hard function $\mathcal{H}$ defined
in eq. 20 for the dominant partonic process $u\bar{u}\to\gamma\gamma\gamma$
with five active massless flavors, i.e. ${N_{c}}=3,$ ${N_{f}}=5$. We observe
overall excellent numerical stability, and the average evaluation time of less
than a second on a single CPU core. We can therefore conclude that our
implementation is suitable for phenomenological applications. Since the
analytic complexity of the subleading-color contributions is notably higher
than that of the leading contributions, it is of interest to compare their
relative numerical stability. In fig. 3(b) we show the $R$-distributions for
each of the four contributions to $\mathcal{H}$ separately. Indeed, we observe
that the subleading functions $H^{(2,1)}$ and $H^{(2,{\tilde{N}_{f}})}$ are
slightly less numerically stable than the leading functions $H^{(2,0)}$,
$H^{(2,{N_{f}})}$. Nevertheless, their numerical behavior is clearly adequate
for the anticipated phenomenological applications.

Figure 4:  Numerical stability of the rational functions contributing to the
finite remainder $R^{(2,0)}_{-++}$ in two different representations: one in
terms of Mandelstam invariants and the other in terms of spinor-helicity
variables. At each phase-space point, we plot the lowest number of correct
digits $R$ over the pentagon-function coefficients $r_{i}$ of eq. 22. See the
caption of fig. 3 for details of the phase-space definition.

Finally, let us recall that we employ a spinor representation of rational
coefficients in our numerical implementation. Nevertheless, as discussed in
section 3.3, we also reconstructed the rational coefficients in terms of
Mandelstam variables. With this data at hand, it is interesting to study the
impact of the choice of variables on the numerical stability of the analytic
expressions. In fig. 4 we show the $R$-distributions for the coefficients
$r_{i}$ of pentagon functions as defined in eq. 22. Specifically, on each
phase-space point we compute the value of $R$ for each $r_{i}$ and plot
$\text{min}_{r_{i}}(R)$. As a representative example, we show this
distribution for the finite remainder $R_{-++}^{(2,0)}$, first using the
compact spinor representation of the ancillary files, and then in a
representation in terms of Mandelstam variables that has been put into a
partial fractions representation using MultivariateApart [92]. We see that the
tail of the distribution of $R$ in the spinor case is improved by about two
digits with respect to the Mandelstam case, a behavior that we find is
consistent across all helicity amplitudes. It is therefore evident that the
spinor representation improves the numerical stability compared to the
representation through Mandelstam invariants. In spite of this, for the
contribution $H^{(2,0)}$ we observe overall similar behavior as in ref. [48].
This implies that the numerical accuracy of $\mathcal{H}^{(2)}$ is limited by
the accuracy of the pentagon-function evaluation.

##  5 Subleading Color Corrections To Hard Function

In all phenomenological studies of triphoton production at hadron colliders to
date, the double-virtual NNLO QCD corrections have been included in the
leading-color approximation. Given our results for the subleading color
contributions to the double-virtual corrections, it is interesting to consider
the numerical impact of these corrections. Specifically, in refs. [12, 13] the
hard function was taken to be

$\mathcal{H}^{(2)}\;\to\;\mathcal{H}^{(2)}_{\text{l.c.}}\coloneqq\frac{{N_{c}}^{2}}{4}H^{(2,0)}\,,$
(71)

instead of the full result in eq. 20. In this section, we study the impact of
the subleading-color contributions that we have calculated in this work
relative to $\mathcal{H}^{(2)}_{\text{l.c.}}$. We focus on the dominant
partonic channel $u\bar{u}\to\gamma\gamma\gamma$, and we consider
distributions of the relative sizes of corrections with respect to the
leading-color result,

$\delta\mathcal{H}_{x}^{(2)}=\frac{\Delta\mathcal{H}_{x}^{(2)}}{\mathcal{H}^{(2)}_{\text{l.c.}}}\,,$
(72)

over the same phase-space points used in section 4.4. Here
$\Delta\mathcal{H}_{x}^{(2)}$ is one of the three corrections in eq. 20:

$\displaystyle\begin{split}\Delta\mathcal{H}_{N_{c}}^{(2)}\,&\coloneqq\,-\frac{1}{4}(H^{(2,0)}+H^{(2,1)})+\frac{1}{4{N_{c}}^{2}}H^{(2,1)}\,,\\\
\Delta\mathcal{H}_{N_{f}}^{(2)}\,&\coloneqq\,C_{F}T_{F}{N_{f}}H^{(2,{N_{f}})}\,,\\\
\Delta\mathcal{H}_{Q_{f}^{2}}^{(2)}\,&\coloneqq\,C_{F}T_{F}\left(\sum_{f=1}^{{N_{f}}}Q_{f}^{2}\right)\,H^{(2,{\tilde{N}_{f}})}\,,\end{split}$
(73)

such that
$\mathcal{H}^{(2)}=\mathcal{H}^{(2)}_{\text{l.c.}}+\Delta\mathcal{H}_{N_{c}}^{(2)}+\Delta\mathcal{H}_{N_{f}}^{(2)}+\Delta\mathcal{H}_{Q_{f}^{2}}^{(2)}$,
and
$\delta{\mathcal{H}}^{(2)}\coloneqq\mathcal{H}^{(2)}/\mathcal{H}^{(2)}_{\text{l.c.}}-1$.
We note that the $\delta\mathcal{H}_{x}$ are scheme dependent, and can vary
substantially between different schemes. In fig. 5, we plot the distribution
of the $\delta\mathcal{H}_{x}$ over phase space in the $q_{T}$,
$\overline{\text{MS}}$ and Catani schemes (see section A.2 for our scheme
definitions). The average correction size in fig. 5 should provide a
reasonable proxy for the subleading-color effects in fiducial cross sections,
whereas the shape demonstrates how the corrections vary over phase space. Our
phase-space sample and conventions are the same as that used in the studies of
numerical stability. Specifically, our 100k sample points are taken from the
phase space used in the Monte-Carlo cross-section computation of ref. [13], as
generated by MATRIX v2 [90]. We work with five active flavors, i.e. $N_{f}=5$,
and set $N_{c}=3$ and the renormalization scale to
$\mu=m_{\gamma\gamma\gamma}$.

((a)) Distribution of subleading-color corrections in the $q_{T}$ and
$\overline{\text{MS}}$ schemes.

((b)) Distribution of subleading-color corrections in the Catani scheme.

Figure 5: Distribution of separate subleading-color corrections
$\delta\mathcal{H}_{x}$ defined in eq. 73 to the the hard function
$\mathcal{H}$ alongside the total correction $\delta\mathcal{H}$. The phase-
space sample of $100k$ phase-space points is defined in ref. [13] and
generated by MATRIX v2 [90]. The renormalization scale is set to
$\mu=m_{\gamma\gamma\gamma}$.

Let us now discuss key features of these distributions. First, we consider the
impact of the choice of IR subtraction scheme and see in fig. 5 that the
$q_{T}$ and $\overline{\text{MS}}$ schemes behave similarly, and both differ
noticeably from the Catani scheme. Next, we note that a highly peaked
distribution implies a simple correction by a factor in all differential cross
sections, while a smeared distribution implies that the corrections are
observable-dependent. Focusing specifically on fig. 5(a) (the schemes used in
refs. [12, 13]), we see that all corrections are negative, but behave
differently over phase space. The correction from the ${N_{c}}$-suppressed
terms is sharply peaked at about $-10\%$, in good agreement with the common
subleading-color behavior. The $H^{(2,{\tilde{N}_{f}})}$ corrections on the
contrary demonstrate significant phase-space dependence. On average, the
combined correction is about $-35\%$, and the shape of the distribution
suggests that for some observables the correction could reach up to about
$-50\%$. Nevertheless, taking into account the smallness of the overall
contribution of $\mathcal{H}^{(2)}$ to NNLO cross sections, we expect that the
approximation employed in refs. [12, 13] should be largely adequate. Finally,
it is worth noting that in the $d\bar{d}$ partonic channel the contribution
from $H^{(2,{\tilde{N}_{f}})}$ to $\mathcal{H}^{(2)}$ is four times larger
than in $u\bar{u}$, while all other contributions are unchanged. Therefore,
$H^{(2,{\tilde{N}_{f}})}$ becomes numerically as important as
$\mathcal{H}^{(2)}_{\text{l.c.}}$ in this channel and the approximation of eq.
71 is no longer valid. However, the overall contribution of the $d\bar{d}$
channel to the cross section is highly suppressed by the ratio of quark
electric charges, so this effect is invisible.

##  6 Conclusion

In this work, we have computed the complete set of two-loop helicity
amplitudes that contribute to the NNLO QCD corrections for triphoton
production at hadron colliders, including all previously unknown non-planar
contributions. We have derived compact analytic expressions for the two-loop
finite remainders that we make available as supplementary material. In
addition, we have implemented our analytic results in a C++ library, providing
a stable and efficient numerical evaluation that is suitable for
phenomenological applications.

To carry out our computation, we employed the numerical unitarity approach as
implemented in Caravel, which we have suitably extended to take care of the
reduction of non-planar five-point integral topologies. To handle the ensuing
increase in complexity, we have proposed a new approach for the construction
of surface terms that capitalizes on formulating the associated syzygy problem
in embedding space, and on the numerical nature of the computational
framework. To reconstruct the analytic form of the results from numerical
evaluations of the amplitude over finite fields we have constructed an
optimized ansatz in spinor-helicity variables. This resulted in a more
efficient reconstruction procedure, with the most complicated ansatz having
only around 4000 free parameters. We anticipate that the advancements in both
of these technical aspects will greatly assist in computations of two-loop
amplitudes with more kinematic scales that currently pose significant
computational challenges.

Using our results, we have estimated the effect of including the subleading-
color corrections in the double-virtual contributions to the cross section. We
have found that, on average, over a representative sample of phase space, the
two-loop hard function receives an extra negative contribution of around
25-35%. The exact effect is observable and scheme dependent. Due to an overall
small contribution from double-virtual corrections to the NNLO cross sections
[12, 13], we conclude that the approximation employed in ref. [12, 13, 14] is
likely sufficient, though in practice it is important to verify this
conclusion for each observable. Nevertheless, it is clear that in view of the
excellent numerical stability and efficiency of our results, including the
subleading-color corrections to the double-virtual contributions in future
phenomenological studies of triphoton production is straightforward.

## Acknowledgements

VS gratefully acknowledges the computing resources provided by the Max Planck
Institute for Physics.

##### Funding information

The work of BP was supported by the European Union’s Horizon 2020 research and
innovation program under the Marie Sklodowska-Curie grant agreement No.896690
(LoopAnsatz). MK’s work was funded by the German Research Foundation (DFG)
within the Research Training Group GRK 2044. VS has received funding from the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme grant agreement 101019620 (ERC Advanced
Grant TOPUP).

##  Appendix A IR Renormalization

###  A.1 Conventions In Catani Scheme

The remainders $\mathcal{R}$ are defined in eq. 18 using the functions
$\mathbf{I}^{(1)}$ and $\mathbf{I}^{(2)}$ [51, 52] which we quote here for
convenience,

$\displaystyle\begin{split}\mathbf{I}^{(1)}&(\epsilon)=C_{F}\frac{e^{\gamma_{E}\epsilon}}{\Gamma(1-\epsilon)}\left(\frac{1}{\epsilon^{2}}+\frac{3}{2\epsilon}\right)\left(-\frac{s_{12}}{\mu^{2}}-\,\mathrm{i}\,0\right)^{-\epsilon},\\\
\mathbf{I}^{(2)}&(\epsilon)=\frac{1}{2}\mathbf{I}^{(1)}(\epsilon)\mathbf{I}^{(1)}(\epsilon)-\frac{\beta_{0}}{\epsilon}\mathbf{I}^{(1)}(\epsilon)+\frac{e^{-\gamma_{E}\epsilon}\Gamma(1-2\epsilon)}{\Gamma(1-\epsilon)}\left(\frac{\beta_{0}}{\epsilon}+K\right)\mathbf{I}^{(1)}(2\epsilon)-\mathbf{H}(\epsilon).\end{split}$
(74)

In $\mathbf{I}^{(2)}(\epsilon)$, we have introduced the functions

$\displaystyle\begin{split}K&=\left(\frac{67}{18}-\frac{\pi^{2}}{6}\right)C_{A}-\frac{10}{9}T_{F}N_{F},\qquad\mathbf{H}(\epsilon)=\frac{e^{\gamma_{E}\epsilon}}{2\epsilon\Gamma(1-\epsilon)}H_{q},\\\
H_{q}&=\left(\frac{\pi^{2}}{2}-6\zeta_{3}-\frac{3}{8}\right)C_{F}^{2}+\left(\frac{13}{2}\zeta_{3}+\frac{245}{216}-\frac{23}{48}\pi^{2}\right)C_{A}C_{F}+\left(\frac{\pi^{2}}{12}-\frac{25}{54}\right)T_{F}C_{F}{N_{f}}.\end{split}$
(75)

###  A.2 IR-Subtraction Scheme Change

Consider two different IR subtraction schemes, in which finite remainders are

$\mathcal{R}=\mathbf{I}\mathcal{A},\qquad\tilde{\mathcal{R}}=\tilde{\mathbf{I}}\mathcal{A}\,,$
(76)

where we suppress the helicity labels. We write the difference between the
squared finite remainders as

$\Delta\coloneqq\absolutevalue*{\tilde{\mathcal{R}}}^{2}-\absolutevalue*{\mathcal{R}}^{2}=\left(\absolutevalue*{\tilde{\mathbf{I}}}^{2}-\absolutevalue*{\mathbf{I}}^{2}\right)\absolutevalue*{\mathcal{A}}^{2}\,,$
(77)

where the absolute values do not include helicity sums. Overall colour factors
have been removed as in eq. 8 and no colour sums are implied.

Upon expansion to the second order in $\alpha_{s}$ we find

$\displaystyle\begin{split}\Delta^{(1)}&=\delta^{(1)}\absolutevalue*{\mathcal{A}^{(0)}}^{2},\qquad\delta^{(1)}\coloneqq
2\Re{\tilde{\mathbf{I}}^{(1)}-\mathbf{I}^{(1)}}\,,\\\
\Delta^{(2)}&=\bar{\delta}^{(2)}\absolutevalue*{\mathcal{A}^{(0)}}^{2}+\delta^{(1)}2\Re{\mathcal{A}^{(0)\star}\mathcal{A}^{(1)}},\quad\bar{\delta}^{(2)}\coloneqq
2\Re{\tilde{\mathbf{I}}^{(2)}-\mathbf{I}^{(2)}}-\absolutevalue*{\tilde{\mathbf{I}}^{(1)}}^{2}+\absolutevalue*{\mathbf{I}^{(1)}}^{2}\,.\end{split}$
(78)

We can rewrite $\Delta^{(2)}$ such that it is defined explicitly through
finite quantities as

$\Delta^{(2)}=\delta^{(2)}\absolutevalue*{\mathcal{A}^{(0)}}^{2}+\delta^{(1)}2\Re{\mathcal{A}^{(0)\star}\mathcal{R}^{(1)}},\qquad\delta^{(2)}\coloneqq\bar{\delta}^{(2)}+\delta^{(1)}2\Re{\mathbf{I}^{(1)}}\,.$
(79)

Using the definitions in eq. 19, and taking advantage of the factorization of
$\absolutevalue*{\mathbf{I}}^{2}$ in the sum, we can therefore write the
scheme shift for the hard functions as

$\displaystyle\begin{split}\tilde{\mathcal{H}}^{(1)}&=\mathcal{H}^{(1)}+\delta^{(1)},\\\
\tilde{\mathcal{H}}^{(2)}&=\mathcal{H}^{(2)}+\delta^{(1)}\mathcal{H}^{(1)}+\delta^{(2)}.\end{split}$
(80)

dependence of the remainder definition in eq. 76 enters in a linear way, the
conversion formulas apply equally to helicity-summed hard functions
$\mathcal{H}^{(L)}$.

Let us consider converting the finite remainders in the Catani scheme, that is
employed in this work, to the ones defined in the $\overline{\text{MS}}$
scheme. Taking the definitions of
$\tilde{\mathbf{I}}=\mathbf{I}^{\overline{\text{MS}}}$ from ref.
[53],444Explicitly, $\tilde{\mathbf{I}}^{(1)}=-\mathbf{Z}^{(1)}/2$,
$\tilde{\mathbf{I}}^{(2)}=((\mathbf{Z}^{(1)})^{2}-\mathbf{Z}^{(2)})/4$, with
the $\mathbf{Z}^{(i)}$ defined in ref. [53]. we obtain

$\displaystyle\delta^{(1)}_{\overline{\text{MS}}}$
$\displaystyle=C_{F}\left(\frac{7\pi^{2}}{6}+3l_{\mu}-l_{\mu}^{2}\right)\,,$
$\displaystyle\delta^{(2)}_{\overline{\text{MS}}}$
$\displaystyle=C_{F}^{2}\left(\frac{49\pi^{4}}{72}+\frac{7\pi^{2}}{2}l_{\mu}+\left(\frac{9}{2}-\frac{7\pi^{2}}{6}\right)l_{\mu}^{2}-3l_{\mu}^{3}+\frac{1}{2}l_{\mu}^{4}\right)$
(81)
$\displaystyle\quad+C_{F}C_{A}\left(\left(\frac{691\pi^{2}}{108}-\frac{25\pi^{4}}{144}-\frac{11\zeta_{3}}{12}\right)+\left(\frac{67}{6}-\frac{157\pi^{2}}{72}\right)l_{\mu}+\left(-\frac{233}{36}+\frac{\pi^{2}}{6}\right)l_{\mu}^{2}+\frac{11}{18}l_{\mu}^{3}\right)$
$\displaystyle\quad+C_{F}T_{F}{N_{f}}\left(-\frac{56\pi^{2}}{27}+\frac{\zeta_{3}}{3}+\left(-\frac{10}{3}+\frac{11\pi^{2}}{18}\right)l_{\mu}+\frac{19l_{\mu}^{2}}{9}-\frac{2l_{\mu}^{3}}{9}\right)\,,$

where $l_{\mu}\coloneqq\log({s_{12}}/{\mu^{2}})$.

Similarly, taking the definitions from eqs. (50-60) from ref. [93] for the
$q_{T}$ scheme we get

$\displaystyle\begin{split}\delta^{(1)}_{q_{T}}&=C_{F}\pi^{2}\,,\\\
\delta^{(2)}_{q_{T}}&=C_{F}^{2}\left(\left(-12\zeta_{3}-\frac{3}{4}+\pi^{2}\right)l_{\mu}+\frac{\pi^{4}}{2}\right)\\\
&\quad+C_{F}C_{A}\left(\left(13\zeta_{3}+\frac{245}{108}-\frac{67\pi^{2}}{24}\right)l_{\mu}-\frac{187\zeta_{3}}{36}-\frac{7\pi^{4}}{48}+\frac{1181\pi^{2}}{216}+\frac{607}{81}\right)\\\
&\quad+C_{F}T_{F}{N_{f}}\left(\left(\frac{5\pi^{2}}{6}-\frac{25}{27}\right)l_{\mu}+\frac{17\zeta_{3}}{9}-\frac{97\pi^{2}}{54}-\frac{164}{81}\right)\,.\end{split}$
(82)

##  Appendix B Reference Values

To facilitate comparison with our results, we present values for the finite
remainders evaluated at a specific kinematic point. The phase-space point we
use is:

$\displaystyle\begin{split}p_{1}^{\mu}&=\\{-0.575,-0.575,0,0\\}\,,\\\
p_{2}^{\mu}&=\\{-0.575,0.575,0,0\\}\,,\\\
p_{3}^{\mu}&=\\{0.458858239,0.405584802,0.207778343,-0.053665747\\}\,,\\\
p_{4}^{\mu}&=\\{0.231129408,-0.097079562,0.009377939,-0.209543351\\}\,,\\\
p_{5}^{\mu}&=\\{0.460012351,-0.308505239,-0.217156282,0.263209099\\}\,.\end{split}$
(83)

This phase-space point was chosen to match the values given in ref. [48] for
the kinematic invariants.

For the spinors, we follow the conventions used, for example, in ref. [85],
i.e. we take

$\lambda_{\alpha}=\begin{pmatrix}\sqrt{p^{0}+p^{3}},\\\
\frac{p^{1}+ip^{2}}{\sqrt{p^{0}+p^{3}}}\end{pmatrix}\,,\quad\text{and}\quad\tilde{\lambda}_{\dot{\alpha}}=\begin{pmatrix}\sqrt{p^{0}+p^{3}},\frac{p^{1}-ip^{2}}{\sqrt{p^{0}+p^{3}}}\end{pmatrix}\,.$
(84)

This implies that with real kinematics we have
$\lambda_{\alpha}^{*}=\tilde{\lambda}_{\dot{\alpha}}$, when the energy is
positive, and $\lambda_{\alpha}^{*}=-\tilde{\lambda}_{\dot{\alpha}}$, when the
energy is negative. In table 3 we provide the evaluations of the finite
remainders at the point of eq. (83). In addition, in tables 4, 5 and 6 we show
the evaluations of the tree, one-loop and two-loop bare amplitudes, which can
be used to derive the results of table 3 following the definitions in eqs. 13,
14, 15 and 18. The one- and two-loop bare amplitudes are normalized by the
following factors

$\Phi_{+++}=\frac{[31]\langle 12\rangle^{3}\langle 13\rangle}{\langle
14\rangle^{2}\langle 15\rangle^{2}\langle
23\rangle^{2}}\quad\text{and}\quad\Phi_{-++}=A^{(0)}_{-++}=-\frac{\langle
12\rangle\langle 23\rangle^{2}}{\langle 14\rangle\langle 15\rangle\langle
24\rangle\langle 25\rangle}\,.$ (85)

The one-loop and two-loop planar remainders reproduce the values of ref. [48].
In this work the finite remainders were not normalized by any spinor weight,
hence a little care is needed when comparing with table 6 of ref. [48]. The
last four lines are the new subleading-color, non-planar contributions.

Finite Remainder | Numerical Evaluation
---|---
$R^{(1)}_{-++}$ | $31.76842068-98.20723767\,i$
$R^{(1)}_{+++}$ | $67.16227913+20.80380252\,i$
$R^{(2,0)}_{-++}$ | $726.0944727-748.8429540\,i$
$R^{(2,0)}_{+++}$ | $1085.896384+310.5673949\,i$
$R^{(2,{N_{f}})}_{-++}$ | $-198.2921242+257.4649652\,i$
$R^{(2,{N_{f}})}_{+++}$ | $-233.4239202-112.3516487\,i$
$R^{(2,1)}_{-++}$ | $-548.9464331-12.71176526\,i$
$R^{(2,1)}_{+++}$ | $-197.0006137-822.9175634\,i$
$R^{(2,{\tilde{N}_{f}})}_{-++}$ | $35.98395348-188.6772927\,i$
$R^{(2,{\tilde{N}_{f}})}_{+++}$ | $321.0180068+179.1161201\,i$
Table 3: Finite remainder evaluations at the point of eq. (83) with the spinors defined as in eq. (84). The remainders are evaluated in the Catani scheme as defined in the text. | $\epsilon^{0}$
---|---
$A^{(0)}_{-++}$ | $6.73642828+10.02454999\,i$
$A^{(0)}_{+++}$ | 0
Table 4: Tree amplitude evaluated at the point of eq. (83) with the spinors defined as in eq. (84). | $\epsilon^{-2}$ | $\epsilon^{-1}$ | $\epsilon^{0}$ | $\epsilon^{1}$ | $\epsilon^{2}$
---|---|---|---|---|---
$A^{(1)}_{-++}\big{/}\Phi_{-++}$ |  $-1$ |  $-3.17428470$ $-3.14159265\,i$ |  $-3.43768120$ $-16.69077767\,i$ |  $-4.54236420$ $-48.29215997\,i$ |  $-28.34154957$ $-104.73071157\,i$
$A^{(1)}_{+++}\big{/}\Phi_{+++}$ | 0 | 0 |  $-122.48761401$ $-218.20999082\,i$ |  $-613.16200128$ $-1772.24966259\,i$ |  $-1264.78147357$ $-6727.58375625\,i$
Table 5: Bare, normalized one-loop amplitudes evaluated at the point of eq. (83). | $\epsilon^{-4}$ | $\epsilon^{-3}$ | $\epsilon^{-2}$ | $\epsilon^{-1}$ | $\epsilon^{0}$
---|---|---|---|---|---
$A^{(2,0)}_{-++}\big{/}\Phi_{-++}$ |  $0.5$ |  $2.25761803$ $+3.14159265\,i$ |  $-3.31724534$ $+20.90350062\,i$ |  $-55.54942677$ $+44.34772278\,i$ |  $-248.76993460$ $-87.79211642\,i$
$A^{(2,0)}_{+++}\big{/}\Phi_{+++}$ |  $0$ |  $0$ |  $122.48761401$ $+218.20999082\,i$ |  $-132.67559542$ $+2049.61318591\,i$ |  $-9927.84571218$ $+3575.60761772\,i$
$A^{(2,{N_{f}})}_{-++}\big{/}\Phi_{-++}$ |  $0$ |  $0.16666667$ |  $1.33587268$ $+1.04719755\,i$ |  $4.64626451$ $+12.87251436\,i$ |  $10.33373683$ $+83.15472522\,i$
$A^{(2,{N_{f}})}_{+++}\big{/}\Phi_{+++}$ |  $0$ |  $0$ |  $0$ |  $81.65840934$ $+145.47332721\,i$ |  $895.94750003$ $+2327.53809534\,i$
$A^{(2,1)}_{-++}\big{/}\Phi_{-++}$ |  $0.5$ |  $3.17428469$ $+3.14159265\,i$ |  $3.54092067$ $+26.66308714\,i$ |  $-38.30735107$ $+112.07323410\,i$ |  $-265.12342759$ $+331.32292317\,i$
$A^{(2,1)}_{+++}\big{/}\Phi_{+++}$ |  $0$ |  $0$ |  $122.48761401$ $+218.20999082\,i$ |  $316.44565596$ $+2849.71648558\,i$ |  $-6265.18093537$ $+17706.12335722\,i$
$A^{(2,{\tilde{N}_{f}})}_{-++}\big{/}\Phi_{-++}$ |  $0$ |  $0$ |  $0$ |  $0$ |  $-11.30451464$ $-11.18613860\,i$
$A^{(2,{\tilde{N}_{f}})}_{+++}\big{/}\Phi_{+++}$ |  $0$ |  $0$ |  $0$ |  $0$ |  $-390.31606513$ $-1248.74833374\,i$
Table 6: Bare, normalized two-loop amplitudes evaluated at the point of eq.
(83).

## References

  * [1] G. Aad _et al._ ,  _Search for new phenomena in events with at least three photons collected in $pp$ collisions at $\sqrt{s}$ = 8 TeV with the ATLAS detector_,  Eur. Phys. J. C 76(4), 210 (2016),  10.1140/epjc/s10052-016-4034-8,  [1509.05051].
  * [2] H. Denizli, K. Y. Oyulmaz and A. Senol,  _Testing for observability of Higgs effective couplings in triphoton production at FCC-hh_ ,  J. Phys. G 46(10), 105007 (2019),  10.1088/1361-6471/ab2bae,  [1901.04784].
  * [3] H. Denizli and A. Senol,  _Sensitivity Reach on Anomalous Higgs Couplings via Triphoton Production for the Post-LHC Circular High-energy Hadron Colliders_ ,  Acta Phys. Polon. B 52(12), 1377 (2021),  10.5506/APhysPolB.52.1377,  [2005.08760].
  * [4] J. A. Aguilar-Saavedra, J. M. Cano and J. M. No,  _More light on Higgs flavor at the LHC: Higgs boson couplings to light quarks through $h+\gamma$ production_,  Phys. Rev. D 103(9), 095023 (2021),  10.1103/PhysRevD.103.095023,  [2008.12538].
  * [5] A. Senol, H. Denizli and C. Helveci,  _Sensitivity of anomalous quartic couplings via tri-photon production at FCC-hh_ (2023),  [2303.14805].
  * [6] A. Zerwekh, C. Dib and R. Rosenfeld,  _Triple photon production at the Tevatron in technicolor models_ ,  Phys. Lett. B 549, 154 (2002),  10.1016/S0370-2693(02)02896-4,  [hep-ph/0207270].
  * [7] N. Toro and I. Yavin,  _Multiphotons and photon jets from new heavy vector bosons_ ,  Phys. Rev. D 86, 055005 (2012),  10.1103/PhysRevD.86.055005,  [1202.6377].
  * [8] G. Das and P. Mathews,  _Neutral Triple Vector Boson Production in Randall-Sundrum Model at the LHC_ ,  Phys. Rev. D 92(9), 094034 (2015),  10.1103/PhysRevD.92.094034,  [1507.08857].
  * [9] S. Catani, L. Cieri, D. de Florian, G. Ferrera and M. Grazzini,  _Diphoton production at hadron colliders: a fully-differential QCD calculation at NNLO_ ,  Phys. Rev. Lett. 108, 072001 (2012),  10.1103/PhysRevLett.108.072001,  [Erratum: Phys.Rev.Lett. 117, 089901 (2016)],  [1110.2375].
  * [10] J. M. Campbell, R. K. Ellis, Y. Li and C. Williams,  _Predictions for diphoton production at the LHC through NNLO in QCD_ ,  JHEP 07, 148 (2016),  10.1007/JHEP07(2016)148,  [1603.02663].
  * [11] S. Catani, L. Cieri, D. de Florian, G. Ferrera and M. Grazzini,  _Diphoton production at the LHC: a QCD study up to NNLO_ ,  JHEP 04, 142 (2018),  10.1007/JHEP04(2018)142,  [1802.02095].
  * [12] H. A. Chawdhry, M. L. Czakon, A. Mitov and R. Poncelet,  _NNLO QCD corrections to three-photon production at the LHC_ ,  JHEP 02, 057 (2020),  10.1007/JHEP02(2020)057,  [1911.00479].
  * [13] S. Kallweit, V. Sotnikov and M. Wiesemann,  _Triphoton production at hadron colliders in NNLO QCD_ ,  Phys. Lett. B 812, 136013 (2021),  10.1016/j.physletb.2020.136013,  [2010.04681].
  * [14] _Modelling of isolated multi-photon production in Monte Carlo generators in ATLAS_ ,  ATL-PHYS-PUB-2021-001 (2021).
  * [15] M. Aaboud _et al._ ,  _Measurement of the production cross section of three isolated photons in $pp$ collisions at $\sqrt{s}$ = 8 TeV using the ATLAS detector_,  Phys. Lett. B 781, 55 (2018),  10.1016/j.physletb.2018.03.057,  [1712.07291].
  * [16] R. K. Valeshabadi, M. Modarres and S. Rezaie,  _Three-photon productions within the $k_{t}$-factorization at the LHC_,  Eur. Phys. J. C 81(11), 961 (2021),  10.1140/epjc/s10052-021-09771-9,  [2107.14607].
  * [17] A. Karpishkov and V. Saleev,  _Production of three isolated photons in the parton Reggeization approach at high energies_ ,  Phys. Rev. D 106(5), 054036 (2022),  10.1103/PhysRevD.106.054036,  [2205.12773].
  * [18] M. Czakon, A. Mitov and R. Poncelet,  _Next-to-Next-to-Leading Order Study of Three-Jet Production at the LHC_ ,  Phys. Rev. Lett. 127(15), 152001 (2021),  10.1103/PhysRevLett.127.152001,  [Erratum: Phys.Rev.Lett. 129, 119901 (2022), Erratum: Phys.Rev.Lett. 129, 119901 (2022)],  [2106.05331].
  * [19] X. Chen, T. Gehrmann, E. W. N. Glover, A. Huss and M. Marcoli,  _Automation of antenna subtraction in colour space: gluonic processes_ ,  JHEP 10, 099 (2022),  10.1007/JHEP10(2022)099,  [2203.13531].
  * [20] S. Catani, S. Devoto, M. Grazzini, S. Kallweit, J. Mazzitelli and C. Savoini,  _Higgs Boson Production in Association with a Top-Antitop Quark Pair in Next-to-Next-to-Leading Order QCD_ ,  Phys. Rev. Lett. 130(11), 111902 (2023),  10.1103/PhysRevLett.130.111902,  [2210.07846].
  * [21] H. A. Chawdhry, M. Czakon, A. Mitov and R. Poncelet,  _NNLO QCD corrections to diphoton production with an additional jet at the LHC_ ,  JHEP 09, 093 (2021),  10.1007/JHEP09(2021)093,  [2105.06940].
  * [22] S. Badger, T. Gehrmann, M. Marcoli and R. Moodie,  _Next-to-leading order QCD corrections to diphoton-plus-jet production through gluon fusion at the LHC_ ,  Phys. Lett. B 824, 136802 (2022),  10.1016/j.physletb.2021.136802,  [2109.12003].
  * [23] S. Badger, C. Brønnum-Hansen, D. Chicherin, T. Gehrmann, H. B. Hartanto, J. Henn, M. Marcoli, R. Moodie, T. Peraro and S. Zoia,  _Virtual QCD corrections to gluon-initiated diphoton plus jet production at hadron colliders_ ,  JHEP 11, 083 (2021),  10.1007/JHEP11(2021)083,  [2106.08664].
  * [24] S. Badger, M. Czakon, H. B. Hartanto, R. Moodie, T. Peraro, R. Poncelet and S. Zoia,  _Isolated photon production in association with a jet pair through next-to-next-to-leading order in QCD_ (2023),  [2304.06682].
  * [25] B. Agarwal, F. Buccioni, A. von Manteuffel and L. Tancredi,  _Two-Loop Helicity Amplitudes for Diphoton Plus Jet Production in Full Color_ ,  Phys. Rev. Lett. 127(26), 262001 (2021),  10.1103/PhysRevLett.127.262001,  [2105.04585].
  * [26] G. ’t Hooft,  _A Planar Diagram Theory for Strong Interactions_ ,  Nucl. Phys. B 72, 461 (1974),  10.1016/0550-3213(74)90154-0.
  * [27] H. Ita,  _Two-loop Integrand Decomposition into Master Integrals and Surface Terms_ ,  Phys. Rev. D94(11), 116015 (2016),  10.1103/PhysRevD.94.116015,  [1510.05626].
  * [28] S. Abreu, F. Febres Cordero, H. Ita, M. Jaquier, B. Page and M. Zeng,  _Two-Loop Four-Gluon Amplitudes from Numerical Unitarity_ ,  Phys. Rev. Lett. 119(14), 142001 (2017),  10.1103/PhysRevLett.119.142001,  [1703.05273].
  * [29] S. Abreu, F. Febres Cordero, H. Ita, B. Page and M. Zeng,  _Planar Two-Loop Five-Gluon Amplitudes from Numerical Unitarity_ ,  Phys. Rev. D 97(11), 116014 (2018),  10.1103/PhysRevD.97.116014,  [1712.03946].
  * [30] S. Abreu, J. Dormans, F. Febres Cordero, H. Ita and B. Page,  _Analytic Form of Planar Two-Loop Five-Gluon Scattering Amplitudes in QCD_ ,  Phys. Rev. Lett. 122(8), 082002 (2019),  10.1103/PhysRevLett.122.082002,  [1812.04586].
  * [31] S. Abreu, J. Dormans, F. Febres Cordero, H. Ita, M. Kraus, B. Page, E. Pascual, M. S. Ruf and V. Sotnikov,  _Caravel: A C++ framework for the computation of multi-loop amplitudes with numerical unitarity_ ,  Comput. Phys. Commun. 267, 108069 (2021),  10.1016/j.cpc.2021.108069,  [2009.11957].
  * [32] A. von Manteuffel and R. M. Schabinger,  _A novel approach to integration by parts reduction_ ,  Phys. Lett. B 744, 101 (2015),  10.1016/j.physletb.2015.03.029,  [1406.4513].
  * [33] T. Peraro,  _Scattering amplitudes over finite fields and multivariate functional reconstruction_ ,  JHEP 12, 030 (2016),  10.1007/JHEP12(2016)030,  [1608.01902].
  * [34] J. Klappert and F. Lange,  _Reconstructing rational functions with FireFly_ ,  Comput. Phys. Commun. 247, 106951 (2020),  10.1016/j.cpc.2019.106951,  [1904.00009].
  * [35] G. De Laurentis and D. Maître,  _Extracting analytical one-loop amplitudes from numerical evaluations_ ,  JHEP 07, 123 (2019),  10.1007/JHEP07(2019)123,  [1904.04067].
  * [36] G. De Laurentis and B. Page,  _Ansätze for scattering amplitudes from p-adic numbers and algebraic geometry_ ,  JHEP 12, 140 (2022),  10.1007/JHEP12(2022)140,  [2203.04269].
  * [37] V. Magerya,  _Rational Tracer: a Tool for Faster Rational Function Reconstruction_ (2022),  [2211.03572].
  * [38] A. V. Belitsky, A. V. Smirnov and R. V. Yakovlev,  _Balancing act: multivariate rational reconstruction for IBP_ (2023),  [2303.02511].
  * [39] F. V. Tkachov,  _A Theorem on Analytical Calculability of Four Loop Renormalization Group Functions_ ,  Phys. Lett. B 100, 65 (1981),  10.1016/0370-2693(81)90288-4.
  * [40] K. G. Chetyrkin and F. V. Tkachov,  _Integration by Parts: The Algorithm to Calculate beta Functions in 4 Loops_ ,  Nucl. Phys. B 192, 159 (1981),  10.1016/0550-3213(81)90199-1.
  * [41] J. Gluza, K. Kajda and D. A. Kosower,  _Towards a Basis for Planar Two-Loop Integrals_ ,  Phys. Rev. D83, 045012 (2011),  10.1103/PhysRevD.83.045012,  [1009.0472].
  * [42] K. J. Larsen and Y. Zhang,  _Integration-by-parts reductions from unitarity cuts and algebraic geometry_ ,  Phys. Rev. D93(4), 041701 (2016),  10.1103/PhysRevD.93.041701,  [1511.01071].
  * [43] S. Caron-Huot and J. M. Henn,  _Iterative structure of finite loop integrals_ ,  JHEP 06, 114 (2014),  10.1007/JHEP06(2014)114,  [1404.2922].
  * [44] Z. Bern, M. Enciso, H. Ita and M. Zeng,  _Dual Conformal Symmetry, Integration-by-Parts Reduction, Differential Equations and the Nonplanar Sector_ ,  Phys. Rev. D 96(9), 096017 (2017),  10.1103/PhysRevD.96.096017,  [1709.06055].
  * [45] S. Abreu, R. Britto, C. Duhr and E. Gardi,  _Cuts from residues: the one-loop case_ ,  JHEP 06, 114 (2017),  10.1007/JHEP06(2017)114,  [1702.03163].
  * [46] S. Abreu, F. Febres Cordero, H. Ita, B. Page and V. Sotnikov,  _Planar Two-Loop Five-Parton Amplitudes from Numerical Unitarity_ ,  JHEP 11, 116 (2018),  10.1007/JHEP11(2018)116,  [1809.09067].
  * [47] S. Abreu, F. Febres Cordero, H. Ita, M. Klinkert, B. Page and V. Sotnikov,  _Leading-color two-loop amplitudes for four partons and a W boson in QCD_ ,  JHEP 04, 042 (2022),  10.1007/JHEP04(2022)042,  [2110.07541].
  * [48] S. Abreu, B. Page, E. Pascual and V. Sotnikov,  _Leading-Color Two-Loop QCD Corrections for Three-Photon Production at Hadron Colliders_ ,  JHEP 01, 078 (2021),  10.1007/JHEP01(2021)078,  [2010.15834].
  * [49] https://gitlab.com/five-point-amplitudes/FivePointAmplitudes-cpp.git.
  * [50] W. H. Furry,  _A Symmetry Theorem in the Positron Theory_ ,  Phys. Rev. 51, 125 (1937),  10.1103/PhysRev.51.125.
  * [51] C. Anastasiou, E. W. N. Glover and M. E. Tejeda-Yeomans,  _Two loop QED and QCD corrections to massless fermion boson scattering_ ,  Nucl. Phys. B 629, 255 (2002),  10.1016/S0550-3213(02)00140-2,  [hep-ph/0201274].
  * [52] E. W. N. Glover and M. E. Tejeda-Yeomans,  _Two loop QCD helicity amplitudes for massless quark massless gauge boson scattering_ ,  JHEP 06, 033 (2003),  10.1088/1126-6708/2003/06/033,  [hep-ph/0304169].
  * [53] H. A. Chawdhry, M. Czakon, A. Mitov and R. Poncelet,  _Two-loop leading-color helicity amplitudes for three-photon production at the LHC_ ,  JHEP 06, 150 (2021),  10.1007/JHEP06(2021)150,  [2012.13553].
  * [54] S. Catani,  _The Singular behavior of QCD amplitudes at two loop order_ ,  Phys. Lett. B 427, 161 (1998),  10.1016/S0370-2693(98)00332-3,  [hep-ph/9802439].
  * [55] D. Chicherin, T. Gehrmann, J. M. Henn, P. Wasser, Y. Zhang and S. Zoia,  _All Master Integrals for Three-Jet Production at Next-to-Next-to-Leading Order_ ,  Phys. Rev. Lett. 123(4), 041603 (2019),  10.1103/PhysRevLett.123.041603,  [1812.11160].
  * [56] D. Chicherin and V. Sotnikov,  _Pentagon Functions for Scattering of Five Massless Particles_ ,  JHEP 12, 167 (2020),  10.1007/JHEP12(2020)167,  [2009.07803].
  * [57] S. Abreu, L. J. Dixon, E. Herrmann, B. Page and M. Zeng,  _The two-loop five-point amplitude in $\mathcal{N}=4$ super-Yang-Mills theory_,  Phys. Rev. Lett. 122(12), 121603 (2019),  10.1103/PhysRevLett.122.121603,  [1812.08941].
  * [58] J. Klappert, F. Lange, P. Maierhöfer and J. Usovitsch,  _Integral reduction with Kira 2.0 and finite field methods_ ,  Comput. Phys. Commun. 266, 108024 (2021),  10.1016/j.cpc.2021.108024,  [2008.06494].
  * [59] T. Peraro,  _FiniteFlow: multivariate functional reconstruction using finite fields and dataflow graphs_ ,  JHEP 07, 031 (2019),  10.1007/JHEP07(2019)031,  [1905.08019].
  * [60] P. Nogueira,  _Automatic Feynman graph generation_ ,  J. Comput. Phys. 105, 279 (1993),  10.1006/jcph.1993.1074.
  * [61] A. Ochirov and B. Page,  _Full Colour for Loop Amplitudes in Yang-Mills Theory_ ,  JHEP 02, 100 (2017),  10.1007/JHEP02(2017)100,  [1612.04366].
  * [62] A. Ochirov and B. Page,  _Multi-Quark Colour Decompositions from Unitarity_ ,  JHEP 10, 058 (2019),  10.1007/JHEP10(2019)058,  [1908.02695].
  * [63] F. R. Anger and V. Sotnikov,  _On the Dimensional Regularization of QCD Helicity Amplitudes With Quarks_ (2018),  [1803.11127].
  * [64] S. Abreu, J. Dormans, F. Febres Cordero, H. Ita, B. Page and V. Sotnikov,  _Analytic Form of the Planar Two-Loop Five-Parton Scattering Amplitudes in QCD_ ,  JHEP 05, 084 (2019),  10.1007/JHEP05(2019)084,  [1904.00945].
  * [65] V. Sotnikov,  _Scattering amplitudes with the multi-loop numerical unitarity method_ ,  Ph.D. thesis, Freiburg U.,  10.6094/UNIFR/151540 (2019).
  * [66] S. Laporta,  _High precision calculation of multiloop Feynman integrals by difference equations_ ,  Int. J. Mod. Phys. A 15, 5087 (2000),  10.1142/S0217751X00002159,  [hep-ph/0102033].
  * [67] R. M. Schabinger,  _A New Algorithm For The Generation Of Unitarity-Compatible Integration By Parts Relations_ ,  JHEP 01, 077 (2012),  10.1007/JHEP01(2012)077,  [1111.4220].
  * [68] D. Cabarcas and J. Ding,  _Linear algebra to compute syzygies and gröbner bases_ ,  In _Proceedings of the 36th International Symposium on Symbolic and Algebraic Computation_ , ISSAC ’11, p. 67–74. Association for Computing Machinery, New York, NY, USA,  ISBN 9781450306751,  10.1145/1993886.1993902 (2011).
  * [69] G. Chen, J. Liu, R. Xie, H. Zhang and Y. Zhou,  _Syzygies Probing Scattering Amplitudes_ ,  JHEP 09, 075 (2016),  10.1007/JHEP09(2016)075,  [1511.01058].
  * [70] J. Böhm, A. Georgoudis, K. J. Larsen, M. Schulze and Y. Zhang,  _Complete sets of logarithmic vector fields for integration-by-parts identities of Feynman integrals_ ,  Phys. Rev. D 98(2), 025023 (2018),  10.1103/PhysRevD.98.025023,  [1712.09737].
  * [71] B. Agarwal, S. P. Jones and A. von Manteuffel,  _Two-loop helicity amplitudes for $gg\to ZZ$ with full top-quark mass effects_,  JHEP 05, 256 (2021),  10.1007/JHEP05(2021)256,  [2011.15113].
  * [72] Z. Wu, J. Boehm, R. Ma, H. Xu and Y. Zhang,  _NeatIBP 1.0, A package generating small-size integration-by-parts relations for Feynman integrals_ (2023),  [2305.08783].
  * [73] D. Simmons-Duffin,  _Projectors, Shadows, and Conformal Blocks_ ,  JHEP 04, 146 (2014),  10.1007/JHEP04(2014)146,  [1204.3894].
  * [74] W. Decker, G.-M. Greuel, G. Pfister and H. Schönemann,  _Singular 4-3-0 — A computer algebra system for polynomial computations_,  http://www.singular.uni-kl.de (2022).
  * [75] G. De Laurentis and D. Maître,  _Two-Loop Five-Parton Leading-Colour Finite Remainders in the Spinor-Helicity Formalism_ ,  JHEP 02, 016 (2021),  10.1007/JHEP02(2021)016,  [2010.14525].
  * [76] D. Chicherin, T. Gehrmann, J. M. Henn, P. Wasser, Y. Zhang and S. Zoia,  _Analytic result for a two-loop five-particle amplitude_ ,  Phys. Rev. Lett. 122(12), 121602 (2019),  10.1103/PhysRevLett.122.121602,  [1812.11057].
  * [77] D. Chicherin, J. M. Henn and G. Papathanasiou,  _Cluster algebras for Feynman integrals_ ,  Phys. Rev. Lett. 126(9), 091603 (2021),  10.1103/PhysRevLett.126.091603,  [2012.12285].
  * [78] B. Page,  _Sagex mathematica and maple schools: Lectures on finite fields and large ansätze_ ,  https://indico.desy.de/event/28075/ (2021).
  * [79] R. Britto, F. Cachazo, B. Feng and E. Witten,  _Direct proof of tree-level recursion relation in Yang-Mills theory_ ,  Phys. Rev. Lett. 94, 181602 (2005),  10.1103/PhysRevLett.94.181602,  [hep-th/0501052].
  * [80] K. Risager,  _A Direct proof of the CSW rules_ ,  JHEP 12, 003 (2005),  10.1088/1126-6708/2005/12/003,  [hep-th/0508206].
  * [81] H. Elvang, D. Z. Freedman and M. Kiermaier,  _Proof of the MHV vertex expansion for all tree amplitudes in N=4 SYM theory_ ,  JHEP 06, 068 (2009),  10.1088/1126-6708/2009/06/068,  [0811.3624].
  * [82] G. De Laurentis and B. Page,  _Constructing Compact Ansätze for Scattering Amplitudes_ ,  PoS LL2022, 038 (2022),  10.22323/1.416.0038,  [2207.10125].
  * [83] Z. Bern, L. J. Dixon and D. A. Kosower,  _One loop corrections to five gluon amplitudes_ ,  Phys. Rev. Lett. 70, 2677 (1993),  10.1103/PhysRevLett.70.2677,  [hep-ph/9302280].
  * [84] Z. Bern, L. J. Dixon and D. A. Kosower,  _One loop amplitudes for e+ e- to four partons_ ,  Nucl. Phys. B 513, 3 (1998),  10.1016/S0550-3213(97)00703-7,  [hep-ph/9708239].
  * [85] G. De Laurentis,  _Numerical techniques for analytical high-multiplicity scattering amplitudes_ ,  Ph.D. thesis, Durham U.,  URL http://etheses.dur.ac.uk/13705/ (2020).
  * [86] G. De Laurentis,  _Lips: p-adic and singular phase space_ (2023), [2305.14075].
  * [87] L. Perron and V. Furnon,  _OR-Tools_ ,  https://developers.google.com/optimization/,  Google (2023).
  * [88] A. V. Smirnov and F. S. Chuharev,  _FIRE6: Feynman Integral REduction with Modular Arithmetic_ ,  Comput. Phys. Commun. 247, 106877 (2020),  10.1016/j.cpc.2019.106877,  [1901.07808].
  * [89] H. Krekel, B. Oliveira, R. Pfannschmidt, F. Bruynooghe, B. Laugher and F. Bruhin,  _pytest 7.1_ ,  https://github.com/pytest-dev/pytest (2004).
  * [90] M. Grazzini, S. Kallweit and M. Wiesemann,  _Fully differential NNLO computations with MATRIX_ ,  Eur. Phys. J. C78(7), 537 (2018),  10.1140/epjc/s10052-018-5771-7,  [1711.06631].
  * [91] S. Abreu, F. Febres Cordero, H. Ita, B. Page and V. Sotnikov,  _Leading-color two-loop QCD corrections for three-jet production at hadron colliders_ ,  JHEP 07, 095 (2021),  10.1007/JHEP07(2021)095,  [2102.13609].
  * [92] M. Heller and A. von Manteuffel,  _MultivariateApart: Generalized partial fractions_ ,  Comput. Phys. Commun. 271, 108174 (2022),  10.1016/j.cpc.2021.108174,  [2101.08283].
  * [93] S. Catani, L. Cieri, D. de Florian, G. Ferrera and M. Grazzini,  _Universality of transverse-momentum resummation and hard factors at the NNLO_ ,  Nucl. Phys. B 881, 414 (2014),  10.1016/j.nuclphysb.2014.02.011,  [1311.1654].