aainstitutetext:  Dipartimento di Fisica and Arnold-Regge Center, Università
di Torino, and INFN, Sezione di Torino, Via P. Giuria 1, I-10125 Torino, Italy
bbinstitutetext:  Institute for Theoretical Particle Physics, KIT, Karlsruhe,
Germany ccinstitutetext:  LAPTh, CNRS – USMB, BP 110, F-74941 Annecy-le-Vieux,
France ddinstitutetext:  Physik-Institut, Universität Zürich,
Wintherturerstrasse 190, CH-8057 Zürich, Switzerland eeinstitutetext:
Cavendish Laboratory, University of Cambridge, Cambridge CB3 0HE, United
Kingdom ffinstitutetext:  Max-Planck-Institut für Physik, Werner-Heisenberg-
Institut, D-80805 München, Germany gginstitutetext:  Institute for Particle
Physics Phenomenology, Department of Physics, Durham University, Durham DH1
3LE, United Kingdom hhinstitutetext:  Dipartimento di Fisica e Astronomia,
Università di Bologna e INFN, Sezione di Bologna, via Irnerio 46, I-40126
Bologna, Italy

# Virtual QCD corrections to gluon-initiated diphoton plus jet production at
hadron colliders

Simon Badger b Christian Brønnum-Hansen c Dmitry Chicherin d Thomas Gehrmann e
Heribertus Bayu Hartanto f Johannes Henn d Matteo Marcoli g Ryan Moodie h
Tiziano Peraro a Simone Zoia<EMAIL_ADDRESS>
###### Abstract

We present an analytic computation of the gluon-initiated contribution to
diphoton plus jet production at hadron colliders up to two loops in QCD. We
reconstruct the analytic form of the finite remainders from numerical
evaluations over finite fields including all colour contributions. Compact
expressions are found using the pentagon function basis. We provide a fast and
stable implementation for the colour- and helicity-summed interference between
the one-loop and two-loop finite remainders in C++ as part of the NJet
library.

††preprint: CAVENDISH-HEP-21/10, IPPP/20/115,
MPP-2021-86, TTP21-019, ZU-TH 26/21

##  1 Introduction

Precise theoretical predictions are in high demand for the current Large
Hadron Collider (LHC) experiments which are aiming to look for tiny deviations
from the Standard Model (SM). Due to the relatively large size of the strong
coupling constant, next-to-next-to-leading order (NNLO) corrections in quantum
chromo-dynamics (QCD) are desirable for a wide variety of final state
processes. In particular, a class of $2\to 3$ scattering processes with many
kinematic scales have presented a considerable challenge to the theoretical
community and there has been a good deal of activity leading to new methods
able of overcoming their algebraic and analytic complexity Kosower:2011ty ;
Mastrolia:2011pr ; Badger:2012dp ; Zhang:2012ce ; Mastrolia:2012an ;
Mastrolia:2012wf ; Ita:2015tya ; Badger:2013gxa ; Badger:2015lda ;
Abreu:2017xsl ; Abreu:2020xvt .

The production of a pair of high energy photons is an important experimental
signature at hadron colliders and can be used for example to study the Higgs
boson through its decay to photons. The SM backgrounds are dominated by QCD
corrections and a precise description of the kinematics of these observables
requires the theoretical predictions to include perturbative information from
the production in association with additional jets. NNLO corrections to the
process $pp\to\gamma\gamma+j$, which is initiated at LO by quark-antiquark and
quark-gluon processes, have been considered a high priority for current and
future experiments for several years Badger:2016bpw ; Bendavid:2018nar ;
Amoroso:2020lgh , and were computed most recently Chawdhry:2021hkp . The Born-
level amplitude for gluon-initiated diphoton plus jet production contains a
closed quark loop coupling to both photons. Consequently, this type of process
starts to contribute to the cross section only from NNLO onwards. Owing to the
large gluon luminosity, it yields a dominant contribution to the NNLO
corrections and dominates their scale uncertainty Chawdhry:2021hkp . To
improve upon this uncertainty requires the NLO corrections to the closed
quark-loop contributions, which amount to the two-loop virtual amplitudes for
$gg\to\gamma\gamma g$ that we derive in this article. Curiously, the gluon
channel has the opposite structure to the conventional expansion in the number
of colour charges, $N_{c}$. The dominant, leading colour, contributions to the
quark-initiated process contain only planar diagrams, while in the gluon-
initiated case the leading-colour limit contains both planar and non-planar
graphs at two loops. Graphs with the highest complexity are thus unavoidable.

The last few years have seen rapid progress in our ability to compute two-loop
$2\to 3$ scattering processes in QCD which had been intractable for a long
time. The analytic computation of the scattering amplitudes in a form suitable
for phenomenological applications requires a number of major technical
bottlenecks to be overcome. A basis of special functions must be identified
that can be evaluated efficiently over the full phase space. For massless
five-point scattering, such a basis has been identified Chicherin:2017dob ;
Papadopoulos:2015jft ; Gehrmann:2018yef ; Chicherin:2018mue ;
Chicherin:2018old  and became recently available as a fast and stable
implementation in C++ valid in the physical scattering region
Chicherin:2020oor . Secondly, the amplitude must be reduced from tensor
Feynman integrals onto a basis of master integrals that can subsequently be
expanded in terms of special functions. Currently, the only viable approach to
this task is through the solution of enormous systems of integration-by-parts
(IBP) identities Tkachov:1981wb ; Chetyrkin:1981qh ; Laporta:2001dd  of which
many public implementations now exist Anastasiou:2004vj ; Studerus:2009ye ;
vonManteuffel:2012np ; Lee:2012cn ; Smirnov:2019qkx ; Klappert:2020nbg . There
has been success in simplifying this problem using syzygy relations
Gluza:2010ws ; Schabinger:2011dz ; Ita:2015tya ; Larsen:2015ped ;
Boehm:2017wjc , module intersection Boehm:2018fpv ; Boehm:2020ijp ,
intersection theory Mastrolia:2018uzb ; Frellesvig:2019uqt ;
Frellesvig:2019kgj ; Frellesvig:2020qot , $\eta$ expansion Liu:2017jxz ;
Liu:2018dmc ; Guan:2019bcx ; Zhang:2018mlo ; Wang:2019mnn , direct solution of
IBPs through recursive relations Kosower:2018obg , multivariate partial
fractioning Boehm:2020ijp , and by-passing complicated algebraic steps through
finite field arithmetic vonManteuffel:2014ixa ; Peraro:2016wsq ;
Klappert:2019emp ; Peraro:2019svx ; Klappert:2020aqs . The latter method can
be applied more broadly Peraro:2016wsq ; Peraro:2019svx , in particular to a
complete reduction of the amplitudes into a representation using special
functions. In this article, we approach the problem through a direct analytic
reconstruction of the amplitudes at the level of the pentagon functions
performing all intermediate steps numerically over finite fields. This
technique has been applied successfully to leading-colour (planar) five-parton
amplitudes first numerically Badger:2017jhb ; Abreu:2017hqn ; Badger:2018gip ;
Abreu:2018jgq  and then analytically Gehrmann:2015bfy ; Badger:2018enw ;
Abreu:2018zmy ; Abreu:2019odu ; Abreu:2021fuk . Leading-colour three-photon
production has also been completed and cross checked by two independent groups
both at the level of the amplitudes Abreu:2020cwb ; Chawdhry:2020for  and of
differential cross sections Chawdhry:2019bji ; Kallweit:2020gcp . Very
recently, NNLO QCD predictions for a number of three-jet observables and
differential three-to-two jet ratios have been computed at leading colour as
well Czakon:2021mjy . The process $gg\to g\gamma\gamma$ contains the most
complicated non-planar topologies with up to rank five tensor numerators even
at leading colour.

Diphoton production has been known at NNLO for some time Catani:2011qz ;
Campbell:2016yrh  and the two-loop scattering amplitudes were among the first
complete $2\to 2$ process to be calculated Anastasiou:2002zn ; Bern:2001df .
The first results for the amplitudes for $pp\to\gamma\gamma j$ appeared in the
last few months both for the amplitudes Agarwal:2021grm ; Chawdhry:2021mkw
and NNLO differential cross section Chawdhry:2021hkp  in the leading-colour
approximation. Very recently, the full-colour two-loop QCD corrections for the
quark-initiated channels to $pp\to\gamma\gamma j$ were presented
Agarwal:2021vdh .

We obtain sufficiently compact analytic expressions for the complete set of
helicity amplitudes for which the ultraviolet (UV) and infrared (IR) poles
have been subtracted, and implement them into an efficient and stable C++ code
as part of the NJet library Badger:2012pg . These expressions take the form of
rational coefficients multiplied by pentagon functions. The code provides
colour- and helicity-summed expressions for the two-loop amplitudes interfered
with the one-loop amplitudes, which can be used directly in phenomenological
applications.

Our paper is organised as follows. We first introduce the notation and
describe the colour decomposition of the amplitudes. We then describe the
methodology used to perform the integration-by-parts reduction and
reconstruction of the finite remainders over finite fields. In particular we
describe a method for performing a univariate partial fractioning of the
rational coefficients of the special functions on the fly. This approach can
be used inside the finite field workflow, reducing significantly the number of
sample points required to complete the analytic reconstruction and yielding
compact analytic expressions. In particular, we show explicitly some
remarkably simple analytic forms we obtained for the all-plus helicity
amplitude, which highlight its conformal properties. Finally, we present the
implementation in the NJet library Badger:2012pg  and the performance of the
code using a realistic set of phase-space points before concluding with a few
remarks on future applications of the results and methods. We also include an
appendix with some details of the momentum twistor formalism used to provide a
rational parametrisation of the kinematics.

##  2 Kinematics and amplitude conventions

We consider the production of a pair of photons in association with a gluon
from gluon fusion,

$\displaystyle g(-p_{1})+g(-p_{2})\to g(p_{3})+\gamma(p_{4})+\gamma(p_{5})\,,$
(1)

up to two-loop order in QCD. All particles are massless, $p_{i}^{2}=0$, and we
take all momenta as outgoing, so that

$\displaystyle\sum_{i=1}^{5}p_{i}=0\,.$ (2)

Without loss of generality, we assume that the external momenta $p_{i}$ live
in a four-dimensional Minkowski space-time, whereas the Feynman loop
integrations are done in $d=4-2\epsilon$ to regulate the divergences. The
kinematics are described by five independent scalar invariants, which can be
chosen as $\\{s_{12},s_{23},s_{34},s_{45},s_{15}\\}$ with
$s_{ij}=(p_{i}+p_{j})^{2}$, and a pseudo-scalar invariant,

$\displaystyle\text{tr}_{5}=4i\epsilon_{\mu\nu\rho\sigma}p_{1}^{\mu}p_{2}^{\mu}p_{3}^{\mu}p_{4}^{\mu}=[12]\langle
23\rangle[34]\langle 41\rangle-\langle 12\rangle[23]\langle 34\rangle[41]\,.$
(3)

The square of $\text{tr}_{5}$ can be expressed in terms of the scalar
invariants through the Gram determinant of the external momenta,

$\displaystyle\text{tr}_{5}^{2}=\Delta:=\text{det}\left(2p_{i}\cdot
p_{j}\right)_{i,j=1,\ldots,4}\,,$ (4)

which is a degree-4 polynomial in the $s_{ij}$. The pseudo-scalar invariant
$\text{tr}_{5}$ therefore introduces an algebraic dependence on the
kinematics, since $\text{tr}_{5}=\pm\sqrt{\Delta}$. We emphasise that the sign
of $\text{tr}_{5}$ changes under parity conjugation, which acts by flipping
the sign of the spatial momentum components,

$\displaystyle
P:\left(p_{i}^{0},\vec{p}_{i}\right)\longrightarrow\left(p_{i}^{0},-\vec{p}_{i}\right)\,,$
(5)

and under odd-signature permutations of the external momenta.

We work in the $s_{12}$ physical scattering region, which is delimited by the
requirements that all $s$-channel invariants are positive and all $t$-channel
invariants are negative,

$\displaystyle s_{12},s_{34},s_{35},s_{45}>0\,,$ (6) $\displaystyle
s_{13},s_{14},s_{15},s_{23},s_{24},s_{25}<0\,,$ (7)

together with the negativity of the Gram determinant, $\Delta<0$, which
follows from the real-valuedness of the momenta Gehrmann:2018yef .

The scattering of gluons and photons is a one-loop process at leading order.
We decompose the scattering amplitude as

