aainstitutetext: Theoretical Physics Department, CERN, 1211 Geneva 23,
Switzerlandbbinstitutetext: Mani L. Bhaumik Institute for Theoretical Physics,
Department of Physics and Astronomy, UCLA, Los Angeles, CA 90095,
USAccinstitutetext: Physics Department, Florida State University, 77 Chieftan
Way, Tallahassee, FL 32306, U.S.A.ddinstitutetext: Physikalisches Institut,
Albert-Ludwigs-Universität Freiburg, Hermann-Herder-Str. 3, D-79104 Freiburg,
Germanyeeinstitutetext: Laboratoire de physique de l’Ecole normale supérieure,
ENS, Université PSL, CNRS, Sorbonne Université, Université Paris-Diderot,
Sorbonne Paris Cité, 24 rue Lhomond, 75005 Paris, Franceffinstitutetext: Max
Planck Institute for Physics (Werner Heisenberg Institute), D–80805 Munich,
Germany

#  Leading-Color Two-Loop QCD Corrections for Three-Jet Production at Hadron
Colliders

S. Abreu c F. Febres Cordero d H. Ita e B. Page f and V. Sotnikov
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We present the complete set of leading-color two-loop contributions required
to obtain next-to-next-to-leading-order (NNLO) QCD corrections to three-jet
production at hadron colliders. We obtain analytic expressions for a
generating set of finite remainders, valid in the physical region for three-
jet production. The analytic continuation of the known Euclidean-region
results is determined from a small set of numerical evaluations of the
amplitudes. We obtain analytic expressions that are suitable for
phenomenological applications and we present a C++ library for their efficient
and stable numerical evaluation.

††preprint:  CERN-TH-2021-023 FR-PHENO-2021-06 MPP-2021-20

##  1 Introduction

Multi-jet events are ubiquitous at high-energy hadron colliders such as the
Large Hadron Collider (LHC). Their precise theoretical description including
higher-order quantum corrections offers a robust path to test the Standard
Model of particle physics. This has been demonstrated, for example, by
comparing measurements of the strong coupling $\alpha_{s}$ extracted from
ratios between three- and two-jet production cross sections with theoretical
predictions computed at next-to-leading order (NLO) in QCD Abazov:2012lua ;
Chatrchyan:2013txa ; CMS:2014mna ; ATLAS:2015yaa . Given the amount of
experimental data that has been collected since these studies were conducted
and the amount of data forecast to be collected in the near future at the LHC,
there is great interest in obtaining theoretical predictions that will allow
these types of studies to be carried out at a new level of precision. In this
paper, we provide a key missing ingredient required to reach this goal, namely
the leading-color two-loop contributions to the three-jet production cross
section at hadron colliders at next-to-next-to-leading order (NNLO) in QCD.
These are the so-called double-virtual $2\rightarrow 3$ contributions, which
together with the real-virtual $2\rightarrow 4$ and the double-real
$2\rightarrow 5$ contributions complete the $\mathcal{O}(\alpha_{s}^{2})$
corrections to the production of three jets at hadron colliders.

The calculation of NNLO QCD corrections to two-jet production at hadron
colliders has been a major success of the field in recent years. Results were
obtained first at leading color Ridder:2013mf ; Currie:2016bfm ;
Currie:2017eqf  and later subleading-color terms were incorporated
Czakon:2019tmo  (see also Currie:2013dwa ). Three-jet production is a more
complex process. Results are currently known at NLO in QCD Nagy:2001fj ;
Nagy:2003tz  (see also Kilgore:1996sq ; Kilgore:2000dr ) and also including
electroweak corrections Frederix:2018nkq ; Reyer:2019obz . Work on the NNLO
QCD corrections to three-jet production has so far been reported only for the
leading-color two-loop scattering amplitudes, first through numerical
evaluations Badger:2013gxa ; Badger:2017jhb ; Abreu:2017hqn ; Badger:2018gip ;
Abreu:2018jgq  and then in analytic form Gehrmann:2015bfy ; Badger:2018enw ;
Abreu:2018zmy ; Abreu:2019odu ; DeLaurentis:2020qle . Progress in including
subleading-color contributions has been made for related supersymmetric
theories Abreu:2018aqd ; Chicherin:2018yne ; Caron-Huot:2020vlo  as well as
for special external gluon states Badger:2019djh . In this work we present the
leading-color double-virtual contributions to the scattering cross section,
starting from the known analytic amplitudes Abreu:2019odu  which we
analytically continue from the (non-physical) Euclidean region to the physical
region of phase space. Furthermore, we implement our expressions in a publicly
available C++ library that allows one to evaluate these contributions
efficiently and in a stable manner, making our results readily available for a
future full calculation of the corresponding NNLO QCD corrections.111The
subleading-color contributions that we neglect are expected to have a minor
impact in phenomenological applications Berger:2009ep ; Ita:2011ar ;
Currie:2013dwa ; Czakon:2019tmo .

Our work builds on a number of recent technical advances that have paved the
way to the calculation of two-loop scattering amplitudes for $2\rightarrow 3$
processes with phenomenological applications Abreu:2020cwb ; Chawdhry:2020for
; Chawdhry:2019bji ; Kallweit:2020gcp ; Agarwal:2021grm ; Badger:2021nhg . We
employ the known pure two-loop master integrals Papadopoulos:2015jft ;
Gehrmann:2015bfy ; Abreu:2018aqd ; Chicherin:2018old , which can be written in
terms of a special set of multi-valued transcendental functions—the so-called
pentagon functions Gehrmann:2018yef ; Chicherin:2020oor . In ref.
Chicherin:2020oor , the five-point integrals are expressed in terms of a basis
of functions for all assignments of the incoming momenta required for the
double-virtual corrections we compute, and they can be efficiently evaluated
in the physical region corresponding to three-jet production.

We obtain the analytic form of a generating set of partial amplitudes in the
physical region by analytically continuing the known amplitudes from the
Euclidean region Abreu:2019odu . Based on considerations on the analytic
continuation of the pentagon functions, the continuation of each helicity
amplitude is reduced to a linear-algebra problem that can be efficiently
solved using evaluations of the amplitudes in a finite field
vonManteuffel:2014ixa ; Peraro:2016wsq ; Abreu:2017hqn . These evaluations are
performed with the program Caravel Abreu:2020xvt . This is a new generic
procedure that allows to upgrade results for amplitudes in the Euclidean
region to results valid in other regions. Finally, we assemble all partial
amplitudes into squared color- and helicity-summed finite remainders, which
can be numerically evaluated over phase space with the provided C++ library
FivePointAmplitudes . We perform several checks on our results and verify the
efficiency and numerical stability of our implementation.

The paper is organized as follows. In section 2 we establish our conventions
and define the objects we will be computing. In section 3 we describe our
procedure to compute the partial amplitudes in the physical region and the
different permutations that are required to obtain squared finite remainders.
In section 4 we present our results in the form of ancillary files and
showcase the efficiency and numerical stability of our code for the numerical
evaluation of the double-virtual NNLO contributions to the production of three
jets at hadron colliders, and discuss the checks that we have performed on our
results. We summarize our results in section 5.

##  2 Notation and Conventions

###  2.1 Helicity Amplitudes

We consider all the channels contributing to the production of three partons
at hadron colliders. We label the initial-state particles with indices 1 and
2, and the final-state particles with indices 3, 4 and 5. There is a single
channel with five gluons,

$\mathrm{g}^{h_{1}}_{p_{1}}\mathrm{g}^{h_{2}}_{p_{2}}\to\mathrm{g}^{h_{3}}_{p_{3}}\mathrm{g}^{h_{4}}_{p_{4}}\mathrm{g}^{h_{5}}_{p_{5}}\,,$
(1)

three channels with a pair of quark and anti-quark,

$\bar{q}^{h_{1}}_{p_{1}}q^{h_{2}}_{p_{2}}\to\mathrm{g}^{h_{3}}_{p_{3}}\mathrm{g}^{h_{4}}_{p_{4}}\mathrm{g}^{h_{5}}_{p_{5}}\,,\quad\bar{q}^{h_{1}}_{p_{1}}\mathrm{g}^{h_{2}}_{p_{2}}\to\bar{q}^{h_{3}}_{p_{3}}\mathrm{g}^{h_{4}}_{p_{4}}\mathrm{g}^{h_{5}}_{p_{5}}\,,\quad\mathrm{g}^{h_{1}}_{p_{1}}\mathrm{g}^{h_{2}}_{p_{2}}\to
q^{h_{3}}_{p_{3}}\bar{q}^{h_{4}}_{p_{4}}\mathrm{g}^{h_{5}}_{p_{5}}\,,$ (2)

and four channels with two distinct pairs of quark and anti-quark

$\displaystyle\begin{split}&\bar{q}^{h_{1}}_{p_{1}}q^{h_{2}}_{p_{2}}\to
Q^{h_{3}}_{p_{3}}\bar{Q}^{h_{4}}_{p_{4}}\mathrm{g}^{h_{5}}_{p_{5}}\,,\quad\bar{q}^{h_{1}}_{p_{1}}Q^{h_{2}}_{p_{2}}\to
Q^{h_{3}}_{p_{3}}\bar{q}^{h_{4}}_{p_{4}}\mathrm{g}^{h_{5}}_{p_{5}}\,,\\\
&\bar{q}^{h_{1}}_{p_{1}}\bar{Q}^{h_{2}}_{p_{2}}\to\bar{q}^{h_{3}}_{p_{3}}\bar{Q}^{h_{4}}_{p_{4}}\mathrm{g}^{h_{5}}_{p_{5}}\,,\quad\bar{q}^{h_{1}}_{p_{1}}\mathrm{g}^{h_{2}}_{p_{2}}\to\bar{q}^{h_{3}}_{p_{3}}Q^{h_{4}}_{p_{4}}\bar{Q}^{h_{5}}_{p_{5}}.\end{split}$
(3)

In eqs. 1, 2 and LABEL:eq:qQProc, we identify each particle by its type
($\mathrm{g}$ for gluon, $q$ or $Q$ for distinct quarks, and $\bar{q}$ and
$\bar{Q}$ for the corresponding anti-quarks), a superscript $h_{i}$ for the
helicity state and a subscript $p_{i}$ for the four momentum. All other
channels involving quarks are related to the ones above by charge conjugation
and permutation of incoming or outgoing labels. While in LABEL:eq:qQProc, we
assumed the two quark pairs to be of distinct flavours, we should also
consider the case with two quark lines of identical flavour. There are three
such channels which, schematically, can be obtained from

