aainstitutetext: Dipartimento di Fisica and Arnold-Regge Center, Università di
Torino, and INFN, Sezione di Torino, Via P. Giuria 1, I-10125 Torino,
Italybbinstitutetext: Institut für Theoretische Teilchenphysik und Kosmologie,
RWTH Aachen University,
D-52056 Aachen, Germanyccinstitutetext: Cavendish Laboratory, University of
Cambridge, Cambridge CB3 0HE, United Kingdomddinstitutetext: Dipartimento di
Fisica e Astronomia, Università di Bologna e INFN, Sezione di Bologna, via
Irnerio 46, I-40126 Bologna, Italy

# Isolated photon production in association with a jet pair through next-to-
next-to-leading order in QCD

Simon Badger b Michał Czakon c Heribertus Bayu Hartanto a Ryan Moodie d
Tiziano Peraro c Rene Poncelet a Simone Zoia<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this work, we provide a comprehensive set of differential cross-section
distributions for photon + di-jet production in proton-proton collisions with
next-to-next-to-leading order precision in massless QCD. The event selection
corresponds to recent measurements by the ATLAS collaboration. We observe an
improved description of data in comparison to lower-order calculations in the
case of observables that are expected to be well described by perturbation
theory. The results also show better agreement with data than parton-shower-
matched and multi-jet-merged predictions generated for the ATLAS analysis
using the Sherpa Monte Carlo. A particular highlight of our study is the use
of exact five-point two-loop virtual amplitudes. This is the first calculation
of a complete two-to-three hadron-collider process at next-to-next-to-leading
order in QCD that does not rely on the leading-colour approximation at two
loops. We demonstrate, nevertheless, that the sub-leading-colour effects
present in the infrared- and ultraviolet-finite double-virtual contributions
are negligible in view of the remaining scale uncertainties.

††preprint: Cavendish-HEP-23/02, P3H-23-022, TTK-23-09

##  1 Introduction

The production of photons in high-energy collisions has a rich phenomenology
due to prompt production in the partonic process as well as secondary emission
after the hadronization phase. In fact, the most significant fraction of
photons produced at colliders such as the Large Hadron Collider (LHC)
originate from hadron decays. The smaller fraction of prompt (or primary)
photons, however, is more interesting to us since those photons allow for the
application and testing of perturbative Quantum Chromodynamics (QCD). Prompt
photons can be produced through two mechanisms: hard emission and
fragmentation. Photons produced by the first mechanism are sensitive to the
QCD dynamics of the interaction and are well separated from the other
participating partons by definition. In contrast, photons produced through
fragmentation of the partonic final state are found in proximity to hadrons,
since they are generated by collinear emissions from partons. The separation
is not unique, since it requires to define the exact meaning of collinear
emissions.

In order to separate the primary from the secondary photons, isolation
criteria are implemented experimentally, restricting the hadronic activity
around a photon candidate by requiring the (transverse) hadronic energy,
$E_{\perp}$, in a cone of size $R$ to be below (a possibly dynamical)
threshold. A non-zero threshold, as practical realities require, implies a
non-vanishing fragmentation component and, therefore, requires a corresponding
treatment of these contributions in theory predictions. The fragmentation
process can be treated in collinear/mass factorization Koller:1978kq ;
Laermann:1982jr . It has been implemented in public codes up to NLO, for
example in JetPhox Catani:2002ny . Recently, first predictions at NNLO QCD,
including photon fragmentation, have become available Gehrmann:2022cih ;
Chen:2022gpk . Within the fragmentation formalism, collinear divergencies are
absorbed into a fragmentation function that contains perturbative and non-
perturbative contributions. The latter must either be modelled or measured
Glover:1993xc ; ALEPH:1995zdi ; Bourhis:1997yu ; OPAL:1997lep ;
Belghobsi:2009hx ; Kaufmann:2017lsd . Alternatively, the so-called smooth-cone
isolation, or Fixione isolation, criterion Frixione:1998jh  can be implemented
to entirely remove the fragmentation component from the theory predictions by
vetoing radiation collinear to the photon. This type of isolation has been
used in many QCD studies Catani:2002ny ; Gehrmann:2013aga ; Chen:2019zmr ;
Gehrmann:2020oec ; Chen:2022gpk . Frixione isolation is handy from the theory
point of view but cannot be implemented in experiment due to the limited
resolution of detectors. A compromise is achieved by the usage of hybrid
isolation Siegert:2016bre ; Chen:2022gpk , which combines a smooth-cone and a
hard-cone prescription. While this is closer to the experimental setup, the
removed fragmentation component can still lead to corrections at the percent
level Chen:2022gpk . It should be mentioned that isolation criteria lead to
logarithmically enhanced contributions for small isolation regions. This may
induce sizeable resummation effects Catani:2013oma ; Balsiger:2018ezi ;
Becher:2022rhu .

The production of isolated photons has been studied experimentally since early
days of collider Athens-Athens-Brookhaven-CERN:1979oxx ; Anassontzis:1982gm ;
CMOR:1989qzc ; UA1:1988zam ; UA2:1991wce  and fixed-target experiments
WA70:1987vvj ; FermilabE706:2004emk . At the Tevatron, for example, cross
sections for photons and photons accompanied by jets have been extensively
measured as demonstrated by Refs. D0:2005ofv ; D0:2008chx ; D0:2013lra ;
CDF:2017cuc . More recently, isolated photon production has become one of the
standard-candle measurements at the LHC ATLAS:2017nah ; ATLAS:2017xqp ;
ATLAS:2019buk ; CMS:2018qao ; CMS:2019jlq ; ALICE:2019rtd ; ATLAS:2019iaa .
Inclusive production of a photon and the production of a photon with a single
resolved jet are particularly sensitive to the gluon content of the proton
since for a non-vanishing transverse momentum of the photon the dominant
partonic-process is $qg\to\gamma q$. Furthermore, these processes can also be
used for the tuning of Monte Carlo generators as well as for jet-energy
calibration Belghobsi:2009hx .

A photon associated with a pair of jets has some unique features, making it
particularly interesting for precision QCD phenomenology Keller:1991xc . Due
to the $2\to 3$ kinematics, it can be used to study angular correlations
between the photon and jets similarly to multi-jet production processes.
However, in contrast to the three-jet case, the photon is identified, and the
phase space can be subdivided into three regions:

  * •

$E_{\perp}(\gamma)>p_{T}(j_{1})$: dominated by diagrams of $pp\to\gamma
q$-type with an additional radiated gluon. This region has enhanced
contributions from photons originating from the hard interaction.

  * •

$p_{T}(j_{1})>E_{\perp}(\gamma)>p_{T}(j_{2})$: dominated by “Bremsstrahlung”
di-jet production $pp\to jj$ where one final state quark radiates a high
energy photon. This region is sensitive to the high $z$ region of the photon
fragmentation process, with $z$ the collinear-momentum fraction of the photon.

  * •

$p_{T}(j_{2})>E_{\perp}(\gamma)$: the “Bremsstrahlung” configuration with
further soft-photon enhancement. This region is sensitive to the lower $z$
region of the photon fragmentation.

Besides providing a test of perturbative QCD and the fragmentation formalism,
the process we wish to study allows us to assess and tune Monte-Carlo event
generators. In this respect it is similar to processes with a lower number of
jets. It is, however, more challenging as far as modelling of jet activity is
concerned. As a last application, we also mention that photon + di-jet
production allows for the calibration of data-driven background estimates
Lindert:2017olm .

Higher-order corrections to the production cross-section of an isolated photon
with two jets have been computed through NLO QCD and matched to parton-shower
simulations Gleisberg:2008ta ; Sherpa:2019gpd . NNLO QCD cross sections for
$2\to 3$ processes Chawdhry:2019bji ; Kallweit:2020gcp ; Chawdhry:2021hkp ;
Czakon:2021mjy ; Badger:2021ohm ; Hartanto:2022qhh ; Hartanto:2022ypo ;
Buonocore:2022pqq ; Alvarez:2023fhi  are the current state-of-the-art in
perturbative QCD. The main bottleneck is the computation of five-point two-
loop amplitudes. In the case of massless final states, all relevant partonic
processes for three-jet Abreu:2019odu ; Abreu:2021oya , three-photon
Chawdhry:2019bji ; Abreu:2020cwb ; Chawdhry:2020for , and di-photon + jet
Agarwal:2021grm ; Chawdhry:2021mkw  production have been computed in leading-
colour approximation. These computations represent the appearing Feynman
integrals in terms of pentagon functions Chicherin:2018old ; Chicherin:2020oor
, which allow for a compact analytical form and a stable numerical evaluation.
The function space for the case of a single external mass has also been worked
out, and the first leading-colour amplitudes involving an external massive
vector boson Hartanto:2019uvl ; Badger:2021imn ; Badger:2021nhg ;
Badger:2021ega ; Badger:2022ncb  have been computed. The pentagon functions
are also available for the non-planar integrals appearing beyond the leading-
colour approximation. Still, only the di-photon + jet amplitude has been
computed at full colour Agarwal:2021vdh . In this article, we evaluate all
partonic amplitudes needed for the photon + di-jet cross sections at full
colour. Equipped with the amplitudes, we perform a first phenomenological
study of the process matching the configuration of recent measurements by
ATLAS ATLAS:2019iaa .

The article is organised as follows. In Section 2, we describe the methods
used for the evaluation of the two-loop contributions. After introducing the
kinematic variables, we provide the helicity- and colour-space decomposition
of the amplitudes. The computation of the latter by means of integration-by-
parts (IBP) identities improved by syzygy equations is described, as well as
the transformations of the occurring special functions necessary to obtain the
amplitudes for all configurations. The section is closed with comments on the
numerical evaluation including methods to improve numerical stability. In
Section 3, we provide differential cross-section distributions for a number of
observables for a setup used in a recent measurement by the ATLAS
collaboration. We comment on the convergence of the perturbative series and
discuss the possible effects of omitted electroweak and fragmentation effects.
We close the main text with conclusions and outlook on future work. In
Appendix A, we provide numerical values of the amplitudes for a benchmark
phase-space point which allows for an easy verification of our results by
independent parties. In Appendix B, we reproduce explicit formulae for the
scale dependence of the amplitudes, which is necessary since our results are
derived for one particular value of the renormalisation scale. Finally, in
Appendix C, we describe the parity transformation and external-momentum
permutations necessary to derive the complete amplitudes from partial results
for specific helicity configurations.

##  2 Two-loop amplitudes

The computation of $pp\to\gamma jj$ production at NNLO QCD accuracy requires
one- and two-loop QCD corrections to the following partonic processes:

  * •

2-quark-2-gluon channel ($0\to\bar{q}qgg\gamma$),

  * •

4-quark channel ($0\to q\bar{q}Q\bar{Q}\gamma$ with either $Q=q$ or $Q\neq
q$).

In addition to these, the $0\to gggg\gamma$ partonic channel also enters
starting at NNLO QCD, since it is a loop-induced process where the leading
order contribution is made of one-loop amplitudes. In this work, we compute
analytically the one- and two-loop QCD corrections to $0\to\bar{q}qgg\gamma$
and $0\to q\bar{q}Q\bar{Q}\gamma$.111We compute the one-loop amplitudes
including the higher-order terms in the dimensional regulator parameter
$\epsilon$. The one-loop amplitude for $0\to gggg\gamma$, as well the
amplitudes contributing to the real-virtual corrections to the NNLO QCD cross
section (i.e. the one-loop amplitudes for $pp\to\gamma jjj$) are obtained from
OpenLoops Buccioni:2019sur .

###  2.1 Kinematics

We take the external momenta $p_{i}$ in the five-particle processes
$0\to\bar{q}qgg\gamma$ and $0\to q\bar{q}Q\bar{Q}\gamma$ to be outgoing. They
satisfy momentum conservation,

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

as well as the on-shell conditions $p_{i}^{2}=0$ for $i=1,\ldots,5$. The five-
particle phase space is described by five independent scalar invariants,

$\vec{s}=\left\\{s_{12},s_{23},s_{34},s_{45},s_{15}\right\\}\,,$ (2)

with $s_{ij}=(p_{i}+p_{j})^{2}$, and a pseudo-scalar quantity,

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

The square of the latter is a scalar quantity, and can thus be expressed in
terms of $\vec{s}$. We have $\operatorname{\rm tr}_{5}^{2}=\Delta_{5}$, where
$\Delta_{5}$ is the Gram determinant of the external momenta,

$\displaystyle\Delta_{5}=\text{det}\left(s_{ij}\right)_{i,j=1,\ldots,4}\,,$
(4)

which is a degree-4 polynomial in the $\vec{s}$. The parameterisation of the
five-particle kinematics therefore contains a square root:
$\sqrt{\Delta_{5}}$. This is important in view of the computation with finite
field arithmetic discussed in Section 2.4, which requires a rational
parameterisation of the kinematics.

We work in dimensional regularisation, with $d=4-2\epsilon$ spacetime
dimensions, and keep the external momenta four-dimensional. This allows us to
use the spinor-helicity formalism, as shown in Eq. 3. We set the
renormalisation scale $\mu_{R}$ to $1$ throughout the computation of the
amplitudes. We restore the dependence on $\mu_{R}$ in our analytic results as
discussed in Appendix B.

###  2.2 Structure of the $0\to\bar{q}qgg\gamma$ amplitude

We label the scattering process as

$0\rightarrow\bar{q}(p_{1},h_{1})+q(p_{2},h_{2})+g(p_{3},h_{3})+g(p_{4},h_{4})+\gamma(p_{5},h_{5})\,,$
(5)

where $h_{i}$ is the helicity of the particle with momentum $p_{i}$. The
$L$-loop amplitude $\mathcal{M}^{(L)}$ has the following colour decomposition:

$\displaystyle\begin{aligned}
\mathcal{M}^{(L)}(1_{\bar{q}},2_{q},3_{g},4_{g},5_{\gamma})&=\sqrt{2}\,e\,g_{s}^{2}\,n^{L}\,\bigg{\\{}(t^{a_{3}}t^{a_{4}})_{i_{2}}^{\;\;\bar{i}_{1}}\mathcal{A}^{(L)}_{34}(1_{\bar{q}},2_{q},3_{g},4_{g},5_{\gamma})\\\
&+(t^{a_{4}}t^{a_{3}})_{i_{2}}^{\;\;\bar{i}_{1}}\mathcal{A}^{(L)}_{43}(1_{\bar{q}},2_{q},3_{g},4_{g},5_{\gamma})+\delta_{i_{2}}^{\;\;\bar{i}_{1}}\delta^{a_{3}a_{4}}\mathcal{A}^{(L)}_{\delta}(1_{\bar{q}},2_{q},4_{g},3_{g},5_{\gamma})\bigg{\\}}\,,\end{aligned}$
(6)

where $g_{s}$ and $e$ are the QCD and QED coupling constants, respectively,
$n=m_{\epsilon}\alpha_{s}/(4\pi)$, $\alpha_{s}=g_{s}^{2}/(4\pi)$,
$m_{\epsilon}=\mathrm{i}(4\pi)^{\epsilon}e^{-\epsilon\gamma_{E}}$, and $t^{a}$
are the generators of $SU(N_{c})$ in the fundamental representation,
normalised according to $\operatorname{\rm tr}(t^{a}t^{b})=\delta^{ab}/2$. We
note that

$\displaystyle\mathcal{A}^{(L)}_{43}(1_{\bar{q}},2_{q},3_{g},4_{g},5_{\gamma})$
$\displaystyle=\mathcal{A}^{(L)}_{34}(1_{\bar{q}},2_{q},4_{g},3_{g},5_{\gamma})\,,$
$\displaystyle\mathcal{A}^{(0)}_{\delta}(1_{\bar{q}},2_{q},4_{g},3_{g},5_{\gamma})$
$\displaystyle=0\,.$ (7)