$\displaystyle\mathcal{A}(1_{g},2_{g},3_{g},4_{\gamma},5_{\gamma})=g_{s}g_{e}^{2}\left(Q_{u}^{2}N_{u}+Q_{d}^{2}N_{d}\right)f^{a_{1}a_{2}a_{3}}\sum_{\ell=1}^{\infty}\left(n_{\epsilon}\frac{\alpha_{s}}{4\pi}\right)^{\ell}A^{(\ell)}(1_{g},2_{g},3_{g},4_{\gamma},5_{\gamma})\,,$
(8)

where $n_{\epsilon}=i(4\pi/\mu_{R}^{2})^{\epsilon}e^{-\epsilon\gamma_{E}}$
with $\mu_{R}$ being the renormalisation scale. In Eq. (8), $g_{s}$ and
$g_{e}$ are the strong and electromagnetic coupling constants,
$\alpha_{s}=g_{s}^{2}/(4\pi)$, $N_{q}$ and $Q_{q}$ are the number of quarks of
type $q$ and their electric charge in units of the electron charge, and
$a_{i}$ is the adjoint $\mathrm{SU}(N_{c})$ colour index of the
$i^{\text{th}}$ gluon. The one-loop amplitude can be obtained from
permutations of pure gluon scattering Dicus:1987fk ; deFlorian:1999tp .

We further expand the loop amplitudes in powers of $N_{c}$ and $n_{f}$ (the
number of light flavour fermions),

$\displaystyle A^{(1)}(1_{g},2_{g},3_{g},4_{\gamma},5_{\gamma})$
$\displaystyle=A^{(1)}_{1}(1_{g},2_{g},3_{g},4_{\gamma},5_{\gamma})\,,$ (9)
$\displaystyle A^{(2)}(1_{g},2_{g},3_{g},4_{\gamma},5_{\gamma})$
$\displaystyle=N_{c}A^{(2)}_{1}(1_{g},2_{g},3_{g},4_{\gamma},5_{\gamma})$
$\displaystyle+\frac{1}{N_{c}}A^{(2)}_{2}(1_{g},2_{g},3_{g},4_{\gamma},5_{\gamma})+n_{f}A^{(2)}_{3}(1_{g},2_{g},3_{g},4_{\gamma},5_{\gamma})\,.$

Surprisingly, the subleading-colour two-loop amplitudes contain only planar
integrals, while the leading colour contains all of the four independent
families shown in Figure 1. This pattern is the opposite to that of the quark-
initiated channels computed in Refs. Agarwal:2021grm ; Chawdhry:2021mkw ;
Agarwal:2021vdh , for which the leading-colour contributions involve only the
planar integrals and are therefore simpler to compute. Providing a prediction
for the gluon-initiated channel necessarily requires handling the most
complicated integral families. A simple analysis of the colour factors of each
of the three-gluon vertex diagrams shown in Figure 2 illustrates how this
pattern arises. Photons couple to any of the fermion propagators, and the
colour factors remain the same. It can then be seen that non-planar
contributions can come from the diagrams (a)–(c) only. Diagrams (d)–(e), which
contribute to the subleading colour, remain planar (allowing for permutations
of the external momenta).

Figure 1: Independent integral families for the $gg\to g\gamma\gamma$
amplitude. The non-planar topologies appear only in the leading-colour
amplitude.

(a) $N_{c}$

(b) $N_{c}$

(c) $N_{c}$

(d) $\frac{1}{N_{c}}$

(e) $N_{c}-\frac{1}{N_{c}}$

Figure 2: The colour factor of each diagram in the $gg\to g\gamma\gamma$
follows from the representative three-gluon, two-loop diagrams with a closed
fermion loop shown here.

In our setup, we reduce directly to the finite remainder where the UV and IR
poles have been subtracted analytically. The poles take a particularly simple
form since there is no tree-level process and the one-loop amplitudes are
finite in $\epsilon$. The one- and two-loop finite remainders are given in
terms of the bare amplitudes by Catani:1998bh ; Becher:2009qa ; Becher:2009cu
; Gardi:2009qi ; Gardi:2009zv ,

$\displaystyle F^{(1)}$
$\displaystyle=A^{(1)}(1_{g},2_{g},3_{g},4_{\gamma},5_{\gamma})\,,$ (10)
$\displaystyle F^{(2)}$
$\displaystyle=A^{(2)}(1_{g},2_{g},3_{g},4_{\gamma},5_{\gamma})-\left(I^{(1)}+\frac{3}{2}\frac{\beta_{0}}{\epsilon}\right)A^{(1)}(1_{g},2_{g},3_{g},4_{\gamma},5_{\gamma})\,,$

where $\beta_{0}=11N_{c}/3-2n_{f}/3$ and

$\displaystyle
I^{(1)}=-n_{\Gamma}(\epsilon)\,\left\\{\frac{N_{c}}{\epsilon^{2}}\left[\left(\frac{\mu_{R}^{2}}{-s_{12}}\right)^{\epsilon}+\left(\frac{\mu_{R}^{2}}{-s_{23}}\right)^{\epsilon}+\left(\frac{\mu_{R}^{2}}{-s_{31}}\right)^{\epsilon}\right]+3\frac{\gamma_{g}}{\epsilon}\right\\}\,,$
(11)

with $n_{\Gamma}(\epsilon)=e^{\epsilon\gamma_{E}}/\Gamma(1-\epsilon)$ and
$\gamma_{g}=\beta_{0}/2$ in the ’t Hooft-Veltman scheme. The logarithms
arising from the $\epsilon$-expansion of $I^{(1)}$ can be analytically
continued to the $s_{12}$ channel by adding a small positive imaginary part to
each $s_{ij}$. The $\beta_{0}$ term in the definition of the two-loop finite
remainder accounts for the strong coupling renormalisation. The finite
remainders inherit from the amplitudes the decomposition in powers of $N_{c}$
and $n_{f}$ given by Eq. (9),

$\displaystyle F^{(1)}(1_{g},2_{g},3_{g},4_{\gamma},5_{\gamma})$
$\displaystyle=F^{(1)}_{1}(1_{g},2_{g},3_{g},4_{\gamma},5_{\gamma})\,,$ (12)
$\displaystyle F^{(2)}(1_{g},2_{g},3_{g},4_{\gamma},5_{\gamma})$
$\displaystyle=N_{c}F^{(2)}_{1}(1_{g},2_{g},3_{g},4_{\gamma},5_{\gamma})$
$\displaystyle+\frac{1}{N_{c}}F^{(2)}_{2}(1_{g},2_{g},3_{g},4_{\gamma},5_{\gamma})+n_{f}F^{(2)}_{3}(1_{g},2_{g},3_{g},4_{\gamma},5_{\gamma})\,.$

Our final results are presented in the ’t Hooft-Veltman scheme, although we
make the distinction between the dimension $d$ of the loop integration and the
dimension $d_{s}=g^{\mu}_{\ \mu}$ arising from the numerator algebra.
Amplitudes with $d_{s}=2$ have a much simpler algebraic structure and contain
information that can then be used to reduce the complexity of the more
difficult $d_{s}-2$ component (see e.g. Section 4.1). The one-loop finite
remainder has only the $d_{s}=2$ term, and we expand the two-loop finite
remainder around $d_{s}=2$ as

$\displaystyle F^{(2)}_{k}$
$\displaystyle=F^{(2)}_{k;0}+F^{(2)}_{k;1}\,(d_{s}-2)\,.$ (13)

The ’t Hooft-Veltman scheme is obtained by setting $d_{s}=d=4-2\epsilon$.

##  3 Computational setup and amplitude reduction

We take a diagrammatic approach to the calculation of the amplitude along the
lines of previous work Hartanto:2019uvl ; Badger:2021owl . Here we briefly
summarise the steps and refer the reader to Ref. Badger:2021owl  for details.
All Feynman diagrams are generated using QGRAF Nogueira:1991ex  and
subsequently processed using a combination of in-house Mathematica and FORM
Kuipers:2012rf ; Ruijl:2017dtg  scripts. In total, including contributions
from ghost diagrams, we find 50 diagrams at one loop and 1527 at two loops.
Aided by the Spinney Cullen:2010jv  package to perform the ’t Hooft algebra,
the numerators are written for each independent helicity configuration. From
the loop denominator structure we assign an integral topology to each diagram.
At this point, the diagram numerators are linear combinations of monomials in
loop-momenta dependent scalar and spinor products with coefficients depending
only on external momenta. These coefficients are loaded into a dataflow graph
using FiniteFlow Peraro:2019svx . This enables numerical sampling over finite
fields, thus sidestepping analytically complicated intermediate expressions in
further steps. We rewrite loop-momenta dependent monomials into inverse
propagator denominators and a choice of irreducible scalar products (ISPs).
The required mapping of the coefficients is performed numerically within the
dataflow framework. After summing all diagrams and dropping scaleless
integrals, we arrive at an expression ready for integration-by-parts (IBP)
reduction.

The reduction to master integrals has been obtained using an improved version
of the Laporta algorithm Laporta:2001dd . For most integral families we
generated identities containing no higher power of propagators with respect to
those appearing in the amplitude, following ideas proposed in Gluza:2010ws ;
Ita:2015tya ; Larsen:2015ped . These identities have been found using the
Baikov representation of loop integrals, for which identities (i) without
higher powers of propagators and (ii) without dimension-shifted integrals can
be found by solving polynomial equations called syzygy equations. Closed form
solutions to both these constraints are separately known. Indeed, the solution
of (i) is almost trivial and the solution for (ii) has been found in Ref.
Boehm:2017wjc . The two syzygy solutions need to be combined for generating
identities that satisfy both constraints. For this purpose we used a custom
syzygy solver that implements the algorithm in Ref. Schabinger:2011dz  using
FiniteFlow Peraro:2019svx . More details on this method can be found in Refs.
Gluza:2010ws ; Ita:2015tya ; Larsen:2015ped ; Boehm:2017wjc .

For each integral family, we generated integral identities only for one
permutation of the external legs. Numerical solutions for all the permutations
contributing to an amplitude have been found by solving the systems of
equations several times, with different numerical inputs for the invariants.
Mappings between master integrals with different permutations of external legs
are applied afterwards to obtain a result in terms of a minimal set of them.

As an additional improvement, for each phase-space point evaluated on a finite
field, we reconstruct the full dependence on the dimensional regulator
$\epsilon$ of the amplitude reduced to master integrals before substituting
their expressions in terms of special functions and computing the Laurent
expansion in $\epsilon$. With this setup, fewer numerical solutions of the
integration-by-parts identities are needed in order to reconstruct analytic
results for the amplitude. This is due to the fact that the expansion of the
integrals into pentagon functions, before performing the Laurent expansion in
$\epsilon$ of the final coefficients, complicates the dependence on $\epsilon$
of the result in this intermediate stage.

To make use of the finite field arithmetic we must have a rational
parametrisation of the external kinematics. We parameterise the kinematics
using momentum twistors Hodges:2009hk ; Badger:2013gxa  where,

$\displaystyle s_{12}=x_{1}\,,$ (14) $\displaystyle s_{23}=x_{1}x_{4}\,,$
$\displaystyle
s_{34}=\frac{x_{1}}{x_{2}}(x_{4}+x_{3}x_{4}+x_{2}x_{3}x_{5}-x_{2}x_{3})\,,$
$\displaystyle s_{45}=x_{1}x_{5}\,,$ $\displaystyle
s_{15}=x_{1}x_{3}(x_{2}-x_{4}+x_{5})\,,$
$\displaystyle\text{tr}_{5}=-\frac{x_{1}^{2}}{x_{2}}\left[x_{2}x_{4}(1+2x_{3})-x_{4}(1+x_{3})(x_{4}-x_{5})+x_{2}^{2}x_{3}(-1+x_{5})\right]\,.$

We stress that the pseudo-scalar invariant $\text{tr}_{5}$, and hence the
square root of the Gram determinant $\Delta$, is a rational function of the
$x_{i}$ variables. Moreover, since $x_{1}$ is the only dimensionful variable,
we can set it to $1$ and recover the dependence on it after the reconstruction
by dimensional analysis. Further details on the momentum twistor
parameterisation are presented in Appendix A. In the following sections, we
will consider all coefficients of the special functions to be rational
functions of the variables $x_{i}$.

##  4 Analytic reconstruction over finite fields

In this section, we present three general strategies to optimise the
reconstruction over finite fields of the rational coefficients in the finite
remainders. At this stage, each component $F(x)$ of the two-loop finite
remainder is expressed as

$\displaystyle F(x)=\sum_{i}r_{i}(x)\,\text{mon}_{i}\left(f\right)\,,$ (15)