$\displaystyle\begin{split}\bar{q}^{h_{1}}_{p_{1}}q^{h_{2}}_{p_{2}}\to
q^{h_{3}}_{p_{3}}\bar{q}^{h_{4}}_{p_{4}}\mathrm{g}^{h_{5}}_{p_{5}}=\left(\bar{q}^{h_{1}}_{p_{1}}q^{h_{2}}_{p_{2}}\to
Q^{h_{3}}_{p_{3}}\bar{Q}^{h_{4}}_{p_{4}}\mathrm{g}^{h_{5}}_{p_{5}}\right)-\left(\bar{q}^{h_{1}}_{p_{1}}Q^{h_{2}}_{p_{2}}\to
Q^{h_{3}}_{p_{3}}\bar{q}^{h_{4}}_{p_{4}}\mathrm{g}^{h_{5}}_{p_{5}}\right)\,,\\\
\bar{q}^{h_{1}}_{p_{1}}\bar{q}^{h_{2}}_{p_{2}}\to\bar{q}^{h_{3}}_{p_{3}}\bar{q}^{h_{4}}_{p_{4}}\mathrm{g}^{h_{5}}_{p_{5}}=\left(\bar{q}^{h_{1}}_{p_{1}}\bar{Q}^{h_{2}}_{p_{2}}\to\bar{q}^{h_{3}}_{p_{3}}\bar{Q}^{h_{4}}_{p_{4}}\mathrm{g}^{h_{5}}_{p_{5}}\right)-\left(\bar{q}^{h_{1}}_{p_{1}}\bar{Q}^{h_{2}}_{p_{2}}\to\bar{Q}^{h_{3}}_{p_{3}}\bar{q}^{h_{4}}_{p_{4}}\mathrm{g}^{h_{5}}_{p_{5}}\right)\,,\\\
\bar{q}^{h_{1}}_{p_{1}}\mathrm{g}^{h_{2}}_{p_{2}}\to\bar{q}^{h_{3}}_{p_{3}}q^{h_{4}}_{p_{4}}\bar{q}^{h_{5}}_{p_{5}}=\left(\bar{q}^{h_{1}}_{p_{1}}\mathrm{g}^{h_{2}}_{p_{2}}\to\bar{q}^{h_{3}}_{p_{3}}Q^{h_{4}}_{p_{4}}\bar{Q}^{h_{5}}_{p_{5}}\right)-\left(\bar{q}^{h_{1}}_{p_{1}}\mathrm{g}^{h_{2}}_{p_{2}}\to\bar{Q}^{h_{3}}_{p_{3}}Q^{h_{4}}_{p_{4}}\bar{q}^{h_{5}}_{p_{5}}\right)\,.\end{split}$
(4)

In order to obtain them we thus also compute the (redundant) channels

$\bar{q}^{h_{1}}_{p_{1}}\bar{Q}^{h_{2}}_{p_{2}}\to\bar{Q}^{h_{3}}_{p_{3}}\bar{q}^{h_{4}}_{p_{4}}\mathrm{g}^{h_{5}}_{p_{5}}\,,\quad\bar{q}^{h_{1}}_{p_{1}}\mathrm{g}^{h_{2}}_{p_{2}}\to\bar{Q}^{h_{3}}_{p_{3}}Q^{h_{4}}_{p_{4}}\bar{q}^{h_{5}}_{p_{5}}.$
(5)

We remark that, when labelling channels, such as in eqs. 5, 4,
LABEL:eq:qQProc, 2 and 1, the momenta $p_{1}$ and $p_{2}$ are taken to be
incoming and $p_{3},p_{4}$ and $p_{5}$ are taken to be outgoing. In all other
contexts we use the symmetric “all-outgoing” convention of refs. Abreu:2018jgq
; Abreu:2019odu , such that momentum conservation is
$p_{1}+p_{2}+p_{3}+p_{4}+p_{5}=0$.

We will consider each of the channels above in the leading-color
approximation. Let us denote the set of helicities by a vector ${\bf\it h}$
and the set of momenta by a vector ${\bf\it p}$. The amplitudes
$\mathcal{M}({\bf\it h},{\bf\it p})$ describing each of the channels above,
which are vectors in color space, can be expanded in powers of the bare strong
coupling $\alpha_{s}^{0}=(g_{s}^{0})^{2}/(4\pi)$,

$\mathcal{M}=\left(g_{s}^{0}\right)^{3}\left(\mathcal{M}^{(0)}+\frac{\alpha^{0}_{s}}{2\pi}\mathcal{M}^{(1)}+\left(\frac{\alpha^{0}_{s}}{2\pi}\right)^{2}\mathcal{M}^{(2)}+\mathcal{O}\left((\alpha^{0}_{s})^{3}\right)\right)\,,$
(6)

where we drop the dependence on ${\bf\it h}$ and ${\bf\it p}$ for compactness.
The order of $\alpha^{0}_{s}$ is aligned with the number of loops of the
Feynman diagrams contributing to the associated coefficient, so we refer to
$\mathcal{M}^{(1)}$ and $\mathcal{M}^{(2)}$ as the one-loop and two-loop
amplitudes respectively, and to $\mathcal{M}^{(0)}$ as the tree-level
amplitude. All calculations will be done in the ’t Hooft-Veltman (HV) scheme
of dimensional regularization with $D=4-2\epsilon$ space-time dimensions and
we mostly keep the dependence on $\epsilon$ implicit. For our definition of
helicity amplitudes with external quarks, we refer the reader to the detailed
discussion in ref. Abreu:2018jgq .

For each channel, the coefficients $\mathcal{M}^{(i)}$ can be decomposed in
terms of different color structures, and at leading color this decomposition
is independent of the loop order. We write

$\mathcal{M}^{(i)}({\bf\it h},{\bf\it
p})=\left(\frac{N_{c}}{2}\right)^{i}\sum_{\sigma\in\Sigma}\mathcal{C}\left(\sigma\right)\Phi\left(\sigma({\bf\it
h}),\sigma({\bf\it p})\right)\mathcal{A}^{(i)}\left(\sigma({\bf\it
h}),\sigma({\bf\it p})\right)\,,$ (7)

where $N_{c}$ is the number of colors, $\mathcal{C}\left(\sigma\right)$
denotes a vector in color space and the sum is over a channel-dependent set of
permutations $\Sigma$. We defer a more detailed discussion of the sets
$\Sigma$ and the vectors $\mathcal{C}\left(\sigma\right)$ to section 3.2. For
now, we simply define the action of a permutation
$\sigma=\\{\sigma_{1},\sigma_{2},\sigma_{3},\sigma_{4},\sigma_{5}\\}\in S_{5}$
on momenta and helicity labels as $\sigma(p_{i})\coloneqq p_{\sigma_{i}}$ and
$\sigma(h_{i})\coloneqq h_{\sigma_{i}}$. The $\Phi\left(\sigma({\bf\it
h}),\sigma({\bf\it p})\right)$ are normalization factors, and the
$\mathcal{A}^{(i)}\left(\sigma({\bf\it h}),\sigma({\bf\it p})\right)$ are
partial amplitudes.

For each term in the sum, the coefficient of each $\mathcal{C}(\sigma)$ has a
fixed ordering of the external momenta determined by the permutation $\sigma$.
The $\Phi\left(\sigma({\bf\it h}),\sigma({\bf\it p})\right)$ are defined so
that $\mathcal{A}^{(0)}\left(\sigma({\bf\it h}),\sigma({\bf\it
p})\right)\equiv 1$ whenever the tree-level amplitude is non-vanishing, and
for the other cases we choose222Here we use the same notation for spinor
brackets as in Abreu:2020xvt .

$\displaystyle\Phi({\bf\it h}=\\{+1,+1,+1,+1,+1\\},\,{\bf\it p})$
$\displaystyle\coloneqq\frac{[21]^{2}\langle 12\rangle^{3}\langle
13\rangle}{\langle 14\rangle^{2}\langle 15\rangle^{2}\langle 23\rangle^{3}},$
(8a) $\displaystyle\Phi({\bf\it h}=\\{-1,+1,+1,+1,+1\\},\,{\bf\it p})$
$\displaystyle\coloneqq\frac{[21]\langle 12\rangle^{4}\langle
13\rangle^{3}}{\langle 14\rangle^{2}\langle 15\rangle^{2}\langle
23\rangle^{5}},$ (8b) $\displaystyle\Phi({\bf\it
h}=\\{+\frac{1}{2},-\frac{1}{2},+1,+1,+1\\},\,{\bf\it p})$
$\displaystyle\coloneqq\frac{[31]\langle 12\rangle^{3}\langle
13\rangle}{\langle 14\rangle^{2}\langle 15\rangle^{2}\langle 23\rangle^{2}},$
(8c)

where we set ${\bf\it p}=\\{p_{1},p_{2},p_{3},p_{4},p_{5}\\}$. With this
choice, the partial amplitudes $\mathcal{A}\left(\sigma({\bf\it
h}),\sigma({\bf\it p})\right)$ only depend on the kinematics through the
Mandelstam invariants

$\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}\,,\quad\\\
&s_{45}=(p_{4}+p_{5})^{2}\,,\quad s_{15}=(p_{1}+p_{5})^{2}\,,\quad\end{split}$
(9)

which we collectively denote as
$\vec{s}=\\{s_{12},s_{23},s_{34},s_{45},s_{15}\\}$, and the parity-odd
contraction of four momenta

${\rm tr}_{5}\coloneqq 4\,\mathrm{i}\,\varepsilon(p_{1},p_{2},p_{3},p_{4}),$
(10)

where
$\varepsilon(\boldsymbol{\cdot},\boldsymbol{\cdot},\boldsymbol{\cdot},\boldsymbol{\cdot})$
is the fully anti-symmetric Levi-Civita tensor, whose sign is fixed such that
$\imaginary({\rm tr}_{5})>0$ corresponds to the choice of the basis helicity
amplitudes in ref. Abreu:2019odu . The action of $\sigma$ on quantities that
depend on momenta is inherited from the action on momentum labels. For
instance, the action on the Mandelstam invariants $s_{ij}=(p_{i}+p_{j})^{2}$
is given by $\sigma(s_{ij})=s_{\sigma_{i}\sigma_{j}}$, and on the parity-odd
invariant in eq. 10 by $\sigma({\rm tr}_{5})=\mathrm{sgn}(\sigma)\,{\rm
tr}_{5}$.

We work in the leading-color approximation where the number of colors $N_{c}$
is large and the ratio ${N_{f}}/{N_{c}}$ is kept fixed, where ${N_{f}}$
denotes the number of massless quark flavors. In this limit, each of the
$\mathcal{A}^{(i)}\left({\bf\it h},{\bf\it p}\right)$ can be further
decomposed in powers of ${N_{f}}$. We write

$\mathcal{A}^{(i)}\left({\bf\it h},{\bf\it
p}\right)=\sum_{j=0}^{i}\left(\frac{N_{f}}{N_{c}}\right)^{j}\mathcal{A}^{(i)[j]}\left({\bf\it
h},{\bf\it p}\right)\,.$ (11)

In fig. 1 we depict characteristic diagrams that contribute to
$\mathcal{A}^{(2)[j]}$ for representative channels with no external quark
lines, a single external quark line or two external quark lines.

$\displaystyle\mathcal{A}^{(2)}_{\mathrm{g}\mathrm{g}\rightarrow\mathrm{g}\mathrm{g}\mathrm{g}}\sim\vbox{\hbox{
\includegraphics[scale={0.38}]{figures/5g}}}+\vbox{\hbox{
\includegraphics[scale={0.38}]{figures/5gnf}}}+\vbox{\hbox{
\includegraphics[scale={0.38}]{figures/5gnf2}}}+\ldots,$
$\displaystyle\mathcal{A}^{(2)}_{\bar{q}q\rightarrow\mathrm{g}\mathrm{g}\mathrm{g}}\sim\vbox{\hbox{
\includegraphics[scale={0.38}]{figures/2q3g}}}+\vbox{\hbox{
\includegraphics[scale={0.38}]{figures/2q3gnf}}}+\vbox{\hbox{
\includegraphics[scale={0.38}]{figures/2q3gnf2}}}+\ldots,$
$\displaystyle\mathcal{A}^{(2)}_{\bar{q}q\rightarrow
Q\bar{Q}\mathrm{g}}\sim\vbox{\hbox{
\includegraphics[scale={0.38}]{figures/4q1g}}}+\vbox{\hbox{
\raisebox{14.22636pt}{\includegraphics[scale={0.38}]{figures/4q1gnf}}}}+\vbox{\hbox{
\raisebox{11.38109pt}{\includegraphics[scale={0.38}]{figures/4q1gnf2}}}}+\ldots.$
Figure 1: Characteristic Feynman diagrams which contribute to
$\mathcal{A}^{(2)[j]}$ for representative channels with no external quark
lines, a single quark line and two external quark lines, for $j=0$ (left
column), $j=1$ (middle column) and $j=2$ (right column).