The partial amplitudes can be further decomposed into gauge invariant sub-
amplitudes according to the closed fermion loop contribution, to the power of
$N_{c}$, and to whether the photon is attached to an internal or external
quark. At one loop we have

$\displaystyle\mathcal{A}^{(1)}_{34}$
$\displaystyle=\mathcal{Q}_{q}N_{c}A_{34;q}^{(1),N_{c}}+\mathcal{Q}_{q}\frac{1}{N_{c}}A_{34;q}^{(1),1/N_{c}}+\mathcal{Q}_{q}n_{f}A_{34;q}^{(1),n_{f}}+\bigg{(}\sum_{l}\mathcal{Q}_{l}\bigg{)}A_{34;l}^{(1),1}\,,$
(8a) $\displaystyle\mathcal{A}^{(1)}_{\delta}$
$\displaystyle=\mathcal{Q}_{q}A_{\delta;q}^{(1),1}+\bigg{(}\sum_{l}\mathcal{Q}_{l}\bigg{)}\frac{1}{N_{c}}A_{\delta;l}^{(1),1/N_{c}}\,,$
(8b)

while at two loops we have

$\displaystyle\mathcal{A}^{(2)}_{34}$
$\displaystyle=\mathcal{Q}_{q}N_{c}^{2}A_{34;q}^{(2),N_{c}^{2}}+\mathcal{Q}_{q}A_{34;q}^{(2),1}+\mathcal{Q}_{q}\frac{1}{N_{c}^{2}}A_{34;q}^{(1),1/N_{c}^{2}}+\mathcal{Q}_{q}N_{c}n_{f}A_{34;q}^{(2),N_{c}n_{f}}+\mathcal{Q}_{q}\frac{n_{f}}{N_{c}}A_{34;q}^{(2),n_{f}/N_{c}}$
$\displaystyle+\mathcal{Q}_{q}n_{f}^{2}A_{34;q}^{(2),n_{f}^{2}}+\bigg{(}\sum_{l}\mathcal{Q}_{l}\bigg{)}N_{c}A_{34;l}^{(2),N_{c}}+\bigg{(}\sum_{l}\mathcal{Q}_{l}\bigg{)}\frac{1}{N_{c}}A_{34;l}^{(2),1/N_{c}}+\bigg{(}\sum_{l}\mathcal{Q}_{l}\bigg{)}n_{f}A_{34;l}^{(2),n_{f}},$
(9a) $\displaystyle\mathcal{A}^{(2)}_{\delta}$
$\displaystyle=\mathcal{Q}_{q}N_{c}A_{\delta;q}^{(2),N_{c}}+\mathcal{Q}_{q}\frac{1}{N_{c}}A_{\delta;q}^{(1),1/N_{c}}+\mathcal{Q}_{q}n_{f}A_{\delta;q}^{(2),n_{f}}+\mathcal{Q}_{q}\frac{n_{f}}{N_{c}^{2}}A_{\delta;q}^{(2),n_{f}/N_{c}^{2}}$
$\displaystyle+\bigg{(}\sum_{l}\mathcal{Q}_{l}\bigg{)}A_{\delta;l}^{(2),1}+\bigg{(}\sum_{l}\mathcal{Q}_{l}\bigg{)}\frac{1}{N_{c}^{2}}A_{\delta;l}^{(2),1/N_{c}^{2}}+\bigg{(}\sum_{l}\mathcal{Q}_{l}\bigg{)}\frac{n_{f}}{N_{c}}A_{\delta;l}^{(2),n_{f}/N_{c}}\,.$
(9b)

We denote by $n_{f}$ the number of light quark flavours, and by
$\mathcal{Q}_{q}$ ($\mathcal{Q}_{l}$) the charge of a flavour-$q$ ($l$) quark
in units of $e$. We show representative Feynman diagrams contributing to the
$0\to\bar{q}qgg\gamma$ process at two loops in Fig. 1. The right-hand sides of
Eqs. 8 and 9 involve the following classes of partial amplitudes:

  * •

$A^{(L),j}_{i;q}$: partial amplitude where the photon is attached to the
external quark line of flavour $q$ (the first and second diagram in Fig. 1).
This type of amplitude always comes with a factor of $\mathcal{Q}_{q}$.

  * •

$A^{(L),j}_{i;l}$: partial amplitude where the photon is attached to quarks
running in the loop (the last diagram in Fig. 1). This type of amplitude
always comes with a factor of $\sum_{l}\mathcal{Q}_{l}$, where we sum over the
internal quarks of flavour $l$.

Figure 1: Representative two-loop Feynman diagrams for the
$0\to\bar{q}qgg\gamma$ partonic process, together with the partial amplitudes
they contribute to.

We derive the analytic expressions of the $A^{(L),j}_{34;i}$ and
$A^{(L),j}_{\delta;i}$ partial amplitudes for the following independent
helicity configurations:

$\displaystyle\texttt{$-++++$},$ $\displaystyle\texttt{$-+++-$},$
$\displaystyle\texttt{$-++-+$},$ $\displaystyle\texttt{$-+-++$}.$ (10)

The remaining helicity configurations and partial amplitudes can be obtained
by performing parity transformation and/or by permuting the external momenta.

It is convenient to factor out a combination of spinor products such that the
helicity amplitude is free from the spinor phase. This is particularly
important when we work with the rational parameterisation of the external
kinematics in terms of momentum twistor variables discussed in Section 2.4:
conjugation and permutation of external momenta must in fact be performed on
phase-free quantities only. For the $0\to\bar{q}qgg\gamma$ partonic process we
choose the following spinor phase factors:

$\displaystyle\Phi_{\bar{q}qgg\gamma}^{-++++}$ $\displaystyle=\frac{\langle
12\rangle\langle 31\rangle[23]}{\langle 23\rangle\langle 34\rangle\langle
45\rangle\langle 51\rangle}\,,$ (11a)
$\displaystyle\Phi_{\bar{q}qgg\gamma}^{-+++-}$ $\displaystyle=\frac{\langle
15\rangle^{3}\langle 25\rangle}{\langle 12\rangle\langle 23\rangle\langle
34\rangle\langle 45\rangle\langle 51\rangle}\,,$ (11b)
$\displaystyle\Phi_{\bar{q}qgg\gamma}^{-++-+}$ $\displaystyle=\frac{\langle
14\rangle^{3}\langle 24\rangle}{\langle 12\rangle\langle 23\rangle\langle
34\rangle\langle 45\rangle\langle 51\rangle}\,,$ (11c)
$\displaystyle\Phi_{\bar{q}qgg\gamma}^{-+-++}$ $\displaystyle=\frac{\langle
13\rangle^{3}}{\langle 12\rangle\langle 34\rangle\langle 45\rangle\langle
51\rangle}\,.$ (11d)

To obtain the squared matrix element we construct the colour-summed
$L_{1}$-loop amplitude interfered with the $L_{2}$-loop amplitude as follows

$\displaystyle\sum_{\mathrm{colour}}\mathcal{M}^{(L_{1})*}\mathcal{M}^{(L_{2})}$
$\displaystyle=2e^{2}g_{s}^{4}\,n^{L_{1}}n^{L_{2}}(N_{c}^{2}-1)\bigg{\\{}\frac{(N_{c}^{2}-1)}{4N_{c}}\left[\mathcal{A}_{34}^{(L_{1})*}\mathcal{A}_{34}^{(L_{2})}+\mathcal{A}_{43}^{(L_{1})*}\mathcal{A}_{43}^{(L_{2})}\right]$
$\displaystyle\quad-\frac{1}{4N_{c}}\left[\mathcal{A}_{34}^{(L_{1})*}\mathcal{A}_{43}^{(L_{2})}+\mathcal{A}_{43}^{(L_{1})*}\mathcal{A}_{34}^{(L_{2})}\right]+N_{c}\;\mathcal{A}_{\delta}^{(L_{1})*}\mathcal{A}_{\delta}^{(L_{2})}$
(12)
$\displaystyle\quad+\frac{1}{2}\left[\mathcal{A}_{34}^{(L_{1})*}\mathcal{A}_{\delta}^{(L_{2})}+\mathcal{A}_{43}^{(L_{1})*}\mathcal{A}_{\delta}^{(L_{2})}+\mathcal{A}_{\delta}^{(L_{1})*}\mathcal{A}_{34}^{(L_{2})}+\mathcal{A}_{\delta}^{(L_{1})*}\mathcal{A}_{43}^{(L_{2})}\right]\bigg{\\}}\,.$

The bare $L$-loop amplitude $\mathcal{M}^{(L)}$ contains both ultraviolet (UV)
and infrared (IR) singularities. The UV singularities are removed by means of
renormalisation in the $\overline{\mathrm{MS}}$ scheme, while the IR ones can
be predicted from the universal IR behaviour of QCD amplitudes Catani:1998bh ;
Becher:2009qa ; Becher:2009cu ; Gardi:2009qi . We define the finite remainder
of an amplitude, $\mathcal{R}^{(L)}$, by subtracting the IR and UV singular
parts from the bare amplitude as

$\displaystyle|\mathcal{R}^{(1)}\rangle$
$\displaystyle=\left[|\mathcal{M}^{(1)}\rangle-\frac{\beta_{0}}{\epsilon}|\mathcal{M}^{(0)}\rangle\right]-\mathbf{Z}^{(1)}|\mathcal{M}^{(0)}\rangle\,,$
(13a) $\displaystyle|\mathcal{R}^{(2)}\rangle$
$\displaystyle=\left[|\mathcal{M}^{(2)}\rangle-2\,\frac{\beta_{0}}{\epsilon}|\mathcal{M}^{(1)}\rangle-\left(\frac{\beta_{1}}{2\epsilon}-\frac{\beta_{0}^{2}}{\epsilon^{2}}\right)|\mathcal{M}^{(0)}\rangle\right]-\mathbf{Z}^{(2)}|\mathcal{M}^{(0)}\rangle-\mathbf{Z}^{(1)}|\mathcal{R}^{(1)}\rangle\,,$
(13b)

where $|\mathcal{M}^{(L)}\rangle$ is a vector in colour space built out of the
partial amplitudes given in Eq. 6. The square brackets on the right-hand sides
of Eqs. 13a and 13b separate the subtraction of the UV poles from that of the
IR ones. The IR pole operator Z,

$\displaystyle\textbf{Z}=\textbf{I}+\frac{\alpha_{s}}{4\pi}\,\textbf{Z}^{(1)}+\left(\frac{\alpha_{s}}{4\pi}\right)^{2}\textbf{Z}^{(2)}+\mathcal{O}\left(\alpha_{s}^{3}\right)\,,$
(14)

where I is the identity, is a process-dependent matrix in colour space that
encodes the IR singularities of the amplitude. We adopt the IR subtraction
scheme of Refs. Becher:2009cu ; Becher:2009qa ; the $\mathbf{Z}^{(L)}$’s can
be read off from Eq. (12) of Ref. Becher:2009qa . We give explicit expressions
of $\mathbf{Z}^{(1)}$ and $\mathbf{Z}^{(2)}$ in the ancillary files
ppajjampsite . The colour decomposition of the colour-dressed finite remainder
$\mathcal{R}^{(L)}$ in terms of partial finite remainders $\mathcal{F}^{(L)}$
for the $0\to\bar{q}qgg\gamma$ partonic process is obtained from Eq. 6 by
substituting

$\mathcal{M}^{(L)}\to\mathcal{R}^{(L)}\qquad\mathrm{and}\qquad\mathcal{A}^{(L)}_{i}\to\mathcal{F}^{(L)}_{i}\,.$
(15)

Similarly, the $(N_{c},n_{f})$-decomposition of the $0\to\bar{q}qgg\gamma$
partial finite remainders is derived from the bare amplitudes decomposition
specified in Eqs. 8 and 9 by replacing

$\mathcal{A}^{(L)}_{i}\to\mathcal{F}^{(L)}_{i}\qquad\mathrm{and}\qquad
A^{(L),k}_{i;j}\to F^{(L),k}_{i;j}\,.$ (16)

###  2.3 Structure of the $0\to q\bar{q}Q\bar{Q}\gamma$ amplitude

For the 4-quark channel, we label the scattering process as follows,

$0\rightarrow
q(p_{1},h_{1})+\bar{q}(p_{2},h_{2})+Q(p_{3},h_{3})+\bar{Q}(p_{4},h_{4})+\gamma(p_{5},h_{5})\,,$
(17)

where $q$ and $Q$ are massless quarks of distinct flavour. The identical-
flavour case is discussed below. Representative Feynman diagrams contributing
to the $0\to q\bar{q}Q\bar{Q}\gamma$ process at two loops are shown in Fig. 2.
The colour decomposition of the $L$-loop amplitude $\mathcal{M}^{(L)}$ is

$\displaystyle\mathcal{M}^{(L)}(1_{q},2_{\bar{q}},3_{Q},4_{\bar{Q}},5_{\gamma})$
$\displaystyle=\sqrt{2}\,e\,g_{s}^{2}\,n^{L}\,\bigg{\\{}\delta_{i_{1}}^{\;\;\bar{i}_{4}}\delta_{i_{3}}^{\;\;\bar{i}_{2}}\mathcal{A}^{(L)}_{1}(1_{q},2_{\bar{q}},3_{Q},4_{\bar{Q}},5_{\gamma})$
$\displaystyle\phantom{=\
}+\frac{1}{N_{c}}\delta_{i_{1}}^{\;\;\bar{i}_{2}}\delta_{i_{3}}^{\;\;\bar{i}_{4}}\mathcal{A}^{(L)}_{2}(1_{q},2_{\bar{q}},3_{Q},4_{\bar{Q}},5_{\gamma})\bigg{\\}}\,.$
(18)

We further decompose the partial amplitudes $\mathcal{A}^{(L)}_{i}$ according
the source of photon radiation,

$\displaystyle\begin{aligned}
\mathcal{A}^{(L)}_{i}(1_{q},2_{\bar{q}},3_{Q},4_{\bar{Q}},5_{\gamma})=\
&\mathcal{Q}_{q}\,\mathcal{A}^{(L)}_{i;q}(1_{q},2_{\bar{q}},3_{Q},4_{\bar{Q}},5_{\gamma})+\mathcal{Q}_{Q}\,\mathcal{A}^{(L)}_{i;Q}(1_{q},2_{\bar{q}},3_{Q},4_{\bar{Q}},5_{\gamma})\\\
&+\bigg{(}\sum_{l}\mathcal{Q}_{l}\bigg{)}\mathcal{A}^{(L)}_{i;l}(1_{q},2_{\bar{q}},3_{Q},4_{\bar{Q}},5_{\gamma})\,,\end{aligned}$
(19)

where we introduced the following sub-amplitudes:

  * •

$\mathcal{A}^{(L)}_{i;q}$: the photon is attached to the external $q$ quark
line (the first diagram in Fig. 2);

  * •

$\mathcal{A}^{(L)}_{i;Q}$: the photon is attached to the external $Q$ quark
line (the second diagram in Fig. 2);

  * •

$\mathcal{A}^{(L)}_{i,l}$: the photon is attached to an internal quark line
(the third diagram in Fig. 2). This class of sub-amplitudes vanishes both at
tree level and one loop for the $0\to q\bar{q}Q\bar{Q}\gamma$ process.

Furthermore, the $\mathcal{A}^{(L)}_{i;Q}$ sub-amplitude can be obtained from
the $\mathcal{A}^{(L)}_{i;q}$ sub-amplitude via

$\displaystyle\mathcal{A}^{(L)}_{i;Q}(1_{q},2_{\bar{q}},3_{Q},4_{\bar{Q}},5_{\gamma})$
$\displaystyle=\mathcal{A}^{(L)}_{i;q}(3_{Q},4_{\bar{Q}},1_{q},2_{\bar{q}},5_{\gamma}),$
$\displaystyle i$ $\displaystyle=1,2\;.$ (20)