where $r_{i}$ are rational functions of the variables $x$ which parameterise
the momentum twistors, and $\text{mon}_{i}(f)$ are linearly independent
monomials of the pentagon functions. The entire chain of operations is
implemented over finite fields in the framework FiniteFlow. We therefore have
a numerical algorithm which evaluates the rational coefficients $r_{i}(x)$
modulo some prime number. The final step consists in reconstructing the
analytic expression of the rational coefficients from a sufficient number of
numerical evaluations. We employ FiniteFlow’s multi-variate functional
reconstruction algorithms, supplemented with three strategies to reduce the
number of required sample points: we determine the linear relations among the
rational coefficients and an ansatz, use univariate slices to identify the
factors belonging to another ansatz, and perform a univariate partial fraction
decomposition on the fly. In the following subsections we discuss thoroughly
each of these procedures, and their application to two-loop diphoton finite
remainders.

###  4.1 Linear relations among the rational coefficients

The representation of the finite remainders in terms of rational coefficients
and special function monomials given by Eq. (15) is in a sense not optimal.
The special function monomials in fact do not all appear independently. They
are present only in a number of independent combinations that is typically
much smaller than the total number of monomials. As a result, the rational
coefficients $r_{i}$ in the finite remainders are not linearly independent.
Expressing the finite remainders in terms of a set of linearly independent
rational coefficients not only leads to more compact expressions, but may also
simplify their reconstruction.

We can determine the linear relations among the rational coefficients
$\\{r_{i}(x)\\}$ of the special function monomials by solving a linear fit
problem,

$\displaystyle\sum_{i}a_{i}\,r_{i}(x)=0\,.$ (16)

Since the coefficients of the linear relations $a_{i}$ are rational numbers,
they require substantially fewer sample points to be reconstructed with
respect to the rational coefficients themselves. We can then use these
relations to express the rational coefficients in terms of a set of linearly
independent ones, which remain to be reconstructed. Choosing the latter to be
the simplest — i.e. the ones with the lowest polynomial degrees — may reduce
the number of sample points required for the reconstruction.

This strategy can be further refined by supplying an ansatz for the rational
coefficients. We then fit the linear relations among the rational coefficients
of the finite remainders and the coefficients of the ansatz, which we denote
by $\\{e_{j}(x)\\}$,

$\displaystyle\sum_{i}a_{i}\,r_{i}(x)+\sum_{j}b_{j}\,e_{j}(x)=0\,,$ (17)

with $a_{i},b_{j}\in\mathbb{Q}$. In the best case scenario, all the rational
coefficients $r_{i}$ can be expressed in terms of the ansatz coefficients
$e_{j}$ and no further reconstruction needs to be performed. Even when the
ansatz does not entirely cover the rational coefficients, it may still lower
the degrees of the linearly independent coefficients which have to be
reconstructed. The ansatz can be constructed from the tree-level amplitude and
the rational coefficients of the one-loop amplitudes up to order
$\epsilon^{2}$ from the analysis of the leading singularities eden2002analytic
; Britto:2004nc ; Cachazo:2008vp ; ArkaniHamed:2010gh  or from other related
amplitudes. In the diphoton case, we can use the two-loop five-gluon
amplitudes. At one loop, the $3g2\gamma$ amplitudes can be expressed in terms
of permutations of the five-gluon ones Dicus:1987fk ; deFlorian:1999tp . While
this is no longer true at two loops, we find there is an important overlap
between the rational coefficients of the $3g2\gamma$ amplitudes and those of
the five-gluon ones. We use as ansatz in the linear relations the rational
coefficients of the leading-colour two-loop five-gluon amplitudes (all two-
loop five-parton amplitudes are available analytically at leading colour
Gehrmann:2015bfy ; Dunbar:2016aux ; Badger:2018enw ; Abreu:2018zmy ;
Badger:2019djh ; Dunbar:2019fcq ; Abreu:2021fuk ; we made use of independent
results, which are being prepared for publication).

###  4.2 Matching factors on univariate slices

The pole structure of the pentagon functions is determined by the letters of
the pentagon alphabet Chicherin:2017dob . The pentagon functions (or their
discontinuities) may in fact have logarithmic singularities in the phase-space
points where one of the letters vanishes. For this reason, it is natural to
expect that the poles of the rational coefficients should be similarly linked
to the pentagon alphabet. Indeed, we observe that the denominators of the
rational coefficients in front of the pentagon functions factorise into a
product of letters of the pentagon alphabet. In other words, each rational
coefficient $r(x)$ has the form

$\displaystyle r(x)=\frac{n(x)}{\prod_{k}{\ell_{k}}^{e_{k}}(x)}\,,$ (18)

where $e_{k}$ are integers, $n(x)$ is a polynomial in the variables $x$, and
$\\{\ell_{k}\\}$ is an ansatz of factors from the pentagon alphabet. The
exponents $e_{k}$ in Eq. (18) may in general be negative, corresponding to
factors in the numerator. We use the following ansatz for the factors,

$\displaystyle\left\\{\ell_{k}(x)\right\\}=\bigl{\\{}$ $\displaystyle\langle
12\rangle\,,\langle 13\rangle\,,\langle 14\rangle\,,\langle
15\rangle\,,\langle 23\rangle\,,\langle 24\rangle\,,\langle
25\rangle\,,\langle 34\rangle\,,\langle 35\rangle\,,\langle
45\rangle\,,[12]\,,[13]\,,[14]\,,[15]\,,$ (19)
$\displaystyle[23]\,,[24]\,,[25]\,,[34]\,,[35]\,,[45]\,,s_{12}-s_{34}\,,s_{12}-s_{35}\,,s_{12}-s_{45}\,,s_{13}-s_{24}\,,$
$\displaystyle
s_{13}-s_{25}\,,s_{13}-s_{45}\,,s_{14}-s_{23}\,,s_{14}-s_{25}\,,s_{14}-s_{35}\,,s_{15}-s_{23}\,,s_{15}-s_{24}\,,$
$\displaystyle
s_{15}-s_{34}\,,s_{23}-s_{45}\,,s_{24}-s_{35}\,,s_{25}-s_{34}\,,\text{tr}_{5}\bigr{\\}}\,.$

The exponents $e_{k}$ in the ansatz (18) can be determined by reconstructing
$r(x)$ on a univariate slice modulo some prime number Abreu:2018zmy . The
univariate slice is defined by parameterising the variables by a single
parameter $t$,

$\displaystyle\\{x_{i}(t)=a_{i}+b_{i}t\\}\,,$ (20)

for constant $a_{i}$ and $b_{i}$. The latter are chosen randomly in the finite
field to avoid artificial simplifications. The dependence on $t$ is chosen to
be linear so that the degrees of the numerator and denominator of
$r(t):=r\left(x(t)\right)$ correspond to the total degrees of $r$ in $x$.
Matching the reconstructed $r(t)$ with the ansatz (18) evaluated on the same
slice allows to determine the exponents $e_{k}$ straightforwardly. With a
univariate reconstruction on just one prime field we can thus infer a lot of
information about the analytic form of the rational coefficients: the
denominators are entirely fixed, and typically some factors of the numerators
are determined as well. What remains to be reconstructed therefore requires
fewer sample points.

###  4.3 Univariate partial fraction decomposition over finite fields

Partial fraction decomposition is a standard and powerful tool for the
simplification of rational functions. The decomposition in partial fractions
is however not unique in the multivariate case. Its application to the
multivariate rational functions in scattering amplitudes is therefore not
straightforward. The necessity to simplify the rational coefficients of two-
loop five-particle scattering amplitudes has recently spurred several
approaches to handle the multivariate case efficiently Abreu:2019odu ;
Boehm:2020ijp ; Heller:2021qkz , based upon Leinartas’ algorithm
Leinartas:1978 ; Raichev:2012 . These algorithms rely on algebraic geometry
techniques, such as multivariate polynomial division and Gröbner bases, and
require the arbitrary choice of a monomial ordering.

Our main goal in this work is actually to simplify the reconstruction of the
rational coefficients over finite fields. In other words, we want to
reconstruct the rational coefficients on the fly, directly in a form which is
decomposed in partial fractions. The simplification of the ensuing analytic
expressions comes as a welcome by-product. We observe that a univariate
partial fraction decomposition is sufficient for this purpose. The advantage
is that it can be straightforwardly implemented over finite fields, avoiding
all algebraic geometry complications. The only arbitrary choice that remains
to be done is which variable to partial fraction with respect to. The latter
can be chosen by observing the impact of the partial fraction decomposition
with respect to each variable separately on the lower order amplitudes. With
the parameterisation of the kinematics in terms of momentum twistors, Eq.
(39), we find it most convenient to partial fraction with respect to $x_{4}$.

We now discuss our algorithm to reconstruct the univariate partial fraction
decomposition of a multivariate rational function $r$ from its numerical
evaluations over finite fields. The algorithm requires as input an ansatz for
the factors which may appear in the denominator of $r$. Only those factors
which depend on the variable with respect to which the partial fraction
decomposition is being performed are strictly necessary. Guessing other
factors may further simplify the reconstruction. In the application to
massless two-loop five-particle scattering amplitudes, the factor ansatz can
be inferred from the letters of the pentagon alphabet Chicherin:2017dob . We
use the factors in Eq. (19).

Let $r$ be a rational function of the variables $x=\\{x_{i}\\}_{i=1}^{n}$. In
this work the $x_{i}$’s are the momentum twistor variables defined by Eq.
(14), so $n=5$, but we outline the algorithm in general. The goal is to
decompose $r$ in partial fractions with respect to one of the variables, say
$x_{k}$. To simplify the notation, we denote the latter by $y=x_{k}$, and the
remaining variables by $\bar{x}=\\{x_{i}\\}_{i=1}^{n}\backslash x_{k}$. We may
not know the analytic expression of $r$, but we must be able to evaluate it
numerically modulo some prime number through some algorithm. Let
$\\{\ell_{i}(\bar{x},y)\\}_{i=1}^{m}$ be an ansatz for the factors which may
appear in the denominator of $r$. Without loss of generality, we assume that
the $\ell_{i}$’s are irreducible polynomials over $\mathbb{Q}$. In other
words, we assume that $r$ has the form

$\displaystyle
r\left(\bar{x},y\right)=\frac{N\left(\bar{x},y\right)}{\prod_{i=1}^{m}\ell_{i}^{e_{i}}\left(\bar{x},y\right)}\,,$
(21)

where $e_{i}\in\mathbb{Z}$, and $N\left(\bar{x},y\right)$ is a function which
depends polynomially on $y$ and rationally on $\bar{x}$. The ansatz
$\\{\ell_{i}\left(\bar{x},y\right)\\}_{i=1}^{m}$ may catch some of the factors
in the numerator of $r\left(\bar{x},y\right)$, corresponding to negative
values of the exponents $e_{i}$ in Eq. (21). This lowers the total degrees of
$N\left(\bar{x},y\right)$ and eventually simplifies its reconstruction, but is
not necessary for the partial fraction decomposition with respect to $y$.
Similarly, the ansatz may cover all the factors in the denominator of $r$, so
that $N\left(\bar{x},y\right)$ is a polynomial in $\bar{x}$ and $y$. What is
necessary for the partial fraction algorithm to work is that the ansatz
contains all the factors in the denominator of $r$ which depend on $y$. We
denote this subset by

$\displaystyle\Lambda_{y}=\bigl{\\{}i\in\\{1,\ldots,m\\}:e_{i}>0\land\text{deg}_{y}\left[\ell_{i}\left(\bar{x},y\right)\right]>0\bigr{\\}}\,,$
(22)

where $\text{deg}_{y}\left[h\right]$ is the degree in $y$ of the polynomial
$h$.

The first step consists of fixing the exponents $e_{i}$ in the ansatz (21). We
do this through the procedure discussed in Section 4.2. In the second step we
determine the degree in $y$ of the numerator $N\left(\bar{x},y\right)$ in the
ansatz (21). We recall that $N\left(\bar{x},y\right)$ is by construction
polynomial in $y$. We compute its degree in $y$ by reconstructing it on
another univariate slice, this time where only $y$ varies,

$\displaystyle\\{x_{i}(t)=a_{i}\ \forall i\neq k\,,\ y(t)=t\\}\,,$ (23)

with $a_{i}$ chosen randomly in the finite field. Clearly, the degree in $t$
of $N(t):=N\left(\bar{x}=\bar{a},y=t\right)$ gives the degree in $y$ of
$N\left(\bar{x},y\right)$. We introduce the short-hand notation