The renormalized amplitudes $\mathcal{M}_{R}$ are obtained from the bare
amplitudes $\mathcal{M}$ by considering an expansion in powers of the
renormalized coupling $\alpha_{s}(\mu)$. In the $\overline{\text{MS}}$-scheme,
the later is related to the bare coupling $\alpha_{s}^{0}$ through

$\alpha_{0}\mu_{0}^{2\epsilon}S_{\epsilon}=\alpha_{s}(\mu)\mu^{2\epsilon}\left(1-\frac{\beta_{0}}{\epsilon}\frac{\alpha_{s}(\mu)}{2\pi}+\left(\frac{\beta_{0}^{2}}{\epsilon^{2}}-\frac{\beta_{1}}{2\epsilon}\right)\left(\frac{\alpha_{s}(\mu)}{2\pi}\right)^{2}+\mathcal{O}\left(\alpha_{s}^{3}(\mu)\right)\right),$
(12)

where $S_{\epsilon}=(4\pi)^{\epsilon}e^{-\epsilon\gamma_{E}}\,$, and $\mu_{0}$
and $\mu$ are the dimensional regularization and renormalization scales
respectively, which from now on we assume to be equal. We will also suppress
the $\mu$ dependence in the coupling and write
$\alpha_{s}\equiv\alpha_{s}(\mu)$. The $\beta_{i}$ are the coefficients of the
QCD $\beta$-function,

$\displaystyle\beta_{0}$
$\displaystyle=\frac{11C_{A}-4T_{F}{N_{f}}}{6}=\frac{11N_{c}-2{N_{f}}}{6}\,,\quad$
(13) $\displaystyle\beta_{1}$
$\displaystyle=\frac{17C_{A}^{2}-6C_{F}T_{F}{N_{f}}-10C_{A}T_{F}{N_{f}}}{6}=\frac{1}{6}\left(17N_{c}^{2}-\frac{13}{2}N_{c}{N_{f}}\right)+\mathcal{O}(N_{c}^{-1})\,,$

where we have set $T_{F}=1/2$, $C_{A}=N_{c}$ and
$C_{F}=(N_{c}^{2}-1)/(2N_{c})$, and only kept terms that contribute at leading
color. The renormalized partial amplitudes are related to their bare
counterparts as

$\displaystyle\begin{split}&\mathcal{A}_{R}^{(0)}=\mathcal{A}^{(0)},\\\
&\mathcal{A}_{R}^{(1)}=S_{\epsilon}^{-1}\mathcal{A}^{(1)}-\frac{3\beta_{0}}{\epsilon
N_{c}}\mathcal{A}^{(0)}\,,\\\
&\mathcal{A}_{R}^{(2)}=S_{\epsilon}^{-2}\mathcal{A}^{(2)}-\frac{5\beta_{0}}{\epsilon
N_{c}}S_{\epsilon}^{-1}\mathcal{A}^{(1)}+\left(\frac{15\beta_{0}^{2}}{2\epsilon^{2}N_{c}^{2}}-\frac{3\beta_{1}}{\epsilon
N_{c}^{2}}\right)\mathcal{A}^{(0)}\,.\end{split}$ (14)

###  2.2 Finite Remainders

Amplitudes contribute to physical observables only through so-called finite
remainders (see e.g. Weinzierl:2011uz ), which we denote
$\mathcal{R}\left({\bf\it h},{\bf\it p}\right)$. We define a remainder for
each partial amplitude $\mathcal{A}\left({\bf\it h},{\bf\it p}\right)$. They
are obtained from the renormalized amplitudes by removing the remaining
singularities of infrared origin that can be determined from lower-loop
amplitudes and well-known universal factors Catani:1998bh ; Becher:2009cu ;
Gardi:2009qi . The finite remainders can be expanded in powers of
$\alpha_{s}$,

$\mathcal{R}=\mathcal{R}^{(0)}+\frac{\alpha_{s}}{2\pi}\mathcal{R}^{(1)}+\left(\frac{\alpha_{s}}{2\pi}\right)^{2}\mathcal{R}^{(2)}+\mathcal{O}(\alpha_{s}^{3})\,,$
(15)

with the $\mathcal{R}^{(i)}$ defined as

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

The operators ${\bf I}^{(i)}$ are channel specific, and we collect them in
appendix A for all the channels in eqs. 1, 2 and LABEL:eq:qQProc. In these
expressions, we extend the expansion of ${\bf I}^{(1)}\mathcal{A}_{R}^{(1)}$
and ${\bf I}^{(2)}\mathcal{A}_{R}^{(0)}$ to also include terms of order
$\epsilon^{0}$. This subtracts non-trivial contributions from the finite term
of $\mathcal{A}_{R}^{(1)}$ and $\mathcal{A}_{R}^{(2)}$ that are related to the
lower-loop amplitudes. It is clear that the finite remainders
$\mathcal{R}^{(k)}$ contain the genuine new information at order $k$. We
stress that even though our calculations are performed in the HV scheme of
dimensional regularization, the remainders computed in this scheme agree with
the ones computed in the conventional dimensional regularization (CDR) scheme,
see e.g. Broggio:2015dga . Finally, the finite remainders can also be
decomposed into powers of ${N_{f}}$, similarly to eq. 11,

$\mathcal{R}^{(i)}\left({\bf\it h},{\bf\it
p}\right)=\sum_{j=0}^{i}\left(\frac{N_{f}}{N_{c}}\right)^{j}\mathcal{R}^{(i)[j]}\left({\bf\it
h},{\bf\it p}\right)\,.$ (17)

Since the functions $\mathcal{R}^{(i)[j]}$ are dimensionless by construction,
their scale dependence can be restored at any point by rescaling the momenta
as ${\bf\it p}\to{\bf\it p}/\mu$ (we set the factorization scale $\mu_{F}$ to
be equal to the renormalization scale, $\mu_{F}=\mu$). Unless it happens to be
important for the point we wish to make, we will keep this dependence
implicit.

The color-dressed remainder $\mathcal{R}^{(i)}_{\mathcal{M}}\left({\bf\it
h},{\bf\it p}\right)$ is trivially obtained from the decomposition of the
color-dressed amplitude $\mathcal{M}^{(i)}({\bf\it h},{\bf\it p})$ into
partial amplitudes, see eq. 7:

$\displaystyle\begin{split}\mathcal{R}_{\mathcal{M}}({\bf\it h},{\bf\it
p})&=g_{s}^{3}\left(\mathcal{R}_{\mathcal{M}}^{(0)}({\bf\it h},{\bf\it
p})+\frac{\alpha_{s}}{2\pi}\mathcal{R}_{\mathcal{M}}^{(1)}({\bf\it h},{\bf\it
p})+\left(\frac{\alpha_{s}}{2\pi}\right)^{2}\mathcal{R}_{\mathcal{M}}^{(2)}({\bf\it
h},{\bf\it p})+\mathcal{O}(\alpha_{s}^{3})\right)\,,\\\
\mathcal{R}_{\mathcal{M}}^{(i)}({\bf\it h},{\bf\it
p})&=\left(\frac{N_{c}}{2}\right)^{i}\sum_{\sigma\in\Sigma}\mathcal{C}\left(\sigma\right)\Phi\left(\sigma({\bf\it
h}),\sigma({\bf\it p})\right)\mathcal{R}^{(i)}\left(\sigma({\bf\it
h}),\sigma({\bf\it p})\right)\,.\end{split}$ (18)

Whenever the distinction is important, we will call the
$\mathcal{R}^{(i)}\left(\sigma({\bf\it h}),\sigma({\bf\it p})\right)$ partial
remainders, in analogy with the concept of partial amplitudes. We also recall
that, even though this is kept implicit, the color-dressed remainders are
vectors in color space (in general, the $\mathcal{C}\left(\sigma\right)$ have
open (anti-)fundamental and adjoint color indices, see section 3.2 for their
explicit form).

###  2.3 Squared Finite Remainders

When computing physical observables, we are interested in the squared
remainder, summed over color states and, in the case of unpolarized
observables, over helicity states. We denote this object as $H$. It admits an
expansion in powers of $\alpha_{s}$ as

$H=H^{(0)}+\frac{\alpha_{s}}{2\pi}H^{(1)}+\left(\frac{\alpha_{s}}{2\pi}\right)^{2}H^{(2)}+\mathcal{O}(\alpha_{s}^{3}).$
(19)

In terms of the quantities in eq. 18, $H$ is defined as

$H=\frac{1}{\mathcal{N}_{H}}\,\sum_{{\bf\it
h}}\mathcal{R}^{\dagger}_{\mathcal{M}}({\bf\it h},{\bf\it
p})\mathcal{R}_{\mathcal{M}}({\bf\it h},{\bf\it p})\,,$ (20)

where the normalization $\mathcal{N}_{H}$ is fixed such that $H^{(0)}=1$,

$\mathcal{N}_{H}\coloneqq(4\pi\alpha_{s})^{3}\sum_{{\bf\it
h}}\absolutevalue{\mathcal{R}^{(0)}_{\mathcal{M}}({\bf\it h},{\bf\it
p})}^{2}\,.$ (21)

To keep the notation light, in this equation and the following we do not write
explicitly the sum over contracted colour indices which is understood to be
performed in the squaring procedure. At order $\alpha_{s}$ we have

$H^{(1)}=\frac{1}{\sum_{{\bf\it
h}}\absolutevalue{\mathcal{R}^{(0)}_{\mathcal{M}}({\bf\it h},{\bf\it
p})}^{2}}\sum_{{\bf\it
h}}\,2\real\left[\mathcal{R}^{(0)\dagger}_{\mathcal{M}}({\bf\it h},{\bf\it
p})\mathcal{R}^{(1)}_{\mathcal{M}}({\bf\it h},{\bf\it p})\right]\,,$ (22)

and at order $\alpha^{2}_{s}$

$\displaystyle H^{(2)}=\frac{1}{\sum_{{\bf\it
h}}\absolutevalue{\mathcal{R}^{(0)}_{\mathcal{M}}({\bf\it h},{\bf\it
p})}^{2}}\left(\sum_{{\bf\it
h}}\absolutevalue{\mathcal{R}^{(1)}_{\mathcal{M}}({\bf\it h},{\bf\it
p})}^{2}+\sum_{{\bf\it
h}}2\real\left[\mathcal{R}^{(0)\dagger}_{\mathcal{M}}({\bf\it h},{\bf\it
p})\mathcal{R}^{(2)}_{\mathcal{M}}({\bf\it h},{\bf\it p})\right]\right)\,,$
(23)

with each of the two contributions depicted in fig. 2. We note that, when
summing over contracted color indices, we have
$\mathcal{C}^{\dagger}(\sigma^{\prime})~{}\mathcal{C}(\sigma)\sim\delta_{\sigma\sigma^{\prime}}$
in the leading-color approximation we adopt. This greatly simplifies the
construction of the squared remainders.

Figure 2: Schematic representation for the two kinds of contributions to
$H^{(2)}$: the square of one-loop amplitudes, and the interference between
tree-level and two-loop amplitudes.

To be clear about which terms we keep in the leading-color approximation, we
explicitly write the decomposition of our remainders in powers of
${N_{f}}/{N_{c}}$. For $H^{(1)}$, we have

$H^{(1)}=\frac{{N_{c}}}{2}\left(H^{(1)[0]}+\frac{{N_{f}}}{{N_{c}}}H^{(1)[1]}\right)\,.$
(24)

For $H^{(2)}$,

$H^{(2)}=\left(\frac{{N_{c}}}{2}\right)^{2}\left(H^{(2)[0]}+\frac{{N_{f}}}{{N_{c}}}H^{(2)[1]}+\left(\frac{{N_{f}}}{{N_{c}}}\right)^{2}H^{(2)[2]}\right)\,.$
(25)

