<figure><img src="image_1.png" /></figure>

# Leading-color two-loop QCD corrections for three-jet production at hadron colliders

S. Abreu, F. Febres Cordero, H. Ita, B. Page, V. Sotnikov

► To cite this version:

S. Abreu, F. Febres Cordero, H. Ita, B. Page, V. Sotnikov. Leading-color two-loop QCD corrections for three-jet production at hadron colliders. Journal of High Energy Physics, 2021, 2021 (07), pp.095. 10.1007/JHEP07(2021)095. hal-03171417

HAL Id: hal-03171417

https://hal.science/hal-03171417v1

Submitted on 10 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# Leading-color two-loop QCD corrections for three-jet production at hadron colliders

S. Abreu,$^{a,b}$ F. Febres Cordero,$^{c}$ H. Ita,$^{d}$ B. Page$^{e}$ and V. Sotnikov$^{f}$

$^{a}$Theoretical Physics Department, CERN,
1211 Geneva 23, Switzerland

$^{b}$Mani L. Bhaumik Institute for Theoretical Physics,
Department of Physics and Astronomy, UCLA,
Los Angeles, CA 90095, U.S.A.

$^{c}$Physics Department, Florida State University,
77 Chieftan Way, Tallahassee, FL 32306, U.S.A.

$^{d}$Physikalisches Institut, Albert-Ludwigs-Universität Freiburg,
Hermann-Herder-Str. 3, D-79104 Freiburg, Germany

$^{e}$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, France

$^{f}$Max Planck Institute for Physics (Werner Heisenberg Institute),
D-80805 Munich, Germany

E-mail: samuel.abreu@cern.ch, ffebres@hep.fsu.edu,
harald.ita@physik.uni-freiburg.de, bpage@ipht.fr, sotnikov@mpp.mpg.de

**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.

**KEYWORDS:** NLO Computations, QCD Phenomenology

ARXIV ePRINT: 2102.13609
# Contents

1 **Introduction** 1

2 **Notation and conventions** 3

2.1 Helicity amplitudes 3

2.2 Finite remainders 6

2.3 Squared finite remainders 7

3 **Physical region helicity amplitudes** 8

3.1 Partial amplitudes in the physical region 8

3.2 From partial amplitudes to amplitudes 11

4 **Results and validation** 13

4.1 Analytic results 13

4.2 Numerical evaluation 14

4.3 Validation 17

5 **Conclusion** 19

A **Infrared subtraction** 19

B **Reference evaluations of squared finite remainders** 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 [1–4]. 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 \to 3$ contributions, which together with the *real-virtual* $2 \to 4$ and the *double-real* $2 \to 5$ contributions complete the $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 [5–7] and later subleading-color terms were incorporated [8] (see also [9]). Three-jet production is a more complex process. Results are currently known at NLO in QCD [10, 11] (see also [12, 13]) and also including electroweak corrections [14, 15]. 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 [16–20] and then in analytic form [21–25]. Progress in including subleading-color contributions has been made for related supersymmetric theories [26–28] as well as for special external gluon states [29]. In this work we present the leading-color double-virtual contributions to the scattering cross section, starting from the known analytic amplitudes [24] 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.<fnref n="1" />

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 \to 3$ processes with phenomenological applications [32–37]. We employ the known pure two-loop master integrals [21, 26, 38, 39], which can be written in terms of a special set of multi-valued transcendental functions — the so-called pentagon functions [40, 41]. In ref. [41], 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 [24]. 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 [18, 42, 43]. These evaluations are performed with the program CARAVEL [44]. 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 [45]. 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 provided as supplementary material of this work and showcase the efficiency and numerical stability of

<footnote n="1">¹The subleading-color contributions that we neglect are expected to have a minor impact in phenomenological applications [8, 9, 30, 31].</footnote>
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,

$$ g_{p_1}^{h_1} g_{p_2}^{h_2} \rightarrow g_{p_3}^{h_3} g_{p_4}^{h_4} g_{p_5}^{h_5}, \qquad (2.1) $$

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

$$ \bar{q}_{p_1}^{h_1} q_{p_2}^{h_2} \rightarrow g_{p_3}^{h_3} g_{p_4}^{h_4} g_{p_5}^{h_5}, \quad \bar{q}_{p_1}^{h_1} g_{p_2}^{h_2} \rightarrow \bar{q}_{p_3}^{h_3} g_{p_4}^{h_4} g_{p_5}^{h_5}, \quad g_{p_1}^{h_1} g_{p_2}^{h_2} \rightarrow q_{p_3}^{h_3} \bar{q}_{p_4}^{h_4} g_{p_5}^{h_5}, \qquad (2.2) $$

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

$$ \begin{aligned} \bar{q}_{p_1}^{h_1} q_{p_2}^{h_2} &\rightarrow Q_{p_3}^{h_3} \bar{Q}_{p_4}^{h_4} g_{p_5}^{h_5}, & \bar{q}_{p_1}^{h_1} Q_{p_2}^{h_2} &\rightarrow Q_{p_3}^{h_3} \bar{q}_{p_4}^{h_4} g_{p_5}^{h_5}, \\ \bar{q}_{p_1}^{h_1} \bar{Q}_{p_2}^{h_2} &\rightarrow \bar{q}_{p_3}^{h_3} \bar{Q}_{p_4}^{h_4} g_{p_5}^{h_5}, & \bar{q}_{p_1}^{h_1} g_{p_2}^{h_2} &\rightarrow \bar{q}_{p_3}^{h_3} Q_{p_4}^{h_4} \bar{Q}_{p_5}^{h_5}. \end{aligned} \qquad (2.3) $$