$\displaystyle
d_{N}=\text{deg}_{y}\left[N\left(\bar{x},y\right)\right]\,,\quad\qquad
d_{i}=\text{deg}_{y}\left[\ell_{i}\left(\bar{x},y\right)\right]\,,\quad\qquad
d_{\Lambda_{y}}=\sum_{i\in\Lambda_{y}}e_{i}d_{i}$ (24)

for the degrees of $N\left(\bar{x},y\right)$ and of the denominator factors
$\ell_{i}\left(\bar{x},y\right)$ in $y$.

Using the information about the factors in the denominator of $r$ and the
degree in $y$ of its numerator, we construct the following ansatz for the
partial fraction decomposition of $r$ with respect to $y$:

$\displaystyle
r\left(\bar{x},y\right)=\sum_{i\in\Lambda_{y}}\sum_{j=1}^{e_{i}}\sum_{k=0}^{d_{i}-1}\frac{U_{ijk}\left(\bar{x}\right)\,y^{k}}{\ell_{i}^{j}\left(\bar{x},y\right)}+R\left(\bar{x}\right)+\sum_{h=1}^{d_{N}-d_{\Lambda_{y}}}V_{h}\left(\bar{x}\right)\,y^{h}\,,$
(25)

where $U_{ijk}\left(\bar{x}\right)$, $R\left(\bar{x}\right)$ and
$V_{h}\left(\bar{x}\right)$ are unknown rational functions of $\bar{x}$. The
right-most term in Eq. (25) is required only if $d_{N}>d_{\Lambda_{y}}$, i.e.
only if the numerator of $r$ has a higher degree in $y$ than the denominator.

The last step of the algorithm consists of reconstructing the analytic
dependence on $\bar{x}$ of the unknown coefficients in the ansatz (25) from
the numerical evaluations of $r\left(\bar{x},y\right)$. To solve this linear
fit problem, we use the algorithm implemented in the FiniteFlow framework
Peraro:2019svx . The solution comes in the form of an algorithm which
numerically evaluates $U_{ijk}\left(\bar{x}\right)$, $R\left(\bar{x}\right)$
and $V_{h}\left(\bar{x}\right)$. The rational reconstruction may be simplified
by first reconstructing the coefficients on a univariate slice where all the
remaining variables $\bar{x}$ vary, and using that to match them with those
factors in the ansatz $\\{\ell_{i}\left(\bar{x},y\right)\\}_{i=1}^{m}$ which
depend only on $\bar{x}$. This may lower the total degrees of the functions
that need to be reconstructed.

In addition to the factors in the original ansatz
$\\{\ell_{i}\left(\bar{x},y\right)\\}_{i=1}^{m}$, the coefficients of the
partial fraction decomposition (25) may also contain spurious factors.
Consider for instance the toy example

$\displaystyle\frac{1}{(y-a)(y-b)}=\frac{1}{(a-b)(y-a)}-\frac{1}{(a-b)(y-b)}\,,$
(26)

where $a$ and $b$ are arbitrary constants such that $a\neq b$. In this
example, the inspection of the left-hand side indicates $\\{y-a,y-b\\}$ as
ansatz for the irreducible denominator factors. The partial fraction
decomposition however contains a factor of $a-b$ in the denominator, which
arises from the residue of the function at the zero of either of the
denominator factors. Clearly $a=b$ is a spurious singularity, manifestly
absent on the left-hand side and produced by the partial fraction
decomposition. In general, we can determine the potential spurious factors by
evaluating the factors in the ansatz $\ell_{i}\left(\bar{x},y\right)$ which
depend on $y$ at their zeros,

$\displaystyle\left\\{\ell_{i}\left(\bar{x},y^{*}_{k}\right)\right\\}_{i\in\Lambda_{y},\,k\in\Lambda_{y}^{1},\,i\neq
k}\,,$ (27)

where $y^{*}_{k}$ is the zero of $\ell_{k}\left(\bar{x},y\right)$,

$\displaystyle\ell_{k}\left(\bar{x},y^{*}_{k}\right)=0\,,$ (28)

and $\Lambda_{y}^{1}$ is the subset of factors which depend linearly on $y$,

$\displaystyle\Lambda_{y}^{1}=\left\\{i\in\Lambda_{y}:\text{deg}_{y}\left[\ell_{i}\left(\bar{x},y\right)\right]=1\right\\}\,.$
(29)

The restriction to zeros of linear functions of $y$ is due to the facts that
the $\ell_{i}$’s are irreducible polynomials over $\mathbb{Q}$ and that we are
factoring over $\mathbb{Q}$. The zeros of higher-degree irreducible
polynomials would introduce algebraic and/or complex dependence.

In practice, we observe that determining the spurious factors does not
simplify the reconstruction. The greatest part of the denominators of the
coefficients in the partial fraction decomposition (25) is in fact determined
by the original ansatz $\\{\ell_{i}\left(\bar{x},y\right)\\}_{i=1}^{m}$. What
remains after they are multiplied away has a total degree which is typically
lower than that of the numerators, which therefore dominates the determination
of the number of sample points required for the reconstruction. While it is
possible to determine entirely the denominators of the coefficients in Eq.
(25), it would not reduce the number of required sample points substantially,
and for this reason we refrain from doing it.

Having determined as many factors as possible in the coefficients of the
partial fraction decomposition, we multiply them away and reconstruct the
remainder using the multivariate rational reconstruction algorithms
implemented in FiniteFlow. It is important to stress that the algorithm which
evaluates the coefficients of the partial fraction decomposition contains the
solution of a linear fit. For each numerical value of $\bar{x}$, Eq. (25) is
sampled for several numerical values of $y$, roughly as many times as the
number of unknowns. This generates a linear system of equations for the
unknowns evaluated at the chosen value of $\bar{x}$. The redundant equations
are removed after the learning phase. The reconstruction on the univariate
slices in the intermediate steps of the algorithm, because it requires several
evaluations of the original functions, obviously has a higher computational
cost with respect to directly evaluating $r$. On the other hand, the
coefficients of the partial fraction decomposition depend on one fewer
variable than the original function $r$, and may have substantially lower
degrees. As a result of all these aspects, the partial fraction decomposition
may be outperformed by a direct reconstruction for simple functions, but
becomes more and more convenient as the complexity of the functions increases.

###  4.4 Summary and impact of the reconstruction strategy

The techniques discussed in the previous sections are general and can be
applied to any rational reconstruction problem, in combination or separately.
In order to reconstruct the rational coefficients of the two-loop diphoton
finite remainders we apply them consecutively as follows.

Stage 1.

We fit the linear relations among the rational coefficients with an ansatz, as
discussed in Section 4.1. We begin with the $(d_{s}-2)^{1}$ components and use
the coefficients of the two-loop leading-colour five-gluon finite remainders
as ansatz. For the $(d_{s}-2)^{0}$ components, which are more complicated, we
add to the ansatz the $(d_{s}-2)^{1}$-coefficients already reconstructed.

Stage 2.

We guess the factors from the ansatz (19) by reconstructing a univariate slice
and multiply them away, as explained in Section 4.2.

Stage 3.

We partial fraction on the fly with respect to $x_{4}$, applying the algorithm
presented in Section 4.3. The coefficients to be reconstructed after this
stage are those in the ansatz for the partial fraction decomposition (25), and
depend on one fewer variable.

Stage 4.

We reconstruct another univariate slice and perform an additional factor
guessing, as in the second stage.

The drop in the complexity of the rational coefficients after each stage for
the most complicated two-loop diphoton finite remainders, which are in the
Maximally-Helicity-Violating (MHV) configurations, is illustrated in Table 1.
As proxy for the complexity of the coefficients we use the maximal
numerator/denominator polynomial degrees, which can be evaluated by
reconstructing univariate slices as discussed in Section 4.3.

finite remainder | original | stage 1 | stage 2 | stage 3* | stage 4*
---|---|---|---|---|---
$F^{(2)}_{1;1}(1^{-},2^{-},3^{+},4^{+},5^{+})$ | $69/60$ | $28/20$ | $24/0$ | 19/10 | 11/5
$F^{(2)}_{1;0}(1^{-},2^{-},3^{+},4^{+},5^{+})$ | $78/69$ | $44/35$ | $43/0$ | $21/10$ | $16/9$
$F^{(2)}_{1;1}(1^{-},2^{+},3^{+},4^{-},5^{+})$ | $59/55$ | $30/27$ | $29/0$ | $18/15$ | $17/4$
$F^{(2)}_{1;0}(1^{-},2^{+},3^{+},4^{-},5^{+})$ | $89/86$ | $38/36$ | $38/0$ | $20/16$ | $17/3$
$F^{(2)}_{1;1}(1^{+},2^{+},3^{+},4^{-},5^{-})$ | $40/42$ | $25/27$ | $25/0$ | $15/18$ | $15/0$
$F^{(2)}_{1;0}(1^{+},2^{+},3^{+},4^{-},5^{-})$ | $66/66$ | $32/33$ | $32/0$ | $13/13$ | $12/3$
Table 1: Maximal numerator/denominator polynomial degrees of the rational
coefficients of the most complicated finite remainders at each stage of our
reconstruction strategy. The column “original” refers to the rational
coefficients prior to any optimisation. The asterisk * highlights that, after
the partial fraction decomposition in stage 3, the coefficients to be
reconstructed depend on one fewer variable.

Interestingly, we observe that the coefficients of the subleading-colour
$3g2\gamma$ two-loop finite remainders $F^{(2)}_{2}$ can be expressed in terms
of those of the leading-colour two-loop five-gluon finite remainders. The
coefficients of the leading-colour $3g2\gamma$ two-loop remainders
$F^{(2)}_{1}$ instead are not entirely fixed by the five-gluon ones, but using
the latter as ansatz in the linear relations reduces significantly the maximal
polynomial degrees of the coefficients which remain to be reconstructed.

As can be appreciated in Table 1, our strategy leads to a substantial drop in
the polynomial degrees. Furthermore, the coefficients to be reconstructed
after the partial fraction decomposition (stage 3) depend on one fewer
variable. This makes the decrease in the number of sample points required for
the reconstruction even more pronounced. The price to pay for this is that
performing the partial fraction decomposition increases the evaluation time
per point, as discussed at the end of Section 4.3. With our setup we observe
that, for the most complicated finite remainders, the evaluation times grows
roughly by one order of magnitude, while the number of sample points required
for the reconstruction decreases by two orders of magnitude. This leads to an
overall gain of roughly one order of magnitude in the reconstruction time. We
stress that the evaluation time relevant here is that of the algorithm which
evaluates the rational coefficients over finite fields, not the final
evaluation time of the finite remainders. Once the reconstruction is
completed, in fact, the rational coefficients are evaluated from their
analytic expressions. For the evaluation time of the finite remainders, we
refer to Section 6.

Our approach therefore leads to an important simplification in the
reconstruction of the rational coefficients. Moreover, the ensuing analytic
expressions are dramatically more compact. This makes them suitable for
compilation in a C++ library, an essential step for future phenomenological
applications.

##  5 Compact analytic expressions for the all-plus configuration

Prior to discussing the numerical implementation of all two-loop helicity
amplitudes, we would like to comment on the all-plus amplitude, which displays
a particularly simple analytic form. We find that the structures appearing are
closely related to those appearing in the five-gluon all-plus amplitudes at
one Bern:1993sx ; Mahlon:1993si ; Bern:1993qk ; Henn:2019mvc  and two loops
Gehrmann:2015bfy ; Dunbar:2016aux ; Badger:2019djh ; Dunbar:2019fcq . We
present the finite remainders in the expansion around $d_{s}=2$.

The all-plus amplitude is finite and rational at one loop. The finite
remainder can be written as

$\displaystyle
F_{1;0}^{(1)}(1^{+},2^{+},3^{+},4^{+}_{\gamma},5^{+}_{\gamma})=-2\frac{[45]^{2}}{\langle
12\rangle\langle 23\rangle\langle 31\rangle}\,.$ (30)

Remarkably, this amplitude is invariant under conformal transformations, and
the expression given here exhibits this property in a manifest way
Henn:2019mvc . If all masses are neglected, the SM Lagrangian is conformally
invariant. This symmetry is obscured at loop level by the appearance of scales
associated with the divergences and it is therefore rather surprising to
observe it in a one-loop amplitude. One might naïvely suppose that this is a
consequence of the finiteness of the all-plus one-loop amplitudes. Yet, the
single-minus one-loop amplitudes are equally finite, but they are not
conformally invariant. This phenomenon still calls for an explanation. These
properties are discussed in detail in Ref. Henn:2019mvc , where the authors
prove that the $n$-gluon amplitudes in QCD are conformally invariant at one
loop. Since the diphoton amplitudes can be expressed as permutations of pure
gluon scattering Dicus:1987fk ; deFlorian:1999tp  and the conformal generators
commute with permutations, all considerations regarding conformal symmetry
trivially extend to the diphoton case.