For most channels the subleading-color corrections are suppressed by a factor
of $1/N_{c}^{2}$. However, for the channels in eq. 4 with two identical quark
lines the suppression is by a single power of $N_{c}$.

##  3 Physical Region Helicity Amplitudes

Our goal in this paper is to obtain the double-virtual contributions for
three-jet production at hadron colliders at NNLO. The central object we need
to compute is the squared finite remainder defined in section 2.3. In this
section, we discuss how we construct the different contributions starting from
the Euclidean-region partial amplitudes obtained in ref. Abreu:2019odu . We
recall the comment below eq. 5 regarding the change of conventions for
denoting each channel between this paper and ref. Abreu:2019odu .

###  3.1 Partial Amplitudes in the Physical Region

The first step towards the construction of the squared finite remainder in the
physical region is to obtain the amplitudes computed in ref. Abreu:2019odu  in
such region. Here we consider only the basis set of partial helicity
amplitudes of ref. Abreu:2019odu , and we will return to their permutations
required to assemble eq. 23 in section 3.2. The amplitudes of ref.
Abreu:2019odu  were computed in the (non-physical) Euclidean region where

$s_{12},s_{23},s_{34},s_{45},s_{15}<0\,,\qquad{\rm tr}_{5}^{2}>0\,.$ (26)

It is our aim to construct results in the physical region associated with the
channels in eqs. 1, 2, LABEL:eq:qQProc, 5 and 4, where $p_{1}$ and $p_{2}$ are
initial-state momenta and $p_{3}$, $p_{4}$ and $p_{5}$ are final-state
momenta. This region is characterized by

$s_{12},s_{34},s_{45}>0,\qquad s_{23},s_{15}<0,\qquad{\rm tr}_{5}^{2}<0.$ (27)

The latter condition is trivially satisfied by real-valued momenta. In this
section, for clarity, we will denote a phase-space point in the physical
region as ${\bf\it p}_{\varphi}$ and a phase-space point in the Euclidean
region as ${\bf\it p}_{E}$. For region-agnostic statements we simply denote a
phase-space point as ${\bf\it p}$.

In order to explain how we convert the results of ref. Abreu:2019odu  to the
kinematic region in eq. 27, let us briefly summarize how they were obtained.
The discussion below holds for each power of ${N_{f}}$, that is for each
$\mathcal{A}^{(i)[j]}$, but for simplicity of the expressions we will suppress
this dependence in intermediate steps. The first step is a computation of the
decomposition of the bare partial amplitudes at one and two loops into a set
of master integrals:

$\mathcal{A}^{(i)}({\bf\it h},{\bf\it p};\epsilon)=\sum_{k}c^{(i)}_{k}({\bf\it
h},{\bf\it p};\epsilon)\,m^{(i)}_{k}({\bf\it p};\epsilon)\,,\qquad i=1,2\,,$
(28)

where we make the dependence on the dimensional regulator $\epsilon$ explicit.
This decomposition was obtained using the two-loop numerical unitarity
approach Ita:2015tya ; Abreu:2017xsl ; Abreu:2017hqn ; Abreu:2018jgq . At this
stage, the master integrals were replaced by their expansion in terms of
Euclidean-region pentagon functions.333These functions are a basis for the
transcendental functions appearing in five-point massless one- and two-loop
master integrals. We refer the reader to ref. Gehrmann:2018yef  for more
details on their definition. Denoting the vector space spanned by relevant
pentagon-function monomials in this region as $\vec{f}({\bf\it
p}_{E})=\\{f_{l}({\bf\it p}_{E})\\}_{l\in B}$, with $B$ denoting the set of
multi-indices $l$, we write

$m^{(i)}_{k}({\bf\it p}_{E};\epsilon)=\sum_{k=0}^{4}\sum_{l\in
B}\epsilon^{k}d^{(i)}_{k,l}f_{l}({\bf\it p}_{E})+\mathcal{O}(\epsilon^{5})\,.$
(29)

In this decomposition, we assume that we have a pure basis of master integrals
ArkaniHamed:2010gh  at one and two loops, in which case the coefficients
$d^{(i)}_{k,l}$ are rational numbers. Furthermore, we assume all integrals are
normalized so that their Laurent expansion around $\epsilon=0$ has no poles,
and stress that the set of pentagon functions is common to one and two loops.
By inserting the decomposition (29) into eq. 28 for both one- and two-loop
amplitudes, and then using these expressions in the definitions for the finite
remainders given in eq. 16, we obtain a decomposition of the two-loop
remainders into pentagon functions:

$\mathcal{R}^{(i)[j]}({\bf\it h},{\bf\it p}_{E})=\sum_{k\in
B}r^{(i)[j]}_{k}({\bf\it h},{\bf\it p}_{E})\,f_{k}({\bf\it p}_{E})\,,$ (30)

where we have reintroduced the superscript $[j]$ to stress that the
decomposition holds for each power of ${N_{f}}$. Obtaining the analytic form
of the coefficients $r^{(i)[j]}_{k}({\bf\it h},{\bf\it p})$ was the main
result of ref. Abreu:2019odu . Let us finish this summary by commenting on the
structure of the coefficients $r^{(i)[j]}_{k}({\bf\it h},{\bf\it p})$. They
can be written as

$r_{k}({\bf\it h},{\bf\it p})=r^{+}_{k}({\bf\it h},\vec{s})+r^{-}_{k}({\bf\it
h},\vec{s})\,{\rm tr}_{5}\,,$ (31)

where we briefly suppress the $i$ and $j$ indices, we recall
$\vec{s}=\\{s_{12},s_{23},s_{34},s_{45},s_{15}\\}$, and $r^{\pm}_{k}({\bf\it
h},\vec{s})$ are rational functions of the Mandelstam variables. The analytic
continuation of the $r_{k}$ from Euclidean to physical region is therefore
trivial. For a given $i$ and $j$, the $r_{k}({\bf\it h},{\bf\it p})$ are not
all linearly independent of each other (over $\mathbb{Q}$) for different
values of $k$, and we can define a basis of rational functions
$\vec{v}({\bf\it h},{\bf\it p})=\\{v_{1}({\bf\it h},{\bf\it p}),\ldots\\}$
such that all $r_{k}({\bf\it h},{\bf\it p})$ live in the $\mathbb{Q}$-linear
span of $\vec{v}({\bf\it h},{\bf\it p})$. We can then write

$\mathcal{R}^{(i)[j]}({\bf\it h},{\bf\it p}_{E})=\vec{v}^{(i)[j]}({\bf\it
h},{\bf\it p}_{E})\cdot M^{(i)[j]}_{E}({\bf\it h})\cdot\vec{f}({\bf\it
p}_{E})\,,$ (32)

where $M_{E}^{(i)[j]}$ is a matrix of rational numbers.

To understand our approach to the analytic continuation of eq. (32), let us
begin by considering what would happen if we had chosen a different basis of
pentagon functions $\vec{f}^{\prime}({\bf\it p}_{E})$, such that
$\vec{f}^{\prime}({\bf\it p}_{E})=N\cdot\vec{f}({\bf\it p}_{E})$ where $N$ is
a matrix of rational numbers. It is clear that to write eq. 32 in this new
basis of functions we do not need to determine the basis $\vec{v}({\bf\it
h},{\bf\it p}_{E})$ of rational functions. The remainders would have exactly
the same form, but with the matrix $M_{E}$ replaced by
$M^{\prime(i)[j]}_{E}=M^{(i)[j]}_{E}\cdot N^{-1}$, which is just a different
matrix of rational numbers.

Let us now discuss what happens with the pentagon functions under analytic
continuation, which we will do through an example. Consider the simplest non-
trivial pentagon function, $\log(-s_{12})$. The choice to write it in this
form follows from the fact that in the Euclidean region $s_{12}<0$ and the
master integrals are real valued. We note that, in this region, we could also
have written $\log(\absolutevalue{s_{12}})$. Let us now continue this function
to the physical region in eq. 27, where $s_{12}>0$. We find

$\log(-s_{12})=\log(\absolutevalue{s_{12}})-\,\mathrm{i}\,\pi\,$ (33)

where the sign is determined by the usual Feynman
$\,\mathrm{i}\,\varepsilon$-prescription of the propagators. We see that the
analytic continuation into the physical region introduces a new pentagon
“function”, the transcendental constant $\,\mathrm{i}\,\pi$. Crucially,
however, the coefficient of this function is generated by the analytic
continuation of $\log(-s_{12})$. This implies that the coefficient of
$\,\mathrm{i}\,\pi$ is trivially related to the coefficient of $\log(-s_{12})$
when analytically continuing eq. 29 or eq. 30, or more precisely that it is in
the $\mathbb{Q}$-linear span of $\vec{v}({\bf\it h},{\bf\it p}_{\varphi})$ (we
recall that ${\bf\it p}_{\varphi}$ is used to denote a point in the physical
region). We assume that this is a general feature of the analytic continuation
of pentagon functions: while the dimension of the vector space in eq. 30 might
increase, the coefficients of the new pentagon functions are not independent
from the coefficients of the pentagon functions in the Euclidean region.444 In
general, the two will be related by the (integer) winding numbers that
naturally appear when analytically continuing multi-valued functions.  Even
though we do not prove this statement, as we shall see below we can explicitly
check that it holds for the remainders in eq. 32.

To analytically continue the remainders computed in ref. Abreu:2019odu , we
use the set of pentagon functions defined in ref. Chicherin:2020oor , and we
denote the vector space of their monomials as $\vec{g}({\bf\it p}_{\varphi})$.
These are valid precisely in the kinematic region of eq. 27, and are related
to the set $\vec{f}({\bf\it p}_{E})$ used in ref. Abreu:2019odu  by a change
of basis and by analytic continuation. Given the points raised above about
these two operations, this means that remainders in the physical region can be
written in terms of the same bases of rational functions that were found for
the Euclidean region. Hence we can write the partial remainder in the physical
region as

$\mathcal{R}^{(i)[j]}_{\varphi}({\bf\it h},{\bf\it
p}_{\varphi};\mu)=\vec{v}^{(i)[j]}({\bf\it h},{\bf\it p}_{\varphi})\cdot
M^{(i)[j]}_{\varphi}({\bf\it h})\cdot\vec{g}({\bf\it p}_{\varphi}/\mu)\,,$
(34)

where for completeness we made the dependence on $\mu$ explicit. We stress
that this representation is valid only when $p_{1}$ and $p_{2}$ are in the
initial state and $p_{3}$, $p_{4}$ and $p_{5}$ are in the final state, because
that is the region of validity of the pentagon functions $\vec{g}({\bf\it
p}_{\varphi})$. We discuss the generalization to other configurations in
section 3.2. The $M^{(i)[j]}_{\varphi}({\bf\it h})$ can be determined from a
small number of numerical evaluations of the amplitudes. To do this, we take
equation eq. 34 as an ansatz for the form of the remainder when expressed in
terms of the physical region pentagon functions $\vec{g}({\bf\it
p}_{\varphi})$, with the entries of the $M^{(i)[j]}_{\varphi}({\bf\it h})$ as
unknowns. We then constrain the unknowns by performing the analogous procedure
of eqs. (29) and (30) at a sufficient number of phase-space points.555This
number is dictated by the number of independent functions in
$\vec{v}^{(i)[j]}({\bf\it h},{\bf\it p})$. The most complicated case, the
five-gluon remainder
$\mathcal{R}^{(2)[0]}_{\varphi}(\\{-1,+1,-1,+1,+1\\},{\bf\it p})$, requires 57
evaluations. We note that a single numerical evaluation gives more than one
constraint. Linear algebra then allows us to extract the values of
$M^{(i)[j]}_{\varphi}({\bf\it h})$.