Figure 2: Representative two-loop Feynman diagrams for the $0\to
q\bar{q}Q\bar{Q}\gamma$ partonic process, together with the partial amplitudes
they contribute to.

In the case where the flavours of the scattering quark pairs are identical
($Q=q$),

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

the $L$-loop colour decomposition is

$\displaystyle\mathcal{M}^{(L)}(1_{q},2_{\bar{q}},3_{q},4_{\bar{q}},5_{\gamma})=\
$
$\displaystyle\sqrt{2}\,e\,g_{s}^{2}\,n^{L}\,\bigg{\\{}\delta_{i_{1}}^{\;\;\bar{i}_{4}}\delta_{i_{3}}^{\;\;\bar{i}_{2}}\mathcal{A}^{(L)}_{1}(1_{q},2_{\bar{q}},3_{Q},4_{\bar{Q}},5_{\gamma})$
$\displaystyle+\frac{1}{N_{c}}\delta_{i_{1}}^{\;\;\bar{i}_{2}}\delta_{i_{3}}^{\;\;\bar{i}_{4}}\mathcal{A}^{(L)}_{2}(1_{q},2_{\bar{q}},3_{Q},4_{\bar{Q}},5_{\gamma})$
$\displaystyle+\delta_{i_{1}}^{\;\;\bar{i}_{2}}\delta_{i_{3}}^{\;\;\bar{i}_{4}}\mathcal{A}^{(L)}_{3}(1_{q},2_{\bar{q}},3_{Q},4_{\bar{Q}},5_{\gamma})$
$\displaystyle+\frac{1}{N_{c}}\delta_{i_{1}}^{\;\;\bar{i}_{4}}\delta_{i_{3}}^{\;\;\bar{i}_{2}}\mathcal{A}^{(L)}_{4}(1_{q},2_{\bar{q}},3_{Q},4_{\bar{Q}},5_{\gamma})\bigg{\\}}\,.$
(21)

The $\mathcal{A}^{(L)}_{3}$ and $\mathcal{A}^{(L)}_{4}$ partial amplitudes are
derived from $\mathcal{A}^{(L)}_{1}$ and $\mathcal{A}^{(L)}_{2}$ via

$\displaystyle\mathcal{A}^{(L)}_{3}(1_{q},2_{\bar{q}},3_{Q},4_{\bar{Q}},5_{\gamma})=-\mathcal{A}^{(L)}_{1}(1_{q},4_{\bar{Q}},3_{Q},2_{\bar{q}},5_{\gamma})\,,$
(22a)
$\displaystyle\mathcal{A}^{(L)}_{4}(1_{q},2_{\bar{q}},3_{Q},4_{\bar{Q}},5_{\gamma})=-\mathcal{A}^{(L)}_{2}(1_{q},4_{\bar{Q}},3_{Q},2_{\bar{q}},5_{\gamma})\,.$
(22b)

We further decompose the partial sub-amplitudes $\mathcal{A}^{(L)}_{i;q}$ and
$\mathcal{A}^{(L)}_{i;l}$ ($i=1,\dots,4$) according to the closed fermion loop
contribution and factor of $N_{c}$. The $(N_{c},n_{f})$ decomposition at one
loop is

$\mathcal{A}^{(1)}_{i;q}=N_{c}A_{i;q}^{(1),N_{c}}+\frac{1}{N_{c}}A_{i;q}^{(1),1/N_{c}}+n_{f}A_{i;q}^{(1),n_{f}}\,,$
(23)

while at two loops we have

$\displaystyle\mathcal{A}^{(2)}_{i;q}$
$\displaystyle=N_{c}^{2}A_{i;q}^{(2),N_{c}^{2}}+A_{i;q}^{(2),1}+\frac{1}{N_{c}^{2}}A_{i;q}^{(2),1/N_{c}^{2}}+N_{c}\,n_{f}A_{i;q}^{(2),N_{c}n_{f}}+\frac{n_{f}}{N_{c}}A_{i;q}^{(2),n_{f}/N_{c}}+n_{f}^{2}A_{i;q}^{(2),n_{f}^{2}}\,,$
(24a) $\displaystyle\mathcal{A}^{(2)}_{i;l}$
$\displaystyle=N_{c}\,A_{i;l}^{(2),N_{c}}+\frac{1}{N_{c}}A_{i;l}^{(2),1/N_{c}}\,.$
(24b)

We derive the analytic expressions of the partial sub-amplitudes
$A^{(L),i}_{1;q}$, $A^{(L),i}_{2;q}$, $A^{(L),i}_{1;l}$, and $A^{(L),i}_{2;l}$
for the following independent helicity configurations:

$\displaystyle\texttt{$-+-++$}\,,$ $\displaystyle\texttt{$-++-+$}\,,$
$\displaystyle\texttt{$+-+-+$}\,,$ $\displaystyle\texttt{$+--++$}\,.$ (25)

We obtain the remaining helicity configurations and partial amplitudes by
performing conjugation and/or permuting the external momenta. We choose the
following spinor phase factors:

$\displaystyle\Phi_{q\bar{q}Q\bar{Q}\gamma}^{-+-++}$
$\displaystyle=\frac{\langle 13\rangle^{3}\langle 24\rangle}{\langle
12\rangle\langle 23\rangle\langle 34\rangle\langle 45\rangle\langle
51\rangle}\,,$ (26a) $\displaystyle\Phi_{q\bar{q}Q\bar{Q}\gamma}^{-++-+}$
$\displaystyle=\frac{\langle 14\rangle^{3}\langle 23\rangle}{\langle
12\rangle\langle 23\rangle\langle 34\rangle\langle 45\rangle\langle
51\rangle}\,,$ (26b) $\displaystyle\Phi_{q\bar{q}Q\bar{Q}\gamma}^{+-+-+}$
$\displaystyle=\frac{\langle 24\rangle^{3}\langle 13\rangle}{\langle
12\rangle\langle 23\rangle\langle 34\rangle\langle 45\rangle\langle
51\rangle}\,,$ (26c) $\displaystyle\Phi_{q\bar{q}Q\bar{Q}\gamma}^{+--++}$
$\displaystyle=\frac{\langle 23\rangle^{3}\langle 14\rangle}{\langle
12\rangle\langle 23\rangle\langle 34\rangle\langle 45\rangle\langle
51\rangle}\,.$ (26d)

We define the colour-summed $L_{1}$-loop amplitude interfered with the
$L_{2}$-loop amplitude to build the squared matrix element for the $0\to
q\bar{q}q\bar{q}\gamma$ process as

$\displaystyle\sum_{\mathrm{colour}}\mathcal{M}^{(L_{1})*}\mathcal{M}^{(L_{2})}=\
$ $\displaystyle
2\,e^{2}\,g_{s}^{4}n^{L_{1}}n^{L_{2}}\bigg{\\{}N_{c}^{2}\left(\mathcal{A}^{(L_{1})*}_{1}+\frac{1}{N_{c}}\mathcal{A}^{(L_{1})*}_{4}\right)\left(\mathcal{A}^{(L_{2})}_{1}+\frac{1}{N_{c}}\mathcal{A}^{(L_{2})}_{4}\right)$
$\displaystyle+N_{c}^{2}\left(\mathcal{A}^{(L_{1})*}_{3}+\frac{1}{N_{c}}\mathcal{A}^{(L_{1})*}_{2}\right)\left(\mathcal{A}^{(L_{2})}_{3}+\frac{1}{N_{c}}\mathcal{A}^{(L_{2})}_{2}\right)$
$\displaystyle+N_{c}\left(\mathcal{A}^{(L_{1})*}_{1}+\frac{1}{N_{c}}\mathcal{A}^{(L_{1})*}_{4}\right)\left(\mathcal{A}^{(L_{2})}_{3}+\frac{1}{N_{c}}\mathcal{A}^{(L_{2})}_{2}\right)$
$\displaystyle+N_{c}\left(\mathcal{A}^{(L_{1})*}_{3}+\frac{1}{N_{c}}\mathcal{A}^{(L_{1})*}_{2}\right)\left(\mathcal{A}^{(L_{2})}_{1}+\frac{1}{N_{c}}\mathcal{A}^{(L_{2})}_{4}\right)\bigg{\\}}\,.$
(27)

To obtain the squared matrix element for $0\to q\bar{q}Q\bar{Q}\gamma$ we set
$\mathcal{A}^{(L_{2})}_{3}=\mathcal{A}^{(L_{2})}_{4}=0$ in Section 2.3.

We define the finite remainders for the $0\to q\bar{q}Q\bar{Q}\gamma$ partonic
process according to Eq. 13. For $0\to q\bar{q}Q\bar{Q}\gamma$,
$|\mathcal{M}^{(L)}\rangle$ is a vector in colour space made of the partial
amplitudes specified in Eq. 18. Similarly to the $0\to\bar{q}qgg\gamma$
process, the colour decomposition of the finite remainder is obtained by
applying to Eq. 18 the replacements in Eq. 15, while the
$(N_{c},n_{f})$-expansion of $0\to q\bar{q}Q\bar{Q}\gamma$ partial finite
remainders is derived by applying the substitutions specified in Eq. 16 to the
bare amplitude decomposition in Eqs. 23 and 24.

###  2.4 Helicity amplitude computation

We compute analytically the one- and two-loop helicity amplitudes contributing
to $pp\to\gamma jj$ production using a workflow that employs Feynman diagrams
in conjunction with numerical sampling over finite fields and functional
reconstruction techniques vonManteuffel:2014ixa ; Peraro:2016wsq  within the
FiniteFlow Peraro:2019svx  framework. This approach has already been employed
in the analytic computation of several two-loop amplitudes Hartanto:2019uvl ;
Badger:2021owl ; Badger:2021imn ; Badger:2021ega .

The computational tool-chain starts with the generation of Feynman diagrams
using Qgraf Nogueira:1991ex , followed by colour decomposition and filtering
the partial amplitudes. We extract only the minimal set of independent partial
amplitudes to be further processed. We then construct the numerators of the
loop amplitudes for all the independent helicity configurations. These loop
numerators are linear combinations of monomials involving loop-momentum
dependent scalar products and spinor strings. The coefficients of these
monomials are functions of spinor products of external momenta. In order to
perform the computation within a finite field framework, we need to use a
rational parametrisation of the external kinematics. To this end, we use the
following momentum-twistor variables Hodges:2009hk ; Badger:2013gxa ,

$x_{1}=s_{12},\qquad x_{2}=-\frac{\operatorname{\rm
tr}_{+}(1234)}{s_{12}s_{34}},\qquad x_{3}=-\frac{\operatorname{\rm
tr}_{+}(1345)}{s_{13}s_{45}},\qquad x_{4}=\frac{s_{23}}{s_{12}},\qquad
x_{5}=\frac{s_{45}}{s_{12}},$ (28)

where $\operatorname{\rm tr}_{+}(ijkl)=\operatorname{\rm
tr}((1+\gamma_{5})\not{p}_{i}\not{p}_{j}\not{p}_{k}\not{p}_{l})/2$. We denote
them cumulatively by $\vec{x}=\\{x_{1},\ldots,x_{5}\\}$. The momentum-twistor
parametrisation in Eq. 28 may be used for the other massless five-point
processes as well ($pp\to jjj$, $pp\to\gamma\gamma j$ and
$pp\to\gamma\gamma\gamma$). The only dimensionful variable in this
parameterisation is $x_{1}$. We can thus set $x_{1}=1$, and recover its
dependence at the end of the calculation. The diagram processing and
filtering, as well as the construction of helicity-dependent loop numerators,
are done using Mathematica and Form Kuipers:2012rf ; Ruijl:2017dtg  scripts,
together with the library Spinney Cullen:2010jv  to carry out the Dirac
algebra and spinor manipulations. We present our results in the ’t Hooft-
Veltman scheme.

We introduce the top-level integral families, which are a set of Feynman
integral families that have the maximum number of loop propagators (in the
two-loop five-point case there are at most 8 propagators, with 3 irreducible
scalar products). The loop-momentum dependent terms in the loop numerators can
then be written in terms of propagators of the top-level integral families,
allowing us to express the loop scattering amplitude as a linear combination
of scalar Feynman integrals. The coefficients of these scalar Feynman
integrals depend only on external kinematics in the form of momentum twistor
variables.

We then reduce the scalar Feynman integrals to a set of pure Henn:2013pwa
master integrals Papadopoulos:2015jft ; Gehrmann:2018yef ; Abreu:2018rcw ;
Chicherin:2018mue ; Chicherin:2018old  by solving the integration-by-parts
relations Tkachov:1981wb ; Chetyrkin:1981qh . We generate the IBP relations
among the integrals of the top-level families for a single ordering of the
external legs, and solve them numerically over finite fields. We obtain the
reduction of the integrals with other orderings of the external legs by
evaluating the IBP solution at permuted points.

For all families, we generated identities among integrals without higher
powers of propagators Gluza:2010ws ; Ita:2015tya ; Larsen:2015ped  with
respect to those appearing in the amplitude. We obtained them by using the
Baikov representation of Feynman integrals Baikov:1996rk ; Baikov:1996cd ,
where IBP identities generally contain both higher powers of propagators and
dimensionally-shifted integrals. The latter can be eliminated by solving
polynomial equations called syzygy equations. Equations without dimension-
shifted integrals can be obtained as closed-form solutions, as illustrated in
Ref. Bohm:2017qme . Following Ref. vonManteuffel:2019wbj , we put these
solutions into a sparse matrix, and eliminated the higher powers of
denominators through Gaussian elimination. We then reconstructed new template
equations using FiniteFlow’s sparse solver Peraro:2019svx . From these, we
generated a significantly smaller system of equations, which we then solved
similarly to the more traditional Laporta algorithm Laporta:2001dd .

Next, we Laurent expand around $\epsilon=0$ up to the required order. The
Laurent expansion of the master integrals is expressed in terms of a basis of
special functions called pentagon functions and transcendental constants (e.g.
$\pi,\zeta_{3}$…) Chicherin:2020oor . We denote the pentagon functions
cumulatively by $f$, and the associated transcendental constants by
$\mathfrak{c}$. At this stage, we write the partial amplitudes as

$A^{(L),k}_{i;j}=\sum_{s=-2L}^{o(L)}\sum_{r}\epsilon^{s}\,c_{r,s}(\vec{x})\,\mathrm{mon}_{r}(f,\mathfrak{c})\,,$
(29)

where $o(1)=2$ and $o(2)=0$, and $\mathrm{mon}_{r}(f,\mathfrak{c})$ are
monomials of pentagon functions $f$ and transcendental constants
$\mathfrak{c}$.

In previous work Badger:2021imn , it was found that directly reconstructing
the finite remainder decreased the reconstruction time, as it had a simpler
expression than the full amplitude. For $pp\to\gamma jj$ production, we find
that the drop in polynomial degrees from amplitude to finite remainder is
insignificant with respect to the reconstruction time. Therefore, we
reconstruct the bare amplitude and subsequently compute the finite remainder.

In order to optimise the rational reconstruction of the amplitudes, we follow
the strategy of Ref. Badger:2021imn . We refer to that work for a detailed
discussion, and give here just an outline of the steps.

original:

Reconstruction of the coefficients $c_{r,s}(\vec{x})$ in Eq. 29 without any
optimisations.

stage 1:

We fit the $\mathbb{Q}$-linear relations among the coefficients
$c_{r,s}(\vec{x})$, and solve them for the linearly independent coefficients,
which are chosen to have the lowest polynomial degrees.

stage 2:

We determine the denominators of the coefficients $c_{r,s}(\vec{x})$ by
reconstructing them analytically on a univariate slice Abreu:2018zmy , and
matching the result against an ansatz. The latter is a product of letters of
the pentagon alphabet Chicherin:2017dob , and spinor products ($\langle
ij\rangle$, $[ij]$). We then multiply away the identified factors.

stage 3:

We perform a univariate partial fraction decomposition of the coefficients
$c_{r,s}(\vec{x})$ with respect to $x_{4}$ on the fly, and reconstruct the
coefficients of the decomposition. The latter depend on $(x_{2},x_{3},x_{5})$
only.

stage 4:

We perform another factor-matching (stage 2), this time on the coefficients of
the partial fraction decompositions of the $c_{r,s}(\vec{x})$’s. We enlarge
the ansatz used in stage 2 to include also the spurious factors introduced by
the univariate partial fraction decomposition (stage 3).

Information on the degrees at each stage of the reconstruction for the most
complicated reconstructions is shown in Table 1 for $0\to\bar{q}qgg\gamma$ and
in Table 2 for $0\to q\bar{q}Q\bar{Q}\gamma$. We provide the complete analytic
expressions for the one- and two-loop finite remainders and IR/UV poles in the
ancillary files ppajjampsite .

amplitude | helicity | original | stage 1 | stage 2 | stage 3 | stage 4
---|---|---|---|---|---|---
$A^{(2),1}_{34;q}$ | $-++-+$ | 94/91 | 74/71 | 74/0 | 22/18 | 22/0
$A^{(2),1}_{34;q}$ | $-+-++$ | 93/89 | 90/86 | 90/0 | 24/14 | 18/0
$A^{(2),1/N_{c}^{2}}_{34;q}$ | $-++-+$ | 90/88 | 73/71 | 73/0 | 23/18 | 22/0
$A^{(2),1/N_{c}^{2}}_{34;q}$ | $-+-++$ | 90/86 | 86/82 | 86/0 | 24/14 | 19/0
$A^{(2),1/N_{c}}_{34;l}$ | $-+-++$ | 89/82 | 74/67 | 73/0 | 27/14 | 20/0
$A^{(2),1/N_{c}}_{34;l}$ | $-++-+$ | 85/81 | 61/58 | 60/0 | 27/18 | 20/0
$A^{(2),N_{c}^{2}}_{34;q}$ | $-+-++$ | 58/55 | 54/51 | 53/0 | 20/16 | 20/0
Table 1:  Maximal numerator/denominator polynomial degrees of the rational coefficients of the six most complicated $0\to\bar{q}qgg\gamma$ amplitudes ordered by the maximal numerator degree prior to any optimisation. The most complicated leading-colour contribution is also included in the last row of the table for comparison.  amplitude | helicity | original | stage 1 | stage 2 | stage 3 | stage 4
---|---|---|---|---|---|---
$A^{(2),1}_{2;q}$ | $-++-+$ | 83/80 | 53/50 | 53/0 | 20/14 | 16/0
$A^{(2),1/N_{c}^{2}}_{2;q}$ | $-++-+$ | 83/80 | 55/52 | 55/0 | 22/14 | 16/0
$A^{(2),1/N_{c}^{2}}_{2;q}$ | $-+-++$ | 83/80 | 62/59 | 62/0 | 22/13 | 17/0
$A^{(2),1}_{1;q}$ | $-+-++$ | 83/80 | 29/26 | 27/0 | 19/10 | 14/0
$A^{(2),1/N_{c}^{2}}_{1;q}$ | $-++-+$ | 83/80 | 39/35 | 39/0 | 26/16 | 20/0
$A^{(2),1/N_{c}^{2}}_{1;q}$ | $-+-++$ | 83/80 | 42/37 | 42/0 | 26/14 | 20/0
$A^{(2),N_{c}^{2}}_{1;q}$ | $-++-+$ | 45/42 | 55/52 | 55/0 | 22/14 | 16/0
Table 2:  Same as Table 1, but for $0\to q\bar{q}Q\bar{Q}\gamma$.

###  2.5 Permutation and conjugation of the pentagon functions

The analytic expressions of the partial amplitudes are written in terms of
rational coefficients and pentagon functions. The latter can be evaluated
numerically in the $s_{12}$ channel through the library PentagonFunctions++
Chicherin:2020oor . To obtain all the helicity configurations and the full set
of partial amplitudes from the independent set, we need to permute the
external momenta entering the rational coefficients as well as the pentagon
functions. Such permutations may require crossings of initial-state momenta to
the final state, preventing the evaluation with PentagonFunctions++. To
overcome this, we exploit the fact that the complete set of pentagon functions
is closed under permutations of the external momenta. In other words, any
permutation $\sigma$ of the pentagon functions $f(\vec{s})$ can be expressed
in terms of pentagon functions in the original ordering of the momenta as

$\displaystyle\sigma\circ
f(\vec{s}):=f\left(\sigma\circ\vec{s}\right)=M^{(\sigma)}\cdot\text{mon}\left(f(\vec{s}),\mathfrak{c}\right)\,,$
(30)

where $M^{(\sigma)}$ is a matrix of rational numbers, and
$\text{mon}\left(f(\vec{s}),\mathfrak{c}\right)$ is the vector of all
monomials in $f(\vec{s})$ and $\mathfrak{c}$ up to transcendental weight 4. We
computed the matrices $M^{(\sigma)}$ from the expressions of all permutations
of the master integrals in terms of pentagon functions provided in Ref.
Chicherin:2020oor . Thanks to Eq. 30, it suffices to evaluate the pentagon
functions in one phase-space point in the $s_{12}$ channel to obtain the
values of all their permutations as well. This not only removes the problem of
evaluating outside the $s_{12}$ channel, but also reduces to the minimum the
number of calls to PentagonFunctions++. Since the latter are the most
expensive part of the numerical evaluation of the two-loop amplitudes, this
approach is also more efficient.

The pentagon functions (as well as the associated transcendental constants)
have well-defined behaviour under parity conjugation. They are either even or
odd, i.e., they either stay the same or change sign. Conjugating them thus
amounts to flipping the sign of the odd ones. In this regard, note that the
pseudo-scalar invariant $\operatorname{\rm tr}_{5}$ is purely imaginary in any
physical channel. The library PentagonFunctions++ takes $\operatorname{\rm
tr}_{5}$ to have a positive imaginary part. The values of the pentagon
functions for the opposite choice can be obtained by parity conjugation.

We provide the pentagon-function permutation rules required to obtain the full
set of helicity configurations and partial amplitudes for $pp\to\gamma jj$
production in the ancillary files under the folder `amplitudes/permrules`
ppajjampsite .

###  2.6 Validation

We perform the following checks to validate our results for the amplitudes.

Ward-identity test.

We verify the gauge invariance of the bare helicity amplitudes by replacing
the gluon or photon polarisation vector by its momentum, and checking that the
resulting amplitudes vanish.

UV and IR poles subtraction.

The fact that the finite remainders in Eq. 13 are indeed finite at
$\epsilon=0$ proves that the amplitudes have the expected UV and IR poles.

Comparison against OpenLoops.

We compare our results numerically against OpenLoops at the level of the
squared matrix element for the tree-level amplitudes and the one-loop finite
remainders (at $\mathcal{O}(\epsilon^{0})$), for all scattering channels
contributing to $pp\to\gamma jj$ production.

Checks on the symmetry of the amplitude.

We perform several checks to demonstrate that the amplitudes have the expected
symmetries. For example, the $0\to\bar{q}qgg\gamma$ subprocess is symmetric
under the exchanges $1\to 2$ and $3\to 4$. We apply the exchange in two ways:
(1) by permuting rational coefficients and pentagon functions, and (2) by
permuting the external momenta in the numerical input. We check that the two
approaches give the same results.

One-loop amplitude check through $\mathcal{O}(\epsilon^{2})$.

To ensure that the $\operatorname{\rm tr}_{5}$ invariant is correctly
evaluated when performing permutations and conjugation as well as in the
numerical evaluation, we derive the full set of partial amplitudes for all the
contributing helicity configurations in terms of scalar master integrals.
Unlike some of the pure master integrals, the latter are parity even, i.e.,
they do not depend on the sign of $\operatorname{\rm tr}_{5}$. We then
evaluate the scalar master integrals numerically using the package AMFlow
Liu:2022chg  up to the necessary $\epsilon$ order such that the amplitude is
evaluated through $\mathcal{O}(\epsilon^{2})$. The amplitude computed using
this approach is then compared against the amplitude obtained by permutation
and/or conjugation.

###  2.7 Numerical evaluation

We implemented the amplitudes presented above in a C++ code that yields the
hard functions, defined in Appendix A, by performing the colour and helicity
sums for the IR- and UV-finite squared remainders according to Eqs. (12) and
(2.3), as required for the evaluation of double-virtual contributions to
differential cross sections. The code utilises the numerical implementation of
the special functions provided in PentagonFunctions++ Chicherin:2020oor , and
has been checked against an independent implementation in Mathematica for a
number of phase-space points. If needed, higher-precision floating-point
arithmetic is available through the qd library QDlib , which implements the
double-double and quad-double data types. This is indeed necessary, since
cancellations between rational terms and/or special functions for some phase-
space points lead to a loss of precision resulting in values without a
sufficient number of correct digits. In order to capture these problematic
cases, we employ the following algorithm for each phase-space point
$p=\\{p_{i}\\}_{i=1}^{5}$:

  1. 1)

The two-loop hard function $\mathcal{H}^{(2)}(p,\mu_{R})$ is evaluated at $p$
with renormalisation scale $\mu_{R}$ in standard double precision;

  2. 2)

Additionally, the two-loop hard function is evaluated at
$p^{\prime}=\\{a\times p_{i}\\}_{i=1}^{5}$ and renormalisation scale
$\mu^{\prime}_{R}=a\times\mu_{R}$, with $a$ a numerical constant. The quantity

$r^{\text{diff}}=\frac{|\mathcal{H}^{(2)}(p,\mu_{R})-a^{2}\,\mathcal{H}^{(2)}(p^{\prime},\mu^{\prime}_{R})|}{|\mathcal{H}^{(2)}(p,\mu_{R})+a^{2}\,\mathcal{H}^{(2)}(p^{\prime},\mu^{\prime}_{R})|}\;,$
(31)

provides an estimate of the relative numerical precision of the squared matrix
element. Here, we have taken into account the mass dimension of
$\mathcal{H}^{(2)}(p,\mu_{R})$ for a $2\to 3$ scattering process. The
estimated “number of correct digits” $n^{\text{cor}}$ is given by

$n^{\text{cor}}=-\log_{10}(r^{\text{diff}})\;.$ (32)

  3. 3)

If $n^{\text{cor}}>4$ the phase-space point is accepted and the algorithm
restarts from 1) for the next phase-space point. If $n^{\text{cor}}\leq 4$ the
floating point precision is increased by switching to the next larger data
type and the algorithm restarts from 1). Instead of simply using the same data
type for every operation during amplitude evaluation, we implement the
following precision levels:

L1:

pentagon and rational functions are evaluated with double precision ($\approx
16$ digits);

L2:

pentagon functions are evaluated with double precision, while the rational
functions with double-double precision ($\approx 32$ digits);

L3:

pentagon functions are evaluated with double-double precision, while the
rational functions with quad-double precision ($\approx 64$ digits).

In order to obtain the differential distributions presented in Section 3, we
have evaluated the squared matrix elements at $\approx 2.2$ million phase-
space points. The distribution of the estimated numerical precision for the
three different precision levels is shown in Fig. 3. About $1\%$ of the phase-
space points required at least level L2 to achieve four correct digits. Among
those points, $10\%$ were identified as still unstable and had to be evaluated
at L3. About $0.1\%$ of the whole set of phase-space points required an
evaluation at L3 to achieve the precision goal of four correct digits.

Figure 3: Fraction of phase-space points with a given number of correct digits
with respect to the total number of points. The labels “L1: sd,sd”, “L2:
sd,dd”, and “L3: dd,qd” denote the implementation levels L1, L2 and L3. The
labels “sd”, “dd”, and “qd” refer to single-double, double-double, and quad-
double precision. Further details are given in the text.

The total computing time for the double-virtual contributions was $10$ kCPUh
corresponding roughly to $16$ seconds on average for the evaluation at a
single phase-space point. This timing includes stability checks and potential
re-evaluations with higher precision.

##  3 LHC Phenomenology at 13 TeV

###  3.1 Setup

The evaluation of perturbative second-order corrections to the differential
cross sections for photon + di-jet production requires the combination of
double-virtual, real-virtual and double-real contributions. To this end, we
employ Stripper, an implementation of the sector-improved residue subtraction
scheme Czakon:2010td ; Czakon:2014oma ; Czakon:2019tmo , which also provides
the cross sections up to next-to-leading order. This software has been used
previously in various cross-section calculations and has been validated in
numerous comparisons to public results. Tree-level matrix elements are
implemented through the AvH library Bury:2015dla . For the single-virtual
contributions, which require one-loop matrix elements with one photon and up
to three jets, we employ OpenLoops2 Buccioni:2019sur . The evaluation of the
squared matrix elements needed for the double-virtual corrections has been
described in the previous Section.

We work in the five-flavour, $n_{f}=5$, scheme, i.e. we consider QCD with five
massless quarks and ignore the existence of top-quarks. All predictions
presented in this section use the NNPDF31 NNPDF:2017mvq  PDF set as
implemented in LHAPDF Buckley:2014ana . We compute the cross section for two
different nominal choices of the renormalization ($\mu_{R}$) and factorization
scale ($\mu_{F}$):

$\displaystyle\mu_{R}$
$\displaystyle=\mu_{F}=H_{T}=E_{\perp}(\gamma)+p_{T}(j_{1})+p_{T}(j_{2})\;\quad\text{and}$
(33a) $\displaystyle\mu_{R}$ $\displaystyle=\mu_{F}=E_{\perp}(\gamma)\;,$
(33b)

where $E_{\perp}(\gamma)$ is the transverse energy of the photon and
$p_{T}(j_{i})$ denotes the $i$-th leading jet transverse momentum. We estimate
the uncertainty from missing higher orders by performing a conventional
7-point scale variation around the nominal $\mu_{R,F}$ values by factors of
$2$ subject to the constraint $\frac{1}{2}\leq\frac{\mu_{R}}{\mu_{F}}\leq 2$.

In order to showcase the theoretical NNLO QCD predictions for differential
cross sections, we will compare our results to measurement data obtained by
the ATLAS collaboration at 13 TeV ATLAS:2019iaa . As already explained in the
introduction, the final-state photon may either be produced through hard
emission, well separated from other partons, or through collinear emission
described within the fragmentation formalism. Since the main achievement of
the current publication is the evaluation of the two-loop virtual corrections
for the two-to-three processes with four partons and one photon, our attention
is devoted to the phase space region which is well described by fixed-order
perturbative methods without collinear enhancements. Fortunately, the
measurements presented in Ref. ATLAS:2019iaa  are divided into two samples,
one of which is more sensitive to the perturbative and one which is more
sensitive to the fragmentation component. Of course, there is always some
contamination in each region from effects characteristic of the other region.
We shall estimate the size of the unwanted effects below, and argue that they
are well within the uncertainties of the measurement data to which we compare.

Since we adopt the phase-space of Ref. ATLAS:2019iaa , our event-selection
cuts are defined as follows:

  1. 1.

We require at least two jets defined with the anti-$k_{T}$ algorithm
Cacciari:2008gp  for jet radius $R=0.4$ that have minimal transverse momentum
of $p_{T}(j)>100\;\text{GeV}$ and maximal rapidity $|\eta(j)|<2.5$.

  2. 2.

The identified jets must be separated from the photon by $\Delta
R(\gamma,j)>0.8$.

  3. 3.

One isolated photon must be present in the final state with
$E_{\perp}(\gamma)\geq 150\;\text{GeV}$, $|\eta(\gamma)|\leq 2.37$ excluding
$1.37\leq|\eta(\gamma)|\leq 1.56$.

  4. 4.