At two-loop order, the $d_{s}=2$ contribution is the only one involving
transcendental functions. Its expression is remarkably simple,

$\displaystyle
F_{1;0}^{(2)}(1^{+},2^{+},3^{+},4^{+}_{\gamma},5^{+}_{\gamma})=\frac{[45]^{2}}{\langle
12\rangle\langle 23\rangle\langle 31\rangle}\sum_{{\rm cyclic}(123)}F_{\rm
box}(s_{12},s_{23};s_{45})\,,$ (31)

where the sum runs over the cyclic permutations of $(1,2,3)$, and

$\displaystyle F_{\rm box}(s_{12},s_{23};s_{45})={\rm
Li}_{2}\left(1-\frac{s_{12}}{s_{45}}\right)+{\rm
Li}_{2}\left(1-\frac{s_{23}}{s_{45}}\right)+\log^{2}\left(\frac{s_{12}}{s_{23}}\right)+\frac{\pi^{2}}{6}$
(32)

is the finite part of the one-loop box with an off-shell leg. The analytic
continuation of the box functions to any scattering region can be easily
achieved by adding a small positive imaginary part to each two-particle
Mandelstam invariant, $s_{ij}\to s_{ij}+i0^{+}$. The other partial amplitudes
at two loops are rational,

$\displaystyle F_{1;1}^{(2)}(1^{+},2^{+},3^{+},4^{+}_{\gamma},5^{+}_{\gamma})$
$\displaystyle=-\frac{[45]^{2}}{\langle 12\rangle\langle 23\rangle\langle
31\rangle}-\frac{1}{2}F_{3;0}^{(2)}(1^{+},2^{+},3^{+},4^{+}_{\gamma},5^{+}_{\gamma})\,,$
$\displaystyle F_{2;0}^{(2)}(1^{+},2^{+},3^{+},4^{+}_{\gamma},5^{+}_{\gamma})$
$\displaystyle=0\,,$ $\displaystyle
F_{2;1}^{(2)}(1^{+},2^{+},3^{+},4^{+}_{\gamma},5^{+}_{\gamma})$
$\displaystyle=-3\frac{[45]^{2}}{\langle 12\rangle\langle 23\rangle\langle
31\rangle}-\frac{1}{2}\frac{\operatorname{\rm
tr}_{5}\left(p_{1},p_{2},p_{3},p_{4}-p_{5}\right)\langle 45\rangle}{\langle
14\rangle\langle 15\rangle\langle 24\rangle\langle 25\rangle\langle
34\rangle\langle 35\rangle}\,,$ $\displaystyle
F_{3;0}^{(2)}(1^{+},2^{+},3^{+},4^{+}_{\gamma},5^{+}_{\gamma})$
$\displaystyle=\frac{1}{3}\operatorname{\rm
tr}_{5}\left(p_{1},p_{2},p_{3},p_{4}-p_{5}\right)\sum_{{\rm
cyclic}(123)}\frac{1}{\langle 23\rangle^{2}\langle 14\rangle\langle
15\rangle\langle 45\rangle}\,,$ (33)

where $\operatorname{\rm tr}_{5}(p_{i},p_{j},p_{k},p_{l})=\operatorname{\rm
tr}(\gamma_{5}\not{p}_{i}\not{p}_{j}\not{p}_{k}\not{p}_{l})$. The peculiar
simplicity of this amplitude at two loops follows from the fact that it
vanishes at tree level and it is rational in four dimensions at one loop. The
one-loop amplitude can in fact be used as an effective on-shell vertex in
four-dimensional unitarity Bern:1994zx ; Bern:1994cg ; Dunbar:2016aux . In
this way, the cuts of the two-loop amplitude become one-loop cuts with an
insertion of the effective vertex. The one- and two-loop all-plus finite
remainders are thus treated as tree-level and one-loop objects, respectively.
As a result, the special functions appearing in the finite remainder at two
loops can have at most transcendental weight two (up to order $\epsilon^{0}$).
Moreover, the rational coefficients of the transcendental functions can be
shown through four-dimensional unitarity to be given by (permutations of) the
one-loop all-plus finite remainder. They thus inherit the symmetry under
conformal transformations from the one-loop amplitude. These beautiful
properties are manifest in our explicit expressions (31) and (30).
Complementing four-dimensional unitarity with recursion relations for the
rational terms allows to compute the two-loop all-plus finite remainders in
the purely gluonic case avoiding altogether the computation of the two-loop
integrals Dunbar:2016aux ; Dunbar:2019fcq . Some results are available even
for amplitudes involving more than five plus-helicity gluons Dunbar:2016cxp ;
Dunbar:2016gjb ; Badger:2016ozq ; Dunbar:2017nfy ; Dunbar:2020wdh ;
Dalgleish:2020mof .

Amplitudes with a single minus helicity share some of the simplicity of the
all-plus case. They also vanish at tree level, and are finite and rational at
one loop. As a result, they also have maximum transcendental weight two at two
loops. Differently from the all-plus amplitudes, however, they do not have the
structure that $F_{1;0}^{(2)}$ has uniform transcendental weight two with all
other contributions being rational. For the amplitudes with two negative
helicities, instead, the finite remainders have maximum weight two and four at
one and two loops, respectively.

##  6 Implementation and performance

The finite remainders are coded up into the NJet C++ library, which is linked
to the PentagonFunctions++ library Chicherin:2020oor  for the evaluation of
the special functions. The six independent helicity amplitudes (shown in Table
2) are permuted analytically onto the global basis of pentagon functions
defined in the $12\to 345$ scattering region to provide a complete list of 16
“mostly-plus” helicity amplitudes required for the sum. This task is performed
using the permuted coefficients from the six fully reconstructed amplitudes as
an ansatz into the linear relations so additional reconstruction time is
avoided (see Section 4.1). Having identified a global basis of pentagon
functions for the complete colour and helicity sum, we formulate the partial
amplitudes as

$F^{h}=c^{h}_{i}M^{h}_{ij}f^{h}_{j}\,,$ (34)

where $h$ is the helicity configuration, $f^{h}_{j}$ is a list of integers
corresponding to the global list of pentagon function monomials, $M^{h}_{ij}$
are sparse matrices of rational numbers, and $c^{h}_{i}$ are the independent
rational coefficients written in terms of independent polynomials in the
momentum twistor variables $x_{i}$. The pentagon function monomials are split
into parity-odd and -even components, which allows the remaining 16 “mostly-
minus” helicities to be computed by simply flipping the parity of the special
functions and applying complex conjugation to the coefficients. The colour-
and helicity-summed matrix element is constructed numerically from these
ingredients. The sparse matrix multiplication is implemented using the Eigen
library eigenweb . Evaluation with 128-bit and 256-bit floating-point numbers
(f128 and f256) is provided via the QD library qdweb . The code is available
through https://bitbucket.org/njet/njet, where we provide additional
installation instructions and example programs demonstrating the usage.

The C++ code returns the values of the one- and two-loop hard functions,
$\mathcal{H}^{(1)}$ and $\mathcal{H}^{(2)}$, obtained by squaring Eq. (8),
substituting the decomposition in $N_{c}$ and $n_{f}$ from Eq. (9),
subtracting the IR and UV poles, and finally summing over colour and helicity,

$\displaystyle\mathcal{H}$
$\displaystyle=\frac{\alpha^{2}\alpha_{s}^{3}}{(4\pi)^{5}}\left(\mathcal{H}^{(1)}+\frac{\alpha_{s}}{4\pi}\,\mathcal{H}^{(2)}\right)+\mathcal{O}(\alpha_{s}^{5})\,,$
(35) $\displaystyle\mathcal{H}^{(2)}$
$\displaystyle=N_{c}\,\mathcal{H}^{(2)}_{1}+\frac{1}{N_{c}}\,\mathcal{H}^{(2)}_{2}+n_{f}\,\mathcal{H}^{(2)}_{3}\,.$

The sum over colours for each helicity can also be returned if required. We
find the evaluation time is dominated by the special functions, particularly
when higher precision is required. In order to ensure fast and stable
numerical evaluation, we adopt the following evaluation strategy.

  1. 1.

The user-provided phase-space point is checked for the precision of the on-
shell constraints. Points are adjusted in case the precision is not acceptable
for the requested number of digits: 64-bit floating-point numbers (f64)
$\sim$15$$; f128 $\sim$31$$; f256 $\sim$62$$.

  2. 2.

The colour- and helicity-summed amplitude is computed using f64 precision at
two points which differ only by overall dimension scaling factor. After
accounting for the overall dimension of the squared amplitude, the two
evaluations should only differ due to rounding errors at intermediate stages
in the evaluation of the coefficients. This accuracy scaling test has been
used extensively at one loop. We refer to this accuracy as f64/f64 since both
coefficients and special functions use f64 precision.

  3. 3.

If the estimated number of correct digits from the scaling tests falls below a
user-defined threshold, the coefficients only are recomputed using f128
precision after the original point is corrected to f128 precision (as in step
1). We refer to this as f128/f64 precision.

  4. 4.

The scaling test is performed again and if it fails the special functions are
re-evaluated in f128 precision. This is f128/f128 precision.

  5. 5.

These steps can be repeated to obtain up to f256/f256 precision. In practice
these steps are rather expensive and unnecessary for standard phenomenological
applications, so they are omitted from our strategy.

While the dimension scaling test has been used successfully at one loop, we
need to be more careful in our applications when linking the
PentagonFunctions++ library, which also makes use of the dimension re-scaling
internally. To validate the reliability of the scaling test as an estimate of
the error of the result, we evaluate both with a direct f128/f128 computation
and via a scaling test with an error cutoff of four digits at f64/f64 for a
set of $60\,000$ points. To ensure a realistic validation, we use “physical”
points with a phase-space sampling density determined by the one-loop process,
obtained from NNLOJET. We compare the estimated error provided by the f64/f64
scaling test to the relative difference between the f64/f64 and f128/f128
evaluations, which is taken as the true error. In the following, percentages
are always with respect to the entire set of points.

The scaling test returns a negative for $2.8\text{\,}\mathrm{\char 37\relax}$
of the points. According to true error, an additional
$0.2\text{\,}\mathrm{\char 37\relax}$ of the points should be failed and are
missed by the scaling test (false positive). Of these points, almost all have
true error of four digits, the remaining $0.008\text{\,}\mathrm{\char
37\relax}$ with three digits, so the effect on stability is small. The scaling
test also fails some points unnecessarily (false negative), this subset
comprising $0.7\text{\,}\mathrm{\char 37\relax}$ of all points, which incurs a
small performance penalty in the evaluation strategy. The effects of the false
estimates are considered to be allowably small.

We note that the dimension scaling test is statistical and therefore one will
always find anomalies in a sufficiently large sample. Care should be taken
when integrating over extreme regions of phase space.

To assess the stability of our implementation (Figure 3) and measure timings,
we evaluate the amplitude squared over $100\,000$ points of the physical phase
space. We find a single f64/f64 call has a mean time of
$9\text{\,}\mathrm{s}$, with $99\text{\,}\mathrm{\char 37\relax}$ of that time
spent evaluating the pentagon functions. Using the full evaluation strategy
with a target minimum accuracy of three digits, we obtain a mean timing per
phase-space point of $26\text{\,}\mathrm{s}$.

We present a benchmark evaluation at a point taken from the physical phase
space. We choose a generic configuration where the momentum invariants
($$\mathrm{G}\mathrm{e}\mathrm{V}$^{2}$) and $\mathrm{tr}_{5}$
($$\mathrm{G}\mathrm{e}\mathrm{V}$^{4}$) take the values, quoted to four
significant figures,

$\displaystyle s_{12}$
$\displaystyle=$1.411\,625\,116\,335\,087\,737\,723\,620\,426\,490\,824\,229\,866\,747\,189\,324\,633\,618\,076\,609\,842\,3\text{\times}{10}^{04}$\,,$
$\displaystyle s_{23}$
$\displaystyle=$-1.404\,683\,473\,797\,232\,141\,566\,057\,759\,505\,570\,635\,159\,964\,600\,893\,962\,013\,660\,115\,990\,2\text{\times}{10}^{03}$\,,$
$\displaystyle s_{34}$
$\displaystyle=$7.666\,979\,944\,894\,607\,509\,160\,964\,631\,355\,036\,312\,271\,936\,864\,719\,917\,660\,163\,425\,977\,2\text{\times}{10}^{03}$\,,$
(36) $\displaystyle s_{45}$
$\displaystyle=$5.493\,245\,056\,556\,116\,203\,228\,130\,479\,601\,100\,098\,132\,017\,248\,174\,921\,334\,783\,005\,715\,5\text{\times}{10}^{03}$\,,$
$\displaystyle s_{15}$
$\displaystyle=$-4.404\,428\,924\,591\,748\,610\,985\,744\,586\,989\,902\,339\,753\,571\,485\,749\,671\,328\,605\,364\,299\,6\text{\times}{10}^{03}$\,,$
$\displaystyle\mathrm{tr}_{5}$
$\displaystyle=$-1.759\,975\,575\,068\,791\,647\,957\,636\,722\,918\,285\,134\,775\,353\,855\,788\,821\,393\,249\,387\,087\,6\text{\times}{10}^{07}$\mathrm{i}\,.$