In eqs. (2.1) to (2.3), we identify each particle by its type (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 eq. (2.3), 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

$$ \begin{aligned} \bar{q}_{p_1}^{h_1} q_{p_2}^{h_2} &\rightarrow q_{p_3}^{h_3} \bar{q}_{p_4}^{h_4} g_{p_5}^{h_5} = (\bar{q}_{p_1}^{h_1} q_{p_2}^{h_2} \rightarrow Q_{p_3}^{h_3} \bar{Q}_{p_4}^{h_4} g_{p_5}^{h_5}) - (\bar{q}_{p_1}^{h_1} Q_{p_2}^{h_2} \rightarrow Q_{p_3}^{h_3} \bar{q}_{p_4}^{h_4} g_{p_5}^{h_5}), \\ \bar{q}_{p_1}^{h_1} \bar{q}_{p_2}^{h_2} &\rightarrow \bar{q}_{p_3}^{h_3} \bar{q}_{p_4}^{h_4} g_{p_5}^{h_5} = (\bar{q}_{p_1}^{h_1} \bar{Q}_{p_2}^{h_2} \rightarrow \bar{q}_{p_3}^{h_3} \bar{Q}_{p_4}^{h_4} g_{p_5}^{h_5}) - (\bar{q}_{p_1}^{h_1} \bar{Q}_{p_2}^{h_2} \rightarrow \bar{Q}_{p_3}^{h_3} \bar{q}_{p_4}^{h_4} g_{p_5}^{h_5}), \\ \bar{q}_{p_1}^{h_1} g_{p_2}^{h_2} &\rightarrow \bar{q}_{p_3}^{h_3} q_{p_4}^{h_4} \bar{q}_{p_5}^{h_5} = (\bar{q}_{p_1}^{h_1} g_{p_2}^{h_2} \rightarrow \bar{q}_{p_3}^{h_3} Q_{p_4}^{h_4} \bar{Q}_{p_5}^{h_5}) - (\bar{q}_{p_1}^{h_1} g_{p_2}^{h_2} \rightarrow \bar{Q}_{p_3}^{h_3} Q_{p_4}^{h_4} \bar{q}_{p_5}^{h_5}). \end{aligned} \qquad (2.4) $$

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

$$ \bar{q}_{p_1}^{h_1} \bar{Q}_{p_2}^{h_2} \rightarrow \bar{Q}_{p_3}^{h_3} \bar{q}_{p_4}^{h_4} g_{p_5}^{h_5}, \quad \bar{q}_{p_1}^{h_1} g_{p_2}^{h_2} \rightarrow \bar{Q}_{p_3}^{h_3} Q_{p_4}^{h_4} \bar{q}_{p_5}^{h_5}. \qquad (2.5) $$

We remark that, when labelling channels, such as in eqs. (2.1) to (2.5), 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. [20, 24], 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 $h$ and the set of momenta by a vector $p$. The
amplitudes $\mathcal{M}(h, 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} = (g_s^0)^3 \left( \mathcal{M}^{(0)} + \frac{\alpha_s^0}{2\pi} \mathcal{M}^{(1)} + \left( \frac{\alpha_s^0}{2\pi} \right)^2 \mathcal{M}^{(2)} + \mathcal{O}((\alpha_s^0)^3) \right), \quad (2.6) $$

where we drop the dependence on $h$ and $p$ for compactness. The order of $\alpha_s^0$ 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. [20].

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)}(h, p) = \left(\frac{N_c}{2}\right)^i \sum_{\sigma \in \Sigma} \mathcal{C}(\sigma) \Phi(\sigma(h), \sigma(p)) \mathcal{A}^{(i)}(\sigma(h), \sigma(p)), \quad (2.7) $$

where $N_c$ is the number of colors, $\mathcal{C}(\sigma)$ 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}(\sigma)$ 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) := p_{\sigma_i}$ and $\sigma(h_i) := h_{\sigma_i}$. The $\Phi(\sigma(h), \sigma(p))$ are normalization factors, and the $\mathcal{A}^{(i)}(\sigma(h), \sigma(p))$ 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(\sigma(h), \sigma(p))$ are defined so that $\mathcal{A}^{(0)}(\sigma(h), \sigma(p)) \equiv 1$ whenever the tree-level amplitude is non-vanishing, and for the other cases we choose<fnref n="2" />

$$ \Phi(h = \{+1, +1, +1, +1, +1\}, p) := \frac{[21]^2 \langle 12 \rangle^3 \langle 13 \rangle}{\langle 14 \rangle^2 \langle 15 \rangle^2 \langle 23 \rangle^3}, \quad (2.8a) $$

$$ \Phi(h = \{-1, +1, +1, +1, +1\}, p) := \frac{[21] \langle 12 \rangle^4 \langle 13 \rangle^3}{\langle 14 \rangle^2 \langle 15 \rangle^2 \langle 23 \rangle^5}, \quad (2.8b) $$

$$ \Phi(h = \{+\frac{1}{2}, -\frac{1}{2}, +1, +1, +1\}, p) := \frac{[31] \langle 12 \rangle^3 \langle 13 \rangle}{\langle 14 \rangle^2 \langle 15 \rangle^2 \langle 23 \rangle^2}, \quad (2.8c) $$

where we set $p = \{p_1, p_2, p_3, p_4, p_5\}$. With this choice, the partial amplitudes $\mathcal{A}(\sigma(h), \sigma(p))$ only depend on the kinematics through the Mandelstam invariants

$$ \begin{aligned} s_{12} &= (p_1 + p_2)^2, & s_{23} &= (p_2 + p_3)^2, & s_{34} &= (p_3 + p_4)^2, \quad (2.9) \\ s_{45} &= (p_4 + p_5)^2, & s_{15} &= (p_1 + p_5)^2, \end{aligned} $$

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

$$ \mathrm{tr}_5 := 4i\epsilon(p_1, p_2, p_3, p_4), \quad (2.10) $$

<footnote n="2">²Here we use the same notation for spinor brackets as in [44].</footnote>
<figure><img src="image_4.png" /><figcaption>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).</figcaption></figure>

where $\epsilon(\cdot, \cdot, \cdot, \cdot)$ is the fully anti-symmetric Levi-Civita tensor, whose sign is fixed such that $\mathrm{Im}(\mathrm{tr}_5) > 0$ corresponds to the choice of the basis helicity amplitudes in ref. [24]. 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. (2.10) by $\sigma(\mathrm{tr}_5) = \mathrm{sgn}(\sigma) \mathrm{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)}(h, p)$ can be further decomposed in powers of $N_f$. We write

$$ \mathcal{A}^{(i)}(h, p) = \sum_{j=0}^{i} \left(\frac{N_f}{N_c}\right)^j \mathcal{A}^{(i)[j]}(h, p). \qquad (2.11) $$

In figure 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.

The renormalized amplitudes $M_R$ are obtained from the bare amplitudes $M$ by considering an expansion in powers of the renormalized coupling $\alpha_s(\mu)$. In the $\overline{\mathrm{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), \qquad (2.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 β-function,

$$ 
\begin{aligned}
\beta_0 &= \frac{11C_A - 4T_F N_f}{6} = \frac{11N_c - 2N_f}{6}, \\
\beta_1 &= \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}),
\end{aligned} 
\quad (2.13) $$

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

$$ 
\begin{aligned}
\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{aligned} 
\quad (2.14) $$

## 2.2 Finite remainders

Amplitudes contribute to physical observables only through so-called finite remainders (see e.g. [46]), which we denote $\mathcal{R}(h, p)$. We define a remainder for each partial amplitude $\mathcal{A}(h, p)$. 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 [47–49]. 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), \quad (2.15) $$

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

$$ 
\begin{aligned}
\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)} \mathcal{A}_R^{(0)} + \mathcal{O}(\epsilon),
\end{aligned} 
\quad (2.16) $$