In practice, we use the implementation of two-loop numerical unitarity in
Caravel Abreu:2020xvt  to generate all the required numerical data. For
efficiency reasons, the procedure is performed over finite fields and we
construct the $\mathbb{Q}$-valued $M^{(i)[j]}_{\varphi}({\bf\it h})$ via the
Chinese remainder theorem and standard rational reconstruction techniques. We
remark that we used at most three finite fields of cardinality
$\mathcal{O}(2^{32})$. Following these steps, we obtain the analytic
continuation of the partial amplitudes computed in ref. Abreu:2019odu  to the
physical region defined in eq. 27, corresponding to the channels in eqs. 1, 2,
LABEL:eq:qQProc, 5 and 4. Crucially, this linear algebra exercise would not
allow one to compute the matrices $M^{(i)[j]}_{\varphi}({\bf\it h})$ in the
ansatz (34) if the assumption on the analytic continuation of the pentagon
functions did not hold.

###  3.2 From Partial Amplitudes to Amplitudes

For physical applications, it is not sufficient to know a single partial
amplitude for each channel. Indeed, we must be able to evaluate all the terms
in the sum in eq. 7 for each helicity configuration. We start by summarizing
this information for each of the channels listed in eqs. 1, 2,
LABEL:eq:qQProc, 5 and 4. We remind the reader of the comment on momenta
conventions under eq. 5.

#### Five gluon:

There is a single channel, see eq. 1. The set $\Sigma=S_{5}/Z_{5}$ contains 24
elements, and the color factor for $\sigma_{0}=\\{1,2,3,4,5\\}$ is

$\mathcal{C}(\sigma_{0})\coloneqq\tr{T^{a_{1}}T^{a_{2}}T^{a_{3}}T^{a_{4}}T^{a_{5}}}\,.$
(35)

For each of the 24 partial amplitudes, there are 32 different helicity
configurations. Accounting for relations under charge conjugation and parity,
there are 192 partial helicity amplitudes to compute. We can further reduce
this set by considering the $S_{2}\times S_{3}$ permutations of initial state
and final state labels which, crucially, still correspond to the physical
region defined at the start of section 3.1. This leaves us with a generating
set of 32 amplitudes.

#### Single quark line:

There are three channels, see eq. 2. For the channel
$\bar{q}^{h_{1}}_{p_{1}}q^{h_{2}}_{p_{2}}\to\mathrm{g}^{h_{3}}_{p_{3}}\mathrm{g}^{h_{4}}_{p_{4}}\mathrm{g}^{h_{5}}_{p_{5}}$,
the set $\Sigma=\\{1,2\\}\otimes S_{3}(3,4,5)$ contains 6 elements, and the
color factor for $\sigma_{0}=\\{1,2,3,4,5\\}$ is

$\mathcal{C}(\sigma_{0})\coloneqq\quantity(T^{a_{3}}T^{a_{4}}T^{a_{5}})^{\bar{i}_{2}}_{i_{1}}\,.$
(36)

For each of the three channels there are 6 partial amplitudes, each with 16
different helicity configurations. Accounting for relations under charge
conjugation and parity, there are 144 partial helicity amplitudes to compute.
Considering the $S_{2}\times S_{3}$ permutation of initial state and final
state labels further reduces the cardinality of the generating set to 40.

#### Two quark lines:

There are four channels, see LABEL:eq:qQProc, but we also include the two
channels in eq. 5 in order to consider two identical quark lines. For the
channel $\bar{q}^{h_{1}}_{p_{1}}q^{h_{2}}_{p_{2}}\to
Q^{h_{3}}_{p_{3}}\bar{Q}^{h_{4}}_{p_{4}}\mathrm{g}^{h_{5}}_{p_{5}}$, the set
$\Sigma=\\{\\{1,2,3,4,5\\},\\{3,4,1,2,5\\}\\}$ contains 2 elements, and the
color factor for $\sigma_{0}=\\{1,2,3,4,5\\}$ is

$\mathcal{C}(\sigma_{0})\coloneqq\quantity(T^{a_{5}})^{\bar{i}_{2}}_{i_{3}}\delta^{\bar{i}_{4}}_{i_{1}}\,.$
(37)

For each of the six channels there are 2 partial amplitudes, each with 8
different helicity configurations. Accounting for relations under charge
conjugation and parity, there are 48 partial helicity amplitudes to compute.
Considering the $S_{2}\times S_{3}$ permutation of initial state and final
state labels further reduces the cardinality of the generating set to 24.

Having classified a generating set of partial amplitudes for each channel
which respect the region constraints of the pentagon functions, we can now
apply an analogous procedure to that of section 3.1 to each of these partial
amplitudes. Let us consider the remainder of a partial amplitude where the
helicities and momenta are related by a permutation $\sigma$ to those in eq.
34. Clearly, the master integrals in the decomposition (28) of the permuted
amplitude are related to those of the original amplitude by the permutation
$\sigma$. We remind the reader that the pentagon functions of ref.
Chicherin:2020oor  form a basis of the space of transcendental functions with
$p_{1}$ and $p_{2}$ initial state and $p_{3}$, $p_{4}$ and $p_{5}$ final
state. Hence we can express any permutation of the master integrals with this
configuration of the momenta as

$m^{(i)}_{k}(\sigma({\bf\it p});\epsilon)=\sum_{k=0}^{4}\sum_{l\in
B_{\varphi}}\epsilon^{k}d^{(i)}_{k,l}(\sigma)g_{l}({\bf\it
p})+\mathcal{O}(\epsilon^{5})\,,$ (38)

where $B_{\varphi}$ is the set of labels distinguishing the elements of
$\vec{g}({\bf\it p})$. It therefore follows that all permutations of the
partial remainders contributing to the remainder in eq. 18 can be expressed in
the same basis of pentagon functions. That is, a partial finite remainder
evaluated with a permutation of the input helicities $\sigma({\bf\it h})$ and
on a permutation of the input point $\sigma({\bf\it p})$ can be expressed in
the form

$\mathcal{R}^{(i)[j]}_{\varphi}(\sigma({\bf\it h}),\sigma({\bf\it
p});\mu)=\vec{v}^{(i)[j]}(\sigma({\bf\it h}),\sigma({\bf\it p}))\cdot
M^{(i)[j]}_{\varphi}(\sigma({\bf\it h}),\sigma)\cdot\vec{g}({\bf\it
p}/\mu)\,,$ (39)

where, as in eq. 34, we made the dependence on $\mu$ explicit. We recall that
we can easily construct the bases $\vec{v}^{(i)[j]}(\sigma({\bf\it
p}),\sigma({\bf\it h}))$, as the action of $\sigma$ on rational functions is
trivial to evaluate. The matrices $M^{(i)[j]}_{\varphi}(\sigma({\bf\it
h}),\sigma)$ can then be determined by the same ansatz procedure described for
the analytic continuation of the results of ref. Abreu:2019odu . Altogether,
this allows us to obtain all the expressions required to evaluate the squared
finite remainders defined in section 2.3 for the channels in eqs. 1, 2,
LABEL:eq:qQProc, 5 and 4.

##  4 Results and Validation

In this section we describe our main results, which are of two types. The
first is a collection of ancillary files containing analytic expressions for a
generating set of (partial) remainders for each of the channels in eqs. 1, 2,
LABEL:eq:qQProc, 5 and 4. The second is a C++ library that allows for an
efficient numerical evaluation of the finite remainders defined in eq. 34 for
all required permutations, as well as the squared finite remainders defined in
eqs. 22 and 23. Finally, we describe the checks we have performed on our
results.

###  4.1 Analytic Results

In the ancillary files we give all necessary ingredients to assemble analytic
expressions for a generating set of two-loop partial remainders
$\mathcal{R}^{(2)[j]}_{\varphi}(\sigma({\bf\it h}),\sigma({\bf\it p}))$, for
each of the channels in eqs. 1, 2, LABEL:eq:qQProc, 5 and 4. From this
generating set, all the permutations required in eq. 18 can be constructed by
applying (a combination of) charge conjugation, parity transformations,666
Parity transformation of the analytic expressions are performed by changing
the signs of all parity-odd pentagon functions and ${\rm tr}_{5}$.  and the
$S_{2}\crossproduct S_{3}$ permutations of initial and final state labels.777
One may wish to express the $S_{2}\crossproduct S_{3}$ permutations of the
partial remainders from the generating set in terms of the pentagon functions
evaluated on a non-permuted phase-space point. This can be achieved by
following the procedure described in section 3.2.

The expressions we present in the ancillary files allow to explicitly
construct the generating set of remainders in the decomposition of eq. 39. The
basis of transcendental functions $\vec{g}({\bf\it p})$ is universal across
all channels and powers of $N_{f}$, and is given in the file
anc/pentagonFunctionBasis.m, written in the notation of the Mathematica
package provided with ref. Chicherin:2020oor . The basis of rational functions
$\vec{v}^{(2)[j]}(\sigma({\bf\it h}),\sigma({\bf\it p}))$ is channel and
power-of-$N_{f}$ specific but trivially dependent on $\sigma$ (more precisely,
the action of $\sigma$ on $\vec{v}^{(2)[j]}(\sigma({\bf\it h}),\sigma({\bf\it
p}))$ simply amounts to the permutation of the labels in the Mandelstam
invariants appearing in the rational functions). For each channel and power of
$N_{f}$, we simply reuse the bases in the ancillary files of ref.
Abreu:2019odu  (we note again the comment below eq. 5 regarding the
conventions to denote each channel). For completeness, we include them in the
file anc/rationalBases.m in the format

{channel, helicity, Nf} -> {list of rational functions in ($\vec{s}$, ${\rm
tr}_{5}$)}.

We note that one can alternatively use the form of the same bases of rational
functions determined in ref. DeLaurentis:2020qle , which are more compactly
expressed in the spinor-helicity formalism. Finally, the matrices
$M^{(2)[j]}_{\varphi}(\sigma({\bf\it h}),\sigma)$ for each member of the
generating sets of remainders are given explicitly in
anc/twoLoopProjectionMatrices.m. Each is expressed in the form

{channel, helicity, NfPower, permutation} -> matrix.

where matrix is in Mathematica’s SparseArray format.

As described in section 3, these matrices constitute the main analytic result
of this paper. Indeed, they were the only missing information required to
obtain analytic expressions for the five-parton amplitudes in the physical
region corresponding to three-jet production at hadron colliders given in eq.
27.

###  4.2 Numerical Evaluation

Having phenomenological applications in mind, we implemented a C++ library
that allows for the efficient evaluation of finite remainders and squared
finite remainders. In this section we comment on this numerical
implementation.

Let us start by describing the quantities which can be evaluated with our C++
library FivePointAmplitudes . The code performs the numerical evaluation of
one- and two-loop (partial) remainders for each channel, each helicity
configuration and each power of $N_{f}$. That is, we evaluate each of the
$\mathcal{R}^{(i)[j]}\left(\sigma({\bf\it h}),\sigma({\bf\it p})\right)$ for
$i=1,2$ and present each power of $(N_{f}/N_{c})^{j}$ separately. We stress
that we have implemented all permutations, and not only the generating set for
which we provide analytic expressions. For physical applications, one might be
interested in the squared finite remainders defined in eqs. 22 and 23. Our
code also outputs the numerical values for these quantities, once again
presenting each power of $N_{f}/N_{c}$ separately, that is the $H^{(i)[j]}$ in
eqs. 24 and 25.