It is well-known that hard photon isolation criteria lead to infrared
sensitivity beyond the Born approximation, which manifests at fixed
perturbative order by uncompensated singularities. These singularities may be
absorbed into a fragmentation function. However, there is a simpler approach
that allows us to avoid the implementation of the fragmentation formalism.
This approach is based on a smooth photon isolation criterion proposed by
Frixione in Ref. Frixione:1998jh . Upon extension with an additional hard-cone
isolation Siegert:2016bre ; Chen:2022gpk , this method simulates the
experimental setup very closely. This hybrid isolation criterion has been
adopted in Section 4.3 of Ref. ATLAS:2019iaa  for the generation of next-to-
leading order QCD predictions with the Monte-Carlo generator Sherpa
Gleisberg:2008ta ; Sherpa:2019gpd  used for comparisons to data in that
publication. Accordingly, we accept an event if the sum, $E_{\perp}(r)$, of
the transverse energies of all partons separated from the photon by the
angular distance $\Delta R\leq r$ satisfies the smooth-cone (Frixione)
condition:

$E_{\perp}(r)\leq
E_{\perp\rm{max}}(r)=0.1\,E_{\perp}(\gamma)\Bigg{(}\frac{1-\cos(r)}{1-\cos(R_{\rm
max})}\Bigg{)}^{2}\quad\text{for}\quad r\leq R_{\rm max}=0.1\;,$ (34)

as well as the hard-cone condition:

$E_{\perp}(r)\leq
E_{\perp\rm{max}}=0.0042\,E_{\perp}(\gamma)+10\;\text{GeV}\quad\text{for}\quad
r\leq R_{\rm max}=0.4\;.$ (35)

These event-selection cuts define the inclusive sample used in the
determination of the differential cross sections. The aforementioned reduction
of the fragmentation effects is further strengthened by imposing the
requirement ATLAS:2019iaa :

$E_{\perp}(\gamma)>p_{T}(j_{1})\;,$ (36)

which defines the direct-enriched sample. The fragmentation-enriched sample
has been defined in Ref. ATLAS:2019iaa  by the requirement
$E_{\perp}(\gamma)<p_{T}(j_{2})$. Here, $j_{1}$ and $j_{2}$ are the hardest
(leading) and next-to-hardest (sub-leading) jet ordered according to their
transverse momenta.

###  3.2 Selected observables

In this section, we discuss selected results for the inclusive and direct-
enriched phase space defined above. Results for the following distributions in
all three phase spaces are provided in numeric format together with this
publication:

  1. 1.

$E_{\perp}(\gamma)$: photon transverse energy,

  2. 2.

$p_{T}^{\text{jet}}$: jet transverse momentum,

  3. 3.

$y^{\text{jet}}$: jet pseudorapidity,

  4. 4.

$|\Delta y^{\gamma-\text{jet}}|$: absolute value of the pseudorapidity
difference between the photon and the jet,

  5. 5.

$|\Delta\phi^{\gamma-\text{jet}}|$: azimuthal-angle difference between the
photon and the jet,

  6. 6.

$|\Delta y^{j_{1}-j_{2}}|$: absolute value of the pseudorapidity difference
between the leading and sub-leading jet,

  7. 7.

$|\Delta\phi^{j_{1}-j_{2}}|$: azimuthal-angle difference between the leading
and sub-leading jet,

  8. 8.

$m(j_{1}j_{2})$: invariant mass of the leading and sub-leading jet,

  9. 9.

$m(\gamma j_{1}j_{2})$: invariant mass of the photon, leading and sub-leading
jet.

In the case of $p_{T}^{\text{jet}},y^{\text{jet}},|\Delta
y^{\gamma-\text{jet}}|$ and $|\Delta\phi^{\gamma-\text{jet}}|$, every event is
counted twice by filling the given histogram using the value of the respective
observable for both the leading and the sub-leading jet.

Together with the LO, NLO and NNLO QCD predictions, we show not only the
measurement results from Ref. ATLAS:2019iaa  but also Sherpa predictions
provided by ATLAS. The latter are obtained from an NLO-matched QCD parton-
shower merged with LO photon+four-jet samples corresponding to the double-real
radiation contributions in the presented NNLO QCD predictions. The Sherpa
predictions use $E_{\perp}(\gamma)$ for the renormalisation and factorisation
scale. More details can be found in Ref. ATLAS:2019iaa .

#### Transverse photon energy $E_{\perp}(\gamma)$.

Figure 4: Differential cross sections w.r.t. the transverse energy of the
photon $E_{\perp}(\gamma)$ in the inclusive (left plot) and direct-enriched
(right plot) phase space at LO (green), NLO (blue) and NNLO (red) QCD compared
to data (black) and Sherpa (purple) prediction provided by ATLASATLAS:2019iaa
. The top panels show the absolute values for the $H_{T}$ scale choice. The
middle (bottom) panel shows the ratio to NLO QCD using the $H_{T}$
($E_{\perp}(\gamma)$) scale. The coloured bands show scale variation and the
vertical coloured bars indicate statistical uncertainties.

The differential cross section is shown in Fig. 4 for the inclusive (left
plot) and for the direct-enriched phase space (right plot). For the $H_{T}$
scale choice, the NLO QCD corrections are about $80\%$ at small transverse
energies and about $30\%$ at high energies relative to LO for both samples.
The NLO QCD predictions have sizeable estimated uncertainties from missing
higher orders of about $10-20\%$. Nevertheless, they describe the ATLAS data
within those uncertainties. The NNLO QCD corrections are positive and much
smaller than the NLO. They vary from $\approx 1\%$ at low transverse energy to
maximally $10\%$ for $E_{\perp}(\gamma)\approx 400-600\;\text{GeV}$. The scale
dependencies decrease significantly to $1-5\%$ when including the second-order
corrections. Moreover, the NNLO QCD predictions describe the data better than
the NLO as far as the distribution’s shape is concerned. On the other hand,
the normalisation is slightly high albeit still compatible within the
systematic uncertainties. The picture changes when considering results
obtained with the $E_{\perp}(\gamma)$ scale. The LO predictions are
significantly larger than for the $H_{T}$ scale, while the NLO corrections are
smaller. In this case, the NNLO QCD corrections are negative and sizeable (up
to $-20\%$) for small transverse energies. The data is again well described
although with larger remaining scale dependence. The large corrections and
uncertainties indicate that the scale choice $E_{\perp}(\gamma)$ does not
capture the relevant kinematic scales entering this observable. This is
different from the inclusive photon production case $pp\to\gamma+X$. The
reason lies most likely in the presence of new scales introduced through the
jet-selection requirements.

At higher energies, while still agreeing with data within the uncertainties,
the predicted spectrum is harder than the measurement. It has been shown in
Refs. Kuhn:2005gv ; Becher:2013zua ; Becher:2015yea  that electroweak (EW)
corrections are significant in the TeV region, and can reach about
$\approx-10\%$ for $E_{\perp}(\gamma)>1\;\text{TeV}$ in $pp\to\gamma j$.
Similar corrections can be expected for $pp\to\gamma jj$ and could alleviate
the observed discrepancy at high energies. The $E_{\perp}(\gamma)$
distribution is not the only one affected by potentially large EW corrections.
In fact, for large invariant jet pair mass $m(j_{1}j_{2})$ and large rapidity
separation between the two jets, sizeable corrections are also expected. The
inclusion of these corrections is, however, left for future work.

A striking observation is that the description of the data by the purely
perturbative QCD (pQCD) calculation is noticeably better than by the Sherpa
prediction which includes parton-shower effects and higher jet multiplicities
through merging. It appears that the Sherpa predictions systematically
overshoot the measurement at high energies, while pQCD predictions are closer
to data both in size and in shape. The uncertainties of the Sherpa predictions
estimated using scale variation are twice as large as those of the fixed order
calculation. This is the reason why the Sherpa predictions still agree with
data. A possible explanation for the size of the uncertainties despite the
fact that both approaches have at least NLO accuracy could be the multi-jet
merging with the photon+4-jet sample, which is likely to come with a large
scale dependence.

#### Jet transverse momentum $p_{T}^{\text{jet}}$.

Figure 5: The same as Fig. 4 but for the $p_{T}^{\text{jet}}$ observable.

In this case, we make similar observations as for the transverse energy of the
photon. The differential cross section is presented in Fig. 5. For the $H_{T}$
scale, we find small NNLO QCD corrections and small remaining scale
dependence. The description of the data is very clearly improved with respect
to NLO QCD and Sherpa. The observable demonstrates once more the sub-optimal
properties of the $E_{\perp}(\gamma)$ scale. Looking at the inclusive phase
space, the NNLO QCD scale dependence is much larger than in the case of
$H_{T}$, especially at high transverse momentum. This is not the case in the
direct-enriched phase space, where both scales behave similarly. This is
consistent with the fact that the two scale choices are closest (kinematics-
wise per event) when $E_{\perp}(\gamma)>p_{T}(j_{1})$ since in that case,
$E_{\perp}(\gamma)$ dominates the value of $H_{T}$. This suggests that the
region where $E_{\perp}(\gamma)<p_{T}(j_{1})$ is particularly poorly described
with the $E_{\perp}(\gamma)$ scale because this scale is relatively low
(compared to the jet momenta), leading to comparatively large
$\alpha_{S}(\mu_{R})$ and thus large corrections.

#### Invariant mass of the jet pair $m(j_{1}j_{2})$.

Figure 6: The same as Fig. 4 but for the $m(j_{1}j_{2})$ observable.

The differential cross section is presented in Fig. 6. Unsurprisingly, the
pattern of perturbative corrections is similar to the previous two cases. The
NNLO QCD corrections are small and lead to a significant reduction of the
scale dependence. The NNLO QCD predictions agree with the data across the
spectrum. Only for small invariant masses do we observe a clear deviation.
This deviation might be related to the different definition of the photon
isolation between theory and experiment. While it seems that Sherpa describes
the data in the first bin better than fixed-order QCD, we point out the large
fluctuations of the results of the Sherpa simulation, which makes any precise
statement impossible. Furthermore, we expect that resummation effects play no
role in the first bin of the distribution. This assessment is supported in
particular by the fact that the NLO and NNLO predictions with the $H_{T}$
scale are very close to each other. For high invariant mass, the Sherpa
predictions completely fail to describe the spectrum, which is noteworthy
because this observable is of special interest in searches for New Physics,
e.g. in the search for a heavy resonance decaying into jets and recoiling
against the prompt photon. The NNLO QCD predictions correctly describe the
falling spectrum with only small corrections w.r.t. NLO QCD. The scale-
variation uncertainty is smaller than the experimental uncertainties. The good
description of the data, in particular in comparison to the Sherpa spectrum,
motivates the usage of the predicted spectrum instead of relying on data-
driven background estimates. Of course, this would require a substantial
investment of computing resources in order to improve the quality of the
Monte-Carlo integration for invariant masses in the TeV range. In this region,
we observe large cancellations between different contributions to the complete
NNLO result, which impact the quality of the histograms. The cancellations are
enhanced for the $E_{\perp}(\gamma)$ scale. The results in the direct-enriched
phase space suffer further from reduced statistics due to the additional cuts.

#### Azimuthal separation between the jets and photon
$|\Delta\phi^{\gamma-\text{jet}}|$ and $|\Delta\phi^{j_{1}-j_{2}}|$.

Figure 7: The same as Fig. 4 but for the $|\Delta\phi^{\gamma-\text{jet}}|$
and $|\Delta\phi^{j_{1}-j_{2}}|$ observables in the inclusive phase space
only.

These are the last two observables that we discuss in detail. We first focus
on the inclusive phase space, for which the differential cross sections are
shown in Fig. 7. For the $H_{T}$ scale, the NLO QCD predictions agree with the
data reasonably well but the NNLO QCD corrections are essential in order to
precisely describe the shape of the distributions. The second-order
corrections improve the data to theory comparison especially for large
azimuthal separations. Interestingly, Sherpa simulations can hardly describe
the photon-jet azimuthal distance for small separation, and the jet-jet
azimuthal distance for large separation. Both cases correspond to the
situation where a photon and a nearby softer jet are likely recoiling against
a harder jet, which leads simultaneously to a small $|\Delta\phi^{\gamma-
j_{2}}|$ and large $|\Delta\phi^{j_{1}-j_{2}}|$. This is also the region of
large $m(j_{1}j_{2})$.

Figure 8: The same as Fig. 4 but for the $|\Delta\phi^{\gamma-\text{jet}}|$
and $|\Delta\phi^{j_{1}-j_{2}}|$ observables in the direct-enriched phase
space only. Note that, in this case, we show the ratio plot with respect to
data since the perturbative corrections are large. See text for more details.

In the case of the direct-enriched phase space, the considered observables
have some peculiar properties due to the geometric constraints arising from
the requirement $E_{\perp}(\gamma)>p_{T}(j_{1})$. The distance between the
photon and the jets corresponds then to the same observable which is studied
in the context of the azimuthal decorrelation in three jet production
ATLAS:2018sjf . The additional phase space constraint leads to the fact that
$|\Delta\phi^{\gamma-\text{jet}}|>\frac{2\pi}{3}$ when considering LO
kinematics. The presence of a kinematical edge is responsible for large
perturbative corrections, as visible in Fig. 8. Indeed this is the only case
where the $E_{\perp}(\gamma)$ scale choice performs better than $H_{T}$, since
it enhances the radiation of high-$p_{T}$ jets filling the phase space better.

#### Estimate of omitted fragmentation contributions.

The fragmentation contributions to the differential cross sections presented
above may be evaluated by convolving parton-to-photon fragmentation functions
with a three-jet production cross section defined with a corresponding
resolved parton. The hard-cone isolation criterion, Eq. (35), implies that
only hard photons with fragmentation-momentum fraction $z\geq
E_{\perp}(\gamma)/(E_{\perp}(\gamma)+E_{\perp,{\rm max}})>0.93$ contribute.
Fragmentation functions are strongly suppressed for large $z$, which is
responsible for the dominance of well-separated photon emissions. This effect
is further strengthened by the cut requiring the photon to be harder than the
hardest jet.

While it is impossible to provide a back-of-the-envelope estimate of the
omitted fragmentation contributions, which have been replaced by contributions
defined by Frixione’s isolation, we may nevertheless invoke the recent Ref.
Chen:2022gpk . There, the fragmentation contributions have been included for
the photon + jet process with a phase space very similar to ours, since
referring to similar measurements by ATLAS. In particular, Figure 9 in Section
5.1 of Ref. Chen:2022gpk  shows a comparison between the transverse
energy/momentum spectra with fragmentation contributions and with a hybrid
photon-isolation criterion. The difference never exceeds 5%. In view of these
findings, we also estimate that the effects omitted in our study should not
exceed 5% in neither the inclusive nor the direct-enriched phase space.

The effect is by design larger in the fragmentation-enriched part of the phase
space. Hence, we do not show any comparison of our predictions to data for
such configurations. However, for reference, we provide the corresponding
plots in the ancillary files attached to the arXiv preprint.

###  3.3 Size of the sub-leading-colour terms of the double-virtual
contributions

Figure 9: A comparison of different approximations to the double virtual
corrections. See text for details.

The differential cross sections presented in the previous Section are the
first NNLO QCD results for a two-to-three process with exact double-virtual
corrections. Previous calculations for three-photon, di-photon + jet, and
three-jet production processes relied on the leading-colour approximation of
the two-loop finite remainder, possibly augmented with some sub-leading
contributions given by additional planar diagrams. In the case of each of
these processes, it was argued that the approximation should not have a
noticeable impact on the results due to the expected size of the sub-leading-
colour contributions in comparison with the scale uncertainty of the NNLO
calculation itself. The argument was based on the small size of the finite
remainder in comparison with the complete cross section, as well as on the
experience-based expectation that sub-leading-colour contributions should
amount to about 10% of the leading-colour contributions. Here, we have the
opportunity to assess the quality of the leading-colour approximation using
the exact result. To this end, we consider the double-virtual finite remainder
cross section

$\sigma^{VVF}=\int\text{d}\Phi\,\mathcal{H}^{(2)}(\mu_{R}^{2})\;,$ (37)