The operators $\mathbf{I}^{(i)}$ are channel specific, and we collect them in appendix A for all the channels in eqs. (2.1) to (2.3). In these expressions, we extend the expansion of $\mathbf{I}^{(1)}\mathcal{A}_R^{(1)}$ and $\mathbf{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. [50]. Finally, the finite remainders can also be decomposed into powers of $N_f$, similarly to eq. (2.11),

$$ \mathcal{R}^{(i)}(h, p) = \sum_{j=0}^{i} \left(\frac{N_f}{N_c}\right)^j \mathcal{R}^{(i)[j]}(h, p). \quad (2.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 $\mathbf{p} \to \mathbf{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}_{\mathcal{M}}^{(i)}(\mathbf{h}, \mathbf{p})$ is trivially obtained from the decomposition of the color-dressed amplitude $\mathcal{M}^{(i)}(\mathbf{h}, \mathbf{p})$ into partial amplitudes, see eq. (2.7):

$$ \begin{aligned} \mathcal{R}_{\mathcal{M}}(\mathbf{h}, \mathbf{p}) &= g_s^3 \left( \mathcal{R}_{\mathcal{M}}^{(0)}(\mathbf{h}, \mathbf{p}) + \frac{\alpha_s}{2\pi} \mathcal{R}_{\mathcal{M}}^{(1)}(\mathbf{h}, \mathbf{p}) + \left(\frac{\alpha_s}{2\pi}\right)^2 \mathcal{R}_{\mathcal{M}}^{(2)}(\mathbf{h}, \mathbf{p}) + \mathcal{O}(\alpha_s^3) \right), \\ \mathcal{R}_{\mathcal{M}}^{(i)}(\mathbf{h}, \mathbf{p}) &= \left(\frac{N_c}{2}\right)^i \sum_{\sigma \in \Sigma} \mathcal{C}(\sigma) \Phi(\sigma(\mathbf{h}), \sigma(\mathbf{p})) \mathcal{R}^{(i)}(\sigma(\mathbf{h}), \sigma(\mathbf{p})). \end{aligned} \quad (2.18) $$

Whenever the distinction is important, we will call the $\mathcal{R}^{(i)}(\sigma(\mathbf{h}), \sigma(\mathbf{p}))$ 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}(\sigma)$ 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). \quad (2.19) $$

In terms of the quantities in eq. (2.18), $H$ is defined as

$$ H = \frac{1}{\mathcal{N}_H} \sum_{\mathbf{h}} \mathcal{R}_{\mathcal{M}}^{\dagger}(\mathbf{h}, \mathbf{p}) \mathcal{R}_{\mathcal{M}}(\mathbf{h}, \mathbf{p}), \quad (2.20) $$

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

$$ \mathcal{N}_H := (4\pi\alpha_s)^3 \sum_{\mathbf{h}} |\mathcal{R}_{\mathcal{M}}^{(0)}(\mathbf{h}, \mathbf{p})|^2. \quad (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_{\mathbf{h}} |\mathcal{R}_{\mathcal{M}}^{(0)}(\mathbf{h}, \mathbf{p})|^2} \sum_{\mathbf{h}} 2 \operatorname{Re} [\mathcal{R}_{\mathcal{M}}^{(0)\dagger}(\mathbf{h}, \mathbf{p}) \mathcal{R}_{\mathcal{M}}^{(1)}(\mathbf{h}, \mathbf{p})], \quad (2.22) $$

and at order $\alpha_s^2$

$$ H^{(2)} = \frac{1}{\sum_{\mathbf{h}} |\mathcal{R}_{\mathcal{M}}^{(0)}(\mathbf{h}, \mathbf{p})|^2} \left( \sum_{\mathbf{h}} |\mathcal{R}_{\mathcal{M}}^{(1)}(\mathbf{h}, \mathbf{p})|^2 + \sum_{\mathbf{h}} 2 \operatorname{Re} [\mathcal{R}_{\mathcal{M}}^{(0)\dagger}(\mathbf{h}, \mathbf{p}) \mathcal{R}_{\mathcal{M}}^{(2)}(\mathbf{h}, \mathbf{p})] \right), \quad (2.23) $$
<figure><img src="image_5.png" /><figcaption>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.</figcaption></figure>

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

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). \qquad (2.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). \qquad (2.25)$$

For most channels the subleading-color corrections are suppressed by a factor of $1/N_c^2$. However, for the channels in eq. (2.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. [24]. We recall the comment below eq. (2.5) regarding the change of conventions for denoting each channel between this paper and ref. [24].

### 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. [24] in such region. Here we consider only the basis set of partial helicity amplitudes of ref. [24], and we will return to their permutations required to assemble eq. (2.23) in section 3.2. The amplitudes of ref. [24] were computed in the (non-physical) Euclidean region where

$$s_{12}, s_{23}, s_{34}, s_{45}, s_{15} < 0, \quad \text{tr}_5^2 > 0. \qquad (3.1)$$
It is our aim to construct results in the physical region associated with the channels in eqs. (2.1) to (2.5), 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, \quad s_{23}, s_{15} < 0, \quad \text{tr}_5^2 < 0. \qquad (3.2)$$

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 $p_\varphi$ and a phase-space point in the Euclidean region as $p_E$. For region-agnostic statements we simply denote a phase-space point as $p$.

In order to explain how we convert the results of ref. [24] to the kinematic region in eq. (3.2), 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)}(\mathbf{h}, \mathbf{p}; \epsilon) = \sum_k c_k^{(i)}(\mathbf{h}, \mathbf{p}; \epsilon) m_k^{(i)}(\mathbf{p}; \epsilon), \quad i = 1, 2, \qquad (3.3)$$

where we make the dependence on the dimensional regulator $\epsilon$ explicit. This decomposition was obtained using the two-loop numerical unitarity approach [18, 20, 51, 52]. At this stage, the master integrals were replaced by their expansion in terms of Euclidean-region pentagon functions.<fnref n="3" /> Denoting the vector space spanned by relevant pentagon-function monomials in this region as $\tilde{f}(\mathbf{p}_E) = \{f_l(\mathbf{p}_E)\}_{l \in B}$, with $B$ denoting the set of multi-indices $l$, we write

$$m_k^{(i)}(\mathbf{p}_E; \epsilon) = \sum_{k=0}^{4} \sum_{l \in B} \epsilon^k d_{k,l}^{(i)} f_l(\mathbf{p}_E) + \mathcal{O}(\epsilon^5). \qquad (3.4)$$

In this decomposition, we assume that we have a pure basis of master integrals [53] at one and two loops, in which case the coefficients $d_{k,l}^{(i)}$ 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 (3.4) into eq. (3.3) for both one- and two-loop amplitudes, and then using these expressions in the definitions for the finite remainders given in eq. (2.16), we obtain a decomposition of the two-loop remainders into pentagon functions:

$$\mathcal{R}^{(i)[j]}(\mathbf{h}, \mathbf{p}_E) = \sum_{k \in B} r_k^{(i)[j]}(\mathbf{h}, \mathbf{p}_E) f_k(\mathbf{p}_E), \qquad (3.5)$$

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_k^{(i)[j]}(\mathbf{h}, \mathbf{p})$ was the main result of ref. [24]. Let us finish this summary by commenting on the structure of the coefficients $r_k^{(i)[j]}(\mathbf{h}, \mathbf{p})$. They can be written as