High precision f128/f128 evaluations are given in the ancillary files. The
values for the finite remainders and the two-loop hard function, normalised by
the leading order, are shown in Tables 2 and 3 respectively. The subleading-
colour corrections are $645$ times smaller than the leading colour at this
point, while the closed fermion loop corrections are $133$ times smaller.
These ratios do change as we sample different points. Averaging over $100$
physical points, the ratio $\lvert
N_{c}\mathcal{H}^{(2)}_{1}/\mathcal{H}^{(1)}\rvert:\lvert\frac{1}{N_{c}}\mathcal{H}^{(2)}_{2}/\mathcal{H}^{(1)}\rvert:\lvert
n_{f}\mathcal{H}^{(2)}_{3}/\mathcal{H}^{(1)}\rvert$ is $2061$:$1$:$14$.

While the evaluation is considerably more difficult than the massless planar
five-gluon scattering owing to the more complicated set of pentagon functions,
our tests show the amplitudes are clearly ready for phenomenological
applications.

helicity | $N_{c}F^{(2)}_{1}/F^{(1)}$ | $\frac{1}{N_{c}}F^{(2)}_{2}/F^{(1)}$ | $n_{f}F^{(2)}_{3}/F^{(1)}$
---|---|---|---
$+++++$ |  $$-2.775\,829\,930\,465\,834\,276\,051\,425\,152\,352\,85\text{\times}{10}^{01}$-$1.017\,450\,762\,905\,025\,192\,080\,363\,744\,770\,10\text{\times}{10}^{01}$$i |  $$-1.673\,274\,842\,489\,072\,398\,154\,590\,327\,135\,80\text{\times}{10}^{00}$-$2.396\,491\,533\,226\,016\,932\,833\,718\,050\,223\,56\text{\times}{10}^{-01}$$i |  $$-5.228\,374\,268\,069\,778\,382\,134\,616\,922\,552\,51\text{\times}{10}^{00}$-$4.034\,281\,977\,353\,245\,802\,041\,120\,025\,534\,37\text{\times}{10}^{00}$$i
$-++++$ |  $$-2.575\,805\,326\,232\,805\,030\,925\,407\,712\,940\,86\text{\times}{10}^{01}$+$2.783\,473\,491\,756\,908\,905\,105\,599\,051\,504\,41\text{\times}{10}^{01}$$i |  $$3.570\,816\,316\,335\,509\,993\,142\,358\,649\,040\,20\text{\times}{10}^{-01}$-$3.213\,389\,187\,883\,823\,411\,550\,480\,184\,420\,21\text{\times}{10}^{-01}$$i |  $$3.362\,681\,184\,972\,715\,526\,330\,549\,968\,893\,30\text{\times}{10}^{-01}$-$4.424\,013\,780\,023\,894\,103\,113\,875\,257\,740\,25\text{\times}{10}^{00}$$i
$+++-+$ |  $$-2.415\,664\,603\,584\,738\,625\,575\,975\,658\,080\,62\text{\times}{10}^{01}$+$1.458\,802\,752\,491\,222\,975\,396\,660\,991\,299\,78\text{\times}{10}^{01}$$i |  $$3.697\,561\,404\,990\,543\,595\,202\,827\,276\,397\,44\text{\times}{10}^{-01}$-$5.539\,324\,451\,257\,584\,084\,765\,220\,175\,645\,49\text{\times}{10}^{-01}$$i |  $$-4.951\,234\,418\,255\,203\,478\,678\,412\,726\,159\,52\text{\times}{10}^{00}$+$6.671\,676\,894\,421\,552\,783\,622\,279\,399\,206\,01\text{\times}{10}^{-01}$$i
$--+++$ |  $$-2.023\,425\,310\,625\,796\,712\,228\,608\,541\,292\,19\text{\times}{10}^{01}$+$8.203\,846\,865\,583\,062\,635\,654\,005\,830\,035\,10\text{\times}{10}^{-01}$$i |  $$-4.054\,874\,106\,792\,913\,778\,381\,461\,703\,511\,51\text{\times}{10}^{-01}$-$3.548\,677\,881\,563\,614\,513\,811\,900\,483\,997\,81\text{\times}{10}^{-01}$$i |  $$5.355\,181\,880\,401\,925\,548\,939\,638\,552\,932\,05\text{\times}{10}^{-02}$+$2.478\,302\,301\,936\,155\,849\,230\,729\,853\,229\,98\text{\times}{10}^{-04}$$i
$-++-+$ |  $$-2.858\,441\,610\,953\,695\,080\,046\,267\,507\,973\,06\text{\times}{10}^{01}$+$3.290\,079\,647\,186\,568\,540\,302\,543\,739\,896\,09\text{\times}{10}^{01}$$i |  $$3.917\,329\,058\,578\,149\,481\,442\,443\,629\,259\,73\text{\times}{10}^{-01}$-$5.488\,698\,247\,220\,500\,419\,957\,286\,870\,064\,86\text{\times}{10}^{-04}$$i |  $$3.021\,564\,000\,286\,861\,279\,021\,431\,513\,096\,53\text{\times}{10}^{00}$+$1.475\,497\,728\,853\,030\,429\,000\,359\,527\,990\,06\text{\times}{10}^{00}$$i
$+++--$ |  $$-2.094\,177\,077\,703\,654\,792\,102\,444\,295\,354\,50\text{\times}{10}^{01}$-$1.534\,273\,031\,454\,048\,296\,629\,736\,311\,868\,51\text{\times}{10}^{01}$$i |  $$-3.079\,681\,646\,858\,099\,804\,886\,346\,025\,152\,62\text{\times}{10}^{-01}$-$4.557\,629\,378\,166\,883\,217\,550\,496\,235\,983\,66\text{\times}{10}^{-01}$$i |  $$-4.879\,848\,569\,201\,212\,005\,975\,476\,436\,907\,81\text{\times}{10}^{00}$-$5.862\,226\,288\,798\,904\,158\,825\,134\,416\,686\,61\text{\times}{10}^{-03}$$i
Table 2: Numerical values of the partial amplitudes for the six independent helicities at the benchmark point in Eq. (36). Values are quoted with $N_{c}=3$ and $n_{f}=5$, to four significant figures. $N_{c}\mathcal{H}^{(2)}_{1}/\mathcal{H}^{(1)}$ | $\frac{1}{N_{c}}\mathcal{H}^{(2)}_{2}/\mathcal{H}^{(1)}$ | $n_{f}\mathcal{H}^{(2)}_{3}/\mathcal{H}^{(1)}$
---|---|---
$52.745\,939\,559\,081\,587\,670\,147\,724\,550\,04$ | $0.081\,762\,809\,080\,571\,194\,970\,251\,870\,265\,06$ | $0.395\,641\,911\,642\,710\,567\,484\,205\,068\,631\,81$
Table 3:  Numerical values for the components of the two-loop hard function
normalised to the one-loop hard function defined in Eq. (35) at the benchmark
point of Eq. (36). Values are quoted with $N_{c}=3$ and $n_{f}=5$, to four
significant figures. Figure 3:  Histogram of the error estimate on the two-
loop evaluations as given by the scaling test. We use the evaluation strategy
with a target accuracy of three digits and show errors for all precision
levels. We see $1.8\text{\,}\mathrm{\char 37\relax}$ of points failing f64/f64
evaluation, with $1.2\text{\,}\mathrm{\char 37\relax}$ passing at f128/f64 and
$0.6\text{\,}\mathrm{\char 37\relax}$ passing at f128/f128. The evaluation
strategy achieves target accuracy for all of the $100\,000$ physical phase-
space points tested.

##  7 Conclusions

In this paper we have presented a complete, full colour, five-point amplitude
at two loops in QCD. All helicity configurations have been implemented into
the NJet C++ library, which provides an efficient and stable evaluation over
the physical scattering region. Though the algebraic complexity of the
amplitude is considerable, the direct analytic reconstruction of the finite
remainders was possible by making use of linear relation amongst the
coefficients and partial fractioning in one variable, which could be done
without any analytic knowledge of the intermediate steps in the reduction. We
expect these techniques will have applications to other important high-
multiplicity two-loop calculations with more external scales such as five-
particle scattering with an off-shell leg, for which there has also been
recent progress Hartanto:2019uvl ; Papadopoulos:2015jft ; Papadopoulos:2019iam
; Abreu:2020jxa ; Canko:2020ylt ; Syrrakos:2020kba ; Badger:2021nhg . We have
found a form that is suitable for phenomenological applications and look
forward to new precision predictions for diphoton production at hadron
colliders including the dominant N3LO corrections we have computed here.

###### Acknowledgements.

We are grateful to Vasily Sotnikov for useful discussions. This project
received funding from the European Union’s Horizon 2020 research and
innovation programmes New level of theoretical precision for LHC Run 2 and
beyond (grant agreement No 683211), High precision multi-jet dynamics at the
LHC (grant agreement No 772009), and Novel structures in scattering amplitudes
(grant agreement No 725110), and from the Swiss National Science Foundation
(SNF) under contract 200020-175595. HBH was partially supported by STFC
consolidated HEP theory grant ST/T000694/1. DC is supported by the French
National Research Agency in the framework of the “Investissements d’avenir”
program (ANR-15-IDEX-02). RM was supported by STFC ST/S505365/1 and
ST/P001246/1. SZ gratefully acknowledges the computing resources provided by
the Max Planck Institute for Physics and by the Max Planck Computing & Data
Facility.

##  Appendix A Momentum twistor parametrisation

Following Hodges:2009hk ; Badger:2013gxa ; Badger:2017jhb , the construction
begins with

$\displaystyle Z_{i}=\begin{pmatrix}\lambda_{i}\\\ \mu_{i}\end{pmatrix}\,,$
(37)

where $\lambda_{i}$ is the negative-helicity spinor, and $\mu_{i}$ is related
to the positive-helicity spinor $\tilde{\lambda}_{i}$ via

$\displaystyle\tilde{\lambda}_{i}=\frac{\langle
i,i+1\rangle\,\mu_{i-1}+\langle i+1,i-1\rangle\,\mu_{i}+\langle
i-1,i\rangle\,\mu_{i+1}}{\langle i,i+1\rangle\langle i-1,i\rangle}\,,$ (38)

with the indices defined modulo $5$. Using the Poincaré and $U(1)$ symmetries
it is possible to fix all but $5$ of the entries of the momentum twistor
matrix $Z=(Z_{i})_{i=1,\ldots,5}$. Explicitly we choose the form,

$\displaystyle Z=\begin{pmatrix}\lambda_{i}\\\
\mu_{i}\end{pmatrix}_{i=1,\ldots,5}\,=\begin{pmatrix}1&0&\frac{1}{x_{1}}&\frac{1+x_{2}}{x_{1}x_{2}}&\frac{1+x_{3}(1+x_{2})}{x_{1}x_{2}x_{3}}\\\
0&1&1&1&1\\\ 0&0&0&\frac{x_{4}}{x_{2}}&1\\\
0&0&1&1&\frac{x_{4}-x_{5}}{x_{4}}\end{pmatrix}\,,$ (39)

The parameterisation used in this work has some benefits: the only
dimensionful quantity is $x_{1}$ and all holomorphic quantities are described
using only $x_{1},x_{2},x_{3}$. For real kinematics only $x_{2}$ and $x_{3}$
are complex. Notice that the conversion between the momentum twistor
coordinates and spinor-helicity expressions is only invertible for phase-free
quantities. For this purpose we may use the following relations,

$\displaystyle x_{1}=s_{12}\,,$ (40) $\displaystyle
x_{2}=-\frac{\text{tr}_{+}(p_{1},p_{2},p_{3},p_{4})}{s_{12}s_{34}}\,,$
$\displaystyle
x_{3}=-\frac{\text{tr}_{+}(p_{1},p_{3},p_{4},p_{5})}{s_{13}s_{45}}\,,$
$\displaystyle x_{4}=\frac{s_{23}}{s_{12}}\,,$ $\displaystyle
x_{5}=\frac{s_{45}}{s_{12}}\,,$