We employ PentagonFunctions++ Chicherin:2020oor  for the numerical evaluation
of the pentagon functions. For the efficient evaluation of the rational
functions, we generate optimized code with FORM Ruijl:2017dtg ;
Kuipers:2013pba . In any of the helicity remainders
$\mathcal{R}^{(2)[j]}\left(\sigma({\bf\it h}),\sigma({\bf\it p})\right)$ that
we compute numerically, the dominant contribution to the evaluation time is
that of the pentagon functions. When computing $H^{(1)}$ and $H^{(2)}$, we
have thus organized our numerical code such that the pentagon functions are
evaluated only once per phase-space point, independently of the number of
permutations contributing to the sums in eqs. 22 and 23. We note nevertheless
that, compared to e.g. the case of triphoton production Abreu:2020cwb , the
evaluation of $H^{(2)}$ for five-gluon channels receives contributions from a
large number of partial amplitudes, and therefore the time spent in the
evaluation of the other components in eq. 39 becomes noticeable. We find that
the typical double-precision evaluation time for any $H^{(2)}$ is
$\mathcal{O}(1\text{s})$ on a single CPU, while the evaluation time for
$H^{(1)}$ is negligible at $\mathcal{O}(10\text{ms})$.

To be suitable for phenomenological applications the numerical evaluations
should not only be fast but also stable across phase space. To guarantee the
latter, we develop a precision-rescue strategy that detects and corrects
unstable numerical evaluations. The main source of numerical instabilities is
the presence of large cancellations between different contributions in eq. 39,
when a phase-space point gets close to the zero sets of the rational
functions’ denominators. As we have already argued in refs. Abreu:2018zmy ;
Abreu:2019odu ; Abreu:2020cwb , these are determined by a set of 25 quantities
$\\{W_{i}(\vec{s})\\}$, with $i=1,\ldots,25$, which are a subset of the so-
called alphabet associated with the five-point two-loop master integrals (see
e.g. ref. Chicherin:2017dob  for the definition of these quantities). This
subset is rational in the Mandelstam variables $\vec{s}$ defined in eq. 9 and
all elements have the same dimensions. We find that potentially unstable
points can be efficiently identified by computing a quantity that
characterizes the spread of scales for a given phase-space point
$\vec{s}$,888We attempted to further reduce the subset of letters to only
include those that correspond to spurious singularities, as these had been
identified in ref. Abreu:2020cwb  as being sufficient to characterize the
points that were less stable. We find that for the five-parton channels
considered in this paper this subset covers most of the unstable points, but
we must also include non-spurious singularities to detect all unstable points.

$\kappa(\vec{s})\coloneqq\frac{\min_{i}\\{\absolutevalue{W_{i}(\vec{s})}\\}}{\max_{i}\\{\absolutevalue{W_{i}(\vec{s})}\\}}\,.$
(40)

We introduce a threshold $\kappa_{\text{thr}}$, and if
$\kappa(\vec{s})>\kappa_{\text{thr}}$ the point is considered stable.
Otherwise, the point is potentially unstable and we perform a second
evaluation on an infinitesimally perturbed point $\vec{s}_{\delta}$ defined as

$\vec{s}_{\delta}-\vec{s}\simeq{\varepsilon_{\text{double}}}\,\vec{s}\,,$ (41)

where $\varepsilon_{\text{double}}\simeq 10^{-16}$ is the machine epsilon of
double-precision floating point numbers.999In practice, each of the components
$(\vec{s}_{\delta})_{i}$ is taken to be the next (or previous) representable
floating-point number after $(\vec{s})_{i}$.  We estimate the accuracy of
potentially unstable points as

$\Delta(\vec{s})\coloneqq\absolutevalue{1-\frac{H(\vec{s}_{\delta})}{H(\vec{s})}}\,.$
(42)

We introduce another threshold value $\Delta_{\text{thr}}$, and if
$\Delta(\vec{s})<\Delta_{\text{thr}}$ the point is considered stable. If
$\Delta(\vec{s})>\Delta_{\text{thr}}$, the point $\vec{s}$ is considered
unstable and evaluated in quadruple precision. The performance of this rescue
strategy and the values we choose for $\kappa_{\text{thr}}$ and
$\Delta_{\text{thr}}$ will be discussed below.

Figure 3: Logarithmic distribution of correct digits (see eq. 44) on samples
of 110k points for each of the channels. The blue line corresponds to
evaluation in double precision. The orange line represents evaluation with the
rescue system enabled. The gray dashed line is the cumulative distribution of
the latter. The rescue system’s thresholds are the same as in table 1. Average
core time for evaluation on a single phase-space point with rescue system
enabled as measured when evaluating all sample points in parallel using all 32
threads of Intel(R) Xeon(R) Silver 4216 CPU @ 2.10GHz.

To demonstrate the performance of our numerical C++ implementation, we employ
Sherpa 2.2 Bothmann:2019yzt  to sample phase-space points from a Monte-Carlo
integration grid optimized on the Born matrix elements for three-jet
production. We sample a sufficient number of points to draw 110k events for
each of the channels $\mathrm{g}\mathrm{g}\to\mathrm{g}\mathrm{g}\mathrm{g}$,
$\bar{q}\mathrm{g}\to\bar{q}\mathrm{g}\mathrm{g}$ and $\bar{q}q\to
Q\bar{Q}\mathrm{g}$, which are representative of all channels in eqs. 1, 2 and
LABEL:eq:qQProc. We use the same phase-space definition that was used in the 7
TeV ATLAS analysis Aad:2014rma , which we summarize here for completeness. We
require three hard jets defined with the anti-$k_{t}$ algorithm, with the
radius parameter of $R=0.4$, as implemented in FastJet Cacciari:2011ma . We
require the jets to be within the rapidity range $\absolutevalue{y}<3$ and
that they all have transverse momenta larger than $50$ GeV. In addition, the
leading and subleading jets are required to have transverse momenta larger
than $150$ GeV and $100$ GeV respectively. We then evaluate $H^{(2)}$ at the
dynamical scale

$\mu=\frac{1}{2}\sum_{j=1}^{3}p_{Tj}\,,$ (43)

and sum the contributions of different powers of ${N_{f}}/{N_{c}}$, setting
${N_{c}}=3$ and ${N_{f}}=5$. The accuracy of double-precision evaluations is
determined by comparing them to an evaluation in quadruple precision as

$d\coloneqq-\log_{10}\absolutevalue{1-\frac{H^{(2)}_{\text{double}}}{H^{(2)}_{\text{quad}}}}\,,$
(44)

where we keep the dependence on the phase-space point implicit. If the rescue
system is turned on, then the accuracy of the rescued points that were
evaluated in quadruple precision is estimated with a second quadruple-
precision evaluation as in eq. 42. In fig. 3 we present the distribution of
correct digits for the three representative channels, obtained over the 110k
phase-space points. We plot the curve with the rescue system turned off (blue
lines), which shows that a rescue system is indeed required for a few points.
In orange we plot the curve with the rescue system turned on, with
$\kappa_{\text{thr}}=10^{-4}$ and $\Delta_{\text{thr}}=10^{-2}$, and we see
that all unstable points have been correctly rescued and the numerical
stability is adequate. The gray dashed line is the cumulative distribution
with the rescue system enabled. On the plots, we also overlay the average core
time for a single phase-space point evaluation with the rescue system turned
on.

To assess the impact of the rescue system on the timings of the numerical
evaluation, we collected in table 1 some characterizing data for the choice of
thresholds $\kappa_{\text{thr}}=10^{-4}$ and $\Delta_{\text{thr}}=10^{-2}$.
The test in eq. 40 catches too many points, but crucially catches all unstable
points as can be seen from fig. 3. The second test, in eq. 42, is much more
precise but also more expensive in terms of evaluation time. Overall, we see
that our rescue system catches all necessary points with only a mild slowdown
factor. The values $\kappa_{\text{thr}}=10^{-4}$ and
$\Delta_{\text{thr}}=10^{-2}$ can be modified by the user according to their
needs. It would be interesting to investigate if using the different form of
the bases of rational functions derived in ref. DeLaurentis:2020qle  reduces
the fraction of unstable points, but we find that our implementation is
already satisfactorily fast and stable for phenomenological applications.

Channels |  % recomputed in double |  % recomputed in quad |  slowdown factor
---|---|---|---
$H^{(2)}_{\mathrm{\mathrm{g}}\,\mathrm{\mathrm{g}}~{}\to~{}\mathrm{\mathrm{g}}\,\mathrm{\mathrm{g}}\,\mathrm{\mathrm{g}}}$ |  2.52 |  0.27 |  1.12
$H^{(2)}_{\bar{\mathrm{q}}\,\mathrm{\mathrm{g}}~{}\to~{}\bar{\mathrm{q}}\,\mathrm{\mathrm{g}}\,\mathrm{\mathrm{g}}}$ |  2.57 |  0.08 |  1.06
$H^{(2)}_{\bar{\mathrm{q}}\,\bar{\mathrm{Q}}~{}\to~{}\bar{\mathrm{q}}\,\bar{\mathrm{Q}}\,\mathrm{\mathrm{g}}}$ |  2.52 |  0.01 |  1.03
Table 1:  Characteristic performance of the rescue system on the samples of
110k points with $\kappa_{\text{thr}}=10^{-4}$ and
$\Delta_{\text{thr}}=10^{-2}$. The first column shows the fraction of points
with $\kappa(\vec{s})<\kappa_{\text{thr}}$, which are recomputed in double
precision at $\vec{s}_{\delta}$. The second column shows the fraction of
points with $\Delta(\vec{s})>\Delta_{\text{thr}}$, which are recomputed in
quadruple precision at $\vec{s}$. The average slowdown due to the usage of the
rescue system w.r.t. only double-precision evaluations is shown in the last
column.

In order to facilitate the comparison with our results, we present in appendix
B the numerical values of squared finite remainders at each power of
${N_{f}}$, see eqs. 24 and 25, at a randomly chosen phase-space point.

###  4.3 Validation

The main analytic results of our paper are the matrices described in section
4.1. These are obtained by building on the analytic results of ref.
Abreu:2019odu , which have undergone very stringent checks. Compared to ref.
Abreu:2019odu , our current results are written in terms of a different basis
of transcendental functions Chicherin:2020oor  which have also undergone
several direct checks (e.g., comparison with the functions of ref.
Gehrmann:2018yef ) and indirect checks Kallweit:2020gcp ; Abreu:2020cwb ;
Chawdhry:2020for ; Agarwal:2021grm . Furthermore, we have verified that we
reproduce the results independently obtained in ref. Badger:2018gip  for the
numerical evaluation of five-gluon amplitudes in the physical region defined
in eq. 27. Finally, we compared our results for the function
$H^{(2)[0]}_{\mathrm{g}\mathrm{g}\to\mathrm{g}\mathrm{g}\mathrm{g}}$ in the
purely gluonic channel against an independent computation
Badger:correspondence  and found agreement.

We have also performed several consistency checks on our numerical results.
Since we obtain all the helicity partial remainders required for the functions
$H^{(1)}$ and $H^{(2)}$ from sums over the permutations of the generating sets
of remainders introduced in section 3.2, it is important to check that these
operations are correctly implemented in the C++ library we have developed. We
have made three types of consistency checks to validate this procedure.

#### Symmetries.

We have checked that the remainders and squared remainders satisfy all
expected charge, parity conjugation and permutation symmetries.

#### Validation of $H^{(1)}$.

We have compared the one-loop squared remainders $H^{(1)}$ of all the channels
in eqs. 1, 2, LABEL:eq:qQProc, 5 and 4 against the evaluations by the BlackHat
library Berger:2008sj  and found full agreement.101010 After adjusting the
definitions of the leading-color squared matrix elements in BlackHat such that
all subleading-color contributions are discarded.  The color-dressed
remainders $\mathcal{R}^{(1)}_{\mathcal{M}}$ in eq. 18 are constructed from
generating sets of partial remainders, as described in section 3.2, and then
are directly fed into the definition of $H^{(1)}$, see eq. 22. The procedure
for assembling the squared remainders is independent of the number of loops,
hence this check validates the corresponding assembly of the contributions for
the calculation of $H^{(2)}$.