$$r_k(\mathbf{h}, \mathbf{p}) = r_k^+(\mathbf{h}, \vec{\mathbf{s}}) + r_k^-(\mathbf{h}, \vec{\mathbf{s}}) \text{tr}_5, \qquad (3.6)$$

<footnote n="3">³These 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. [40] for more details on their definition.</footnote>
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_k^{\pm}(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(h, 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}(h, p) = \{v_1(h, p), \dots\}$ such that all $r_k(h, p)$ live in the $\mathbb{Q}$-linear span of $\vec{v}(h, p)$. We can then write

$$ \mathcal{R}^{(i)[j]}(h, p_E) = \vec{v}^{(i)[j]}(h, p_E) \cdot M_E^{(i)[j]}(h) \cdot \vec{f}(p_E), \qquad (3.7) $$

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

To understand our approach to the analytic continuation of eq. (3.7), let us begin by considering what would happen if we had chosen a different basis of pentagon functions $\vec{f}'(p_E)$, such that $\vec{f}'(p_E) = N \cdot \vec{f}(p_E)$ where $N$ is a matrix of rational numbers. It is clear that to write eq. (3.7) in this new basis of functions we do not need to determine the basis $\vec{v}(h, p_E)$ of rational functions. The remainders would have exactly the same form, but with the matrix $M_E$ replaced by $M_E^{(i)[j]} = M_E^{(i)[j]} \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(|s_{12}|)$. Let us now continue this function to the physical region in eq. (3.2), where $s_{12} > 0$. We find

$$ \log(-s_{12}) = \log(|s_{12}|) - i\pi \qquad (3.8) $$

where the sign is determined by the usual Feynman $i\epsilon$-prescription of the propagators. We see that the analytic continuation into the physical region introduces a new pentagon “function”, the transcendental constant $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 $i\pi$ is trivially related to the coefficient of $\log(-s_{12})$ when analytically continuing eq. (3.4) or eq. (3.5), or more precisely that it is in the $\mathbb{Q}$-linear span of $\vec{v}(h, p_\phi)$ (we recall that $p_\phi$ 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. (3.5) might increase, the coefficients of the new pentagon functions are not independent from the coefficients of the pentagon functions in the Euclidean region.<fnref n="4" /> 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. (3.7).

To analytically continue the remainders computed in ref. [24], we use the set of pentagon functions defined in ref. [41], and we denote the vector space of their monomials as $\vec{g}(p_\phi)$. These are valid precisely in the kinematic region of eq. (3.2), and are related to the set $\vec{f}(p_E)$ used in ref. [24] by a change of basis and by analytic continuation. Given the

<footnote n="4">⁴In general, the two will be related by the (integer) winding numbers that naturally appear when analytically continuing multi-valued functions.</footnote>
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}_{\varphi}^{(i)[j]}(\mathbf{h}, \mathbf{p}_{\varphi}; \mu) = \bar{v}^{(i)[j]}(\mathbf{h}, \mathbf{p}_{\varphi}) \cdot M_{\varphi}^{(i)[j]}(\mathbf{h}) \cdot \bar{g}(\mathbf{p}_{\varphi}/\mu), \quad (3.9) $$

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 $\bar{g}(\mathbf{p}_{\varphi})$. We discuss the generalization to other configurations in section 3.2. The $M_{\varphi}^{(i)[j]}(\mathbf{h})$ can be determined from a small number of numerical evaluations of the amplitudes. To do this, we take equation eq. (3.9) as an ansatz for the form of the remainder when expressed in terms of the physical region pentagon functions $\bar{g}(\mathbf{p}_{\varphi})$, with the entries of the $M_{\varphi}^{(i)[j]}(\mathbf{h})$ as unknowns. We then constrain the unknowns by performing the analogous procedure of eqs. (3.4) and (3.5) at a sufficient number of phase-space points.<fnref n="5" /> Linear algebra then allows us to extract the values of $M_{\varphi}^{(i)[j]}(\mathbf{h})$.

In practice, we use the implementation of two-loop numerical unitarity in CARAVEL [44] to generate all the required numerical data. For efficiency reasons, the procedure is performed over finite fields and we construct the Q-valued $M_{\varphi}^{(i)[j]}(\mathbf{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. [24] to the physical region defined in eq. (3.2), corresponding to the channels in eqs. (2.1) to (2.5). Crucially, this linear algebra exercise would not allow one to compute the matrices $M_{\varphi}^{(i)[j]}(\mathbf{h})$ in the ansatz (3.9) 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. (2.7) for each helicity configuration. We start by summarizing this information for each of the channels listed in eqs. (2.1) to (2.5). We remind the reader of the comment on momenta conventions under eq. (2.5).

**Five gluon.** There is a single channel, see eq. (2.1). The set $\Sigma = S_5/Z_5$ contains 24 elements, and the color factor for $\sigma_0 = \{1, 2, 3, 4, 5\}$ is

$$ C(\sigma_0) := \mathrm{tr}\{T^{a_1}T^{a_2}T^{a_3}T^{a_4}T^{a_5}\}. \quad (3.10) $$

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

<footnote n="5">⁵This number is dictated by the number of independent functions in $\bar{v}^{(i)[j]}(\mathbf{h}, \mathbf{p})$. The most complicated case, the five-gluon remainder $\mathcal{R}_{\varphi}^{(2)[0]}(\{-1, +1, -1, +1, +1\}, \mathbf{p})$, requires 57 evaluations. We note that a single numerical evaluation gives more than one constraint.</footnote>
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.2). For the channel $\bar{q}_{p1}^{h_1} q_{p2}^{h_2} \to g_{p3}^{h_3} g_{p4}^{h_4} g_{p5}^{h_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

$$ C(\sigma_0) := (T^{a_3} T^{a_4} T^{a_5})_{i_1}^{\bar{i}_2}. \qquad (3.11) $$

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 eq. (2.3), but we also include the two channels in eq. (2.5) in order to consider two identical quark lines. For the channel $\bar{q}_{p1}^{h_1} q_{p2}^{h_2} \to Q_{p3}^{h_3} \bar{Q}_{p4}^{h_4} g_{p5}^{h_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

$$ C(\sigma_0) := (T^{a_5})_{i_3}^{\bar{i}_2} \delta_{i_1}^{\bar{i}_4}. \qquad (3.12) $$

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. (3.9). Clearly, the master integrals in the decomposition (3.3) 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. [41] 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_k^{(i)}(\sigma(\mathbf{p}); \epsilon) = \sum_{k=0}^{4} \sum_{l \in B_{\varphi}} \epsilon^k d_{k,l}^{(i)}(\sigma) g_l(\mathbf{p}) + \mathcal{O}(\epsilon^5), \qquad (3.13) $$

where $B_{\varphi}$ is the set of labels distinguishing the elements of $\vec{g}(\mathbf{p})$. It therefore follows that all permutations of the partial remainders contributing to the remainder in eq. (2.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(\mathbf{h})$ and on a permutation of the input point $\sigma(\mathbf{p})$ can be expressed in the form

$$ \mathcal{R}_{\varphi}^{(i)[j]}(\sigma(\mathbf{h}), \sigma(\mathbf{p}); \mu) = \vec{v}^{(i)[j]}(\sigma(\mathbf{h}), \sigma(\mathbf{p})) \cdot M_{\varphi}^{(i)[j]}(\sigma(\mathbf{h}), \sigma) \cdot \vec{g}(\mathbf{p}/\mu), \qquad (3.14) $$

where, as in eq. (3.9), we made the dependence on $\mu$ explicit. We recall that we can easily construct the bases $\vec{v}^{(i)[j]}(\sigma(\mathbf{p}), \sigma(\mathbf{h}))$, as the action of $\sigma$ on rational functions is trivial
to evaluate. The matrices $M_{\phi}^{(i)[j]}(\sigma(\mathbf{h}), \sigma)$ can then be determined by the same ansatz procedure described for the analytic continuation of the results of ref. [24]. 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. (2.1) to (2.5).

# 4 Results and validation

In this section we describe our main results, which are of two types. The first is a collection of files, provided as supplementary material of this work, containing analytic expressions for a generating set of (partial) remainders for each of the channels in eqs. (2.1) to (2.5). The second is a C++ library that allows for an efficient numerical evaluation of the finite remainders defined in eq. (3.9) for all required permutations, as well as the squared finite remainders defined in eqs. (2.22) and (2.23). Finally, we describe the checks we have performed on our results.

## 4.1 Analytic results

In the supplementary material we give all necessary ingredients to assemble analytic expressions for a generating set of two-loop partial remainders $\mathcal{R}_{\phi}^{(2)[j]}(\sigma(\mathbf{h}), \sigma(\mathbf{p}))$, for each of the channels in eqs. (2.1) to (2.5). From this generating set, all the permutations required in eq. (2.18) can be constructed by applying (a combination of) charge conjugation, parity transformations,<fnref n="6" /> and the $S_2 \times S_3$ permutations of initial and final state labels.<fnref n="7" />

The expressions we present in the supplementary material allow to explicitly construct the generating set of remainders in the decomposition of eq. (3.14). The basis of transcendental functions $\tilde{g}(\mathbf{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. [41]. The basis of rational functions $\tilde{v}^{(2)[j]}(\sigma(\mathbf{h}), \sigma(\mathbf{p}))$ is channel and power-of-$N_f$ specific but trivially dependent on $\sigma$ (more precisely, the action of $\sigma$ on $\tilde{v}^{(2)[j]}(\sigma(\mathbf{h}), \sigma(\mathbf{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 supplementary material of ref. [24] (we note again the comment below eq. (2.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 (s, tr5)}.`

We note that one can alternatively use the form of the same bases of rational functions determined in ref. [25], which are more compactly expressed in the spinor-helicity formalism. Finally, the matrices $M_{\phi}^{(2)[j]}(\sigma(\mathbf{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.

<footnote n="6">⁶Parity transformation of the analytic expressions are performed by changing the signs of all parity-odd pentagon functions and tr5.</footnote>

<footnote n="7">⁷One may wish to express the $S_2 \times 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.</footnote>
As described in section 3, these matrices constitute the main analytic result of this pa-
per. 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. (3.2).

## 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++ li-
brary [45]. The code performs the numerical evaluation of one- and two-loop (partial)
remainders for each channel, each helicity configuration and each power of N<sub>f</sub>. That is,
we evaluate each of the R<sup>(i)[j]</sup>(&sigma;<span style="text-decoration: overline">h</span>), &sigma;<span style="text-decoration: overline">p</span>) for i = 1, 2 and present each power of (N<sub>f</sub>/N<sub>c</sub>)<sup>j</sup>
separately. We stress that we have implemented all permutations, and not only the gener-
ating set for which we provide analytic expressions. For physical applications, one might
be interested in the squared finite remainders defined in eqs. (2.22) and (2.23). Our code
also outputs the numerical values for these quantities, once again presenting each power of
N<sub>f</sub>/N<sub>c</sub> separately, that is the H<sup>(i)[j]</sup> in eqs. (2.24) and (2.25).

We employ PentagonFunctions++ [41] for the numerical evaluation of the pentagon
functions. For the efficient evaluation of the rational functions, we generate optimized code
with FORM [54, 55]. In any of the helicity remainders $\mathcal{R}^{(2)[j]}(\sigma(\mathbf{h}), \sigma(\mathbf{p}))$ that we compute
numerically, the dominant contribution to the evaluation time is that of the pentagon func-
tions. 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. (2.22) and (2.23). We note
nevertheless that, compared to e.g. the case of triphoton production [32], the evaluation of
$H^{(2)}$ for five-gluon channels receives contributions from a large number of partial ampli-
tudes, and therefore the time spent in the evaluation of the other components in eq. (3.14)
becomes noticeable. We find that the typical double-precision evaluation time for any $H^{(2)}$
is $O(1s)$ on a single CPU, while the evaluation time for $H^{(1)}$ is negligible at $O(10ms)$.

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. (3.14), when a phase-space point gets close to the zero sets of the
rational functions' denominators. As we have already argued in refs. [23, 24, 32], these are
determined by a set of 25 quantities {$W_i(\vec{s})$}, with $i = 1, \dots, 25$, which are a subset of the
so-called alphabet associated with the five-point two-loop master integrals (see e.g. ref. [56]
for the definition of these quantities). This subset is rational in the Mandelstam variables
$\vec{s}$ defined in eq. (2.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}$,<fnref n="8" />

$$ \kappa(\vec{s}) := \frac{\min_i \{|W_i(\vec{s})|\}}{\max_i \{|W_i(\vec{s})|\}}. \qquad (4.1) $$

We introduce a threshold $\kappa_{thr}$, and if $\kappa(\vec{s}) > \kappa_{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} \approx \varepsilon_{\text{double}} \vec{s}, \qquad (4.2) $$

where $\varepsilon_{\text{double}} \approx 10^{-16}$ is the machine epsilon of double-precision floating point numbers.<fnref n="9" /> We estimate the accuracy of potentially unstable points as

$$ \Delta(\vec{s}) := \left| 1 - \frac{H(\vec{s}_{\delta})}{H(\vec{s})} \right|. \qquad (4.3) $$

We introduce another threshold value $\Delta_{thr}$, and if $\Delta(\vec{s}) < \Delta_{thr}$ the point is considered stable. If $\Delta(\vec{s}) > \Delta_{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_{thr}$ and $\Delta_{thr}$ will be discussed below.

To demonstrate the performance of our numerical C++ implementation, we employ SHERPA 2.2 [57] 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 $gg \to ggg$, $\bar{q}g \to \bar{q}gg$ and $\bar{q}q \to Q\bar{Q}g$, which are representative of all channels in eqs. (2.1) to (2.3). We use the same phase-space definition that was used in the 7 TeV ATLAS analysis [58], 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 [59]. We require the jets to be within the rapidity range $|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}, \qquad (4.4) $$

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 := - \log_{10} \left| 1 - \frac{H_{\text{double}}^{(2)}}{H_{\text{quad}}^{(2)}} \right|, \qquad (4.5) $$

<footnote n="8">⁸We attempted to further reduce the subset of letters to only include those that correspond to spurious singularities, as these had been identified in ref. [32] 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.</footnote>

<footnote n="9">⁹In 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$.</footnote>
<figure><img src="image_6.png" /><figcaption>Figure 3. Logarithmic distribution of correct digits (see eq. (4.5)) 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.</figcaption></figure>
<table><caption>**Table 1.** Characteristic performance of the rescue system on the samples of 110k points with $\kappa_{thr} = 10^{-4}$ and $\Delta_{thr} = 10^{-2}$. The first column shows the fraction of points with $\kappa(\vec{s}) < \kappa_{thr}$, which are recomputed in double precision at $\vec{s}_\delta$. The second column shows the fraction of points with $\Delta(\vec{s}) > \Delta_{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.</caption><thead><tr><th>Channels</th><th>% recomputed in double</th><th>% recomputed in quad</th><th>slowdown factor</th></tr></thead><tbody><tr><td>H<sup>(2)</sup><sub>gg</sub> → ggg</td><td>2.52</td><td>0.27</td><td>1.12</td></tr><tr><td>H<sup>(2)</sup><sub>qg</sub> → qgg</td><td>2.57</td><td>0.08</td><td>1.06</td></tr><tr><td>H<sup>(2)</sup><sub>qQ</sub> → qQg</td><td>2.52</td><td>0.01</td><td>1.03</td></tr></tbody></table>

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. (4.3). In figure 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_{thr} = 10^{-4}$ and $\Delta_{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_{thr} = 10^{-4}$ and $\Delta_{thr} = 10^{-2}$. The test in eq. (4.1) catches too many points, but crucially catches all unstable points as can be seen from figure 3. The second test, in eq. (4.3), 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_{thr} = 10^{-4}$ and $\Delta_{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. [25] reduces the fraction of unstable points, but we find that our implementation is already satisfactorily fast and stable for phenomenological applications.

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. (2.24) and (2.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. [24], which have undergone very stringent checks. Compared to ref. [24], our current results are written in terms of a
different basis of transcendental functions [41] which have also undergone several direct checks (e.g., comparison with the functions of ref. [40]) and indirect checks [32, 33, 35, 36]. Furthermore, we have verified that we reproduce the results independently obtained in ref. [19] for the numerical evaluation of five-gluon amplitudes in the physical region defined in eq. (3.2). Finally, we compared our results for the function $H_{gg→ggg}^{(2)[0]}$ in the purely gluonic channel against an independent computation [60] 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. (2.1) to (2.5) against the evaluations by the BlackHat library [61] and found full agreement.<fnref n="10" /> The color-dressed remainders $\mathcal{R}_M^{(1)}$ in eq. (2.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. (2.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

$$ 
\begin{aligned}
\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} + \mathbf{H}_{[n]} \right) \mathcal{R}^{(0)} + \frac{2\beta_0(2+\lambda)}{N_c} \mathcal{R}^{(1)}(\mu),
\end{aligned}
\qquad (4.6) $$

where we denote by $\lambda$ the power of $g_s^0$ in the tree-level amplitude (for five-parton amplitudes $\lambda = 3$, see eq. (2.6)) and the $\beta_i$ are the coefficients of the QCD $\beta$ function explicitly given in eq. (2.13). $\mathbf{H}_{[n]}$ is given explicitly in eqs. (A.6) and (A.7) 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. (4.6). 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.

<footnote n="10">¹⁰After adjusting the definitions of the leading-color squared matrix elements in BlackHat such that all subleading-color contributions are discarded.</footnote>
# 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. [24] to the physical region corresponding to three-jet production is performed by rewriting them in terms of the pentagon functions defined in ref. [41]. 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 [44]. We provide analytic expressions for a generating set of partial remainders. For all channels in eqs. (2.1) to (2.5), these allow one to evaluate the color-dressed remainders defined in eq. (2.18) as well as the squared remainders $H^{(1)}$ and $H^{(2)}$ defined in eqs. (2.24) and (2.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 $O(1s)$. 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.

# A Infrared subtraction

The divergences of renormalized two-loop amplitudes obey a universal structure [47–49],

$$ 
\begin{align} 
A_R^{(1)} &= \mathbf{I}_{[n]}^{(1)}(\epsilon) A_R^{(0)} + \mathcal{O}(\epsilon^0), \tag{A.1} \\ 
A_R^{(2)} &= \mathbf{I}_{[n]}^{(2)}(\epsilon) A_R^{(0)} + \mathbf{I}_{[n]}^{(1)}(\epsilon) A_R^{(1)} + \mathcal{O}(\epsilon^0), 
\end{align} 
$$
with the renormalized amplitudes $\mathcal{A}_R^{(i)}$ related to the bare amplitudes $\mathcal{A}^{(i)}$ as in eq. (2.14). For amplitudes in the leading-color approximation the operators $\mathbf{I}_{[n]}^{(1)}$ and $\mathbf{I}_{[n]}^{(2)}$ are diagonal in color space. The operator $\mathbf{I}_{[n]}^{(1)}$ is given by

$$ \mathbf{I}_{[n]}^{(1)}(\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} - i\epsilon\right)^{-\epsilon}, \qquad (A.2) $$

where $s_{i,j} = (p_i + p_j)^2$, with the indices defined cyclically, and we made explicit the $i\epsilon$ associated with Mandelstam variables that follows from the Feynman $i\epsilon$-prescription of the propagators. The index $a_i$ denotes a type of particle, i.e. $a_i \in \{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:

$$ \begin{aligned} \gamma_{g,g} &= \frac{1}{\epsilon^2} + \frac{1}{\epsilon} \frac{\beta_0}{N_c}, & \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_{g,q} &= \gamma_{g,\bar{q}} = \gamma_{g,Q} = \gamma_{g,\bar{Q}} = \frac{\gamma_{g,g} + \gamma_{q,Q}}{2}, & \gamma_{q,\bar{q}} &= \gamma_{Q,\bar{Q}} = 0. \end{aligned} \qquad (A.3) $$

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

$$ \mathbf{I}_{[n]}^{(2)}(\epsilon) = -\frac{1}{2}\mathbf{I}_{[n]}^{(1)}(\epsilon)\mathbf{I}_{[n]}^{(1)}(\epsilon) - \frac{2\beta_0}{N_c\epsilon}\mathbf{I}_{[n]}^{(1)}(\epsilon) + \frac{e^{-\gamma_E\epsilon}\Gamma(1-2\epsilon)}{\Gamma(1-\epsilon)}\left(\frac{2\beta_0}{N_c\epsilon} + K\right)\mathbf{I}_{[n]}^{(1)}(2\epsilon) \qquad (A.4) \\ + \frac{e^{\gamma_E\epsilon}}{\epsilon\Gamma(1-\epsilon)}\mathbf{H}_{[n]}, $$

where

$$ K = \frac{67}{9} - \frac{\pi^2}{3} - \frac{10}{9} \frac{N_f}{N_c}, \qquad (A.5) $$

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

$$ \mathbf{H}_{[n]} = \sum_{i=1}^{n} \left( \delta_{a_i,g} H_g + (\delta_{a_i,q} + \delta_{a_i,\bar{q}} + \delta_{a_i,Q} + \delta_{a_i,\bar{Q}}) H_q \right), \qquad (A.6) $$

with

$$ \begin{aligned} H_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{aligned} \qquad (A.7) $$

## B Reference evaluations of squared finite remainders

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. (2.24) and (2.25), at the randomly chosen phase-space point

$$ \begin{aligned} s_{12} &= 1.322500000, & s_{23} &= -0.994109498, & s_{34} &= 0.264471591, \\ s_{45} &= 0.267126049, & s_{15} &= -0.883795230, & \mu &= 1. \end{aligned} \qquad (B.1) $$

The results are obtained with our C++ library.
<table><caption>**Table 2.** Reference values for the evaluation of squared finite remainders at each power of $N_f$, as defined in eqs. (2.24) and (2.25) on the phase-space point given in eq. (B.1).</caption><thead><tr><th>Channel</th><th><i>H</i><sup>(1)</sup>[0]</th><th><i>H</i><sup>(1)</sup>[1]</th><th><i>H</i><sup>(2)</sup>[0]</th><th><i>H</i><sup>(2)</sup>[1]</th><th><i>H</i><sup>(2)</sup>[2]</th></tr></thead><tbody><tr><td>gg → ggg</td><td>16.135254222</td><td>0.19163044752</td><td>464.47846208</td><td>-58.116292408</td><td>0.60077232705</td></tr><tr><td><span style="text-decoration: overline;">q</span><i>q</i> → ggg</td><td>9.5879406141</td><td>-3.0604943308</td><td>184.44415807</td><td>-61.765802987</td><td>6.3615768297</td></tr><tr><td><span style="text-decoration: overline;">q</span>g → <span style="text-decoration: overline;">q</span>gg</td><td>26.908169290</td><td>-3.6373308269</td><td>867.25232363</td><td>-230.76277359</td><td>12.598811302</td></tr><tr><td>gg → <i>q</i><span style="text-decoration: overline;">q</span>g</td><td>24.495592766</td><td>-2.5939909248</td><td>745.87682394</td><td>-166.84486839</td><td>6.1899943330</td></tr><tr><td><span style="text-decoration: overline;">q</span><i>q</i> → <i>Q</i><span style="text-decoration: overline;">Q</span>g</td><td>10.460907919</td><td>-4.2060557725</td><td>212.42454564</td><td>-80.136400792</td><td>8.2094005806</td></tr><tr><td><span style="text-decoration: overline;">q</span><i>Q</i> → <i>Q</i><span style="text-decoration: overline;">q</span>g</td><td>27.104747640</td><td>-4.0829938180</td><td>705.58902507</td><td>-209.42216177</td><td>12.483148067</td></tr><tr><td><span style="text-decoration: overline;">q</span><span style="text-decoration: overline;">Q</span> → <span style="text-decoration: overline;">q</span><i>Q</i>g</td><td>42.313652168</td><td>-8.0064067852</td><td>1628.2933493</td><td>-562.78735847</td><td>44.198947852</td></tr><tr><td><span style="text-decoration: overline;">q</span>g → <span style="text-decoration: overline;">q</span><i>Q</i><span style="text-decoration: overline;">Q</span></td><td>28.068256507</td><td>-6.3593609865</td><td>935.81439233</td><td>-324.32790785</td><td>29.070926975</td></tr><tr><td><span style="text-decoration: overline;">q</span><i>q</i> → <i>q</i><span style="text-decoration: overline;">q</span>g</td><td>20.846053179</td><td>-4.1292696285</td><td>520.14108472</td><td>-160.80597165</td><td>10.876062192</td></tr><tr><td><span style="text-decoration: overline;">q</span><span style="text-decoration: overline;">q</span> → <span style="text-decoration: overline;">q</span><span style="text-decoration: overline;">q</span>g</td><td>42.259655399</td><td>-7.9918854619</td><td>1624.7163564</td><td>-561.33769564</td><td>44.056509019</td></tr><tr><td><span style="text-decoration: overline;">q</span>g → <span style="text-decoration: overline;">q</span><i>q</i><span style="text-decoration: overline;">q</span></td><td>28.497167934</td><td>-6.2611415380</td><td>947.84964732</td><td>-322.54996102</td><td>28.093290494</td></tr></tbody></table>

**Open Access.** This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.

**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 [arXiv:1207.4957] [INSPIRE].

[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 [arXiv:1304.7498] [INSPIRE].

[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 [arXiv:1412.1633] [INSPIRE].

[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_s(m_Z)$, Phys. Lett. B 750 (2015) 427 [arXiv:1508.01579] [INSPIRE].

[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 [arXiv:1301.7310] [INSPIRE].

[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 [arXiv:1611.01460] [INSPIRE].

[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 [arXiv:1705.10271] [INSPIRE].
[8] M. Czakon, A. van Hameren, A. Mitov and R. Poncelet, Single-jet inclusive rates with exact color at $O(\alpha_s^4)$, JHEP 10 (2019) 262 [arXiv:1907.12911] [INSPIRE].

[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 [arXiv:1310.3993] [INSPIRE].

[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] [INSPIRE].

[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] [INSPIRE].

[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] [INSPIRE].

[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, (2000) [hep-ph/0009193] [INSPIRE].

[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 [arXiv:1804.10017] [INSPIRE].

[15] M. Reyer, M. Schönherr and S. Schumann, Full NLO corrections to 3-jet production and $R_{32}$ at the LHC, Eur. Phys. J. C 79 (2019) 321 [arXiv:1902.01763] [INSPIRE].

[16] S. Badger, H. Frellesvig and Y. Zhang, A Two-Loop Five-Gluon Helicity Amplitude in QCD, JHEP 12 (2013) 045 [arXiv:1310.1051] [INSPIRE].

[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 [arXiv:1712.02229] [INSPIRE].

[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 [arXiv:1712.03946] [INSPIRE].

[19] S. Badger et al., Applications of integrand reduction to two-loop five-point scattering amplitudes in QCD, PoS LL2018 (2018) 006 [arXiv:1807.09709] [INSPIRE].

[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 [arXiv:1809.09067] [INSPIRE].

[21] T. Gehrmann, J.M. Henn and N.A. Lo Presti, Analytic form of the two-loop planar five-gluon all-plus-helicity amplitude in QCD, Phys. Rev. Lett. 116 (2016) 062001 [Erratum ibid. 116 (2016) 189903] [arXiv:1511.05409] [INSPIRE].

[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 [arXiv:1811.11699] [INSPIRE].

[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 [arXiv:1812.04586] [INSPIRE].

[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 [arXiv:1904.00945] [INSPIRE].
[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 [arXiv:2010.14525] [INSPIRE].

[26] S. Abreu, L.J. Dixon, E. Herrmann, B. Page and M. Zeng, *The two-loop five-point amplitude in $N=4$ super-Yang-Mills theory*, Phys. Rev. Lett. 122 (2019) 121603 [arXiv:1812.08941] [INSPIRE].

[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 [arXiv:1812.11057] [INSPIRE].

[28] S. Caron-Huot, D. Chicherin, J. Henn, Y. Zhang and S. Zoia, *Multi-Regge Limit of the Two-Loop Five-Point Amplitudes in $N=4$ Super Yang-Mills and $N=8$ Supergravity*, JHEP 10 (2020) 188 [arXiv:2003.03120] [INSPIRE].

[29] S. Badger et al., *Analytic form of the full two-loop five-gluon all-plus helicity amplitude*, Phys. Rev. Lett. 123 (2019) 071601 [arXiv:1905.03733] [INSPIRE].

[30] C.F. Berger et al., *Next-to-Leading Order QCD Predictions for $W+3$-Jet Distributions at Hadron Colliders*, Phys. Rev. D 80 (2009) 074036 [arXiv:0907.1984] [INSPIRE].

[31] H. Ita and K. Ozeren, *Colour Decompositions of Multi-quark One-loop QCD Amplitudes*, JHEP 02 (2012) 118 [arXiv:1111.4193] [INSPIRE].

[32] S. Abreu, B. Page, E. Pascual and V. Sotnikov, *Leading-Color Two-Loop QCD Corrections for Three-Photon Production at Hadron Colliders*, JHEP 01 (2021) 078 [arXiv:2010.15834] [INSPIRE].

[33] H.A. Chawdhry, M. Czakon, A. Mitov and R. Poncelet, *Two-loop leading-color helicity amplitudes for three-photon production at the LHC*, JHEP 06 (2021) 150 [arXiv:2012.13553] [INSPIRE].

[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 [arXiv:1911.00479] [INSPIRE].

[35] S. Kallweit, V. Sotnikov and M. Wiesemann, *Triphoton production at hadron colliders in NNLO QCD*, Phys. Lett. B 812 (2021) 136013 [arXiv:2010.04681] [INSPIRE].

[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$, JHEP 04 (2021) 201 [arXiv:2102.01820] [INSPIRE].

[37] S. Badger, H.B. Hartanto and S. Zoia, *Two-loop QCD corrections to $Wb\bar{b}$ production at hadron colliders*, Phys. Rev. Lett. 127 (2021) 012001 [arXiv:2102.02516] [INSPIRE].

[38] C.G. Papadopoulos, D. Tommasini and C. Wever, *The Pentabox Master Integrals with the Simplified Differential Equations approach*, JHEP 04 (2016) 078 [arXiv:1511.09404] [INSPIRE].

[39] D. Chicherin, T. Gehrmann, J.M. Henn, P. Wasser, Y. Zhang and S. Zoia, *All Master Integrals for Three-Jet Production at Next-to-Next-to-Leading Order*, Phys. Rev. Lett. 123 (2019) 041603 [arXiv:1812.11160] [INSPIRE].

[40] T. Gehrmann, J.M. Henn and N.A. Lo Presti, *Pentagon functions for massless planar scattering amplitudes*, JHEP 10 (2018) 103 [arXiv:1807.09812] [INSPIRE].

[41] D. Chicherin and V. Sotnikov, *Pentagon Functions for Scattering of Five Massless Particles*, JHEP 12 (2020) 167 [arXiv:2009.07803] [INSPIRE].

[42] A. von Manteuffel and R.M. Schabinger, *A novel approach to integration by parts reduction*, Phys. Lett. B 744 (2015) 101 [arXiv:1406.4513] [INSPIRE].
[43] T. Peraro, *Scattering amplitudes over finite fields and multivariate functional reconstruction*, JHEP 12 (2016) 030 [arXiv:1608.01902] [INSPIRE].

[44] S. Abreu et al., *Caravel: A C++ Framework for the Computation of Multi-Loop Amplitudes with Numerical Unitarity*, Comput. Phys. Commun. 267 (2021) 108069 [arXiv:2009.11957] [INSPIRE].

[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 [arXiv:1107.5131] [INSPIRE].

[47] S. Catani, *The singular behavior of QCD amplitudes at two loop order*, Phys. Lett. B 427 (1998) 161 [hep-ph/9802439] [INSPIRE].

[48] T. Becher and M. Neubert, *Infrared singularities of scattering amplitudes in perturbative QCD*, Phys. Rev. Lett. 102 (2009) 162001 [Erratum ibid. 111 (2013) 199905] [arXiv:0901.0722] [INSPIRE].

[49] E. Gardi and L. Magnea, *Factorization constraints for soft anomalous dimensions in QCD scattering amplitudes*, JHEP 03 (2009) 079 [arXiv:0901.1091] [INSPIRE].

[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 [arXiv:1506.05301] [INSPIRE].

[51] H. Ita, *Two-loop Integrand Decomposition into Master Integrals and Surface Terms*, Phys. Rev. D 94 (2016) 116015 [arXiv:1510.05626] [INSPIRE].

[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 [arXiv:1703.05273] [INSPIRE].

[53] N. Arkani-Hamed, J.L. Bourjaily, F. Cachazo and J. Trnka, *Local Integrals for Planar Scattering Amplitudes*, JHEP 06 (2012) 125 [arXiv:1012.6032] [INSPIRE].

[54] B. Ruijl, T. Ueda and J. Vermaseren, *FORM version 4.2*, arXiv:1707.06453 [INSPIRE].

[55] J. Kuipers, T. Ueda and J.A.M. Vermaseren, *Code Optimization in FORM*, Comput. Phys. Commun. 189 (2015) 1 [arXiv:1310.7007] [INSPIRE].

[56] D. Chicherin, J. Henn and V. Mitev, *Bootstrapping pentagon functions*, JHEP 05 (2018) 164 [arXiv:1712.09610] [INSPIRE].

[57] SHERPA collaboration, *Event Generation with Sherpa 2.2*, SciPost Phys. 7 (2019) 034 [arXiv:1905.09127] [INSPIRE].

[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 [arXiv:1411.1855] [INSPIRE].

[59] M. Cacciari, G.P. Salam and G. Soyez, *FastJet User Manual*, Eur. Phys. J. C 72 (2012) 1896 [arXiv:1111.6097] [INSPIRE].

[60] S. Badger, C. Brønnum-Hansen, H.B. Hartanto, R. Moodie, T. Peraro and S. Zoia, private correspondence (2021).

[61] C.F. Berger et al., *An Automated Implementation of On-Shell Methods for One-Loop Amplitudes*, Phys. Rev. D 78 (2008) 036003 [arXiv:0803.4180] [INSPIRE].