with
$\text{tr}_{+}(p_{i},p_{j},p_{k},p_{l})=\text{tr}[(1+\gamma_{5})\not{p}_{i}\not{p}_{j}\not{p}_{k}\not{p}_{l}]/2=[ij]\langle
jk\rangle[kl]\langle li\rangle$.

In our work we express the helicity amplitudes in terms of the momentum
twistors variables $x_{i}$. The phase information can be restored by
multiplying and dividing by a suitable phase factor,

$\displaystyle\mathcal{A}=\Phi(\lambda_{i},\tilde{\lambda}_{i})\,\left(\frac{\mathcal{A}(x_{i})}{\Phi(x_{i})}\right)\,,$
(41)

where $\mathcal{A}$ is an helicity amplitude – or in general some object with
a non-trivial phase – and $\Phi$ is an arbitrary factor with the same helicity
weights as $\mathcal{A}$. The quantities in the parentheses in Eq. (41) are
both written in terms of momentum twistors. Their ratio is phase-free and can
thus be expressed in terms of the scalar and pseudo-scalar invariants
$s_{ij}$’s and $\text{tr}_{5}$, e.g. through Eqs. (40). The factor outside the
parenthesis is instead written in terms of the spinor helicity variables and
carries all the phase information of $\mathcal{A}$.