#### Scale dependence of remainders.

Given that we evaluate the finite remainders directly, the correct pole
structure of the underlying scattering amplitudes is satisfied by
construction. We note nevertheless that finite remainders retain a memory of
the poles of the amplitudes and of their explicit (scheme-dependent)
definition. This information is encoded in the scale dependence of the finite
remainders. Quite generally, we find that

$\displaystyle\begin{split}\frac{\partial\,\mathcal{R}^{(1)}(\mu)}{\partial\log\mu}~{}=~{}&\frac{2\beta_{0}\lambda}{N_{c}}\mathcal{R}^{(0)}\,,\\\
\frac{\partial\,\mathcal{R}^{(2)}(\mu)}{\partial\log\mu}~{}=~{}&2\left(\frac{\beta_{1}\lambda}{N_{c}^{2}}+{\bf
H}_{[n]}\right)\mathcal{R}^{(0)}+\frac{2\beta_{0}\left(2+\lambda\right)}{N_{c}}\mathcal{R}^{(1)}(\mu)\,,\end{split}$
(45)

where we denote by $\lambda$ the power of $g_{s}^{0}$ in the tree-level
amplitude (for five-parton amplitudes $\lambda=3$, see eq. 6) and the
$\beta_{i}$ are the coefficients of the QCD $\beta$ function explicitly given
in eq. 13. ${\bf H}_{[n]}$ is given explicitly in eqs. 51 and 52 and controls
the collinear singularities of scattering amplitudes. We have verified that
the scale dependence of the remainders obtained from our numerical
implementation is described by eq. 45. Furthermore, we have independently
checked that, starting from this scale dependence, we reproduce the scale
dependence of the functions $H^{(1)}$ and $H^{(2)}$, which again validates the
correct assembly of these functions.

##  5 Conclusion

In this paper we present the NNLO double-virtual contributions to three-jet
production at hadron colliders at leading color. While all required amplitudes
were already available in the literature in analytic form, they were only
available in the non-physical Euclidean region, and as such could not be
directly used for phenomenological applications. Aside from analytically
continuing those expressions to the physical region, we also provide a C++
library that allows to efficiently evaluate these contributions and obtain
stable results across phase space.

The analytic continuation of the results of ref. Abreu:2019odu  to the
physical region corresponding to three-jet production is performed by
rewriting them in terms of the pentagon functions defined in ref.
Chicherin:2020oor . By noticing that this rewriting amounts to computing a
matrix of rational numbers, we obtained the physical-region expressions from a
few numerical evaluations of the amplitudes (57 in the most complicated case).
These evaluations were performed with the implementation of the two-loop
numerical unitarity approach in Caravel Abreu:2020xvt . We provide analytic
expressions for a generating set of partial remainders. For all channels in
eqs. 1, 2, LABEL:eq:qQProc, 5 and 4, these allow one to evaluate the color-
dressed remainders defined in eq. 18 as well as the squared remainders
$H^{(1)}$ and $H^{(2)}$ defined in eqs. 24 and 25 in the physical region.

Our driving motivation was to provide results that are ready to be used for
phenomenological applications. To this end, we developed a C++ library which
efficiently evaluates individual partial remainders as well as the one- and
two-loop squared remainders. The numerical evaluations can be performed in
double or quadruple precision. We find that double precision is sufficient for
the vast majority of phase-space points, and have implemented a rescue system
which detects unstable points and reevaluates them in quadruple precision.
With the rescue system enabled, we find very good overall numerical stability
and our evaluation time for the two-loop squared remainders is
$\mathcal{O}(1\text{s})$. This shows that our results and the numerical
library we provide are ready to be used in phenomenological studies.

## Acknowledgments

We thank Simon Badger and Simone Zoia for the correspondence on an independent
calculation of the finite remainders in the purely gluonic channel. The work
of B.P is supported by the French Agence Nationale pour la Recherche, under
grant ANR–17–CE31–0001–01. V.S. is supported by the European Research Council
(ERC) under the European Union’s Horizon 2020 research and innovation
programme, Novel structures in scattering amplitudes (grant agreement No.
725110). The work of F. F. C. is supported in part by the U.S. Department of
Energy under grant DE-SC0010102.

##  Appendix A Infrared Subtraction

The divergences of renormalized two-loop amplitudes obey a universal structure
Catani:1998bh ; Becher:2009cu ; Gardi:2009qi ,

$\displaystyle\begin{split}\mathcal{A}_{R}^{(1)}&={\bf
I}^{(1)}_{[n]}(\epsilon)\mathcal{A}_{R}^{(0)}+\mathcal{O}(\epsilon^{0})\,,\\\
\mathcal{A}_{R}^{(2)}&={\bf I}^{(2)}_{[n]}(\epsilon)A_{R}^{(0)}+{\bf
I}^{(1)}_{[n]}(\epsilon)\mathcal{A}_{R}^{(1)}+\mathcal{O}(\epsilon^{0})\,,\end{split}$
(46)

with the renormalized amplitudes $\mathcal{A}_{R}^{(i)}$ related to the bare
amplitudes $\mathcal{A}^{(i)}$ as in LABEL:eq:twoLoopUnRenorm. For amplitudes
in the leading-color approximation the operators $\mathbf{I}^{(1)}_{[n]}$ and
$\mathbf{I}^{(2)}_{[n]}$ are diagonal in color space. The operator
$\mathbf{I}^{(1)}_{[n]}$ is given by

${\bf
I}^{(1)}_{[n]}(\epsilon)=-\frac{e^{\gamma_{E}\epsilon}}{\Gamma(1-\epsilon)}\sum_{i=1}^{n}\gamma_{a_{i},a_{i+1}}\left(-\frac{s_{i,i+1}}{\mu^{2}}-\,\mathrm{i}\,\varepsilon\right)^{-\epsilon}\,,$
(47)

where $s_{i,j}=(p_{i}+p_{j})^{2}$, with the indices defined cyclically, and we
made explicit the $\,\mathrm{i}\,\varepsilon$ associated with Mandelstam
variables that follows from the Feynman
$\,\mathrm{i}\,\varepsilon$-prescription of the propagators. The index $a_{i}$
denotes a type of particle, i.e.
$a_{i}\in\\{\mathrm{g},q,\bar{q},Q,\bar{Q}\\}$. The symbols $\gamma_{a,b}$ are
symmetric under exchange of indices, $\gamma_{a,b}=\gamma_{b,a}$, and given
by:

$\displaystyle\begin{split}\gamma_{\mathrm{g},\mathrm{g}}&=\frac{1}{\epsilon^{2}}+\frac{1}{\epsilon}\frac{\beta_{0}}{{N_{c}}}\,,\qquad\gamma_{q,Q}=\gamma_{q,\bar{Q}}=\gamma_{\bar{q},Q}=\gamma_{\bar{q},\bar{Q}}=\frac{1}{\epsilon^{2}}+\frac{3}{2\epsilon}\,,\\\
\gamma_{\mathrm{g},q}&=\gamma_{\mathrm{g},\bar{q}}=\gamma_{\mathrm{g},Q}=\gamma_{\mathrm{g},\bar{Q}}=\frac{\gamma_{\mathrm{g},\mathrm{g}}+\gamma_{q,Q}}{2}\,,\qquad\gamma_{q,\bar{q}}=\gamma_{Q,\bar{Q}}=0\,.\end{split}$
(48)

The operator ${\bf I}^{(2)}_{[n]}$ is

$\displaystyle\begin{split}{\bf I}^{(2)}_{[n]}(\epsilon)=&-\frac{1}{2}{\bf
I}^{(1)}_{[n]}(\epsilon){\bf
I}^{(1)}_{[n]}(\epsilon)-\frac{2\beta_{0}}{{N_{c}}\epsilon}{\bf
I}^{(1)}_{[n]}(\epsilon)+\frac{e^{-\gamma_{E}\epsilon}\Gamma(1-2\epsilon)}{\Gamma(1-\epsilon)}\left(\frac{2\beta_{0}}{{N_{c}}\epsilon}+K\right){\bf
I}^{(1)}_{[n]}(2\epsilon)\\\
&+\frac{e^{\gamma_{E}\epsilon}}{\epsilon\Gamma(1-\epsilon)}{\bf
H}_{[n]}\,,\end{split}$ (49)

where

$K=\frac{67}{9}-\frac{\pi^{2}}{3}-\frac{10}{9}\frac{{N_{f}}}{{N_{c}}}\,,$ (50)

and ${\bf H}_{[n]}$ is a diagonal operator at leading color that depends on
the number of external quarks and gluons,

$\displaystyle\begin{split}{\bf
H}_{[n]}&=\sum_{i=1}^{n}\left(\delta_{a_{i},\mathrm{g}}H_{\mathrm{g}}+(\delta_{a_{i},q}+\delta_{a_{i},\bar{q}}+\delta_{a_{i},Q}+\delta_{a_{i},\bar{Q}})H_{q}\right)\,,\end{split}$
(51)

with

$\displaystyle\begin{split}H_{\mathrm{g}}&=\left(\frac{\zeta_{3}}{2}+\frac{5}{12}+\frac{11\pi^{2}}{144}\right)-\left(\frac{\pi^{2}}{72}+\frac{89}{108}\right)\frac{N_{f}}{{N_{c}}}+\frac{5}{27}\left(\frac{N_{f}}{{N_{c}}}\right)^{2}\,,\\\
H_{q}&=\left(\frac{7\zeta_{3}}{4}+\frac{409}{864}-\frac{11\pi^{2}}{96}\right)+\left(\frac{\pi^{2}}{48}-\frac{25}{216}\right)\frac{N_{f}}{{N_{c}}}\,.\end{split}$
(52)

##  Appendix B Reference Evaluations of Squared Finite Remainders

Channel $H^{(1)[0]}$ $H^{(1)[1]}$ $H^{(2)[0]}$ $H^{(2)[1]}$ $H^{(2)[2]}$
$\mathrm{g}\mathrm{g}\rightarrow\mathrm{g}\mathrm{g}\mathrm{g}$ $16.135254222$
$0.19163044752$ $464.47846208$ $-58.116292408$ $0.60077232705$
${\bar{q}}q\rightarrow\mathrm{g}\mathrm{g}\mathrm{g}$ $9.5879406141$
$-3.0604943308$ $184.44415807$ $-61.765802987$ $6.3615768297$
${\bar{q}}\mathrm{g}\rightarrow{\bar{q}}\mathrm{g}\mathrm{g}$ $26.908169290$
$-3.6373308269$ $867.25232363$ $-230.76277359$ $12.598811302$
$\mathrm{g}\mathrm{g}\rightarrow q{\bar{q}}\mathrm{g}$ $24.495592766$
$-2.5939909248$ $745.87682394$ $-166.84486839$ $6.1899943330$
${\bar{q}}q\rightarrow Q{\bar{Q}}\mathrm{g}$ $10.460907919$ $-4.2060557725$
$212.42454564$ $-80.136400792$ $8.2094005806$ ${\bar{q}}Q\rightarrow
Q{\bar{q}}\mathrm{g}$ $27.104747640$ $-4.0829938180$ $705.58902507$
$-209.42216177$ $12.483148067$
${\bar{q}}{\bar{Q}}\rightarrow{\bar{q}}{\bar{Q}}\mathrm{g}$ $42.313652168$
$-8.0064067852$ $1628.2933493$ $-562.78735847$ $44.198947852$
${\bar{q}}\mathrm{g}\rightarrow{\bar{q}}Q{\bar{Q}}$ $28.068256507$
$-6.3593609865$ $935.81439233$ $-324.32790785$ $29.070926975$
${\bar{q}}q\rightarrow q{\bar{q}}\mathrm{g}$ $20.846053179$ $-4.1292696285$
$520.14108472$ $-160.80597165$ $10.876062192$
${\bar{q}}{\bar{q}}\rightarrow{\bar{q}}{\bar{q}}\mathrm{g}$ $42.259655399$
$-7.9918854619$ $1624.7163564$ $-561.33769564$ $44.056509019$
${\bar{q}}\mathrm{g}\rightarrow{\bar{q}}q{\bar{q}}$ $28.497167934$
$-6.2611415380$ $947.84964732$ $-322.54996102$ $28.093290494$