where $\mathcal{H}^{(2)}(\mu_{R}^{2})$ is defined in Appendix A. Terms
depending on the logarithm of the renormalisation scale,
$\ln(\mu_{R}^{2}/\mu_{R0}^{2})$, occurring in the finite remainder
$\mathcal{H}^{(2)}(\mu_{R}^{2})$ can be derived using the formulae presented
in Refs. Becher:2009cu ; Becher:2009qa , see Appendix B. They require the
evaluation of tree-level and one-loop amplitudes, which can be done at full
colour for any process. Hence, we decompose the finite remainder as follows

$\displaystyle\begin{aligned}
\mathcal{H}^{(2)}(\mu_{R}^{2})&=\mathcal{H}^{(2)}(s_{12})+\sum_{i=1}^{4}c_{i}\ln^{i}\left(\frac{\mu_{R}^{2}}{s_{12}}\right)\\\
&\approx\mathcal{H}^{(2)}_{l.c.}(s_{12})+\sum_{i=1}^{4}c_{i}\ln^{i}\left(\frac{\mu_{R}^{2}}{s_{12}}\right)\;.\end{aligned}$
(38)

Note that the choice of the reference scale does influence the size of the
sub-leading colour effects in $\mathcal{H}^{(2)}(\mu_{R0}^{2})$. Here, just as
in our previous calculations for the aforementioned photon and jet processes,
we take $\mu_{R0}^{2}=s_{12}$. This corresponds to the value used in the
implementation of the sector-improved residue subtraction scheme in the C++
library Stripper. The choice should not lead to unwanted enhancements of
$\mathcal{H}^{(2)}(\mu_{R0}^{2})$ as long as the final states are well
separated. We compare three different approximations:

  1. 1)

$\mathcal{H}^{(2)}_{l.c.}(s_{12})$: the large $N_{c}$ limit assuming that
$n_{f}$ and $Q_{l}$ scale with the same rate. This is more than just the
planar approximation since the $Q_{l}$ terms can have non-planar
contributions.

  2. 2)

$\mathcal{H}^{(2)}_{l.c.;Q_{l}=0}(s_{12})$: same as 1) but with $Q_{l}=0$.

  3. 3)

$\mathcal{H}^{(2)}_{0}(s_{12})$: here we just set
$\mathcal{H}^{(2)}(s_{12})=0$.

A comparison of the differential cross sections (all contributions, not only
the double virtual) assuming each of the three approximations is shown in Fig.
9. We plot the ratio with respect to the full-colour result labelled “NNLO”.
We observe that there is no noticeable difference between the exact and the
approximated results in the case of $\mathcal{H}^{(2)}_{l.c.}(s_{12})$. The
$\mathcal{H}^{(2)}_{l.c.;Q_{l}=0}(s_{12})$ approximation has a small impact
reaching at most $1\%$ at the cross section level.

Finally, we can use the $\mathcal{H}^{(2)}_{0}(s_{12})$ approximation to
estimate the overall contribution of the two-loop hard functions in the NNLO
QCD result. This approximation is clearly distinguishable from the full result
and indicates that the double-virtual contribution is about $10\%$ across the
phase space. In some cases, the double-virtual contribution may be even
larger, as demonstrated by the invariant-mass distribution for large
$m(j_{1}j_{2})$ values.

As a side note, we point out that the loop-induced corrections from the $gg\to
gg\gamma$ one-loop squared matrix elements never exceed $1$ permille relative
to NLO QCD.

##  4 Conclusions

We have derived the exact two-loop amplitudes in the four-quark and two-quark-
two-gluon channels necessary for the complete description of single photon
hadro-production in conjunction with two final-state jets at NNLO in QCD. We
employed finite-field arithmetic throughout the computation, and reconstructed
the analytic expression of the amplitudes using an optimised algorithm
informed by the expected singularity structure, and boosted by a partial
fraction decomposition. We optimised the IBP systems with the help of syzygy
equations. The results are explicit for a minimal set of helicity
configurations and partial amplitudes expressed in terms of pentagon functions
suitable for efficient numerical evaluation. The complete hard functions
required for the squared matrix elements are obtained by a combination of
parity conjugation and permutation of the external momenta. Each of these
operations has been described in detail. We have also paid particular
attention to numerical stability, which is achieved with the help of higher-
precision floating-point arithmetic in our implementation in C++. The two-loop
amplitudes obtained here are made available in a public repository
ppajjampsite .

Equipped with the two-loop hard functions, we have performed a
phenomenological analysis of the photon + di-jet production process in the
setup corresponding to recent ATLAS measurements at the LHC. We have
demonstrated that the NNLO results improve the description of the data in the
case of the photon transverse energy, the jet transverse momentum, and the di-
jet invariant mass. Nevertheless, there is a noticeable difference between
theory and experiment in the case of the photon transverse energy starting
around $1\;\text{TeV}$. This difference increases with energy and can be
attributed to electroweak effects not accounted for in our study. This
explanation is consistent with the negative sign of electroweak Sudakov
corrections typical of high-transverse-momentum processes. In the case of the
azimuthal angle between the photon and the leading jet, and the azimuthal
angle between the two leading jets, the description is satisfactory in the
inclusive case, while perturbation theory fails as expected for exclusive cuts
designed to remove fragmentation effects. Finally, we notice that the
theoretical predictions obtained from parton-shower-matched and multi-jet-
merged simulations generated with Sherpa have larger scale uncertainties and
provide a poorer data description than NLO let alone NNLO QCD predictions. One
possible explanation of this unsatisfactory situation could be the reliance on
a multi-jet-merged sample of leading-order accuracy. We do not exclude the
possibility of improving the description by a better setup of the Sherpa
simulation. Should this be achievable, it would demonstrate the pitfalls of
using the very complicated machinery of Monte-Carlo event generators.

Since our phenomenological study is designed to showcase the two-loop two-to-
three virtual corrections, we have not evaluated the fragmentation
contributions. These contributions may be omitted thanks to the special set of
cuts called here the direct-enriched phase space. While we estimate the
systematic uncertainty implied by the omission as lower than 5% and decreasing
with energy, we intend to include fragmentation in a future study. Finally, we
stress that infrared safety of the observables is guaranteed by the use of the
well-known Fixione isolation condition in conjunction with the hard-cone
isolation used in the experimental setup.

Our study not only completes the set of predictions for two-to-three processes
with photons and jets, but is also the first not to rely on the leading-colour
approximation of double-virtual contributions.

###### Acknowledgements.

The authors would like to thank Alexander Mitov for fruitful discussions. The
work of M.C. was supported by the Deutsche Forschungsgemeinschaft under grant
396021762 – TRR 257. H.B.H. is supported by STFC consolidated HEP theory grant
ST/T000694/1. R.P. acknowledges support from the Leverhulme Trust and the
Isaac Newton Trust. This project received funding from the European Research
Council (ERC) under the European Union’s research and innovation programmes
New level of theoretical precision for LHC Run 2 and beyond (grant agreement
No. 683211), High precision multi-jet dynamics at the LHC (grant agreement No.
772099), and High-precision multi-leg Higgs and top physics with finite fields
(grant agreement No. 101040760). This work was performed using the Cambridge
Service for Data Driven Discovery (CSD3), part of which is operated by the
University of Cambridge Research Computing on behalf of the STFC DiRAC HPC
Facility (www.dirac.ac.uk). The DiRAC component of CSD3 was funded by BEIS
capital funding via STFC capital grants ST/P002307/1 and ST/R002452/1 and STFC
operations grant ST/R00689X/1. DiRAC is part of the National e-Infrastructure.
Simulations were performed with computing resources granted by RWTH Aachen
University under project p0020025.

##  Appendix A Numerical benchmark for the hard functions

In this Appendix we provide numerical values of the hard functions for the
$0\to\bar{q}qgg\gamma$ and $0\to q\bar{q}Q\bar{Q}\gamma$ partonic subprocesses
at a benchmark phase-space point. We define the tree-level, one-loop, and two-
loop hard functions as follows

$\displaystyle\mathcal{H}^{(0)}$
$\displaystyle=\underset{\mathrm{colour}}{\overline{\sum}}\,\underset{\mathrm{helicity}}{\overline{\sum}}\big{|}\mathcal{M}^{(0)}\big{|}^{2}\,,$
(39a) $\displaystyle\mathcal{H}^{(1)}$
$\displaystyle=2\,\mathrm{Re}\,\underset{\mathrm{colour}}{\overline{\sum}}\,\underset{\mathrm{helicity}}{\overline{\sum}}\mathcal{M}^{(0)*}\mathcal{R}^{(1)}\,,$
(39b) $\displaystyle\mathcal{H}^{(2)}$
$\displaystyle=2\,\mathrm{Re}\,\underset{\mathrm{colour}}{\overline{\sum}}\,\underset{\mathrm{helicity}}{\overline{\sum}}\mathcal{M}^{(0)*}\mathcal{R}^{(2)}+\underset{\mathrm{colour}}{\overline{\sum}}\,\underset{\mathrm{helicity}}{\overline{\sum}}\mathcal{R}^{(1)*}\mathcal{R}^{(1)}\,.$
(39c)

The factors required for averaging over initial state helicities and colours,
as well as the relevant symmetry factors in the case of identical final state
partons are included in the normalisation. We define the following short-hand
notation

$f_{1}f_{2}\to f_{3}f_{4}\gamma\qquad\Longrightarrow\qquad
f_{1}(-p_{1})+f_{2}(-p_{2})\to f_{3}(p_{3})+f_{4}(p_{4})+\gamma(p_{5})\,,$
(40)

to represent a scattering process involving two incoming partons ($f_{1}$,
$f_{2}$), and two outgoing partons ($f_{3}$, $f_{4}$) in association with a
photon.

In Table 3, we show the values of the tree-level, one-loop, and two-loop hard
functions for a selection of partonic scattering channels that contribute to
$pp\to\gamma jj$ at the following phase-space point,

$\displaystyle p_{1}$ $\displaystyle=(-500,0,0,-500)\,,$ $\displaystyle p_{2}$
$\displaystyle=(-500,0,0,500)\,,$ $\displaystyle p_{3}$
$\displaystyle=(458.57878788544,169.45322030968,379.65366207819,-193.50247465025)\,,$
(41) $\displaystyle p_{4}$
$\displaystyle=(364.06662073682,-18.329869293192,-347.70430131937,106.34960775871)\,,$
$\displaystyle p_{5}$
$\displaystyle=(177.35459137774,-151.12335101649,-31.949360758832,87.152866891544)\,,$

and using the following parameters,

$\mu_{R}=173.2\,,\qquad\qquad\alpha_{s}=0.118\,,\qquad\qquad
e=0.30795376724\,.$ (42)

The momenta and the renormalisation scale $\mu_{R}$ are given in units of GeV.

$0\to\bar{q}qgg\gamma$ | $\mathcal{H}^{(0)}\,\left[\mathrm{GeV}^{-2}\right]$ | $\mathcal{H}^{(1)}/\mathcal{H}^{(0)}$ | $\mathcal{H}^{(2)}/\mathcal{H}^{(0)}$
---|---|---|---
$u\bar{u}\to gg\gamma$ | $1.42202\cdot 10^{-5}$ | -0.0868074 | 0.0219214
$gg\to u\bar{u}\gamma$ | $2.69842\cdot 10^{-6}$ | -0.246582 | 0.0412084
$ug\to ug\gamma$ | $4.53256\cdot 10^{-5}$ | 0.0299060 | 0.0169971
$gu\to ug\gamma$ | $2.70402\cdot 10^{-5}$ | 0.0155916 | 0.00457076
$0\to q\bar{q}Q\bar{Q}\gamma$ | $\mathcal{H}^{(0)}\,\left[\mathrm{GeV}^{-2}\right]$ | $\mathcal{H}^{(1)}/\mathcal{H}^{(0)}$ | $\mathcal{H}^{(2)}/\mathcal{H}^{(0)}$
$u\bar{u}\to u\bar{u}\gamma$ | $4.29458\cdot 10^{-5}$ | 0.0756652 | -0.0388519
$u\bar{u}\to d\bar{d}\gamma$ | $4.41265\cdot 10^{-6}$ | -0.253477 | 0.0493917
$u\bar{u}\to c\bar{c}\gamma$ | $8.04351\cdot 10^{-6}$ | -0.290635 | 0.0616392
$uu\to uu\gamma$ | $1.88278\cdot 10^{-5}$ | 0.171349 | 0.0140249
$ud\to ud\gamma$ | $2.00840\cdot 10^{-5}$ | 0.233819 | -0.0167785
$uc\to uc\gamma$ | $6.53549\cdot 10^{-6}$ | 0.228568 | -0.0168878
Table 3: Numerical results for the hard functions defined in Eq. 39 for a
selection of partonic reactions contributing to NNLO QCD computation of
$pp\to\gamma jj$ production at the LHC. We employ the phase-space point given
in Eq. 41 and parameters specified in Eq. 42. The notation used for the
partonic scattering process is defined in Eq. 40.

Numerical results for the hard functions for all the scattering channels
contributing to $pp\to\gamma jj$, as well as all partial helicity finite
remainders for all contributing helicity configurations can be accessed by
running the scripts `SquaredFiniteRemainder_2q2gA.wl` and
`SquaredFiniteRemainder_2q2QA.wl` available in the public repository
ppajjampsite . These two files also provide all the necessary rules to perform
parity transformation and permutations to obtain the full set of helicity
configurations, partial amplitudes and scattering channels.

##  Appendix B $\mu_{R}$ dependence of the finite remainder

In the computation of the finite remainders $\mathcal{R}^{(L)}$ we set
$\mu_{R}=1$. We recover the dependence on $\mu_{R}$ a posteriori, by adding
suitable $\mu_{R}$-restoring terms $\delta\mathcal{R}^{(L)}$ to the finite
remainders,

$\displaystyle\mathcal{R}^{(L)}(\vec{s},\mu_{R}^{2})=\mathcal{R}^{(L)}(\vec{s},\mu_{R}^{2}=1)+\delta\mathcal{R}^{(L)}(\vec{s},\mu_{R}^{2})\,.$
(43)

$\vec{s}$ is the set of kinematic invariants defined in Eq. 2. The
$\mu_{R}$-restoring terms depend only on the $\beta$ function, the IR
operators $\textbf{Z}^{(L)}$, and the lower-loop finite remainders evaluated
at $\mu_{R}^{2}=1$. Computing the finite remainders with $\mu_{R}^{2}=1$ is
thus sufficient to recover the full dependence. Explicitly, the one- and two-
loop $\mu_{R}$-restoring terms are given by