## References

  * (1) D. A. Kosower and K. J. Larsen, _Maximal Unitarity at Two Loops_ , _Phys. Rev. D_ 85 (2012) 045017, [1108.1180].
  * (2) P. Mastrolia and G. Ossola, _On the Integrand-Reduction Method for Two-Loop Scattering Amplitudes_ , _JHEP_ 11 (2011) 014, [1107.6041].
  * (3) S. Badger, H. Frellesvig and Y. Zhang, _Hepta-Cuts of Two-Loop Scattering Amplitudes_ , _JHEP_ 04 (2012) 055, [1202.2019].
  * (4) Y. Zhang, _Integrand-Level Reduction of Loop Amplitudes by Computational Algebraic Geometry Methods_ , _JHEP_ 09 (2012) 042, [1205.5707].
  * (5) P. Mastrolia, E. Mirabella, G. Ossola and T. Peraro, _Scattering Amplitudes from Multivariate Polynomial Division_ , _Phys. Lett. B_ 718 (2012) 173–177, [1205.7087].
  * (6) P. Mastrolia, E. Mirabella, G. Ossola and T. Peraro, _Integrand-Reduction for Two-Loop Scattering Amplitudes through Multivariate Polynomial Division_ , _Phys. Rev. D_ 87 (2013) 085026, [1209.4319].
  * (7) H. Ita, _Two-loop Integrand Decomposition into Master Integrals and Surface Terms_ , _Phys. Rev. D_ 94 (2016) 116015, [1510.05626].
  * (8) S. Badger, H. Frellesvig and Y. Zhang, _A Two-Loop Five-Gluon Helicity Amplitude in QCD_ , _JHEP_ 12 (2013) 045, [1310.1051].
  * (9) S. Badger, G. Mogull, A. Ochirov and D. O’Connell, _A Complete Two-Loop, Five-Gluon Helicity Amplitude in Yang-Mills Theory_ , _JHEP_ 10 (2015) 064, [1507.08797].
  * (10) 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 (2017) 142001, [1703.05273].
  * (11) S. Abreu, J. Dormans, F. Febres Cordero, H. Ita, M. Kraus, B. Page et al., _Caravel: A C++ Framework for the Computation of Multi-Loop Amplitudes with Numerical Unitarity_ , 2009.11957.
  * (12) J. R. Andersen et al., _Les Houches 2015: Physics at TeV Colliders Standard Model Working Group Report_ , in _9th Les Houches Workshop on Physics at TeV Colliders_ , 5, 2016.  1605.04692.
  * (13) _Les Houches 2017: Physics at TeV Colliders Standard Model Working Group Report_ , 3, 2018.
  * (14) S. Amoroso et al., _Les Houches 2019: Physics at TeV Colliders: Standard Model Working Group Report_ , in _11th Les Houches Workshop on Physics at TeV Colliders: PhysTeV Les Houches_ , 3, 2020.  2003.01700.
  * (15) H. A. Chawdhry, M. Czakon, A. Mitov and R. Poncelet, _NNLO QCD corrections to diphoton production with an additional jet at the LHC_ , 2105.06940.
  * (16) D. Chicherin, J. Henn and V. Mitev, _Bootstrapping pentagon functions_ , _JHEP_ 05 (2018) 164, [1712.09610].
  * (17) C. G. Papadopoulos, D. Tommasini and C. Wever, _The Pentabox Master Integrals with the Simplified Differential Equations approach_ , _JHEP_ 04 (2016) 078, [1511.09404].
  * (18) T. Gehrmann, J. Henn and N. Lo Presti, _Pentagon functions for massless planar scattering amplitudes_ , _JHEP_ 10 (2018) 103, [1807.09812].
  * (19) D. Chicherin, T. Gehrmann, J. Henn, N. Lo Presti, V. Mitev and P. Wasser, _Analytic result for the nonplanar hexa-box integrals_ , _JHEP_ 03 (2019) 042, [1809.06240].
  * (20) D. Chicherin, T. Gehrmann, J. Henn, P. Wasser, Y. Zhang and S. Zoia, _All Master Integrals for Three-Jet Production at Next-to-Next-to-Leading Order_ , _Phys. Rev. Lett._ 123 (2019) 041603, [1812.11160].
  * (21) D. Chicherin and V. Sotnikov, _Pentagon Functions for Scattering of Five Massless Particles_ , _JHEP_ 12 (2020) 167, [2009.07803].
  * (22) F. V. Tkachov, _A Theorem on Analytical Calculability of Four Loop Renormalization Group Functions_ , _Phys. Lett._ 100B (1981) 65–68.
  * (23) K. G. Chetyrkin and F. V. Tkachov, _Integration by Parts: The Algorithm to Calculate beta Functions in 4 Loops_ , _Nucl. Phys._ B192 (1981) 159–204.
  * (24) S. Laporta, _High precision calculation of multiloop Feynman integrals by difference equations_ , _Int. J. Mod. Phys._ A15 (2000) 5087–5159, [hep-ph/0102033].
  * (25) C. Anastasiou and A. Lazopoulos, _Automatic integral reduction for higher order perturbative calculations_ , _JHEP_ 07 (2004) 046, [hep-ph/0404258].
  * (26) C. Studerus, _Reduze-Feynman Integral Reduction in C++_ , _Comput. Phys. Commun._ 181 (2010) 1293–1300, [0912.2546].
  * (27) A. von Manteuffel and C. Studerus, _Reduze 2 - Distributed Feynman Integral Reduction_ , 1201.4330.
  * (28) R. N. Lee, _Presenting LiteRed: a tool for the Loop InTEgrals REDuction_ , 1212.2685.
  * (29) A. V. Smirnov and F. S. Chuharev, _FIRE6: Feynman Integral REduction with Modular Arithmetic_ , _Comput. Phys. Commun._ 247 (2020) 106877, [1901.07808].
  * (30) J. Klappert, F. Lange, P. Maierhöfer and J. Usovitsch, _Integral Reduction with Kira 2.0 and Finite Field Methods_ , _Comput. Phys. Commun._ 266 (2021) 108024, [2008.06494].
  * (31) J. Gluza, K. Kajda and D. A. Kosower, _Towards a Basis for Planar Two-Loop Integrals_ , _Phys. Rev. D_ 83 (2011) 045012, [1009.0472].
  * (32) R. M. Schabinger, _A New Algorithm For The Generation Of Unitarity-Compatible Integration By Parts Relations_ , _JHEP_ 01 (2012) 077, [1111.4220].
  * (33) K. J. Larsen and Y. Zhang, _Integration-by-parts reductions from unitarity cuts and algebraic geometry_ , _Phys. Rev. D_ 93 (2016) 041701, [1511.01071].
  * (34) 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 (2018) 025023, [1712.09737].
  * (35) J. Böhm, A. Georgoudis, K. J. Larsen, H. Schönemann and Y. Zhang, _Complete integration-by-parts reductions of the non-planar hexagon-box via module intersections_ , _JHEP_ 09 (2018) 024, [1805.01873].
  * (36) J. Boehm, M. Wittmann, Z. Wu, Y. Xu and Y. Zhang, _IBP reduction coefficients made simple_ , _JHEP_ 12 (2020) 054, [2008.13194].
  * (37) P. Mastrolia and S. Mizera, _Feynman Integrals and Intersection Theory_ , _JHEP_ 02 (2019) 139, [1810.03818].
  * (38) H. Frellesvig, F. Gasparotto, M. K. Mandal, P. Mastrolia, L. Mattiazzi and S. Mizera, _Vector Space of Feynman Integrals and Multivariate Intersection Numbers_ , _Phys. Rev. Lett._ 123 (2019) 201602, [1907.02000].
  * (39) H. Frellesvig, F. Gasparotto, S. Laporta, M. K. Mandal, P. Mastrolia, L. Mattiazzi et al., _Decomposition of Feynman Integrals on the Maximal Cut by Intersection Numbers_ , _JHEP_ 05 (2019) 153, [1901.11510].
  * (40) H. Frellesvig, F. Gasparotto, S. Laporta, M. K. Mandal, P. Mastrolia, L. Mattiazzi et al., _Decomposition of Feynman Integrals by Multivariate Intersection Numbers_ , _JHEP_ 03 (2021) 027, [2008.04823].
  * (41) X. Liu, Y.-Q. Ma and C.-Y. Wang, _A Systematic and Efficient Method to Compute Multi-loop Master Integrals_ , _Phys. Lett. B_ 779 (2018) 353–357, [1711.09572].
  * (42) X. Liu and Y.-Q. Ma, _Determining arbitrary Feynman integrals by vacuum integrals_ , _Phys. Rev. D_ 99 (2019) 071501, [1801.10523].
  * (43) X. Guan, X. Liu and Y.-Q. Ma, _Complete reduction of integrals in two-loop five-light-parton scattering amplitudes_ , _Chin. Phys. C_ 44 (2020) 093106, [1912.09294].
  * (44) P. Zhang, C.-Y. Wang, X. Liu, Y.-Q. Ma, C. Meng and K.-T. Chao, _Semi-analytical calculation of gluon fragmentation into 1S${}_{0}^{[1,8]}$ quarkonia at next-to-leading order_, _JHEP_ 04 (2019) 116, [1810.07656].
  * (45) Y. Wang, Z. Li and N. Ul Basat, _Direct reduction of multiloop multiscale scattering amplitudes_ , _Phys. Rev. D_ 101 (2020) 076023, [1901.09390].
  * (46) D. A. Kosower, _Direct Solution of Integration-by-Parts Systems_ , _Phys. Rev. D_ 98 (2018) 025008, [1804.00131].
  * (47) A. von Manteuffel and R. M. Schabinger, _A novel approach to integration by parts reduction_ , _Phys. Lett._ B744 (2015) 101–104, [1406.4513].
  * (48) T. Peraro, _Scattering amplitudes over finite fields and multivariate functional reconstruction_ , _JHEP_ 12 (2016) 030, [1608.01902].
  * (49) J. Klappert and F. Lange, _Reconstructing rational functions with FireFly_ , _Comput. Phys. Commun._ 247 (2020) 106951, [1904.00009].
  * (50) T. Peraro, _FiniteFlow: multivariate functional reconstruction using finite fields and dataflow graphs_ , _JHEP_ 07 (2019) 031, [1905.08019].
  * (51) J. Klappert, S. Y. Klein and F. Lange, _Interpolation of dense and sparse rational functions and other improvements in FireFly_ , _Comput. Phys. Commun._ 264 (2021) 107968, [2004.01463].
  * (52) S. Badger, C. Brønnum-Hansen, H. B. Hartanto and T. Peraro, _First look at two-loop five-gluon scattering in QCD_ , _Phys. Rev. Lett._ 120 (2018) 092001, [1712.02229].
  * (53) 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 (2018) 116014, [1712.03946].
  * (54) S. Badger, C. Brønnum-Hansen, T. Gehrmann, H. B. Hartanto, J. Henn, N. A. Lo Presti et al., _Applications of integrand reduction to two-loop five-point scattering amplitudes in QCD_ , _PoS_ LL2018 (2018) 006, [1807.09709].
  * (55) S. Abreu, F. Febres Cordero, H. Ita, B. Page and V. Sotnikov, _Planar Two-Loop Five-Parton Amplitudes from Numerical Unitarity_ , _JHEP_ 11 (2018) 116, [1809.09067].
  * (56) T. Gehrmann, J. Henn and N. Lo Presti, _Analytic form of the two-loop planar five-gluon all-plus-helicity amplitude in QCD_ , _Phys. Rev. Lett._ 116 (2016) 062001, [1511.05409].
  * (57) S. Badger, C. Brønnum-Hansen, H. B. Hartanto and T. Peraro, _Analytic helicity amplitudes for two-loop five-gluon scattering: the single-minus case_ , _JHEP_ 01 (2019) 186, [1811.11699].
  * (58) S. Abreu, J. Dormans, F. Febres Cordero, H. Ita and B. Page, _Analytic Form of Planar Two-Loop Five-Gluon Scattering Amplitudes in QCD_ , _Phys. Rev. Lett._ 122 (2019) 082002, [1812.04586].
  * (59) S. Abreu, J. Dormans, F. Febres Cordero, H. Ita, B. Page and V. Sotnikov, _Analytic Form of the Planar Two-Loop Five-Parton Scattering Amplitudes in QCD_ , _JHEP_ 05 (2019) 084, [1904.00945].
  * (60) S. Abreu, F. F. Cordero, H. Ita, B. Page and V. Sotnikov, _Leading-Color Two-Loop QCD Corrections for Three-Jet Production at Hadron Colliders_ , 2102.13609.
  * (61) S. Abreu, B. Page, E. Pascual and V. Sotnikov, _Leading-Color Two-Loop QCD Corrections for Three-Photon Production at Hadron Colliders_ , 2010.15834.
  * (62) H. A. Chawdhry, M. Czakon, A. Mitov and R. Poncelet, _Two-loop leading-color helicity amplitudes for three-photon production at the LHC_ , 2012.13553.
  * (63) H. A. Chawdhry, M. L. Czakon, A. Mitov and R. Poncelet, _NNLO QCD corrections to three-photon production at the LHC_ , _JHEP_ 02 (2020) 057, [1911.00479].
  * (64) S. Kallweit, V. Sotnikov and M. Wiesemann, _Triphoton production at hadron colliders in NNLO QCD_ , 2010.04681.
  * (65) M. Czakon, A. Mitov and R. Poncelet, _Tour de force in Quantum Chromodynamics: A first next-to-next-to-leading order study of three-jet production at the LHC_ , 2106.05331.
  * (66) 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 (2012) 072001, [1110.2375].
  * (67) J. M. Campbell, R. K. Ellis, Y. Li and C. Williams, _Predictions for diphoton production at the LHC through NNLO in QCD_ , _JHEP_ 07 (2016) 148, [1603.02663].
  * (68) 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 (2002) 255–289, [hep-ph/0201274].
  * (69) Z. Bern, A. De Freitas and L. J. Dixon, _Two loop amplitudes for gluon fusion into two photons_ , _JHEP_ 09 (2001) 037, [hep-ph/0109078].
  * (70) B. Agarwal, F. Buccioni, A. von Manteuffel and L. Tancredi, _Two-loop leading colour QCD corrections to $q\bar{q}\to\gamma\gamma g$ and $qg\to\gamma\gamma q$_, _JHEP_ 04 (2021) 201, [2102.01820].
  * (71) H. A. Chawdhry, M. Czakon, A. Mitov and R. Poncelet, _Two-loop leading-colour QCD helicity amplitudes for two-photon plus jet production at the LHC_ , 2103.04319.
  * (72) B. Agarwal, F. Buccioni, A. von Manteuffel and L. Tancredi, _Two-loop helicity amplitudes for diphoton plus jet production in full color_ , 2105.04585.
  * (73) S. Badger, B. Biedermann, P. Uwer and V. Yundin, _Numerical evaluation of virtual corrections to multi-jet production in massless QCD_ , _Comput. Phys. Commun._ 184 (2013) 1981–1998, [1209.0100].
  * (74) D. A. Dicus and S. S. D. Willenbrock, _Photon Pair Production and the Intermediate Mass Higgs Boson_ , _Phys. Rev. D_ 37 (1988) 1801.
  * (75) D. de Florian and Z. Kunszt, _Two photons plus jet at LHC: The NNLO contribution from the g g initiated process_ , _Phys. Lett. B_ 460 (1999) 184–188, [hep-ph/9905283].
  * (76) S. Catani, _The Singular behavior of QCD amplitudes at two loop order_ , _Phys. Lett._ B427 (1998) 161–171, [hep-ph/9802439].
  * (77) T. Becher and M. Neubert, _On the Structure of Infrared Singularities of Gauge-Theory Amplitudes_ , _JHEP_ 06 (2009) 081, [0903.1126].
  * (78) T. Becher and M. Neubert, _Infrared singularities of scattering amplitudes in perturbative QCD_ , _Phys. Rev. Lett._ 102 (2009) 162001, [0901.0722].
  * (79) E. Gardi and L. Magnea, _Factorization constraints for soft anomalous dimensions in QCD scattering amplitudes_ , _JHEP_ 03 (2009) 079, [0901.1091].
  * (80) E. Gardi and L. Magnea, _Infrared singularities in QCD amplitudes_ , _Frascati Phys. Ser._ 50 (2010) 137–157, [0908.3273].
  * (81) H. B. Hartanto, S. Badger, C. Brønnum-Hansen and T. Peraro, _A numerical evaluation of planar two-loop helicity amplitudes for a W-boson plus four partons_ , _JHEP_ 09 (2019) 119, [1906.11862].
  * (82) S. Badger, E. Chaubey, H. B. Hartanto and R. Marzucca, _Two-loop leading colour QCD helicity amplitudes for top quark pair production in the gluon fusion channel_ , 2102.13450.
  * (83) P. Nogueira, _Automatic Feynman graph generation_ , _J. Comput. Phys._ 105 (1993) 279–289.
  * (84) J. Kuipers, T. Ueda, J. A. M. Vermaseren and J. Vollinga, _FORM version 4.0_ , _Comput. Phys. Commun._ 184 (2013) 1453–1467, [1203.6543].
  * (85) B. Ruijl, T. Ueda and J. Vermaseren, _FORM version 4.2_ , 1707.06453.
  * (86) G. Cullen, M. Koch-Janusz and T. Reiter, _Spinney: A Form Library for Helicity Spinors_ , _Comput. Phys. Commun._ 182 (2011) 2368–2387, [1008.0803].
  * (87) A. Hodges, _Eliminating spurious poles from gauge-theoretic amplitudes_ , _JHEP_ 05 (2013) 135, [0905.1473].
  * (88) R. Eden, P. Landshoff, D. Olive and J. Polkinghorne, _The Analytic S-Matrix_.  Cambridge University Press, 2002.
  * (89) R. Britto, F. Cachazo and B. Feng, _Generalized unitarity and one-loop amplitudes in N=4 super-Yang-Mills_ , _Nucl. Phys. B_ 725 (2005) 275–305, [hep-th/0412103].
  * (90) F. Cachazo, _Sharpening The Leading Singularity_ , 0803.1988.
  * (91) N. Arkani-Hamed, J. L. Bourjaily, F. Cachazo and J. Trnka, _Local Integrals for Planar Scattering Amplitudes_ , _JHEP_ 06 (2012) 125, [1012.6032].
  * (92) D. C. Dunbar and W. B. Perkins, _Two-loop five-point all plus helicity Yang-Mills amplitude_ , _Phys. Rev. D_ 93 (2016) 085029, [1603.07514].
  * (93) S. Badger, D. Chicherin, T. Gehrmann, G. Heinrich, J. Henn, T. Peraro et al., _Analytic form of the full two-loop five-gluon all-plus helicity amplitude_ , _Phys. Rev. Lett._ 123 (2019) 071601, [1905.03733].
  * (94) D. C. Dunbar, J. H. Godwin, W. B. Perkins and J. M. W. Strong, _Color Dressed Unitarity and Recursion for Yang-Mills Two-Loop All-Plus Amplitudes_ , _Phys. Rev. D_ 101 (2020) 016009, [1911.06547].
  * (95) M. Heller and A. von Manteuffel, _MultivariateApart: Generalized Partial Fractions_ , 2101.08283.
  * (96) L. E. K., _Factorization of rational functions of several variables into partial fractions_ , _Izvestiya Vysshikh Uchebnykh Zavedenii. Matematika_ 47 (1978) .
  * (97) A. Raichev, _Leinartas’s partial fraction decomposition_ , _arXiv e-prints_ (June, 2012) arXiv:1206.4740, [1206.4740].
  * (98) Z. Bern, L. J. Dixon and D. A. Kosower, _New QCD results from string theory_ , in _International Conference on Strings 93_ , 5, 1993.  hep-th/9311026.
  * (99) G. Mahlon, _Multi - gluon helicity amplitudes involving a quark loop_ , _Phys. Rev. D_ 49 (1994) 4438–4453, [hep-ph/9312276].
  * (100) Z. Bern, G. Chalmers, L. J. Dixon and D. A. Kosower, _One loop N gluon amplitudes with maximal helicity violation via collinear limits_ , _Phys. Rev. Lett._ 72 (1994) 2134–2137, [hep-ph/9312333].
  * (101) J. Henn, B. Power and S. Zoia, _Conformal Invariance of the One-Loop All-Plus Helicity Scattering Amplitudes_ , _JHEP_ 02 (2020) 019, [1911.12142].
  * (102) Z. Bern, L. J. Dixon, D. C. Dunbar and D. A. Kosower, _One loop n point gauge theory amplitudes, unitarity and collinear limits_ , _Nucl. Phys. B_ 425 (1994) 217–260, [hep-ph/9403226].
  * (103) Z. Bern, L. J. Dixon, D. C. Dunbar and D. A. Kosower, _Fusing gauge theory tree amplitudes into loop amplitudes_ , _Nucl. Phys. B_ 435 (1995) 59–101, [hep-ph/9409265].
  * (104) D. C. Dunbar, G. R. Jehu and W. B. Perkins, _The two-loop n-point all-plus helicity amplitude_ , _Phys. Rev. D_ 93 (2016) 125006, [1604.06631].
  * (105) D. C. Dunbar, G. R. Jehu and W. B. Perkins, _Two-loop six gluon all plus helicity amplitude_ , _Phys. Rev. Lett._ 117 (2016) 061602, [1605.06351].
  * (106) S. Badger, G. Mogull and T. Peraro, _Local integrands for two-loop all-plus Yang-Mills amplitudes_ , _JHEP_ 08 (2016) 063, [1606.02244].
  * (107) D. C. Dunbar, J. H. Godwin, G. R. Jehu and W. B. Perkins, _Analytic all-plus-helicity gluon amplitudes in QCD_ , _Phys. Rev. D_ 96 (2017) 116013, [1710.10071].
  * (108) D. C. Dunbar, W. B. Perkins and J. M. W. Strong, _$n$ -point QCD two-loop amplitude_, _Phys. Rev. D_ 101 (2020) 076001, [2001.11347].
  * (109) A. R. Dalgleish, D. C. Dunbar, W. B. Perkins and J. M. W. Strong, _Full color two-loop six-gluon all-plus helicity amplitude_ , _Phys. Rev. D_ 101 (2020) 076024, [2003.00897].
  * (110) G. Guennebaud, B. Jacob et al., “Eigen v3.” http://eigen.tuxfamily.org, 2010.
  * (111) Y. Hida, X. S. Li and D. H. Bailey, “libqd: quad-double / double-double computation package.” https://www.davidhbailey.com/dhbsoftware/, 2010.
  * (112) C. G. Papadopoulos and C. Wever, _Internal Reduction method for computing Feynman Integrals_ , _JHEP_ 02 (2020) 112, [1910.06275].
  * (113) S. Abreu, H. Ita, F. Moriello, B. Page, W. Tschernow and M. Zeng, _Two-Loop Integrals for Planar Five-Point One-Mass Processes_ , _JHEP_ 11 (2020) 117, [2005.04195].
  * (114) D. D. Canko, C. G. Papadopoulos and N. Syrrakos, _Analytic representation of all planar two-loop five-point Master Integrals with one off-shell leg_ , _JHEP_ 01 (2021) 199, [2009.13917].
  * (115) N. Syrrakos, _Pentagon integrals to arbitrary order in the dimensional regulator_ , 2012.10635.
  * (116) S. Badger, H. B. Hartanto and S. Zoia, _Two-loop QCD corrections to $Wb\bar{b}$ production at hadron colliders_, 2102.02516.