Table 2:  Reference values for the evaluation of squared finite remainders at
each power of ${N_{f}}$, as defined in eqs. 24 and 25 on the phase-space point
given in eq. 53.

In order to facilitate the comparison with our results, we present in table 2
reference values for the evaluation of squared finite remainders at each power
of ${N_{f}}$, as defined in eqs. 24 and 25, at the randomly chosen phase-space
point

$\displaystyle\begin{split}s_{12}&=1.322500000,\quad\,\,s_{23}=-0.994109498,\quad\,\,s_{34}=0.264471591,\\\
s_{45}&=0.267126049,\quad\,\,s_{15}=-0.883795230,\quad\,\,\mu=1.\end{split}$
(53)

The results are obtained with our C++ library.

## References

  * (1) D0 collaboration, _Measurement of angular correlations of jets at $\sqrt{s}=1.96$ TeV and determination of the strong coupling at high momentum transfers_, _Phys. Lett. B_ 718 (2012) 56 [1207.4957].
  * (2) CMS collaboration, _Measurement of the Ratio of the Inclusive 3-Jet Cross Section to the Inclusive 2-Jet Cross Section in pp Collisions at $\sqrt{s}$ = 7 TeV and First Determination of the Strong Coupling Constant in the TeV Range_, _Eur. Phys. J. C_ 73 (2013) 2604 [1304.7498].
  * (3) CMS collaboration, _Measurement of the inclusive 3-jet production differential cross section in proton–proton collisions at 7 TeV and determination of the strong coupling constant in the TeV range_ , _Eur. Phys. J. C_ 75 (2015) 186 [1412.1633].
  * (4) ATLAS collaboration, _Measurement of transverse energy-energy correlations in multi-jet events in $pp$ collisions at $\sqrt{s}=7$ TeV using the ATLAS detector and determination of the strong coupling constant $\alpha_{\mathrm{s}}(m_{Z})$_, _Phys. Lett. B_ 750 (2015) 427 [1508.01579].
  * (5) A. Gehrmann-De Ridder, T. Gehrmann, E.W.N. Glover and J. Pires, _Second order QCD corrections to jet production at hadron colliders: the all-gluon contribution_ , _Phys. Rev. Lett._ 110 (2013) 162003 [1301.7310].
  * (6) J. Currie, E.W.N. Glover and J. Pires, _Next-to-Next-to Leading Order QCD Predictions for Single Jet Inclusive Production at the LHC_ , _Phys. Rev. Lett._ 118 (2017) 072002 [1611.01460].
  * (7) J. Currie, A. Gehrmann-De Ridder, T. Gehrmann, E.W.N. Glover, A. Huss and J. Pires, _Precise predictions for dijet production at the LHC_ , _Phys. Rev. Lett._ 119 (2017) 152001 [1705.10271].
  * (8) M. Czakon, A. van Hameren, A. Mitov and R. Poncelet, _Single-jet inclusive rates with exact color at $\mathcal{O}$ (${\alpha}_{s}^{4}$)_, _JHEP_ 10 (2019) 262 [1907.12911].
  * (9) J. Currie, A. Gehrmann-De Ridder, E.W.N. Glover and J. Pires, _NNLO QCD corrections to jet production at hadron colliders from gluon scattering_ , _JHEP_ 01 (2014) 110 [1310.3993].
  * (10) Z. Nagy, _Three jet cross-sections in hadron hadron collisions at next-to-leading order_ , _Phys. Rev. Lett._ 88 (2002) 122003 [hep-ph/0110315].
  * (11) Z. Nagy, _Next-to-leading order calculation of three jet observables in hadron hadron collision_ , _Phys. Rev. D_ 68 (2003) 094002 [hep-ph/0307268].
  * (12) W.B. Kilgore and W.T. Giele, _Next-to-leading order gluonic three jet production at hadron colliders_ , _Phys. Rev. D_ 55 (1997) 7183 [hep-ph/9610433].
  * (13) W.B. Kilgore and W.T. Giele, _A Next-to-leading order calculation of hadronic three jet production_ , in _30th International Conference on High-Energy Physics_ , pp. 502–505, 7, 2000 [hep-ph/0009193].
  * (14) R. Frederix, S. Frixione, V. Hirschi, D. Pagani, H.S. Shao and M. Zaro, _The automation of next-to-leading order electroweak calculations_ , _JHEP_ 07 (2018) 185 [1804.10017].
  * (15) M. Reyer, M. Schönherr and S. Schumann, _Full NLO corrections to 3-jet production and $\mathbf{R_{32}}$ at the LHC_, _Eur. Phys. J. C_ 79 (2019) 321 [1902.01763].
  * (16) S. Badger, H. Frellesvig and Y. Zhang, _A Two-Loop Five-Gluon Helicity Amplitude in QCD_ , _JHEP_ 12 (2013) 045 [1310.1051].
  * (17) 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].
  * (18) 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].
  * (19) 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].
  * (20) 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].
  * (21) 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].
  * (22) 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].
  * (23) 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].
  * (24) 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].
  * (25) G. De Laurentis and D. Maître, _Two-Loop Five-Parton Leading-Colour Finite Remainders in the Spinor-Helicity Formalism_ , _JHEP_ 02 (2021) 016 [2010.14525].
  * (26) S. Abreu, L.J. Dixon, E. Herrmann, B. Page and M. Zeng, _The two-loop five-point amplitude in $\mathcal{N}=4$ super-Yang-Mills theory_, _Phys. Rev. Lett._ 122 (2019) 121603 [1812.08941].
  * (27) 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 (2019) 121602 [1812.11057].
  * (28) S. Caron-Huot, D. Chicherin, J. Henn, Y. Zhang and S. Zoia, _Multi-Regge Limit of the Two-Loop Five-Point Amplitudes in $\mathcal{N}=4$ Super Yang-Mills and $\mathcal{N}=8$ Supergravity_, _JHEP_ 10 (2020) 188 [2003.03120].
  * (29) 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].
  * (30) C.F. Berger, Z. Bern, L.J. Dixon, F. Febres Cordero, D. Forde, T. Gleisberg et al., _Next-to-Leading Order QCD Predictions for W+3-Jet Distributions at Hadron Colliders_ , _Phys. Rev. D_ 80 (2009) 074036 [0907.1984].
  * (31) H. Ita and K. Ozeren, _Colour Decompositions of Multi-quark One-loop QCD Amplitudes_ , _JHEP_ 02 (2012) 118 [1111.4193].
  * (32) S. Abreu, B. Page, E. Pascual and V. Sotnikov, _Leading-Color Two-Loop QCD Corrections for Three-Photon Production at Hadron Colliders_ , _JHEP_ 21 (2020) 078 [2010.15834].
  * (33) 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.
  * (34) 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].
  * (35) S. Kallweit, V. Sotnikov and M. Wiesemann, _Triphoton production at hadron colliders in NNLO QCD_ , _Phys. Lett. B_ 812 (2021) 136013 [2010.04681].
  * (36) 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$_, 2102.01820.
  * (37) S. Badger, H.B. Hartanto and S. Zoia, _Two-loop QCD corrections to $Wb\bar{b}$ production at hadron colliders_, 2102.02516.
  * (38) C.G. Papadopoulos, D. Tommasini and C. Wever, _The Pentabox Master Integrals with the Simplified Differential Equations approach_ , _JHEP_ 04 (2016) 078 [1511.09404].
  * (39) 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].
  * (40) T. Gehrmann, J. Henn and N. Lo Presti, _Pentagon functions for massless planar scattering amplitudes_ , _JHEP_ 10 (2018) 103 [1807.09812].
  * (41) D. Chicherin and V. Sotnikov, _Pentagon Functions for Scattering of Five Massless Particles_ , _JHEP_ 12 (2020) 167 [2009.07803].
  * (42) A. von Manteuffel and R.M. Schabinger, _A novel approach to integration by parts reduction_ , _Phys. Lett. B_ 744 (2015) 101 [1406.4513].
  * (43) T. Peraro, _Scattering amplitudes over finite fields and multivariate functional reconstruction_ , _JHEP_ 12 (2016) 030 [1608.01902].
  * (44) 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_ , _Computer Physics Communications_ 267 (2021) 108069 [2009.11957].
  * (45) https://gitlab.com/five-point-amplitudes/FivePointAmplitudes-cpp.git.
  * (46) S. Weinzierl, _Does one need the O(epsilon)- and O(epsilon^2)-terms of one-loop amplitudes in an NNLO calculation ?_ , _Phys. Rev. D_ 84 (2011) 074007 [1107.5131].
  * (47) S. Catani, _The Singular behavior of QCD amplitudes at two loop order_ , _Phys. Lett. B_ 427 (1998) 161 [hep-ph/9802439].
  * (48) T. Becher and M. Neubert, _Infrared singularities of scattering amplitudes in perturbative QCD_ , _Phys. Rev. Lett._ 102 (2009) 162001 [0901.0722].
  * (49) E. Gardi and L. Magnea, _Factorization constraints for soft anomalous dimensions in QCD scattering amplitudes_ , _JHEP_ 03 (2009) 079 [0901.1091].
  * (50) A. Broggio, C. Gnendiger, A. Signer, D. Stöckinger and A. Visconti, _SCET approach to regularization-scheme dependence of QCD amplitudes_ , _JHEP_ 01 (2016) 078 [1506.05301].
  * (51) H. Ita, _Two-loop Integrand Decomposition into Master Integrals and Surface Terms_ , _Phys. Rev. D_ 94 (2016) 116015 [1510.05626].
  * (52) 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].
  * (53) N. Arkani-Hamed, J.L. Bourjaily, F. Cachazo and J. Trnka, _Local Integrals for Planar Scattering Amplitudes_ , _JHEP_ 06 (2012) 125 [1012.6032].
  * (54) B. Ruijl, T. Ueda and J. Vermaseren, _FORM version 4.2_ , 1707.06453.
  * (55) J. Kuipers, T. Ueda and J. Vermaseren, _Code Optimization in FORM_ , _Comput. Phys. Commun._ 189 (2015) 1 [1310.7007].
  * (56) D. Chicherin, J. Henn and V. Mitev, _Bootstrapping pentagon functions_ , _JHEP_ 05 (2018) 164 [1712.09610].
  * (57) Sherpa collaboration, _Event Generation with Sherpa 2.2_ , _SciPost Phys._ 7 (2019) 034 [1905.09127].
  * (58) ATLAS collaboration, _Measurement of three-jet production cross-sections in $pp$ collisions at 7 TeV centre-of-mass energy using the ATLAS detector_, _Eur. Phys. J. C_ 75 (2015) 228 [1411.1855].
  * (59) M. Cacciari, G.P. Salam and G. Soyez, _FastJet User Manual_ , _Eur. Phys. J. C_ 72 (2012) 1896 [1111.6097].
  * (60) S. Badger, C. Brønnum-Hansen, H.B. Hartanto, R. Moodie, T. Peraro and S. Zoia. private correspondence, 2021.
  * (61) C.F. Berger, Z. Bern, L.J. Dixon, F. Febres Cordero, D. Forde, H. Ita et al., _An Automated Implementation of On-Shell Methods for One-Loop Amplitudes_ , _Phys. Rev._ D78 (2008) 036003 [0803.4180].