$\displaystyle|\delta\mathcal{R}^{(1)}(\vec{s},\mu_{R}^{2})\rangle=\left[\frac{1}{2}\log^{2}(\mu_{R}^{2})\,\textbf{Z}^{(1)}_{-2}+\log(\mu_{R}^{2})\left(\beta_{0}\,\textbf{I}+\textbf{Z}^{(1)}_{-1}(\vec{s})\right)\right]|\mathcal{M}^{(0)}(\vec{s})\rangle\,,$
(44)
$\displaystyle|\delta\mathcal{R}^{(2)}(\vec{s},\mu_{R}^{2})\rangle=\left[\frac{1}{2}\log^{2}(\mu_{R}^{2})\,\textbf{Z}^{(1)}_{-2}+\log(\mu_{R}^{2})\left(2\beta_{0}\,\textbf{I}+\textbf{Z}^{(1)}_{-1}(\vec{s})\right)\right]|\mathcal{R}^{(1)}(\vec{s},\mu_{R}^{2}=1)\rangle$
$\displaystyle\phantom{|\delta\mathcal{R}^{(2)}(\vec{s},\mu_{R}^{2})\rangle=}+\biggl{\\{}\log^{4}(\mu_{R}^{2})\left[-\frac{5}{24}\left(\textbf{Z}^{(1)}_{-2}\right)^{2}+\frac{2}{3}\,\textbf{Z}^{(2)}_{-4}\right]$
$\displaystyle\phantom{|\delta\mathcal{R}^{(2)}(\vec{s},\mu_{R}^{2})\rangle=+\biggl{\\{}}+\frac{1}{6}\log^{3}(\mu_{R}^{2})\left(10\beta_{0}\,\textbf{Z}^{(1)}_{-2}-5\,\textbf{Z}^{(1)}_{-2}\cdot\textbf{Z}^{(1)}_{-1}(\vec{s})+8\,\textbf{Z}^{(2)}_{-3}(\vec{s})\right)$
$\displaystyle\phantom{|\delta\mathcal{R}^{(2)}(\vec{s},\mu_{R}^{2})\rangle=+\biggl{\\{}}+\log^{2}(\mu_{R}^{2})\left[\beta_{0}^{2}\,\textbf{I}+\frac{5}{2}\beta_{0}\,\textbf{Z}^{(1)}_{-1}(\vec{s})-\frac{1}{2}\left(\textbf{Z}^{(1)}_{-1}(\vec{s})\right)^{2}+2\,\textbf{Z}^{(2)}_{-2}(\vec{s})\right]+$
$\displaystyle\phantom{|\delta\mathcal{R}^{(2)}(\vec{s},\mu_{R}^{2})\rangle=+\biggl{\\{}}+\log(\mu_{R}^{2})\left(\beta_{1}\,\textbf{I}+2\,\textbf{Z}^{(2)}_{-1}(\vec{s})\right)\biggr{\\}}|\mathcal{M}^{(0)}(\vec{s})\rangle\,,$
(45)

where I is the identity operator in colour space, and $\textbf{Z}^{(L)}_{k}$
is the order-$\epsilon^{k}$ term in the Laurent expansion around $\epsilon=0$
of the IR operators at $\mu_{R}^{2}=1$,

$\displaystyle\textbf{Z}^{(L)}(\vec{s},\mu_{R}^{2}=1)=\sum_{w=-2L}^{-1}\textbf{Z}^{(L)}_{w}(\vec{s})\,\epsilon^{w}\,.$
(46)

Note that — in the scheme we adopted — the IR operators do not contain finite
terms, i.e. they have no terms of order $\epsilon^{0}$ or higher. We further
recall that the IR operators $\textbf{Z}^{(L)}$ are process-dependent, and
that the leading poles — $\textbf{Z}^{(L)}_{-2L}$ — are constant. In Eqs. 44
and B we spelled out the dependence on the kinematic invariants $\vec{s}$, and
all logarithms of $\mu_{R}^{2}$ are explicit.

We verify that the finite remainders have the correct scale dependence by
evaluating them at two phase-space points related by a rescaling by some
factor $a>0$, and checking that they satisfy the relation

$\displaystyle\frac{\mathcal{R}^{(L)}(a\,\vec{s},1)+\delta\mathcal{R}^{(L)}(a\,\vec{s},a)}{\mathcal{R}^{(L)}(\vec{s},1)}=\frac{\mathcal{M}^{(0)}(a\,\vec{s})}{\mathcal{M}^{(0)}(\vec{s})}\,.$
(47)

This follows from the fact that amplitudes and finite remainders are
homogeneous functions of $\vec{s}$ and $\mu_{R}^{2}$.

##  Appendix C Parity transformation and external-momentum permutation for the
amplitudes

We derived analytic expressions for the independent helicity configurations
and the minimal set of partial amplitudes. In order to obtain the expressions
for the complete set of helicity configurations and partial amplitudes
required to evaluate the squared matrix elements in the various scattering
channels, we apply a parity transformation and/or a permutation of the
external momenta to the minimal set of partial amplitudes. In this Appendix we
spell out the algorithms required to perform these operations. We first define
the following short-hand notation for the partial amplitudes:

$\mathcal{A}^{\vec{h}}(g)=\mathcal{A}^{(L)}\left(1^{h_{1}},2^{h_{2}},3^{h_{3}},4^{h_{4}},5^{h_{5}}_{\gamma}\right)\,,$
(48)

where $g=(12345)$ is the default external-momentum ordering specified in Eqs.
5 and 17, and $\vec{h}=(h_{1},h_{2},h_{3},h_{4},h_{5})$. Our analytic
expressions for the loop partial amplitudes have the form

$\mathcal{A}^{\vec{h}}(g)=\sum_{k}c^{\vec{h}}_{k}(\vec{x}(g))\,\mathrm{mon}_{k}(f(g),\mathfrak{c})\,.$
(49)

Here, we extend the notation for the momentum-twistor variables $\vec{x}$ and
pentagon functions $f$ to show explicitly their dependence on the routing of
the external momenta: $\vec{x}(g)$ and $f(g)$. $\vec{x}(g)$ indicates an
evaluation of the right-hand-side of Eq. 28 at a given phase-space point with
external-momentum ordering $g$. Since the spinor-phase information is lost
when using momentum-twistor variables, we restore it in our analytic
expressions for the partial amplitudes through suitable helicity-dependent
spinor-phase factors $\Phi^{\vec{h}}$ (defined in Eqs. 11 and 26 for
$0\to\bar{q}qgg\gamma$ and $0\to q\bar{q}Q\bar{Q}\gamma$, respectively). In
other words, we divide the coefficients $c_{k}^{\vec{h}}$ in Eq. 49 by the
phase factor expressed in terms of momentum-twistor variables
$\Phi^{\vec{h}}(\vec{x})$, and multiply them by the same expressed in terms of
spinor products $\Phi^{\vec{h}}(\langle ij\rangle,[ij])$. We thus have

$\mathcal{A}^{\vec{h}}(g)=\Phi^{\vec{h}}(\langle
ij\rangle,[ij];g)\,\sum_{k}\frac{c^{\vec{h}}_{k}(\vec{x}(g))}{\Phi^{\vec{h}}(\vec{x}(g))}\,\mathrm{mon}_{k}(f(g),\mathfrak{c})\,,$
(50)

where we spell out the dependence on the external-momentum ordering $g$. The
benefit of this approach is that the momentum-twistor dependent coefficients
$c^{\vec{h}}_{k}(\vec{x}(g))/\Phi^{\vec{h}}(\vec{x}(g))$ are phase-free
quantities, and the phase of the partial amplitude is explicitly restored.
This allows us to perform the following two operations.

  * •

Parity transformation
The action of parity $\mathbf{P}$ flips the helicities and the sign of all
pseudo-scalar quantities. In particular, it flips the sign of
$\operatorname{\rm tr}_{5}$ and thus exchanges $\operatorname{\rm tr}_{+}$
with $\operatorname{\rm tr}_{-}$. Thus, there is

$\displaystyle\mathcal{A}^{-\vec{h}}(g)$
$\displaystyle=\mathbf{P}\,\mathcal{A}^{\vec{h}}(g)$
$\displaystyle=\Phi^{\vec{h}}(\langle ij\rangle,[ij];g)\big{|}_{\langle
ij\rangle\leftrightarrow[ji]}\,\sum_{k}\frac{c^{\vec{h}}_{k}(\vec{x}(g))}{\Phi^{\vec{h}}(\vec{x}(g))}\bigg{|}_{\operatorname{\rm
tr}_{+}\leftrightarrow\operatorname{\rm
tr}_{-}}\,\mathrm{mon}_{k}\left(\mathbf{P}f(g),\mathbf{P}\mathfrak{c}\right)\,.$
(51)

We recall from Section 2.5 that parity acts on the pentagon functions $f$ and
the associated transcendental constants $\mathfrak{c}$ by flipping the sign of
all parity-odd functions/constants.

Complex conjugation needed in order to construct the squared matrix elements
(Eqs. 12 and 2.3) requires that — on top of the parity transformation — we
also complex-conjugate the values of the pentagon functions and constants,

$\displaystyle\left(\mathcal{A}^{\vec{h}}(g)\right)^{*}$
$\displaystyle=\Phi^{\vec{h}}(\langle ij\rangle,[ij];g)\big{|}_{\langle
ij\rangle\leftrightarrow[ji]}\,\sum_{k}\frac{c^{\vec{h}}_{k}(\vec{x}(g))}{\Phi^{\vec{h}}(\vec{x}(g))}\bigg{|}_{\operatorname{\rm
tr}_{+}\leftrightarrow\operatorname{\rm
tr}_{-}}\,\mathrm{mon}_{k}\left(\mathbf{P}(f(g))^{*},\mathbf{P}(\mathfrak{c})^{*}\right)\,.$
(52)

  * •

Permutation of the external momenta
Let $g^{\prime}$ be an ordering of the external momenta different from
$g=(12345)$. We denote by $\sigma$ the permutation of the external momenta
which takes from $g$ to $g^{\prime}$, and by $\vec{h}^{\prime}$ the permuted
helicities,

$g^{\prime}=\sigma\circ g\,,\qquad\vec{h}^{\prime}=\sigma\circ\vec{h}\,.$ (53)

The partial amplitude for the ordering $g^{\prime}$ is given by

$\displaystyle\mathcal{A}^{\vec{h^{\prime}}}(g^{\prime})$
$\displaystyle=\Phi^{\vec{h}}(\langle
ij\rangle,[ij];g^{\prime})\,\sum_{k}\frac{c^{\vec{h}}_{k}(\vec{x}(g^{\prime}))}{\Phi^{\vec{h}}(\vec{x}(g^{\prime}))}\,\mathrm{mon}_{k}\left(\sigma\circ
f(g),\sigma\circ\mathfrak{c}\right)\,.$ (54)

We emphasise that some of the constants $\mathfrak{c}$ have odd parity, and
change sign if the permutation $\sigma$ has odd signature. A more thorough
discussion of the permutation of pentagon functions and transcendental
constants can be found in Section 2.5.

## References

  * (1) K. Koller, T. F. Walsh and P. M. Zerwas, _Testing QCD: Direct Photons in $e^{+}e^{-}$ Collisions_, _Z. Phys. C_ 2 (1979) 197.
  * (2) E. Laermann, T. F. Walsh, I. Schmitt and P. M. Zerwas, _Direct Photons in $e^{+}e^{-}$ Annihilation_, _Nucl. Phys. B_ 207 (1982) 205.
  * (3) S. Catani, M. Fontannaz, J. P. Guillet and E. Pilon, _Cross-section of isolated prompt photons in hadron hadron collisions_ , _JHEP_ 05 (2002) 028 [hep-ph/0204023].
  * (4) T. Gehrmann and R. Schürmann, _Photon fragmentation in the antenna subtraction formalism_ , _JHEP_ 04 (2022) 031 [2201.06982].
  * (5) X. Chen, T. Gehrmann, E. W. N. Glover, M. Höfer, A. Huss and R. Schürmann, _Single photon production at hadron colliders at NNLO QCD with realistic photon isolation_ , _JHEP_ 08 (2022) 094 [2205.01516].
  * (6) E. W. N. Glover and A. G. Morgan, _Measuring the photon fragmentation function at LEP_ , _Z. Phys. C_ 62 (1994) 311.
  * (7) ALEPH collaboration, _First measurement of the quark to photon fragmentation function_ , _Z. Phys. C_ 69 (1996) 365.
  * (8) L. Bourhis, M. Fontannaz and J. P. Guillet, _Quarks and gluon fragmentation functions into photons_ , _Eur. Phys. J. C_ 2 (1998) 529 [hep-ph/9704447].
  * (9) OPAL collaboration, _Measurement of the quark to photon fragmentation function through the inclusive production of prompt photons in hadronic Z0 decays_ , _Eur. Phys. J. C_ 2 (1998) 39 [hep-ex/9708020].
  * (10) Z. Belghobsi, M. Fontannaz, J. P. Guillet, G. Heinrich, E. Pilon and M. Werlen, _Photon - Jet Correlations and Constraints on Fragmentation Functions_ , _Phys. Rev. D_ 79 (2009) 114024 [0903.4834].
  * (11) T. Kaufmann, A. Mukherjee and W. Vogelsang, _Recent developments on parton-to-photon fragmentation functions_ , _CERN Proc._ 1 (2018) 211 [1708.06683].
  * (12) S. Frixione, _Isolated photons in perturbative QCD_ , _Phys. Lett. B_ 429 (1998) 369 [hep-ph/9801442].
  * (13) T. Gehrmann, N. Greiner and G. Heinrich, _Photon isolation effects at NLO in $\gamma\gamma$ \+ jet final states in hadronic collisions_, _JHEP_ 06 (2013) 058 [1303.0824].
  * (14) X. Chen, T. Gehrmann, N. Glover, M. Höfer and A. Huss, _Isolated photon and photon+jet production at NNLO QCD accuracy_ , _JHEP_ 04 (2020) 166 [1904.01044].
  * (15) T. Gehrmann, N. Glover, A. Huss and J. Whitehead, _Scale and isolation sensitivity of diphoton distributions at the LHC_ , _JHEP_ 01 (2021) 108 [2009.11310].
  * (16) F. Siegert, _A practical guide to event generation for prompt photon production with Sherpa_ , _J. Phys. G_ 44 (2017) 044007 [1611.07226].
  * (17) S. Catani, M. Fontannaz, J. P. Guillet and E. Pilon, _Isolating Prompt Photons with Narrow Cones_ , _JHEP_ 09 (2013) 007 [1306.6498].
  * (18) M. Balsiger, T. Becher and D. Y. Shao, _Non-global logarithms in jet and isolation cone cross sections_ , _JHEP_ 08 (2018) 104 [1803.07045].
  * (19) T. Becher, S. Favrod and X. Xu, _QCD anatomy of photon isolation_ , _JHEP_ 01 (2023) 005 [2208.01554].
  * (20) Athens-Athens-Brookhaven-CERN collaboration, _Direct Production of High $p_{T}$ Single Photons at the CERN Intersecting Storage Rings_, _Phys. Lett. B_ 87 (1979) 292.
  * (21) E. Anassontzis et al., _High p(t) Direct Photon Production in p p Collisions_ , _Z. Phys. C_ 13 (1982) 277.
  * (22) CMOR collaboration, _Direct Photon Production at the CERN ISR_ , _Nucl. Phys. B_ 327 (1989) 541.
  * (23) UA1 collaboration, _Direct Photon Production at the CERN Proton - anti-Proton Collider_ , _Phys. Lett. B_ 209 (1988) 385.
  * (24) UA2 collaboration, _A Measurement of the direct photon production cross-section at the CERN anti-p p collider_ , _Phys. Lett. B_ 263 (1991) 544.
  * (25) WA70 collaboration, _Production of High Transverse Momentum Prompt Photons and Neutral Pions in Proton Proton Interactions at 280-GeV/c_ , _Z. Phys. C_ 38 (1988) 371.
  * (26) Fermilab E706 collaboration, _Measurement of direct photon production at Tevatron fixed target energies_ , _Phys. Rev. D_ 70 (2004) 092009 [hep-ex/0407011].
  * (27) D0 collaboration, _Measurement of the isolated photon cross section in $p\bar{p}$ collisions at $\sqrt{s}$ = 1.96-TeV_, _Phys. Lett. B_ 639 (2006) 151 [hep-ex/0511054].
  * (28) D0 collaboration, _Measurement of the Differential Cross-Section for the Production of an Isolated Photon with Associated Jet in $p\bar{p}$ Collisions at $\sqrt{s}$ = 1.96-TeV_, _Phys. Lett. B_ 666 (2008) 435 [0804.1107].
  * (29) D0 collaboration, _Measurement of the differential cross section of photon plus jet production in $p\bar{p}$ collisions at $\sqrt{s}=1.96$ TeV_, _Phys. Rev. D_ 88 (2013) 072008 [1308.2708].
  * (30) CDF collaboration, _Measurement of the inclusive-isolated prompt-photon cross section in $p\bar{p}$ collisions using the full CDF data set_, _Phys. Rev. D_ 96 (2017) 092003 [1703.00599].
  * (31) ATLAS collaboration, _Measurement of the cross section for inclusive isolated-photon production in $pp$ collisions at $\sqrt{s}=13$ TeV using the ATLAS detector_, _Phys. Lett. B_ 770 (2017) 473 [1701.06882].
  * (32) ATLAS collaboration, _Measurement of the cross section for isolated-photon plus jet production in $pp$ collisions at $\sqrt{s}=13$ TeV using the ATLAS detector_, _Phys. Lett. B_ 780 (2018) 578 [1801.00112].
  * (33) ATLAS collaboration, _Measurement of the inclusive isolated-photon cross section in $pp$ collisions at $\sqrt{s}=13$ TeV using 36 fb-1 of ATLAS data_, _JHEP_ 10 (2019) 203 [1908.02746].
  * (34) CMS collaboration, _Measurement of differential cross sections for inclusive isolated-photon and photon+jets production in proton-proton collisions at $\sqrt{s}=$ 13 TeV_, _Eur. Phys. J. C_ 79 (2019) 20 [1807.00782].
  * (35) CMS collaboration, _Measurements of triple-differential cross sections for inclusive isolated-photon+jet events in pp collisions at $\sqrt{s}=8\,\text{TeV}$_, _Eur. Phys. J. C_ 79 (2019) 969 [1907.08155].
  * (36) ALICE collaboration, _Measurement of the inclusive isolated photon production cross section in pp collisions at $\sqrt{s}=7$ TeV_, _Eur. Phys. J. C_ 79 (2019) 896 [1906.01371].
  * (37) ATLAS collaboration, _Measurement of isolated-photon plus two-jet production in $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector_, _JHEP_ 03 (2020) 179 [1912.09866].
  * (38) S. Keller and J. F. Owens, _Event structure in photon plus two jet final states_ , _Phys. Lett. B_ 269 (1991) 445.
  * (39) J. M. Lindert et al., _Precise predictions for $V+$ jets dark matter backgrounds_, _Eur. Phys. J. C_ 77 (2017) 829 [1705.04664].
  * (40) T. Gleisberg, S. Hoeche, F. Krauss, M. Schonherr, S. Schumann, F. Siegert et al., _Event generation with SHERPA 1.1_ , _JHEP_ 02 (2009) 007 [0811.4622].
  * (41) Sherpa collaboration, _Event Generation with Sherpa 2.2_ , _SciPost Phys._ 7 (2019) 034 [1905.09127].
  * (42) 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].
  * (43) S. Kallweit, V. Sotnikov and M. Wiesemann, _Triphoton production at hadron colliders in NNLO QCD_ , _Phys. Lett. B_ 812 (2021) 136013 [2010.04681].
  * (44) H. A. Chawdhry, M. Czakon, A. Mitov and R. Poncelet, _NNLO QCD corrections to diphoton production with an additional jet at the LHC_ , _JHEP_ 09 (2021) 093 [2105.06940].
  * (45) M. Czakon, A. Mitov and R. Poncelet, _Next-to-Next-to-Leading Order Study of Three-Jet Production at the LHC_ , _Phys. Rev. Lett._ 127 (2021) 152001 [2106.05331].
  * (46) S. Badger, T. Gehrmann, M. Marcoli and R. Moodie, _Next-to-leading order QCD corrections to diphoton-plus-jet production through gluon fusion at the LHC_ , _Phys. Lett. B_ 824 (2022) 136802 [2109.12003].
  * (47) H. B. Hartanto, R. Poncelet, A. Popescu and S. Zoia, _Next-to-next-to-leading order QCD corrections to Wbb¯ production at the LHC_ , _Phys. Rev. D_ 106 (2022) 074016 [2205.01687].
  * (48) H. B. Hartanto, R. Poncelet, A. Popescu and S. Zoia, _Flavour anti- $k_{\text{T}}$ algorithm applied to $Wb\bar{b}$ production at the LHC_, 2209.03280.
  * (49) L. Buonocore, S. Devoto, S. Kallweit, J. Mazzitelli, L. Rottoli and C. Savoini, _Associated production of a $W$ boson and massive bottom quarks at next-to-next-to-leading order in QCD_, 2212.04954.
  * (50) M. Alvarez, J. Cantero, M. Czakon, J. Llorente, A. Mitov and R. Poncelet, _NNLO QCD corrections to event shapes at the LHC_ , 2301.01086.
  * (51) 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].
  * (52) S. Abreu, F. Febres Cordero, H. Ita, B. Page and V. Sotnikov, _Leading-color two-loop QCD corrections for three-jet production at hadron colliders_ , _JHEP_ 07 (2021) 095 [2102.13609].
  * (53) 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 [2010.15834].
  * (54) 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 [2012.13553].
  * (55) B. Agarwal, F. Buccioni, A. von Manteuffel and L. Tancredi, _Two-loop leading colour QCD corrections to $q\bar{q}\to\gamma\gamma g$ and $qg\to\gamma\gamma q$_, _JHEP_ 04 (2021) 201 [2102.01820].
  * (56) H. A. Chawdhry, M. Czakon, A. Mitov and R. Poncelet, _Two-loop leading-colour QCD helicity amplitudes for two-photon plus jet production at the LHC_ , _JHEP_ 07 (2021) 164 [2103.04319].
  * (57) 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 [1812.11160].
  * (58) D. Chicherin and V. Sotnikov, _Pentagon Functions for Scattering of Five Massless Particles_ , _JHEP_ 20 (2020) 167 [2009.07803].
  * (59) H. B. Hartanto, S. Badger, C. Brønnum-Hansen and T. Peraro, _A numerical evaluation of planar two-loop helicity amplitudes for a W-boson plus four partons_ , _JHEP_ 09 (2019) 119 [1906.11862].
  * (60) S. Badger, C. Brønnum-Hansen, D. Chicherin, T. Gehrmann, H. B. Hartanto, J. Henn et al., _Virtual QCD corrections to gluon-initiated diphoton plus jet production at hadron colliders_ , _JHEP_ 11 (2021) 083 [2106.08664].
  * (61) S. Badger, H. B. Hartanto and S. Zoia, _Two-Loop QCD Corrections to Wbb¯ Production at Hadron Colliders_ , _Phys. Rev. Lett._ 127 (2021) 012001 [2102.02516].
  * (62) S. Badger, H. B. Hartanto, J. Kryś and S. Zoia, _Two-loop leading-colour QCD helicity amplitudes for Higgs boson production in association with a bottom-quark pair at the LHC_ , _JHEP_ 11 (2021) 012 [2107.14733].
  * (63) S. Badger, H. B. Hartanto, J. Kryś and S. Zoia, _Two-loop leading colour helicity amplitudes for W ±$\gamma$ \+ j production at the LHC_, _JHEP_ 05 (2022) 035 [2201.04075].
  * (64) B. Agarwal, F. Buccioni, A. von Manteuffel and L. Tancredi, _Two-Loop Helicity Amplitudes for Diphoton Plus Jet Production in Full Color_ , _Phys. Rev. Lett._ 127 (2021) 262001 [2105.04585].
  * (65) F. Buccioni, J.-N. Lang, J. M. Lindert, P. Maierhöfer, S. Pozzorini, H. Zhang et al., _OpenLoops 2_ , _Eur. Phys. J. C_ 79 (2019) 866 [1907.13071].
  * (66) S. Catani, _The Singular behavior of QCD amplitudes at two loop order_ , _Phys. Lett._ B427 (1998) 161 [hep-ph/9802439].
  * (67) T. Becher and M. Neubert, _On the Structure of Infrared Singularities of Gauge-Theory Amplitudes_ , _JHEP_ 06 (2009) 081 [0903.1126].
  * (68) T. Becher and M. Neubert, _Infrared singularities of scattering amplitudes in perturbative QCD_ , _Phys. Rev. Lett._ 102 (2009) 162001 [0901.0722].
  * (69) E. Gardi and L. Magnea, _Factorization constraints for soft anomalous dimensions in QCD scattering amplitudes_ , _JHEP_ 03 (2009) 079 [0901.1091].
  * (70) S. Badger, M. Czakon, H. B. Hartanto, R. Moodie, T. Peraro, R. Poncelet et al., _Amplitudes for “Isolated photon production in association with a jet pair through next-to-next- to-leading order in QCD”_ , Apr., 2023.  10.5281/zenodo.7824963.
  * (71) A. von Manteuffel and R. M. Schabinger, _A novel approach to integration by parts reduction_ , _Phys. Lett. B_ 744 (2015) 101 [1406.4513].
  * (72) T. Peraro, _Scattering amplitudes over finite fields and multivariate functional reconstruction_ , _JHEP_ 12 (2016) 030 [1608.01902].
  * (73) T. Peraro, _FiniteFlow: multivariate functional reconstruction using finite fields and dataflow graphs_ , _JHEP_ 07 (2019) 031 [1905.08019].
  * (74) S. Badger, E. Chaubey, H. B. Hartanto and R. Marzucca, _Two-loop leading colour QCD helicity amplitudes for top quark pair production in the gluon fusion channel_ , _JHEP_ 06 (2021) 163 [2102.13450].
  * (75) P. Nogueira, _Automatic Feynman graph generation_ , _J. Comput. Phys._ 105 (1993) 279.
  * (76) A. Hodges, _Eliminating spurious poles from gauge-theoretic amplitudes_ , _JHEP_ 05 (2013) 135 [0905.1473].
  * (77) S. Badger, H. Frellesvig and Y. Zhang, _A Two-Loop Five-Gluon Helicity Amplitude in QCD_ , _JHEP_ 12 (2013) 045 [1310.1051].
  * (78) J. Kuipers, T. Ueda, J. A. M. Vermaseren and J. Vollinga, _FORM version 4.0_ , _Comput. Phys. Commun._ 184 (2013) 1453 [1203.6543].
  * (79) B. Ruijl, T. Ueda and J. Vermaseren, _FORM version 4.2_ , 1707.06453.
  * (80) G. Cullen, M. Koch-Janusz and T. Reiter, _Spinney: A Form Library for Helicity Spinors_ , _Comput. Phys. Commun._ 182 (2011) 2368 [1008.0803].
  * (81) J. M. Henn, _Multiloop integrals in dimensional regularization made simple_ , _Phys. Rev. Lett._ 110 (2013) 251601 [1304.1806].
  * (82) C. G. Papadopoulos, D. Tommasini and C. Wever, _The Pentabox Master Integrals with the Simplified Differential Equations approach_ , _JHEP_ 04 (2016) 078 [1511.09404].
  * (83) T. Gehrmann, J. M. Henn and N. A. Lo Presti, _Pentagon functions for massless planar scattering amplitudes_ , _JHEP_ 10 (2018) 103 [1807.09812].
  * (84) S. Abreu, B. Page and M. Zeng, _Differential equations from unitarity cuts: nonplanar hexa-box integrals_ , _JHEP_ 01 (2019) 006 [1807.11522].
  * (85) D. Chicherin, T. Gehrmann, J. M. Henn, N. A. Lo Presti, V. Mitev and P. Wasser, _Analytic result for the nonplanar hexa-box integrals_ , _JHEP_ 03 (2019) 042 [1809.06240].
  * (86) F. V. Tkachov, _A Theorem on Analytical Calculability of Four Loop Renormalization Group Functions_ , _Phys. Lett._ 100B (1981) 65.
  * (87) K. G. Chetyrkin and F. V. Tkachov, _Integration by Parts: The Algorithm to Calculate $\beta$ Functions in 4 Loops_, _Nucl. Phys._ B192 (1981) 159.
  * (88) J. Gluza, K. Kajda and D. A. Kosower, _Towards a Basis for Planar Two-Loop Integrals_ , _Phys. Rev. D_ 83 (2011) 045012 [1009.0472].
  * (89) H. Ita, _Two-loop Integrand Decomposition into Master Integrals and Surface Terms_ , _Phys. Rev. D_ 94 (2016) 116015 [1510.05626].
  * (90) K. J. Larsen and Y. Zhang, _Integration-by-parts reductions from unitarity cuts and algebraic geometry_ , _Phys. Rev. D_ 93 (2016) 041701 [1511.01071].
  * (91) P. A. Baikov, _Explicit solutions of the three loop vacuum integral recurrence relations_ , _Phys. Lett. B_ 385 (1996) 404 [hep-ph/9603267].
  * (92) P. A. Baikov, _Explicit solutions of n loop vacuum integral recurrence relations_ , hep-ph/9604254.
  * (93) J. Böhm, A. Georgoudis, K. J. Larsen, M. Schulze and Y. Zhang, _Complete sets of logarithmic vector fields for integration-by-parts identities of Feynman integrals_ , _Phys. Rev. D_ 98 (2018) 025023 [1712.09737].
  * (94) A. von Manteuffel and R. M. Schabinger, _Quark and gluon form factors in four loop QCD: The $N_{f}^{2}$ and $N_{q\gamma}N_{f}$ contributions_, _Phys. Rev. D_ 99 (2019) 094014 [1902.08208].
  * (95) S. Laporta, _High precision calculation of multiloop Feynman integrals by difference equations_ , _Int. J. Mod. Phys._ A15 (2000) 5087 [hep-ph/0102033].
  * (96) 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].
  * (97) D. Chicherin, J. Henn and V. Mitev, _Bootstrapping pentagon functions_ , _JHEP_ 05 (2018) 164 [1712.09610].
  * (98) X. Liu and Y.-Q. Ma, _AMFlow: A Mathematica package for Feynman integrals computation via auxiliary mass flow_ , _Comput. Phys. Commun._ 283 (2023) 108565 [2201.11669].
  * (99) X. S. L. Y. Hida and D. H. Bailey, “libqd: quad-double / double-double computation package.” www.davidhbailey.com/dhbsoftware, 2010\.
  * (100) M. Czakon, _A novel subtraction scheme for double-real radiation at NNLO_ , _Phys. Lett. B_ 693 (2010) 259 [1005.0274].
  * (101) M. Czakon and D. Heymes, _Four-dimensional formulation of the sector-improved residue subtraction scheme_ , _Nucl. Phys. B_ 890 (2014) 152 [1408.2500].
  * (102) 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].
  * (103) M. Bury and A. van Hameren, _Numerical evaluation of multi-gluon amplitudes for High Energy Factorization_ , _Comput. Phys. Commun._ 196 (2015) 592 [1503.08612].
  * (104) NNPDF collaboration, _Parton distributions from high-precision collider data_ , _Eur. Phys. J. C_ 77 (2017) 663 [1706.00428].
  * (105) A. Buckley, J. Ferrando, S. Lloyd, K. Nordström, B. Page, M. Rüfenacht et al., _LHAPDF6: parton density access in the LHC precision era_ , _Eur. Phys. J. C_ 75 (2015) 132 [1412.7420].
  * (106) M. Cacciari, G. P. Salam and G. Soyez, _The anti- $k_{t}$ jet clustering algorithm_, _JHEP_ 04 (2008) 063 [0802.1189].
  * (107) J. H. Kuhn, A. Kulesza, S. Pozzorini and M. Schulze, _Electroweak corrections to hadronic photon production at large transverse momenta_ , _JHEP_ 03 (2006) 059 [hep-ph/0508253].
  * (108) T. Becher and X. Garcia i Tormo, _Electroweak Sudakov effects in $W,Z$ and $\gamma$ production at large transverse momentum_, _Phys. Rev. D_ 88 (2013) 013009 [1305.4202].
  * (109) T. Becher and X. Garcia i Tormo, _Addendum: Electroweak Sudakov effects in W, Z and gamma production at large transverse momentum_ , 1509.01961.
  * (110) ATLAS collaboration, _Measurement of dijet azimuthal decorrelations in $pp$ collisions at $\sqrt{s}=8$ TeV with the ATLAS detector and determination of the strong coupling_, _Phys. Rev. D_ 98 (2018) 092004 [1805.04691].