aainstitutetext: Physik-Institut, Universität Zurich, Winterthurerstrasse 190,
CH-8057 Zürich, Switzerlandbbinstitutetext: Physik Department, Technische
Universität München, James-Franck-Straße 1, 85748 Garching, Germany

# Two-loop helicity amplitudes for $H+$jet production to higher orders in the
dimensional regulator

Thomas Gehrmann a Petr Jakubčík b Cesare Carlo Mella b Nikolaos Syrrakos b
Lorenzo Tancredi<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In view of the forthcoming High-Luminosity phase of the LHC, next-to-next-to-
next-to-leading (N3LO) calculations for the most phenomenologically relevant
processes become necessary. In this work, we take the first step towards this
goal for H$+$jet production by computing the one- and two-loop helicity
amplitudes for the two contributing processes, $H\to ggg$, $H\to q\bar{q}g$,
in an effective theory with infinite top quark mass, to higher orders in the
dimensional regulator. We decompose the amplitude in scalar form factors
related to the helicity amplitudes and in a new basis of tensorial structures.
The form factors receive contributions from Feynman integrals which were
reduced to a novel canonical basis of master integrals. We derive and solve a
set of differential equations for these integrals in terms of Multiple
Polylogarithms (MPLs) of two variables up to transcendental weight six.

††preprint: ZU-TH 07/23, TUM-HEP 1451/23

##  1 Introduction

A little over a decade ago, the Higgs boson was discovered after analysing the
data collected during Run I at the Large Hadron Collider (LHC) at CERN
ATLAS:2012yve ; CMS:2012qbp . The discovery was achieved while the collider
was running at reduced centre-of-mass energies of 7 and 8 TeV and with only a
small fraction of the total dataset which will be accumulated during its
entire runtime. Indeed, it is expected that the forthcoming High-Luminosity
phase of the LHC (HL-LHC) will yield a dataset corresponding to 3 $ab^{-1}$ of
integrated luminosity for $pp$ collisions at 14 TeV Dainese:2019rgk .

Since its discovery, the Higgs boson has been at the centre of the
experimental effort at the LHC Heinrich:2020ybq . Studying its properties
improves our understanding of electroweak symmetry breaking (EWSB), the
mechanism which is believed to be responsible for the generation of the masses
of fermions and weak gauge bosons. The dominant production channel for the
Higgs boson at the LHC is gluon fusion. In the Standard Model, the Higgs
coupling to two gluons is mainly mediated through a loop of top quarks, making
it a loop-induced process already at Leading Order (LO). For this reason,
computing higher order perturbative corrections to Higgs production in the
full theory quickly becomes prohibitive.

It was realized long ago that radiative corrections can increase the LO Higgs
cross-section in gluon fusion by as much as $\mathcal{O}(100\%)$ Dawson:1990zj
; Spira:1995rr  and describing its production in hadron collisions thus
requires higher-order calculations. Indeed, quite recently the computation of
the fully inclusive Higgs cross-section with full dependence on the top quark
mass has been pushed to next-to-next-to-leading-order (NNLO) Czakon:2021yub
using numerical techniques to handle the required two- and three-loop
scattering amplitudes.

A useful alternative to performing calculations with full dependence on the
top mass is to work in the heavy top quark mass limit $M_{t}\to\infty$. Under
the assumption that the top quark is the largest scale involved in the
calculation, one can integrate out the top mass and formulate an effective
Lagrangian for the $Hgg$ coupling Wilczek:1977 ; Shifman:1978 ; Inami:1982xt .
In this description, the top-quark loop mediating the $Hgg$ interaction
shrinks to a point and calculations start at tree level, involving only
massless partons. In this limit, inclusive Anastasiou:2015vya ;
Mistlberger:2018etf ; Anastasiou:2016cez  as well as fully differential
Chen:2021isd  predictions for Higgs boson production via gluon fusion are
known up to next-to-next-to-next-to-leading-order (N3LO).

One of the most promising observables to study the EWSB mechanism is the Higgs
transverse momentum, see for example Bishara:2016jga . To remain differential
in external radiation, one must study the production of a Higgs boson in
association with (at least) one resolved jet. For this process, results in the
heavy top quark limit are currently known up to NNLO Boughezal:2015dra ;
Boughezal:2015aha ; Caola:2015wna ; Chen:2016zka ; Campbell:2019gmd , and the
residual theoretical uncertainty can be estimated at around
$\mathcal{O}(5\%)$. The heavy top quark limit approximation is valid for
transverse momentum of the Higgs which is lower than two times the top quark
mass, $p_{T}<2\,m_{t}$. For higher $p_{T}$, the heavy quark loop is resolved
and finite mass corrections are needed. Results with finite top mass at NLO
were first obtained for the very high transverse momentum kinematic region,
$p_{T}\gg 2m_{t}$, Lindert:2018iug , and numerically for general kinematics
Jones:2018hbb ; Chen:2021azt . More recently, following the calculation of the
relevant two-loop master integrals Bonciani:2016qxi ; Bonciani:2019jyb , the
full NLO analytical calculation has been completed Bonciani:2022jmb .

Taking into account the wealth of currently available data, as well as the
expected high-luminosity phase, N3LO calculations will become essential
Caola:2022ayt  in order to perform phenomenological studies at the $1\%$ level
at the LHC. A key ingredient to extend the current calculations to N3LO are
the four-point amplitudes for the production of a Higgs boson and a parton in
parton-parton collisions. Working in the effective theory described above, one
needs to compute tree level, one-loop, two-loop and three-loop amplitudes for
the scattering of three massless partons and one massive scalar. On top of
their phenomenological importance, the structure of these amplitudes to higher
loops is also of formal interest and has been the subject of thorough
investigation. In this context it is worth noticing that contrary to the naive
expectation, up to two loops, the finite remainder amplitudes for the decay of
a Higgs boson to three gluons have been shown to be expressible in terms of
just classical polylogarithms Duhr:2012fh .

Starting at one loop, amplitudes exhibit singularities which can be regulated
in the framework of dimensional regularization. In $D=4-2\epsilon$ dimensions,
the loop amplitudes are computed as a Laurent expansion in $\epsilon$. An N3LO
computation requires one-loop amplitudes up to $\mathcal{O}(\epsilon^{4})$,
two-loop amplitudes up to $\mathcal{O}(\epsilon^{2})$ and three-loop
amplitudes up to $\mathcal{O}(\epsilon^{0})$. Our goal in this paper is to
provide the first ingredient for such a calculation, namely the two-loop
amplitudes up to order $\epsilon^{2}$. These amplitudes were previously
obtained in Gehrmann:2012  to order $\epsilon^{0}$, providing a key ingredient
to the calculation of NNLO QCD corrections to Higgs+jet production
Boughezal:2015dra ; Boughezal:2015aha ; Caola:2015wna ; Chen:2016zka ;
Campbell:2019gmd  and the Higgs boson transverse momentum distribution
Chen:2016zka ; Chen:2018pzu .

While our overall approach is relatively standard, we include multiple new
elements which help us organize the calculation more efficiently in view of a
subsequent extension to three loops. First of all, we cast the relevant
amplitudes into a compact tensorial basis and construct new helicity
projectors to extract the corresponding helicity amplitudes directly from the
Feynman diagrams Peraro:2019cjj ; Peraro:2020sfm . We then employ standard
integration-by-parts identities Chetyrkin:1981qh ; Laporta:2000dsw  to express
these amplitudes in terms of so-called master integrals, and evaluate them
using the differential equation method Kotikov:1990kg ; Kotikov:1991hm ;
Kotikov:1991pm ; Gehrmann:1999as . At two loops, the master integrals were
computed up to transcendental weight four more than two decades ago
Gehrmann:2000zt ; Gehrmann:2001ck . Since then, substantial advances were made
in the understanding of mathematical structures underlying Feynman integrals
and their associated differential equations.

Indeed, about a decade ago it was realized that a special class of Feynman
integrals, dubbed _local integrals_ Arkani-Hamed:2010pyv ; Kotikov:2010gf ,
plays a crucial role in representing scattering amplitudes, in particular in
the case of (planar) $N=4$ Super Yang-Mills theory (SYM). These integrals only
feature singularities of the logarithmic type and exhibit uniform maximum
transcendentality Kotikov:2006ts ; Kotikov:2010gf , which has for a long time
been conjectured to characterize scattering amplitudes in $N=4$ SYM
Kotikov:2001sc . While these properties do not translate in an obvious way to
non-supersymmetric theories like QCD, it has been shown that such integrals
can still substantially simplify the calculations of scattering amplitudes
within the Standard Model. Namely, integrals of this type fulfil particularly
simple systems of differential equations, in the so-called canonical form
Henn:2013pwa . Canonical sets of equations are more easily solved and provide
a direct handle on the analytic properties of the corresponding integrals in
various singular regions. Importantly, the solution of a canonical system with
rational coefficients is straightforwardly expressed in terms of a well-
understood class of functions, the Multiple Polylogarithms (MPLs)
Goncharov:1998kja ; Remiddi:1999ew ; Vollinga:2004sn .

Among the initial applications of this formalism were planar ladder-type
integrals in $H+$jet production up to three loops DiVita:2014pza , which
include the planar two-loop integrals up to $\mathcal{O}(\epsilon^{2})$ as a
subset. Here we consider the full set of two-loop integrals: planar and non-
planar. We construct a pure basis of uniform transcendental weight and
demonstrate that it can be solved in terms of MPLs to any order in $\epsilon$.
At variance with Gehrmann:2000zt ; Gehrmann:2001ck , we show that using a
generalized set of regularity conditions on the canonical master integrals,
all boundary conditions required to fix the solution of the differential
equations can be inferred in terms of a small number of one-scale two- and
three-point functions which are known in closed form in the literature. In
view of applications at N3LO, we limit ourselves to perform the calculation
explicitly to order $\epsilon^{2}$, which corresponds to transcendental weight
six.

The structure of the paper is as follows. The effective coupling of the Higgs
boson to light partons and the definition of kinematics is given in Section 2.
In Section 3, we formulate the general structure of the amplitude and use
projectors to obtain tensor coefficients and construct the helicity
amplitudes. In Section 4, we describe in detail the construction of pure bases
for the two-loop integral families and solve their canonical differential
equations analytically. The ultraviolet renormalization and the subtraction of
infrared singularities of our amplitudes are discussed in Section 5. Finally,
the crossing of the helicity amplitudes to all appropriate kinematic
configurations and the necessary analytic continuation of the relevant
multiple polylogarithms is discussed in Section 6. We conclude with a brief
summary in Section 7.

##  2 Notation and kinematics

###  2.1 The effective Lagrangian

The Higgs boson interacts with Standard Model particles with a coupling
strength proportional to their mass, and therefore cannot couple directly to
gluons or massless quarks. Nevertheless, starting at one loop, the Higgs can
interact with gluons through virtual loops of massive quarks. In the limit of
a very heavy quark mass, $m_{q}\rightarrow\infty$, one can show that this
coupling becomes independent of $m_{q}$ and an effective theory can be
formulated by integrating out the corresponding quark from the full theory
Wilczek:1977 ; Shifman:1978 ; Inami:1982xt . While most quarks have relatively
small masses compared to the typical energy scales of scattering processes at
the LHC and can often be considered as massless, the same is not true for the
top quark. Since it is the heaviest of all known Standard Model particles,
this effective field theory (EFT) works extremely well for the top quark, at
least as long as all scales involved are smaller than twice its mass
Ellis:1988 ; Bauri:1990 .

In the following, we will work in the EFT with a matter content of $N_{f}=5$
massless quarks and one very heavy quark, the top quark, integrated out. In
this case, the effective Lagrangian becomes

$\displaystyle\mathcal{L}_{int}=-\,\frac{\lambda}{4}HG_{a}^{\mu\nu}G_{a,\mu\nu},$
(1)

where $G_{a}^{\mu\nu}$ is the field strength tensor of the gluons and $H$ is
the Higgs field. From dimensional analysis, the effective coupling $\lambda$
has inverse mass dimension. It was shown long ago Spira:95 ; Chetyrkin:98  how
to perform the matching of this effective theory to the Standard Model
Lagrangian SpiraHiggs:1993 ; Spira:1995rr ; SpiraHiggs3:1998 .

###  2.2 Kinematics

We are ultimately interested in computing the amplitude for the production of
a Higgs boson and a hadronic jet in parton-parton annihilation at the LHC. For
simplicity, we start considering the problem in the crossed kinematics which
corresponds to the decay of a Higgs boson into three partons. There are two
relevant partonic channels, namely the decay into three gluons

$H(p_{4})\rightarrow g_{1}(p_{1})+g_{2}(p_{2})+g_{3}(p_{3}),$ (2)

and into a quark, anti-quark and a gluon

$H(p_{4})\rightarrow q(p_{1})+\bar{q}(p_{2})+g(p_{3}).$ (3)

The amplitudes for the production processes can then be obtained through an
analytic continuation of the decay kinematics Gehrmann:2002zr , see Section 6
for more details.

The Mandelstam invariants are defined as

$s_{12}=(p_{1}+p_{2})^{2}\,,\quad\quad s_{13}=(p_{1}+p_{3})^{2}\,,\quad\quad
s_{23}=(p_{2}+p_{3})^{2},$ (4)

and satisfy the conservation equation

$s_{12}+s_{13}+s_{23}=M_{H}^{2},$ (5)

where $M_{H}$ is the mass of the Higgs particle. It is more convenient to work
with dimension-less ratios

$x=\frac{s_{12}}{M_{H}^{2}},\quad\quad y=\frac{s_{13}}{M_{H}^{2}},\quad\quad
z=\frac{s_{23}}{M_{H}^{2}},$ (6)

such that (5) implies the relation,

$x+y+z=1.$ (7)

In the decay kinematic region, all these invariants are non-negative. This,
together with (6), defines the corresponding kinematic region

$z\geq 0,\quad\quad 0\leq y\leq 1-z,\quad\quad x=1-y-z.$ (8)

##  3 Tensor Decomposition

Following earlier work on this process Gehrmann:2012 , we decompose the
amplitudes $\mathcal{M}_{ggg}$ and $\mathcal{M}_{q\overline{q}g}$ as

$\displaystyle\mathcal{M}_{ggg}$
$\displaystyle=\mathcal{S}_{\mu\nu\rho}(g_{1},g_{2},g_{3})\,\epsilon_{1}^{\mu}\epsilon_{2}^{\nu}\epsilon_{3}^{\rho},$
(9) $\displaystyle\mathcal{M}_{q\overline{q}g}$
$\displaystyle=\mathcal{T}_{\mu}(q,\bar{q},g)\,\epsilon^{\mu},$

where we used $\epsilon_{i}$ to denote the polarization vectors of external
gluons. The above tensors can be expanded perturbatively in the QCD coupling
constant $\alpha_{s}$ as

$\displaystyle\mathcal{S}_{\mu\nu\rho}(g_{1},g_{2},g_{3})=\,$
$\displaystyle\lambda\sqrt{4\pi\alpha_{s}}f^{a_{1}a_{2}a_{3}}\Big{[}\mathcal{S}^{(0)}_{\mu\nu\rho}(g_{1},g_{2},g_{3})+\Big{(}\frac{\alpha_{s}}{2\pi}\Big{)}\,\mathcal{S}^{(1)}_{\mu\nu\rho}(g_{1},g_{2},g_{3})$
$\displaystyle\phantom{\lambda\sqrt{4\pi\alpha_{s}}f^{a_{1}a_{2}a_{3}}\;}+\Big{(}\frac{\alpha_{s}}{2\pi}\Big{)}^{2}\,\mathcal{S}^{(2)}_{\mu\nu\rho}(g_{1},g_{2},g_{3})+\mathcal{O}(\alpha_{s}^{3})\Big{]}\,,$
(10) $\displaystyle\mathcal{T}_{\mu}(q,\overline{q},g)=\,$
$\displaystyle\lambda\sqrt{4\pi\alpha_{s}}T^{a}_{ij}\Big{[}\mathcal{T}^{(0)}_{\mu}(q,\overline{q},g)+\Big{(}\frac{\alpha_{s}}{2\pi}\Big{)}\,\mathcal{T}^{(1)}_{\mu}(q,\overline{q},g)$
$\displaystyle\phantom{\lambda\sqrt{4\pi\alpha_{s}}T^{a}_{ij}\;}+\Big{(}\frac{\alpha_{s}}{2\pi}\Big{)}^{2}\,\mathcal{T}^{(2)}_{\mu}(q,\overline{q},g)+\mathcal{O}(\alpha_{s}^{3})\Big{]}\,,$
(11)

where the coefficients $\mathcal{S}^{(i)}_{\mu\nu\rho}$ and
$\mathcal{T}_{\mu}^{(i)}$ are the $i$-loop contributions to the amplitude. The
$SU(3)$ group generators are normalised as
$\operatorname{Tr}(T^{a}T^{b})=\delta^{ab}/2$.

Given the external states and their possible helicity and spin configurations,
the amplitudes $\mathcal{S}_{\mu\nu\rho}$ and $\mathcal{T}_{\mu}$ can only
depend on a limited number of tensor structures. These structures can be
further constrained by exploiting symmetries and choosing a gauge or,
following Gehrmann:2012 , by enforcing gauge invariance through the Ward
identities. In this Section, we aim to find such a tensor basis in order to be
able to work with the scalar coefficients of the amplitudes with respect to
this basis, known as the form factors. The form factors, in turn, are obtained
by applying projector operators on the full amplitude expanded in Feynman
diagrams. We derive a basis of 4 tensor structures for $H\rightarrow ggg$ and
a basis of 2 tensor structures for $H\rightarrow q\bar{q}g$, and the
corresponding projectors in Subsections 3.1 and 3.2, respectively.

Since we are ultimately interested in fixing the helicities of external states
and computing the associated helicity amplitudes, we introduce the spinor-
helicity formalism in Subsection 3.3 and derive a set of genuinely independent
helicity amplitudes. These are written as a unique spinor factor times a
scalar coefficient which is a linear combination of the form factors. The same
linear combination of form factor projectors also defines a helicity amplitude
projector.

###  3.1 Tensor decomposition for $H\rightarrow ggg$

Let us start considering the decay of a scalar boson into three massless spin-
one particles. The most general tensor structure one can build using the four-
vectors associated with the external particles is

$\displaystyle\mathcal{S}_{\mu\nu\rho}(g_{1},g_{2},g_{3})\epsilon_{1}^{\mu}\epsilon_{2}^{\nu}\epsilon_{3}^{\rho}$
$\displaystyle=\,\sum_{i,j,k=1}^{3}A_{ijk}\,p_{i}\cdot\epsilon_{1}\,p_{j}\cdot\epsilon_{2}\,p_{k}\cdot\epsilon_{3}+\sum_{i=1}^{3}B_{i}\,p_{i}\cdot\epsilon_{1}\epsilon_{2}\cdot\epsilon_{3}$
$\displaystyle\phantom{=\;}+\sum_{i=1}^{3}C_{i}\,p_{i}\cdot\epsilon_{2}\epsilon_{3}\cdot\epsilon_{1}+\sum_{i=1}^{3}D_{i}\,p_{i}\cdot\epsilon_{3}\epsilon_{1}\cdot\epsilon_{2}\,.$

There are four helicity configurations for this amplitude, so we expect four
independent tensor structures in $D=4$. Indeed, the transversality conditions
$p_{i}\cdot\epsilon_{i}=0$, $i=1,2,3$ and the cyclic gauge choice
$\epsilon_{1}\cdot p_{2}=0$, $\epsilon_{2}\cdot p_{3}=0$, $\epsilon_{3}\cdot
p_{1}=0$ restrict the tensor structures considerably:

$\displaystyle\mathcal{S}_{\mu\nu\rho}(g_{1},g_{2},g_{3})\epsilon_{1}^{\mu}\epsilon_{2}^{\nu}\epsilon_{3}^{\rho}$
$\displaystyle=A_{312}\,p_{3}\cdot\epsilon_{1}\,p_{1}\cdot\epsilon_{2}\,p_{2}\cdot\epsilon_{3}+B_{3}\,\epsilon_{2}\cdot\epsilon_{3}\,p_{3}\cdot\epsilon_{1}$
$\displaystyle\phantom{=\,\,}+C_{1}\,p_{1}\cdot\epsilon_{2}\,\epsilon_{3}\cdot\epsilon_{1}+D_{2}\,p_{2}\cdot\epsilon_{3}\,\epsilon_{1}\cdot\epsilon_{2}\,,$
(12)
$\displaystyle=\mathcal{G}_{1}T_{1}+\mathcal{G}_{2}T_{2}+\mathcal{G}_{3}T_{3}+\mathcal{G}_{4}T_{4}\,,$
(13)

where we relabelled the coefficients $\\{A_{i},...,D_{i}\\}$ as the form
factors $\mathcal{G}_{i}$ and defined the basis

$\displaystyle T_{1}$
$\displaystyle=p_{1}\cdot\epsilon_{2}\,\,\epsilon_{3}\cdot\epsilon_{1},$
$\displaystyle T_{2}$
$\displaystyle=p_{2}\cdot\epsilon_{3}\,\,\epsilon_{1}\cdot\epsilon_{2},$
$\displaystyle T_{3}$
$\displaystyle=\epsilon_{2}\cdot\epsilon_{3}\,\,p_{3}\cdot\epsilon_{1},$
$\displaystyle T_{4}$
$\displaystyle=p_{3}\cdot\epsilon_{1}\,p_{1}\cdot\epsilon_{2}\,p_{2}\cdot\epsilon_{3}\,.$
(14)

The form factors can be obtained from a Feynman diagram decomposition of the
amplitude $\mathcal{S}^{(i)}_{\mu\nu\rho}$ at any loop order by applying
suitable projectors $\mathcal{P}_{i}$ defined as

$\displaystyle\sum_{pol}\mathcal{P}_{i}\,\mathcal{S}_{\mu\nu\rho}(g_{1},g_{2},g_{3})\epsilon_{1}^{\mu}\epsilon_{2}^{\nu}\epsilon_{3}^{\rho}=\mathcal{G}_{i}\,,$
(15)

where polarization vectors satisfy the gauge-fixed polarization sum, with
reference vectors defined above. The projectors can in turn be decomposed in
terms of the dual of the tensor basis in (13):

$\displaystyle\mathcal{P}_{i}=\sum_{j=1}^{4}c_{i}^{(j)}T_{j}^{\dagger}\,.$
(16)

To work out the projectors explicitly, we insert the decompositions (16) and
(13) into the definition (15) and obtain the requirement

$\displaystyle\sum_{j=1}^{4}c_{i}^{(j)}T_{j}^{\dagger}\sum_{k=1}^{4}\mathcal{G}_{k}T_{k}\overset{!}{=}\mathcal{G}_{i}\,,$
(17)

which is satisfied by the coefficients
$c_{i}^{(j)}=(T^{\dagger}_{j}T_{i})^{-1}$. In particular for the amplitude
$H\rightarrow ggg$ with external states in $D$ dimensions Peraro:2020sfm , we
get

$\displaystyle\mathcal{P}_{1}$
$\displaystyle=\frac{1}{s_{12}s_{13}(D-3)}(s_{23}\,T_{1}^{\dagger}-T_{4}^{\dagger}),$
$\displaystyle\mathcal{P}_{2}$
$\displaystyle=\frac{1}{s_{12}s_{23}(D-3)}(s_{13}\,T_{2}^{\dagger}-T_{4}^{\dagger}),$
$\displaystyle\mathcal{P}_{3}$
$\displaystyle=\frac{1}{s_{13}s_{23}(D-3)}(s_{12}\,T_{3}^{\dagger}-T_{4}^{\dagger}),$
$\displaystyle\mathcal{P}_{4}$
$\displaystyle=\frac{1}{s_{12}s_{13}s_{23}(D-3)}(D\,T_{4}^{\dagger}-s_{12}T_{3}^{\dagger}-s_{23}T_{1}^{\dagger}-s_{13}T_{2}^{\dagger})\,.$
(18)

###  3.2 Tensor decomposition for $H\rightarrow q\bar{q}g$

Similarly, it is easy to see that the most general tensor decomposition for
two external spinors and a four-vector is

$\displaystyle\mathcal{T}_{\mu}(q,\overline{q},g)\epsilon_{3}^{\mu}$
$\displaystyle=\sum_{\Gamma}A_{\Gamma}\,\,\bar{u}(p_{1})\,\Gamma_{\mu}\,v(p_{2})\epsilon_{3}^{\mu}$
$\displaystyle+\sum_{j=1,2,3}\,\sum_{\Gamma^{\prime}}B_{\Gamma^{\prime}}\,\bar{u}(p_{1})\,\Gamma^{\prime}\,v(p_{2})\epsilon_{3}\cdot
p_{j}\,,$ (19)

where we sum over odd products of gamma matrices. In the first sum, all
indices but one are contracted amongst each other or with external momenta
while in the second, there are no indices left, and the polarization vector
contracts with an external momentum vector.
In general, one might expect the length of the spinor chains to be bound by
the loop order. By simple enumeration, one can see that at tree level there
can only be one gamma matrix, at one loop up to three gamma matrices and at
two loops up to five. Nonetheless, it is easy to see that with the momenta at
hand, one cannot build any spinor chain with more than one Dirac $\gamma$
matrix. Enforcing the transversality condition $p_{3}\cdot\epsilon_{3}=0$ and
gauge choice $\epsilon_{3}\cdot p_{1}=0$, one is left with the decomposition

$\displaystyle\mathcal{T}_{\mu}(q,\overline{q},g)\epsilon_{3}^{\mu}$
$\displaystyle=\,\mathcal{F}_{1}\,\bar{u}(p_{1})\not{\epsilon_{3}}v(p_{2})+\mathcal{F}_{2}\,\bar{u}(p_{1})\not{p_{3}}v(p_{2})\epsilon_{3}\cdot
p_{2}.$ (20)

Hence we define the two tensor structures

$\displaystyle T_{1}$
$\displaystyle=\bar{u}(p_{1})\not{\epsilon_{3}}v(p_{2}),$ $\displaystyle
T_{2}$ $\displaystyle=\bar{u}(p_{1})\not{p_{3}}v(p_{2})\epsilon_{3}\cdot
p_{2}.$ (21)

The form factors are extracted from the amplitude using the projectors
$\mathcal{P}_{i}$ with

$\displaystyle\sum_{pol}\mathcal{P}_{i}\,\mathcal{T}_{\mu}(q,\bar{q},g)\epsilon_{3}^{\mu}=\mathcal{F}_{i},$
(22)

where polarization vector of the gluon satisfies the gauge-fixed polarization
sum, with reference vector defined above. Following the strategy outlined in
Subsection 3.1, we get for the two projectors

$\displaystyle\mathcal{P}_{1}$
$\displaystyle=\frac{1}{2s_{12}(D-3)}\Big{(}T_{1}^{\dagger}-\frac{1}{s_{23}}T_{2}^{\dagger}\Big{)},$
$\displaystyle\mathcal{P}_{2}$
$\displaystyle=\frac{1}{2s_{12}s_{23}(D-3)}\Big{(}\frac{D-2}{s_{23}}T_{2}^{\dagger}-T_{1}^{\dagger}\Big{)}.$
(23)

###  3.3 Helicity amplitudes

From the tensors found above, we can easily obtain compact expressions for the
relevant helicity amplitudes. We work in the ’t Hooft-Veltman scheme and
consider external states as four-dimensional and assume fixed helicity states.
In massless QCD, both gluons and quarks have two helicity configurations,
$\lambda_{i}=\pm$, and the amplitudes can be written as

$\displaystyle\mathcal{M}^{\lambda_{1}\lambda_{2}\lambda_{3}}_{ggg}$
$\displaystyle=\mathcal{S}_{\mu\nu\rho}(g_{1},g_{2},g_{3})\,\epsilon_{1,\lambda_{1}}^{\mu}(p_{1})\epsilon_{2,\lambda_{1}}^{\nu}(p_{2})\epsilon_{3,\lambda_{3}}^{\rho}(p_{3}),$
(24)
$\displaystyle\mathcal{M}_{q\overline{q}g}^{\lambda_{1}\lambda_{2}\lambda_{3}}$
$\displaystyle=\mathcal{T}_{\mu}(q_{\lambda_{1}},\bar{q}_{\lambda_{2}},g)\,\epsilon^{\mu}_{3,\lambda_{3}}(p_{3})\,.$

The two helicity states of a four-component massless spinor are projected out
through

$\psi_{L}=\frac{1}{2}\big{(}1-\gamma_{5}\big{)}\,\psi,\quad\psi_{R}=\frac{1}{2}\big{(}1+\gamma_{5}\big{)}\,\psi,$
(25)

and we fix the spinor-helicity bracket representation of incoming fermions as

$\displaystyle|p]=\psi_{L}(p)\,,\quad|p\rangle=\psi_{R}(p)\,,$ (26)

and of incoming anti-fermions as

$\displaystyle\langle p|=\bar{\psi}_{L}(p),\quad[p|=\bar{\psi}_{R}(p).$ (27)

For massless vector bosons, incoming states of positive and negative helicity
are given by

$\epsilon_{+}^{\mu}(p;r)=-\frac{[r\,\gamma^{\mu}\,p\rangle}{\sqrt{2}\,[rp]}\,,\qquad\epsilon_{-}^{\mu}(p;r)=\frac{\langle
r\,\gamma^{\mu}\,p]}{\sqrt{2}\,\langle rp\rangle},$ (28)

where the reference momentum $r$ is an arbitrary light-like vector such that
$r\cdot p\neq 0$. Outgoing states have the same representation with switched
helicities.

Let us start by considering the decay $H\rightarrow ggg$. There are two
independent helicity configurations, which we choose to be
$\mathcal{M}^{+++}_{ggg}$ and $\mathcal{M}^{++-}_{ggg}$, while the other
helicity amplitudes are obtained through parity conjugation and by relabelling
of the gluon momenta. Starting from the decomposition (13) with the basis (14)
and applying the definitions (26)–(28), we can cast the independent helicity
amplitudes in terms of spinor products,

$\displaystyle\mathcal{M}^{+++}_{ggg}$
$\displaystyle=\alpha\frac{1}{\sqrt{2}}\frac{M_{H}^{4}}{\langle
1\,2\rangle\langle 2\,3\rangle\langle 3\,1\rangle},$ (29)
$\displaystyle\mathcal{M}^{++-}_{ggg}$
$\displaystyle=\beta\frac{1}{\sqrt{2}}\frac{[1\,2]^{3}}{[2\,3][1\,3]}\,,$ (30)

where the coefficients $\alpha$ and $\beta$ are simple linear combinations of
the original form factors:

$\displaystyle\alpha$
$\displaystyle=-\frac{2s_{12}s_{13}\,\mathcal{G}_{1}+2s_{12}s_{23}\,\mathcal{G}_{2}+2s_{13}s_{23}\,\mathcal{G}_{3}+s_{12}s_{13}s_{23}\,\mathcal{G}_{4}}{2M_{H}^{4}},$
(31) $\displaystyle\beta$
$\displaystyle=-\frac{2s_{23}\,\mathcal{G}_{2}+s_{23}s_{13}\mathcal{G}_{4}}{2s_{12}}.$
(32)

The helicity projectors for $\alpha$ and $\beta$ are built by replacing the
$\mathcal{G}_{i}$ with the corresponding form factor projectors from (18):

$\displaystyle\mathcal{P}_{\alpha}$
$\displaystyle=-\frac{1}{2M_{H}^{4}}\Big{(}2s_{12}s_{13}\mathcal{P}_{1}+2s_{12}s_{23}\,\mathcal{P}_{2}+2s_{13}s_{23}\,\mathcal{P}_{3}+s_{12}s_{13}s_{23}\,\mathcal{P}_{4}\Big{)}$
$\displaystyle=\frac{1}{2(D-3)M_{H}^{4}}\Big{(}(6-D)T_{4}^{\dagger}-s_{23}T_{1}^{\dagger}-s_{12}T_{3}^{\dagger}-s_{13}T_{2}^{\dagger}\Big{)},$
(33) $\displaystyle\mathcal{P}_{\beta}$
$\displaystyle=-\frac{2s_{23}\,\mathcal{P}_{2}+s_{23}s_{13}\mathcal{P}_{4}}{2s_{12}}$
$\displaystyle=\frac{1}{2(D-3)s_{12}^{2}}\Big{(}s_{23}T_{1}^{\dagger}-s_{13}T_{2}^{\dagger}+s_{12}T_{3}^{\dagger}+(2-D)T_{4}^{\dagger}\Big{)}.$
(34)

Let us now consider the other partonic process $H\rightarrow q\bar{q}g$. In
this case, there is only one independent helicity amplitude, which we choose
to be $\mathcal{M}_{q\bar{q}g}^{LR+}$. In this helicity configuration, the
first tensor in (20) vanishes, while the second tensor yields

$\mathcal{M}_{q\bar{q}g}^{LR+}=\gamma\frac{1}{\sqrt{2}}\frac{[2\,3]^{2}}{[1\,2]},$
(35)

where

$\gamma=s_{12}\,\mathcal{F}_{2}.$ (36)

Again, the helicity projector is easily obtained by replacing
$\mathcal{F}_{2}$ with the associated projector from (23),

$\mathcal{P}_{\gamma}=s_{12}\,\mathcal{P}_{2}=\frac{1}{2s_{23}(D-3)}\Big{(}\frac{D-2}{s_{23}}T_{2}^{\dagger}-T_{1}^{\dagger}\Big{)}.$
(37)

Just like full amplitudes in (10) and (11), the helicity amplitude
coefficients $\alpha$, $\beta$, $\gamma$ also have a perturbative expansion,

$\Omega=\lambda\sqrt{4\pi\alpha_{s}}\,T_{\Omega}\,\Big{[}\Omega^{(0)}+\Big{(}\frac{\alpha_{s}}{2\pi}\Big{)}\,\Omega^{(1)}+\Big{(}\frac{\alpha_{s}}{2\pi}\Big{)}^{2}\,\Omega^{(2)}+\mathcal{O}(\alpha_{s}^{3})\Big{]}$
(38)

for $\Omega=\alpha,\beta,\gamma$. The overall colour factors for the two
processes are $T_{\alpha}=T_{\beta}=f^{a_{1}a_{2}a_{3}}$ and
$T_{\gamma}=T^{a_{3}}_{\,\,\,i,j}$. The tree-level helicity amplitudes are
well-known and the coefficients evaluate to

$\displaystyle\alpha^{(0)}$ $\displaystyle=\beta^{(0)}=-1$ (39)
$\displaystyle\gamma^{(0)}$ $\displaystyle=1.$ (40)

##  4 Master integrals

Family PL: $\\{P_{i}\\}$ | Family NPL: $\\{N_{i}\\}$
---|---
$k_{1}$ | $k_{1}$
$k_{2}$ | $k_{2}$
$k_{1}-k_{2}$ | $k_{1}-p_{1}$
$k_{1}-p_{1}$ | $k_{2}-p_{1}$
$k_{2}-p_{1}$ | $k_{1}+p_{3}$
$k_{1}-p_{1}-p_{2}$ | $k_{2}+p_{3}$
$k_{2}-p_{1}-p_{2}$ | $k_{1}+p_{2}+p_{3}$
$k_{1}-p_{1}-p_{2}-p_{3}$ | $k_{1}-k_{2}$
$k_{2}-p_{1}-p_{2}-p_{3}$ | $k_{1}-k_{2}+p_{2}$
Table 1: The two auxiliary topologies of momenta which label the propagators
of every diagram appearing in the amplitudes. PL labels exclusively planar
diagrams, whereas NPL accommodates also all non-planar sectors.

The form factors $\mathcal{F}_{i}$ and $\mathcal{G}_{i}$, which relate to the
helicity amplitudes via (32), (36), receive contributions from all relevant
Feynman diagrams at a given perturbative order. All tree-level, one- and two-
loop diagrams with $H$ and $ggg$ or $q\bar{q}g$ external states are generated
with standard QCD vertices and the effective Higgs interactions using QGRAF
Nogueira:1991ex . We shift each diagram to a kinematic crossing of one of the
auxiliary topologies presented in Table 1 using Reduze2 Studerus:2009ye ;
vonManteuffel:2012np . After inserting Feynman rules and evaluating the Dirac
and Lorentz algebra in FORM Vermaseren:2000nd , the contribution of each
diagram to the form factors in (13) and (20) can be written as a combination
of scalar integrals of the form

$I_{a_{1},...,a_{9}}=\int\bigg{(}\prod_{l=1}^{2}(-M_{H}^{2})^{-\epsilon}e^{\gamma_{E}\epsilon}\frac{\mathrm{d}^{D}k_{l}}{i\pi^{d/2}}\bigg{)}\prod_{i=1}^{9}D_{i}^{-a_{i}}\,,$
(41)

with $k_{l}$ the loop momenta, $D_{i}\in\\{P_{i}\\}\cup\\{N_{i}\\}$ the
internal propagators, and $\gamma_{E}=0.577\dots$ the Euler-Mascheroni
constant. The topology of each scalar integral is uniquely identified by the
propagators which appear in the denominator ($a_{i}\geq 1$). This information
can be used to compute the diagram’s sector ID within one of the auxiliary
topologies in Table 1 as a binary nine-bit number. The scalar integrals
defined in (41) satisfy integration-by-parts (IBP) identities Chetyrkin:1981qh
; Laporta:2000dsw , allowing them to be expressed in terms of a minimal set of
so-called master integrals. Note that the physical measure in the bare
amplitude is $\mathrm{d}^{D}k/(2\pi)^{D}$ per loop compared to the integration
measure in (41), which then requires a simple conversion factor when inserting
the solutions of the master integrals into the amplitude.

To evaluate the master integrals, we use the method of differential equations
Kotikov:1990kg ; Kotikov:1991hm ; Kotikov:1991pm ; Gehrmann:1999as , augmented
with a canonical basis Henn:2013pwa  for the master integrals. In particular,
we reduce all scalar integrals appearing in the amplitudes using Reduze2
vonManteuffel:2012np  and Kira Maierhofer:2017gsa ; Klappert:2020nbg  directly
to a canonical basis which involves 89 master integrals. Our basis is
different from the one considered in Gehrmann:2000zt ; Gehrmann:2001ck , but
relations amongst the two sets of integrals can easily be established. There
are 4 planar and 2 non-planar sectors with the maximum number of propagators,
$t=7$, depicted in Figure 1. 85 of the 89 IBP master integrals are subsectors
of the non-trivial top sectors (d)-(f). They are kinematic crossings of the 16
planar and 8 non-planar topologies depicted in Figures 2 and 3, which will be
explained in more detail after introducing the canonical basis.

(a) PL 253

(b) PL 493

(c) PL 471

(d) PL 247

(e) NPL 487

(f) NPLx23 494

Figure 1: The 4 planar and 2 non-planar sectors with $t=7$, labelled with the
propagators $P_{i}$ and $N_{i}$, respectively, as listed in Table 1. Their
names indicate the integral family and sector ID. The sectors (d)-(f) contain
new master integrals and are the focus of this Section. Sectors (a) and (b)
are entirely reducible to master integrals found in the trees of the former
three top sectors. Sector (c) can be reduced to masters from (d)-(f) and four
additional integrals which are easily computed from one-loop results.

###  4.1 Canonical basis

Multiple public packages exist for the derivation of a candidate canonical
basis like CANONICA Meyer:2017joq , Fuchsia Gituliar:2017vzm , or DlogBasis
Henn:2020lye . We use DLogBasis to provide a list of candidate UV-finite
integrals with unit leading singularities for the planar family. For the non-
planar family, we constructed canonical candidates by studying the leading
singularities of the relevant master integrals, following a loop-by-loop
approach. Note that for all planar and some non-planar integrals, canonical
candidates can be obtained just by rescaling a single integral in the sector
by its maximal cut, i.e. a solution to the homogeneous part of its
differential equation Primo:2016ebd . This can be confirmed by studying the
maximal cut of the specific canonical integral in the Baikov representation
Baikov:1996iu . The canonical bases sufficient for the computation of all
integrals in the families PL and NPL are presented in the supplementary
material. The amplitude, however, contains diagrams which produce integrals in
kinematic crossings of the families in Table 1. To obtain a basis that is
sufficient to represent the physical amplitude, we proceed from the lowest
sectors in applying crossings to the canonical integrals in PL and NPL, and
appending to our basis those canonical integrals whose reduction contains new
masters. This yields a minimal set of crossed and uncrossed canonical master
integrals, sufficient for physical applications, whose generic topologies are
depicted in Figures 2 and 3. Note that we obtain results first in the
Euclidean region, where all invariants are negative $s_{ij}<0$. For
simplicity, we also set $M_{H}^{2}=-1$ in explicit formulas below.

Up to now, we only considered the non-trivial top sectors (d)-(f) in Figure 1,
but the amplitude requires the reduction of integrals from all six topologies.
While integrals in (a) and (b) can be reduced to integrals considered earlier,
(c) contains four new master integrals, $I_{86},\dots,I_{89}$ which need to be
added to the basis. They are the integrals (the numbering refers to the
canonical basis in the supplementary material):

$\displaystyle I_{86}$ $\displaystyle=\epsilon^{2}\left[\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}} \feynman[inline=(a1.base)]
\vertex(x1) at (-1,0.25) {$p_{1}$}; \vertex(a1) at (-1,0) {$p_{2}$};
\vertex(y1) at (-1,-0.25) {$p_{3}$}; \vertex(b1) at (-0.5,0) ; \vertex(b2) at
(0,0); \vertex(b3) at (0.5,0) ; \vertex(a3) at (1,0) {$q$} ; \diagram* { (a1)
-- [] (b1), (x1) -- [] (b1), (y1) -- [] (b1), (b1) -- [ half left] (b2), (b2)
-- [ half left] (b1), (b2) -- [ half left] (b3), (b3) -- [ half left] (b2),
(b3) -- [scalar] (a3), }; \vertex[dot] (i) at ($(b1)!0.5!(b2)+(0,0.21)$) {};
\vertex[dot] (j) at ($(b2)!0.5!(b3)+(0,0.21)$) {};
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\right],$ (42)
$\displaystyle I_{87}$
$\displaystyle=-\epsilon^{3}(1-y-z)z\left[\leavevmode\hbox to0pt{\vbox
to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}} \feynman[inline=(a1.base)]
\vertex(b1) at (-1,1) {$p_{3}$}; \vertex(b2) at (-1,-1) {$p_{2}$}; \vertex(b3)
at (1,1) {$q$}; \vertex(b4) at (1,-1) {$p_{1}$}; \vertex(x1) at (-0.35,0.35) ;
\vertex(x2) at (-0.35,-0.35); \vertex(x3) at (0.35,0.35) ; \vertex(x4) at
(0.35,-0.35) ; \vertex(y3) at (0.55,0.55) ; \diagram* { (b1) -- [] (x1), (b2)
-- [] (x2), (y3) -- [scalar] (b3), (x3) -- [half right] (y3), (x3) -- [half
left] (y3), (b4) -- [] (x4), (x3) -- [] (x1), (x1) -- [] (x2), (x2) -- []
(x4), (x4) -- [] (x3) }; \vertex[dot] (i) at ($(x3)!0.5!(y3)+(-0.09,0.09)$)
{}; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\right]$ (43)

and

$I_{88}=I_{87}(p_{1}\leftrightarrow p_{2})\,,\qquad
I_{89}=I_{87}(p_{2}\leftrightarrow p_{3})\,.$

The full expressions up to $\mathcal{O}(\epsilon^{6})$ for $I_{86}$ as well as
$I_{87}$ and its two crossings are given in the supplementary material. In
this way, one can complete the set of 89 master integrals sufficient to reduce
any integral in this process. The full canonical basis is given in the
supplementary material accompanying this paper.

Figure 2: Planar topologies which appear in the canonical basis. The dashed
line is the massive leg, dots represent squared propagators. A propagator in
brackets denotes an integral with this propagator in the numerator.

###  4.2 Solution of differential equations

For the purpose of computing the master integrals, we consider first the 16
integrals in the tree of top sector (d) in PL and the 36 integrals in the tree
of the non-planar top sectors (e), (f) in NPL. Separately for the two
auxiliary topologies, we compute the derivatives of the candidate canonical
combinations, and insert the IBP reduction to obtain differential equations in
the following form

$\displaystyle\frac{\partial}{\partial
y}\vec{I}(y,z;\epsilon)=\epsilon\left(\frac{1}{y}A_{0}+\frac{1}{y-1}A_{1}+\frac{1}{y+z}A_{z}+\frac{1}{y-(1-z)}A_{1-z}\right)\vec{I}(y,z;\epsilon)\,,$
(44) $\displaystyle\frac{\partial}{\partial
z}\vec{I}(y,z;\epsilon)=\epsilon\left(\frac{1}{z}B_{0}+\frac{1}{z-1}B_{1}+\frac{1}{z+y}B_{y}+\frac{1}{z-(1-y)}B_{1-y}\right)\vec{I}(y,z;\epsilon)\,,$
(45)

where $A_{i},B_{i}$ are sparse matrices of rational numbers. It is obvious
from this form that the solutions for the canonical combinations can be
expressed in terms of MPLs Goncharov:1998kja  with the alphabet
$\\{y,z,y-1,z-1,y+z,1-y-z\\}$, which are usually written out in terms of a
fibration in either $y$ or $z$. We recall here that MPLs are defined as
iterated integrals over rational functions

$\displaystyle
G(l_{1},...,l_{n};x)=\int_{0}^{x}\frac{dt}{t-l_{1}}G(l_{2},...,l_{n};t)\,,\quad
G(\underbrace{0,...,0}_{n};x)=\frac{1}{n!}\log^{n}(x)\,,$ (46)

with $G(x)=1$. In this context, $n$ is referred to as the _transcendental
weight_ of the polylogarithm.

By construction, the method of differential equations cannot be used to
compute purely one-scale integrals directly. These must instead be obtained
from alternative methods and added to the system of differential equations. At
two loops, there are five such two- and three-point functions:

$\displaystyle I_{1}$ $\displaystyle=\epsilon^{2}(1-y-z)\left[\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}} \feynman[inline=(a1.base)]
\vertex(x1) at (-0.5,0.25) {$p_{1}$}; \vertex(y1) at (-0.5,-0.25) {$p_{2}$};
\vertex(b1) at (-0.25,0) ; \vertex(b3) at (0.25,0) ; \vertex(a3) at (0.5,0)
{$p_{12}$} ; \diagram* { (x1) -- [] (b1), (y1) -- [] (b1), (b1) -- [] (b3),
(b3) -- [ half left] (b1), (b1) -- [ half left] (b3), (b3) -- [scalar] (a3),
}; \vertex[dot] (i) at ($(b1)!0.5!(b3)+(0,0.21)$) {}; \vertex[dot] (j) at
($(b1)!0.5!(b3)-(0,0.21)$) {}; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\right]\,,$ (47)
$\displaystyle I_{2}$ $\displaystyle=-\epsilon^{2}\left[\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}} \feynman[inline=(a1.base)]
\vertex(x1) at (-0.5,0.25) {$p_{1}$}; \vertex(a1) at (-0.5,0) {$p_{2}$};
\vertex(y1) at (-0.5,-0.25) {$p_{3}$}; \vertex(b1) at (-0.25,0) ; \vertex(b3)
at (0.25,0) ; \vertex(a3) at (0.5,0) {$p_{4}$} ; \diagram* { (x1) -- [] (b1),
(a1) -- [] (b1), (y1) -- [] (b1), (b1) -- [] (b3), (b3) -- [ half left] (b1),
(b1) -- [ half left] (b3), (b3) -- [scalar] (a3), }; \vertex[dot] (i) at
($(b1)!0.5!(b3)+(0,0.21)$) {}; \vertex[dot] (j) at ($(b1)!0.5!(b3)-(0,0.21)$)
{}; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\right]\,,$ (48)
$\displaystyle I_{6}$ $\displaystyle=\epsilon^{3}(1-y-z)\left[\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}} \feynman[inline=(a1.base)]
\vertex(x1) at (-1,0.25) {$p_{1}$}; \vertex(y1) at (-1,-0.25) {$p_{2}$};
\vertex(b1) at (-0.5,0) ; \vertex(b2) at (0,0.25); \vertex(b3) at (0,-0.25);
\vertex(v2) at (0.5,0.25) {$p_{1}$} ; \vertex(v3) at (0.5,-0.25) {$p_{2}$} ;
\diagram* { (x1) -- [] (b1), (y1) -- [] (b1), (b1) -- [ ] (b2), (b2) --
[quarter left ] (b3), (b2) -- [ quarter right] (b3), (b1) -- [ ] (b3), (b3) --
[ ] (v3), (b2) -- [] (v2), }; \vertex[dot] (i) at ($(b2)!0.5!(b3)+(0.1,0)$)
{}; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\right]\,,$ (49)
$\displaystyle I_{9}$
$\displaystyle=\epsilon^{2}(1-y-z)^{2}\left[\leavevmode\hbox to0pt{\vbox
to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}} \feynman[inline=(a1.base)]
\vertex(x1) at (-1,0.25) {$p_{1}$}; \vertex(y1) at (-1,-0.25) {$p_{2}$};
\vertex(b1) at (-0.5,0) ; \vertex(b2) at (0,0); \vertex(b3) at (0.5,0) ;
\vertex(a3) at (1,0) {$p_{12}$} ; \diagram* { (x1) -- [] (b1), (y1) -- []
(b1), (b1) -- [ half left] (b2), (b2) -- [ half left] (b1), (b2) -- [ half
left] (b3), (b3) -- [ half left] (b2), (b3) -- [scalar] (a3), }; \vertex[dot]
(i) at ($(b1)!0.5!(b2)+(0,0.21)$) {}; \vertex[dot] (j) at
($(b2)!0.5!(b3)+(0,0.21)$) {}; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\right]\,,$ (50)
$\displaystyle I_{68}$ $\displaystyle=\epsilon^{4}z^{2}\left[\leavevmode\hbox
to0pt{\vbox to0pt{\pgfpicture\makeatletter\hbox{\hskip 0.0pt\lower 0.0pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{ }{{}} \feynman[inline=(a1.base)]
\vertex(x1) at (-1,0.25) {$p_{2}$}; \vertex(y1) at (-1,-0.25) {$p_{3}$};
\vertex(b1) at (-0.5,0) ; \vertex(b2) at (0,0.25); \vertex(b3) at (0,-0.25);
\vertex(z2) at (-0.25,0.25); \vertex(z3) at (-0.25,-0.25); \vertex(v2) at
(0.5,0.25) {$p_{2}$} ; \vertex(v3) at (0.5,-0.25) {$p_{3}$} ; \diagram* { (x1)
-- [] (b1), (y1) -- [] (b1), (b1) -- [ ] (z2), (z2) -- [ ] (b3), (z3) -- [ ]
(b2), (b1) -- [ ] (z3), (b3) -- [ ] (v3), (b2) -- [] (v2), (b3) -- [ ] (z3),
(b2) -- [] (z2), }; \pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}\right]\,,$ (51)

and their kinematic crossings. Analytical expressions for all these one-scale
integrals are known in closed form in the dimensional regulator $\epsilon$
Gehrmann:1999as ; Gehrmann:2005pd .

We approach the solution of the differential equations for the remaining
integrals as follows. At each order $n$ in $\epsilon$, we consider the vector
of master integrals $\vec{I}^{(n)}(y,z)$ and by choice integrate first the
equations in the variable $y$. If the set of differential equations in the two
variables are consistent, we obtain a partial solution which differs from the
full solution $\vec{I}^{(n)}(y,z)$ only by a function of the other variable,
$\vec{f}(z)$,

$\vec{I}^{(n)}(y,z)=\vec{I}^{(n)}_{y}(y,z)+\vec{f}(z)\,.$ (52)

The intermediate result is then substituted into the set of equations in $z$

$\displaystyle\frac{\partial}{\partial
z}\vec{I}^{(n)}(y,z)=\frac{\partial}{\partial
z}\vec{I}^{(n)}_{y}(y,z)+\frac{\partial}{\partial
z}\vec{f}(z)\overset{!}{=}B\vec{I}^{(n-1)}(y,z)\,$ (53)

and solved for $\vec{f}(z)$, which fixes the final solution up to a numerical
constant

$\displaystyle\vec{I}^{(n)}(y,z)=\vec{I}^{(n)}_{y}(y,z)+\vec{I}^{(n)}_{z}(z)+\vec{c}\,,$
(54)

with,

$\vec{I}^{(n)}_{z}(z)=\int^{z}dz^{\prime}\big{[}B\vec{I}^{(n-1)}(y,z^{\prime})-\frac{\partial}{\partial
z^{\prime}}\vec{I}^{(n)}_{y}(y,z^{\prime})\big{]}.$ (55)

As stated above, (55) cannot depend on $y$ and the spurious dependence must
cancel from the right-hand side of the equation.

Figure 3: The non-planar topologies in the canonical basis.

###  4.3 Fixing boundary conditions

The remaining numerical constants $\vec{c}$ can in principle be fixed by
evaluating the integrals at special kinematical points. We opted to obtain all
boundary conditions without any additional computations simply by imposing a
set of regularity conditions on the general solution (54). As suggested in
Gehrmann:2000zt ; Gehrmann:2001ck , a large set of boundary conditions can be
obtained exploiting regularity of the master integrals at various pseudo-
thresholds. This can be achieved in practice by multiplying the differential
equations (44) and (45) by those letters which correspond to pseudo-thresholds
of the integrals, and taking the limit as follows:

left-hand side $\displaystyle:\lim_{x\to l_{i}}\left(\frac{\partial}{\partial
x}\vec{I}\right)(x-l_{i})$ $\displaystyle=0\,,$ (56) right-hand side
$\displaystyle:\lim_{x\to l_{i}}\left(\epsilon A\vec{I}\right)(x-l_{i})$
$\displaystyle=\lim_{x\to l_{i}}\epsilon\sum_{j}A_{j}(x-l_{i})\vec{I}=\epsilon
A_{i}\lim_{x\to l_{i}}\vec{I}\,.$ (57)

Requiring the right-hand side to vanish yields non-trivial relations between
the integrals in the limit $x\to l_{i}$. Namely, if the rational factor in
question appears in the homogeneous term of the differential equation for a
given master integral, its value at a regular kinematic point can be typically
related to other integrals in the same sector and its subtopologies. This
approach is sufficient for the planar topology.

In order to impose these regularity conditions, we used PolyLogTools
Duhr:2019tlz  to manipulate multiple polylogarithms up to weight 5, evaluate
the required limits, and perform changes in the fibration basis. Beyond weight
5, we had to carry out the required fibrations ourselves, building on the
implementation of differentiation and integration of MPLs in PolyLogTools. In
particular, we differentiated the integrals with respect to the variable we
intend to fibrate into and obtained linear combinations of MPLs of weight 5,
which can be treated with automated routines. The result can subsequently be
integrated back and expressed in the required form, up to an integration
constant. All constants can be fixed by comparing the original and fibrated
function at a kinematic point, and reconstructing their difference as a number
of the appropriate weight using the PSLQ algorithm pslq . In practice, we find
that our definition of the variables $y$ and $z$, (6), allows us to consider
just the 3 limits $y\to 1$, $y\to 0$, and $y\to-z$.

In the NPL topology, integrals possess branch points when $y\to 0$, $z\to 0$
and $1-y-z\to 0$. For this reason, the strategy outlined above can only be
applied to a small number of letters and one can show that these conditions
are not enough to fix all remaining constants. Taking inspiration from
Henn:2020lye ; Wasser:2022kwg , we also consider the singular limits (i.e.
genuine thresholds of the master integrals). In terms of the Mandelstam
variables, these correspond to the additional limits $s\to 0$, $t\to 0$ and
$u\to 0$. Crucially, our canonical basis consists of UV-finite integrals.
Hence, we only need to regulate IR divergences and can assume that
$\epsilon<0$. Within this condition and keeping $\epsilon$ fixed, if a linear
combination of integrals develops a singular behaviour for one of the singular
limits, this must correspond to a spurious UV-type divergence. Since no new UV
divergences can appear when the kinematical invariants take special values, we
may impose that such spurious divergences do not occur at the kinematic points
where one of the letters vanishes. In this way, we obtain yet more equations
between the boundary constants. Explicitly, the solution to the DE near the
point $y\to l_{i}$ is

$\displaystyle\vec{I}(y,z;\epsilon)=\exp{\\{\epsilon\log{(y-l_{i})}\lim_{y\to
l_{i}}[yA_{y}]\\}}\vec{I}\big{|}_{y=l_{i}}+\mathcal{O}(y)\,.$ (58)

The matrix exponential contains elements with terms of the type
$(y-l_{i})^{a\epsilon}$. If $a>0$, such an expression diverges in the limit
$y\to l_{i}$. Since all our integrals must be finite at this kinematic point,
the constants $\vec{I}\big{|}_{y=l_{i}}$ ought to take specific values such
that these terms cancel.

We also studied the asymptotic behaviour of our master integrals with asy
Jantzen:2012mw , which is implemented in FIESTA Smirnov:2015mct , and verified
the $\epsilon$-dependence in the limit $y\to l_{i}$. Note that in
Gehrmann:2001ck , it was required to integrate one order in $\epsilon$ higher
than necessary to enforce the non-appearance of spurious singularities at the
previous order, which is no longer required. We checked our solutions for all
top sector master integrals numerically against pySecDec Borowka:2017idc  for
several Euclidean points up to weight six and found perfect agreement.

In Appendix B, we present the computation of the one-loop master integrals to
order $\mathcal{O}(\epsilon^{4})$, which also serves as a simple example of
some of the techniques discussed in this section.

##  5 UV renormalisation and IR regularisation

The bare helicity amplitudes (38) contain ultraviolet (UV) as well as infrared
(IR) divergences that manifest as poles in the Laurent expansion in the
dimensional regulator $\epsilon$. The former are treated in the
$\overline{\text{MS}}$ scheme by expressing the amplitudes in terms of the
renormalized couplings, $\alpha_{s}\equiv\alpha_{s}(\mu^{2})$ and
$\lambda\equiv\lambda(\mu^{2})$, evaluated at the renormalization scale
$\mu^{2}$. The resulting amplitudes still contain IR singularities, which will
be cancelled analytically by those occurring in radiative processes of the
same order Kinoshita:1962 ; Lee:1964 . Their structure is universal and it was
originally determined up to two loops by Catani Catani:1997 ; Catani:1998 .
These results were later systematised and extended to general processes and up
to three loops in Sterman:2002qn ; Aybat:2006wq ; Aybat:2006mz ; Becher:2009cu
; Becher:2009qa ; Dixon:2009gx ; Gardi:2009qi ; Gardi:2009zv ; Almelid:2015jia
.

In this Section, we present the necessary steps and formulae to perform the UV
renormalization and subtraction of the IR poles. This allows us to obtain the
one-loop and two-loop finite remainders, which we decompose according to their
colour structure. In contrast to Gehrmann:2012 , where the IR subtraction was
performed in the Catani scheme, we followed a subtraction scheme based on
Soft-Collinear Effective Theory Becher:2009cu ; Becher:2009qa , which can be
more naturally extended to higher loops. In subsection 5.5 we provide
conversion formulas between the two different schemes.

###  5.1 Ultraviolet renormalization

We start by denoting all unrenormalized quantities with a superscript $U$ and
then replace the bare coupling $\alpha^{U}$ with the renormalized strong
coupling $\alpha_{s}\equiv\alpha_{s}(\mu^{2})$, evaluated at the
renormalization scale $\mu^{2}$,

$\alpha^{U}\mu_{0}^{2\epsilon}S_{\epsilon}=\alpha_{s}\mu^{2\epsilon}\bigg{[}1-\frac{\beta_{0}}{\epsilon}\bigg{(}\frac{\alpha_{s}}{2\pi}\bigg{)}+\bigg{(}\frac{\beta_{0}^{2}}{\epsilon^{2}}-\frac{\beta_{1}}{2\epsilon}\bigg{)}\bigg{(}\frac{\alpha_{s}}{2\pi}\bigg{)}^{2}+\mathcal{O}(\alpha_{s}^{3})\bigg{]},$
(59)

where $S_{\epsilon}=(4\pi)^{\epsilon}\mathrm{e}^{-\epsilon\gamma_{E}}$ and
$\mu_{0}^{2}$ is the mass parameter in dimensional regularization introduced
to maintain a dimensionless coupling in the bare QCD Lagrangian density. The
explicit form of the first two $\beta$-function coefficients
$\beta_{0},\,\beta_{1}$ reads

$\displaystyle\beta_{0}$
$\displaystyle=\frac{11C_{A}}{6}-\frac{2T_{R}N_{F}}{3},$ (60)
$\displaystyle\beta_{1}$
$\displaystyle=\frac{17C_{A}^{2}}{6}-\frac{5C_{A}T_{R}N_{F}}{3}-C_{F}T_{R}N_{F},$
(61)

with the QCD colour factors,

$C_{A}=N,\quad C_{F}=\frac{N^{2}-1}{2N},\quad T_{R}=\frac{1}{2}.$ (62)

The effective coupling $\lambda$ is renormalized as follows,

$\lambda^{U}=\lambda\bigg{[}1-\frac{\beta_{0}}{\epsilon}\bigg{(}\frac{\alpha_{S}}{2\pi}\bigg{)}+\bigg{(}\frac{\beta_{0}^{2}}{\epsilon^{2}}-\frac{\beta_{1}}{\epsilon}\bigg{)}\bigg{(}\frac{\alpha_{S}}{2\pi}\bigg{)}^{2}+\mathcal{O}(\alpha_{S}^{3})\bigg{]}.$
(63)

The renormalized coefficients of the UV-finite but IR-divergent amplitudes can
be written in terms of the $i$-loop contribution to the unrenormalized
coefficients $\Omega^{(i),\,U}$ as

$\displaystyle\Omega^{(0)}$ $\displaystyle=\Omega^{(0),\,U},$ (64)
$\displaystyle\Omega^{(1)}$
$\displaystyle=S_{\epsilon}^{-1}\Omega^{(1),\,U}-\frac{3\beta_{0}}{2\epsilon}\Omega^{(0),\,U},$
(65) $\displaystyle\Omega^{(2)}$
$\displaystyle=S_{\epsilon}^{-2}\Omega^{(2),\,U}-\frac{5\beta_{0}}{2\epsilon}S_{\epsilon}^{-1}\Omega^{(1),\,U}-\bigg{(}\frac{5\beta_{1}}{4\epsilon}-\frac{15\beta_{0}^{2}}{8\epsilon^{2}}\bigg{)}\Omega^{(0),\,U}.$
(66)

For the remainder of the paper, we will set $\mu^{2}=\mu_{0}^{2}=M_{H}^{2}$
for simplicity.

###  5.2 Infrared factorization

Since the IR poles of $l$-loop amplitudes in gauge theories factorize in
colour-space in terms of lower loop amplitudes, the IR poles can be subtracted
multiplicatively as

$\mathbf{\Omega}({p},\epsilon)=\mathbfcal{Z}(\\{p\\};\epsilon)\,\,\mathbf{\Omega}^{\text{finite}}(\\{p\\}).$
(67)

Contrary to multiplicative renormalization of UV divergences, the bold
notation in (67) indicates that, in general, $\mathbfcal{Z}$ and
$\mathbf{\Omega}$ are operators and vectors in colour space, respectively. The
all-order nature of (67) makes this approach particularly advantageous for
generalizations to higher orders in perturbation theory.

Solving a renormalization group equation for $\mathbfcal{Z}$, one finds

$\mathbfcal{Z}(\epsilon,\\{p\\},\mu)=\mathbb{P}\,\text{exp}\Big{[}\int_{\mu}^{\infty}\frac{d\mu^{\prime}}{\mu^{\prime}}\mathbf{\Gamma}(\\{p\\},\mu^{\prime})\Big{]}=\sum_{l=0}^{\infty}\Big{(}\frac{\alpha_{s}}{2\pi}\Big{)}^{l}\mathbfcal{Z}^{(l)},$
(68)

where $\mathbb{P}$ is the path-ordering symbol, meaning that the colour
operators are ordered from left to right in decreasing values of
$\mu^{\prime}$. As originally proposed by Catani Catani:1998 , the anomalous-
dimension matrix $\mathbf{\Gamma}$ for amplitudes with $n$ QCD partons up to
two loops is entirely governed by the dipole colour correlations operator

$\mathbf{\Gamma}(\\{p\\},\mu)=\mathbf{\Gamma}_{dipole}(\\{p\\},\mu),$ (69)
$\mathbf{\Gamma}_{dipole}(\\{p\\},\mu)=\sum_{1\leq i<j\leq
n}\mathbf{T}_{i}^{a}\mathbf{T}_{j}^{a}\,\gamma^{\text{cusp}}(\alpha_{s})\log\Big{(}-\frac{\mu^{2}}{s_{ij}+i\eta}\Big{)}\,+\,\sum_{i=1}^{n}\gamma^{i}(\alpha_{s}),$
(70)

where $\gamma^{\text{cusp}}$ is the cusp anomalous dimension and $\gamma^{i}$
is the anomalous dimension of the $i$-th external particle, the latter
depending on the nature of the particle. The cusp anomalous dimension carries
the information about overlapping soft and collinear divergences, while
$\gamma^{i}$ only involves collinear divergences associated to the $i$-th
parton. The perturbative expansions up to two loops for the anomalous
dimensions are listed in Appendix A. The coupling constant is evaluated at the
renormalization scale $\alpha_{s}=\alpha_{s}(\mu)$.

The colour operators $\mathbf{T}^{a}_{i}$ are related to the $SU(N)$
generators and their action on the $i$-th coloured parton is defined following
the convention in Catani:1998  as:

$\displaystyle(\mathbf{T}_{i}^{a})_{b_{i}c_{i}}$
$\displaystyle=i\,f^{a}_{b_{i}c_{i}}\quad\text{for a gluon}\,,$
$\displaystyle(\mathbf{T}_{i}^{a})_{l_{i}k_{i}}$
$\displaystyle=+T^{a}_{l_{i}k_{i}}\quad\text{for a final(initial) state
quark(anti-quark)}\,,$ (71) $\displaystyle(\mathbf{T}_{i}^{a})_{l_{i}k_{i}}$
$\displaystyle=-T^{a}_{k_{i}l_{i}}\quad\text{for a initial(final) state
quark(anti-quark)}\,.$

where $i$ labels the particle on which the operator is acting. Importantly,
colour conservation can be rephrased as

$\displaystyle\sum_{i}\mathbf{T}_{i}^{a}=0.$ (72)

It follows from these definitions that the repeated action of one operator
evaluates to a Casimir,

$(\mathbf{T}_{i}^{a})^{2}=C_{i}\,,$ (73)

where $C_{i}=C_{A}$ if particle $i$ is a gluon and $C_{i}=C_{F}$ in case of
(anti-)quark.

We also define the expansions,

$\mathbf{\Gamma}_{dipole}=\sum_{l=0}^{\infty}\mathbf{\Gamma}_{l}\,\,\Big{(}\frac{\alpha_{s}}{2\pi}\Big{)}^{l+1},\quad\quad\,\Gamma^{\prime}=\frac{\partial\mathbf{\Gamma}_{dipole}}{\partial\log(\mu)}=\sum_{l=0}^{\infty}\Gamma^{\prime}_{l}\,\,\Big{(}\frac{\alpha_{s}}{2\pi}\Big{)}^{l+1},$
(74)

where one can drop the bold notation in the derivative because the resulting
operator is always diagonal in colour space:

$\displaystyle\Gamma^{\prime}=-\gamma^{\text{cusp}}(\alpha_{s})\sum_{i}C_{i}.$
(75)

In our specific case, expanding (67) up to two loops, IR divergences in the
renormalized two loop amplitudes can be expressed in terms of the renormalized
tree and one-loop amplitudes multiplied by appropriate operators

$\displaystyle\mathbf{\Omega}^{(0)}_{\text{finite}}$
$\displaystyle=\mathbf{\Omega}^{(0)},$ (76)
$\displaystyle\mathbf{\Omega}^{(1)}_{\text{finite}}$
$\displaystyle=\mathbf{\Omega}^{(1)}-\mathbf{I}_{\Omega}^{(1)}\,\mathbf{\Omega}^{(0)},$
(77) $\displaystyle\mathbf{\Omega}^{(2)}_{\text{finite}}$
$\displaystyle=\mathbf{\Omega}^{(2)}-\mathbf{I}_{\Omega}^{(2)}\,\mathbf{\Omega}^{(0)}-\mathbf{I}_{\Omega}^{(1)}\,\mathbf{\Omega}^{(1)},$
(78)

With the definitions given in (68)–(74), we can express the subtraction
operators in (77)–(78) as

$\displaystyle\mathbf{I}^{\text{SCET},(1)}$
$\displaystyle=\mathbfcal{Z}^{(1)},$ (79)
$\displaystyle\mathbf{I}^{\text{SCET},(2)}$
$\displaystyle=\mathbfcal{Z}^{(2)}-\big{(}\mathbfcal{Z}^{(1)}\big{)}^{2},$
(80)

where the process dependence is implicit and

$\displaystyle\mathbfcal{Z}^{(1)}$
$\displaystyle=\frac{\Gamma^{\prime}_{0}}{4\epsilon^{2}}+\frac{\mathbf{\Gamma_{0}}}{2\epsilon},$
(81) $\displaystyle\mathbfcal{Z}^{(2)}$
$\displaystyle=\frac{{\Gamma_{0}^{\prime}}^{2}}{32\epsilon^{4}}+\frac{\Gamma^{\prime}_{0}}{8\epsilon^{3}}\Big{(}\mathbf{\Gamma_{0}}-\frac{3}{2}\beta_{0}\Big{)}+\frac{\mathbf{\Gamma_{0}}}{8\epsilon^{2}}\Big{(}\mathbf{\Gamma}_{0}-2\beta_{0}\Big{)}+\frac{\Gamma^{\prime}_{1}}{16\epsilon^{2}}+\frac{\mathbf{\Gamma}_{1}}{4\epsilon}.$
(82)

It is now manifest that the IR operators in SCET scheme have only pole terms,
without any $\mathcal{O}(\epsilon^{0})$ contributions.

The cancellation of IR poles according to (78) is a strong analytic check on
multiloop amplitudes. We shall now follow the SCET approach to derive explicit
expressions for the subtraction operators.

###  5.3 SCET operators for three coloured partons: $ggg$ and $q\bar{q}g$

The definition of the finite parts of a two-loop amplitude is based on (80),
which we now work out for our process. We start from $\Gamma^{\prime}_{n}$,
which can be read off directly from (75) and the dependence on the nature of
the three partons is entirely contained in the sum over the Casimirs. On the
other hand, for the $\mathbf{\Gamma}_{n}$ coefficients, we can first notice
that an additional simplification occurs when only three coloured partons are
involved. In this case, the resulting amplitude is proportional to one single
colour structure, namely $f^{abc}$ is the colour factor for $H\rightarrow
g(a)+g(b)+g(c)$, while $T^{a}_{ij}$ is the factor for $H\rightarrow
q(i)+\bar{q}(j)+g(c)$.

As a consequence Catani:1998 , the dipole operator itself diagonalizes in
colour space. In fact, using colour conservation (72) and the property (73),
one can rewrite the product of two $\mathbf{T}^{a}$ operators as a sum over
Casimirs

$\displaystyle 2\,\mathbf{T}_{1}^{a}\mathbf{T}_{2}^{a}$
$\displaystyle=(\mathbf{T}_{3}^{a})^{2}-(\mathbf{T}_{1}^{a})^{2}-(\mathbf{T}_{1}^{a})^{2}=C_{3}-C_{1}-C_{2}\,,$
$\displaystyle 2\,\mathbf{T}_{2}^{a}\mathbf{T}_{3}^{a}$
$\displaystyle=(\mathbf{T}_{1}^{a})^{2}-(\mathbf{T}_{2}^{a})^{2}-(\mathbf{T}_{3}^{a})^{2}=C_{1}-C_{2}-C_{3}\,,$
$\displaystyle 2\,\mathbf{T}_{1}^{a}\mathbf{T}_{3}^{a}$
$\displaystyle=(\mathbf{T}_{2}^{a})^{2}-(\mathbf{T}_{1}^{a})^{2}-(\mathbf{T}_{3}^{a})^{2}=C_{2}-C_{1}-C_{3}.$
(83)

Accordingly, we can drop the bold notation and we find expressions relevant to
the three helicity coefficients $\alpha,\beta,\gamma$:

$\displaystyle\Gamma^{\alpha}_{n}$
$\displaystyle=\Gamma^{\beta}_{n}=-\frac{C_{A}}{2}\Big{(}L_{12}+L_{23}+L_{13}\Big{)}\gamma^{\text{cusp}}_{n}+3\gamma_{n}^{g},$
(84) $\displaystyle\Gamma^{\gamma}_{n}$
$\displaystyle=-C_{F}\,L_{12}\,\gamma^{\text{cusp}}_{n}-\frac{C_{A}}{2}\Big{(}-L_{12}+L_{23}+L_{13}\Big{)}\gamma^{\text{cusp}}_{n}+2\gamma_{n}^{q}+\gamma_{n}^{g},$
(85)

with

$\displaystyle L_{ij}=\log{\Big{(}-\frac{\mu^{2}}{s_{ij}+i\eta}\Big{)}}$ (86)

and the anomalous dimension coefficients defined as in Appendix A.

We stress here that there is an imaginary part arising from the logarithms
above whenever the corresponding invariant is positive, which depends on the
kinematical region considered.

###  5.4 Results for the helicity amplitudes

We computed the renormalized amplitudes for the decay process up to order
$\epsilon^{4}$ at one loop and up to order $\epsilon^{2}$ at two loops.
Working in the SCET subtraction scheme, we derived the finite remainder for
all the helicity amplitudes. Decomposed according to their colour structure,
they read

$\displaystyle\Omega^{(1)}_{{\rm finite}}$ $\displaystyle=\bigg{(}N\
A^{(1)}_{\Omega}+\frac{1}{N}\ B^{(1)}_{\Omega}+N_{F}C^{(1)}_{\Omega}\bigg{)},$
(87) $\displaystyle\Omega^{(2)}_{{\rm finite}}$
$\displaystyle=\bigg{(}N^{2}A^{(2)}_{\Omega}+N^{0}B^{(2)}_{\Omega}+\frac{1}{N^{2}}C^{(2)}_{\Omega}+\frac{N_{F}}{N}D^{(2)}_{\Omega}+NN_{F}E^{(2)}_{\Omega}+N_{F}^{2}F^{(2)}_{\Omega}\bigg{)}.$
(88)

The same structure holds for the renormalized amplitudes, though the
respective coefficients will still contain poles in $\epsilon$.

In the supplementary material, we provide the coefficients of the renormalized
amplitudes for the decay processes and of the finite remainder for the Higgs
decay kinematics as well as for all crossings Gehrmann:2002zr  relevant to
$H+$jet production processes.

###  5.5 Conversion to the Catani scheme

The two-loop helicity amplitudes for the decay of a Higgs boson into three
partons were first computed in Gehrmann:2012 . The authors obtained the finite
remainder by subtracting the IR singularities according to the original Catani
prescription Catani:1998 . Here we give the conversion rules between the two
subtraction schemes, which also served as a cross-check of our results.

Following closely the notation of Gehrmann:2012 , we write the subtraction
operators $\mathbf{I}_{\Omega}$ of (77) in Catani scheme as

$\displaystyle\mathbf{I}_{\alpha}^{\text{C},(1)}$
$\displaystyle=\mathbf{I}_{\beta}^{\text{C},(1)}=-\frac{\mathrm{e}^{\epsilon\gamma}}{2\Gamma(1-\epsilon)}\bigg{[}N\bigg{(}\frac{1}{\epsilon^{2}}+\frac{\beta_{0}}{N\epsilon}\bigg{)}\big{(}{\tt
S}_{12}+{\tt S}_{23}+{\tt S}_{13}\big{)}\bigg{]},$
$\displaystyle\mathbf{I}_{\gamma}^{\text{C},(1)}$
$\displaystyle=-\frac{\mathrm{e}^{\epsilon\gamma}}{2\Gamma(1-\epsilon)}\bigg{[}N\bigg{(}\frac{1}{\epsilon^{2}}+\frac{3}{4\epsilon}+\frac{\beta_{0}}{2N\epsilon}\bigg{)}\big{(}{\tt
S}_{23}+{\tt
S}_{13}\big{)}-\frac{1}{N}\bigg{(}\frac{1}{\epsilon^{2}}+\frac{3}{2\epsilon}\bigg{)}{\tt
S}_{12}\bigg{]},$ (89)

with,

${\tt S}_{ij}=\bigg{(}-\frac{\mu^{2}}{s_{ij}}\bigg{)}^{\epsilon}.$ (90)

The second-order operator can be built starting from the one-loop operator as

$\displaystyle\mathbf{I}_{\Omega}^{\text{C},(2)}=$
$\displaystyle-\frac{1}{2}\mathbf{I}^{\text{C},(1)}_{\Omega}(\epsilon)\,\mathbf{I}^{\text{C},(1)}_{\Omega}(\epsilon)-\frac{\beta_{0}}{\epsilon}\mathbf{I}^{\text{C},(1)}_{\Omega}(\epsilon)$
$\displaystyle+\mathrm{e}^{-\epsilon\gamma}\frac{\Gamma(1-2\epsilon)}{\Gamma(1-\epsilon)}\bigg{(}\frac{\beta_{0}}{\epsilon}+K\bigg{)}\mathbf{I}^{\text{C},(1)}_{\Omega}(2\epsilon)+\mathbf{H}^{(2)}_{\Omega}(\epsilon)\,,$
(91)

where we introduced the constant

$K=\bigg{(}\frac{67}{18}-\frac{\pi^{2}}{6}\bigg{)}C_{A}-\frac{10}{9}T_{R}N_{F}.$
(92)

The remaining term in (91) involves the operator
$\mathbf{H}^{(2)}_{\Omega}(\epsilon)$ and produces only a single pole in
$\epsilon$. Its explicit form is

$\mathbf{H}^{(2)}_{\Omega}(\epsilon)=\frac{\mathrm{e}^{\epsilon\gamma}}{4\epsilon\Gamma(1-\epsilon)}H^{(2)}_{\Omega}.$
(93)

The constant $H^{(2)}_{\Omega}$ is renormalization scheme and process
dependent and in our case it reads

$\displaystyle H^{(2)}_{\alpha}$ $\displaystyle=H^{(2)}_{\beta}=3H^{(2)}_{g},$
(94) $\displaystyle H^{(2)}_{\gamma}$
$\displaystyle=2H^{(2)}_{q}+H^{(2)}_{g},$ (95)

where in the $\overline{\text{MS}}$ scheme the constants
$H^{(2)}_{q},\,H^{(2)}_{g}$ are

$\displaystyle H^{(2)}_{q}$ $\displaystyle=\left({7\over
4}\zeta_{3}+{\frac{409}{864}}-{\frac{11\pi^{2}}{96}}\right)N^{2}+\left(-{1\over
4}\zeta_{3}-{41\over 108}-{\pi^{2}\over 96}\right)+\left(-{3\over
2}\zeta_{3}-{3\over 32}+{\pi^{2}\over 8}\right){1\over N^{2}}$
$\displaystyle+\left({\pi^{2}\over 48}-{25\over
216}\right){(N^{2}-1)N_{F}\over N}\;,$ (96) $\displaystyle H^{(2)}_{g}$
$\displaystyle=\left(\frac{1}{2}\zeta_{3}+{\frac{5}{12}}+{\frac{11\pi^{2}}{144}}\right)N^{2}+{\frac{5}{27}}\,N_{F}^{2}+\left(-{\frac{{\pi}^{2}}{72}}-{\frac{89}{108}}\right)NN_{F}-\frac{N_{F}}{4N}.$
(97)

When describing the SCET subtraction scheme, we pointed out that the
subtraction operators have no finite $\mathcal{O}(\epsilon^{0})$ contribution.

In contrast, the Catani operators contain coefficients at higher order in the
dimensional regulator, generated by the $\epsilon$ expansion of the resummed
coefficient defined in Eq. (90). The tree-level amplitude is finite. One-loop
amplitudes have poles starting from order $\mathcal{O}(1/\epsilon^{2})$ and
two-loop amplitudes have poles starting from order
$\mathcal{O}(1/\epsilon^{4})$. We indicate with $\mathbf{\Omega}_{n}^{(l)}$,
$\mathbf{I}_{n}^{\text{SCET},(l)}$ and $\mathbf{I}_{n}^{\text{C},(l)}$ the
coefficients of order $\epsilon^{n}$ of the renormalized amplitude, the SCET
operator and the Catani operator, respectively.

From the subtraction formulae (77)–(78), it is easy to obtain the following
conversion rules for the finite remainders

$\displaystyle\mathbf{\Omega}_{\text{finite}}^{\text{SCET},(0)}-\mathbf{\Omega}_{\text{finite}}^{\text{C},(0)}$
$\displaystyle=0,$ (98)
$\displaystyle\mathbf{\Omega}_{\text{finite}}^{\text{SCET},(1)}-\mathbf{\Omega}_{\text{finite}}^{\text{C},(1)}$
$\displaystyle=\mathbf{I}_{0}^{\text{C},(1)}\mathbf{\Omega}^{(0)},$ (99)
$\displaystyle\mathbf{\Omega}_{\text{finite}}^{\text{SCET},(2)}-\mathbf{\Omega}_{\text{finite}}^{\text{C},(2)}$
$\displaystyle=\mathbf{I}_{0}^{\text{C},(2)}\mathbf{\Omega}^{(0)}+\mathbf{I}_{2}^{\text{C},(1)}\mathbf{\Omega}_{-2}^{(1)}+\mathbf{I}_{1}^{\text{C},(1)}\mathbf{\Omega}_{-1}^{(1)}+\mathbf{I}_{0}^{\text{C},(1)}\mathbf{\Omega}_{0}^{(1)}.$
(100)

In deriving the above rules, we made use of the fact that
$\mathbf{I}^{\text{C},(1)}$ and $\mathbf{I}^{\text{SCET},(1)}$ have the same
pole structure. Consequently, terms multiplying
$(\mathbf{I}^{\text{C},(1)}_{-1}-\mathbf{I}^{\text{SCET},(1)}_{-1})$ or
$(\mathbf{I}^{\text{C},(1)}_{-2}-\mathbf{I}^{\text{SCET},(1)}_{-2})$ in (100)
vanish. This can be easily understood by inspecting the origin of the poles in
the one-loop cancellation in (77).

For both the decay and the production kinematics, we verified that the result
given in Gehrmann:2012  is correctly reproduced converting our finite
remainder with the above rules.

##  6 Analytic Continuation

So far, we have described our calculation making explicit reference to decay
processes, for which all the kinematic invariants are positive. We considered
the decay of a Higgs boson into three gluons, $H\rightarrow ggg$, and into a
quark-antiquark pair and gluon, $H\rightarrow q\bar{q}g$. In view of
applications to LHC physics, we are interested in the production regions, in
which a Higgs is produced together with a parton: $gg\rightarrow Hg$,
$qg\rightarrow Hq$, $\bar{q}g\rightarrow H\bar{q}$ and $q\bar{q}\rightarrow
Hg$.

The general strategy for performing the analytic continuation for the MPLs
appearing in $2\to 2$ scattering involving 4-point functions with one external
off-shell leg and massless propagators was outlined in detail in
Gehrmann:2002zr . Our aim is to describe how this strategy can be applied to
the process at hand. Referring to Fig. 1 of reference Gehrmann:2002zr  and
using the same labels for the various kinematic regions, the goal is to find a
procedure to analytically continue MPLs from the decay region (1a) to the
three production regions, (2a), (3a), (4a). Whenever a particle is crossed
from the initial state to the final state (or vice versa), two of the
invariants become negative and one remains positive, representing the centre-
of-mass energy of the incoming partons. We recall here that results in the
decay region (1a) are expressed as MPLs of the variables $y$ and $z$ defined
in (6), which fulfil the constraints $0<z<1\,,\,\,0<y<1-z$. Within these
bounds, our set of MPLs are real and no branching point is crossed.

To describe the analytic continuation to the scattering kinematics, let us
consider the case of region (3a). In (3a), the kinematic constraints are
$z<0\,,\,\,1-z<y<+\infty$ and MPLs must be evaluated across a branch cut,
developing an imaginary part. Physically, the particle with momentum $p_{2}$
is crossed and the process is $p_{1}+p_{3}\rightarrow p_{2}+p_{4}$. Crucially,
there exists a change of variables which maps region (3a) _linearly_ back into
the decay region (1a) and can be implemented analytically on our MPLs. In
fact, by defining

$\displaystyle v\equiv 1/y\,,\quad\quad\quad u\equiv-z/y\,,$ (101)

the new variables $u$ and $v$ satisfy again $0<u<1\,,\,\,0<v<1-u$. By re-
expressing the MPLs in terms of $u$ and $v$, the imaginary part can be made
explicit in terms of multiple zeta values and real-valued MPLs. The same
manipulations can be performed in regions (2a) and (4a) with different
definitions of the $(u,v)$ variables. In general, the $v$ variable is the
reciprocal of the centre-of-mass energy and so its definition depends on which
particle is crossed from the final to the initial state.

Instead of performing the analytic continuation in all the three regions (2a),
(3a) and (4a), we found it simpler to first consider suitable crossings of the
amplitude in the decay region and in a second step continue these crossed
amplitudes only to the region (3a). Explicitly, for the case of three gluons,
the two independent amplitudes in the decay kinematics are

$\displaystyle\mathcal{M}^{+++}:\quad\quad H$ $\displaystyle\rightarrow
g_{+}(p_{1})+g_{+}(p_{2})+g_{+}(p_{3})\,,$ (102)
$\displaystyle\mathcal{M}^{++-}:\quad\quad H$ $\displaystyle\rightarrow
g_{+}(p_{1})+g_{+}(p_{2})+g_{-}(p_{3})\,.$ (103)

In the production region, there are eight different helicity configurations.
Thanks to parity symmetry, we can limit ourselves to consider just four of
them (see the left column of Table 2) and relate them to the other four (right
column).

Process | Parity Related
---|---
$\mathcal{M}^{--,+}_{gg,g}:\,\,g_{-}(p_{1})+g_{-}(p_{3})\rightarrow H+g_{+}(p_{2})\,\quad$ |  $\quad\quad\,g_{+}(p_{1})+g_{+}(p_{3})\rightarrow H+g_{-}(p_{2})$,
$\mathcal{M}^{-+,+}_{gg,g}:\,\,g_{-}(p_{1})+g_{+}(p_{3})\rightarrow H+g_{+}(p_{2})\,\quad$ |  $\quad\quad\,g_{+}(p_{1})+g_{-}(p_{3})\rightarrow H+g_{-}(p_{2})$,
$\mathcal{M}^{--,-}_{gg,g}:\,\,g_{-}(p_{1})+g_{-}(p_{3})\rightarrow H+g_{-}(p_{2})\,\quad$ |  $\quad\quad\,g_{+}(p_{1})+g_{+}(p_{3})\rightarrow H+g_{+}(p_{2})$,
$\mathcal{M}^{+-,+}_{gg,g}:\,\,g_{+}(p_{1})+g_{-}(p_{3})\rightarrow H+g_{+}(p_{2})\,\quad$ |  $\quad\quad\,g_{-}(p_{1})+g_{+}(p_{3})\rightarrow H+g_{-}(p_{2})$.
Table 2: Different kinematical crossings of the $H\to ggg$ amplitudes. The
comma is used to separate the helicities of initial and final state partons.

Note that in continuing to the region (3a), the momenta $p_{1}$ and $p_{3}$
are always in the initial state. We computed the first three amplitudes with a
combination of crossings and analytic continuation as follows (an extra
helicity flip due to time reversal is always understood after continuation to
(3a)):

$\displaystyle\mathcal{M}^{+++}_{ggg}\,$ $\displaystyle\xrightarrow[]{\text{to
(3a)}}\,\mathcal{M}^{--,+}_{gg,g}\,,$ (104)
$\displaystyle\mathcal{M}^{++-}_{ggg}\,$ $\displaystyle\xrightarrow[]{\text{to
(3a)}}\,\mathcal{M}^{-+,+}_{gg,g}\,,$ (105)
$\displaystyle\mathcal{M}^{++-}_{ggg}\,\xrightarrow[]{p_{2}\leftrightarrow
p_{3}}\,$ $\displaystyle\mathcal{M}^{+-+}_{ggg}\,\xrightarrow[]{\text{to
(3a)}}\,\mathcal{M}^{--,-}_{gg,g}\,.$ (106)

The fourth amplitude, $\mathcal{M}^{+-,+}_{gg,g}$, can be derived from
$\mathcal{M}^{-+,+}_{gg,g}$ by the crossing $p_{1}\leftrightarrow p_{3}$,
which implies $v\rightarrow v$ and $u\rightarrow 1-u-v$. Since no branch cut
is crossed under this transformation, no new analytic continuation is needed
and we can conclude that

$\displaystyle\mathcal{M}^{+-,+}_{gg,g}(v,u)=\mathcal{M}^{-+,+}_{gg,g}(v,u)|_{u\rightarrow
1-u-v}\,.$ (107)

Let us consider now the decay into a quark-antiquark pair and a gluon:

$\displaystyle\mathcal{M}^{LR+}:\quad\quad H$ $\displaystyle\rightarrow
q_{L}(p_{1})+\bar{q}_{R}(p_{2})+g_{+}(p_{3})\,.$ (108)

In the production region there are 12 non-zero independent helicity
configurations, but only 3 of them need to be computed, while the others are
related by parity and charge conjugation. Our choice of the independent
configurations is according to the left column of Table 3.

Process | C, P, CP
---|---
$\mathcal{M}^{R-,R}_{qg,q}:\,\,q_{R}(p_{1})+g_{-}(p_{3})\rightarrow H+q_{R}(p_{2})\,\quad$ | $\quad\quad\,\bar{q}_{L}(p_{1})+g_{+}(p_{3})\rightarrow H+\bar{q}_{L}(p_{2})$
| $\quad\quad\,q_{L}(p_{1})+g_{+}(p_{3})\rightarrow H+q_{L}(p_{2})$
| $\quad\quad\,\bar{q}_{R}(p_{1})+g_{-}(p_{3})\rightarrow
H+\bar{q}_{R}(p_{2})$
$\mathcal{M}^{L-,L}_{\bar{q}g,\bar{q}}:\,\,\bar{q}_{L}(p_{1})+g_{-}(p_{3})\rightarrow H+\bar{q}_{L}(p_{2})\,\quad$ | $\quad\quad\,q_{R}(p_{1})+g_{+}(p_{3})\rightarrow H+q_{R}(p_{2})$
| $\quad\quad\,\bar{q}_{R}(p_{1})+g_{+}(p_{3})\rightarrow
H+\bar{q}_{R}(p_{2})$
| $\quad\quad\,q_{L}(p_{1})+g_{-}(p_{3})\rightarrow H+q_{L}(p_{2})$
$\mathcal{M}^{RL,+}_{q\bar{q},g}:\,\,q_{R}(p_{1})+\bar{q}_{L}(p_{3})\rightarrow H+g_{+}(p_{2})\,\quad$ | $\quad\quad\,q_{R}(p_{1})+\bar{q}_{L}(p_{3})\rightarrow H+g_{-}(p_{2})$
| $\quad\quad\,q_{L}(p_{1})+\bar{q}_{R}(p_{3})\rightarrow H+g_{-}(p_{2})$
| $\quad\quad\,\bar{q}_{R}(p_{1})+q_{L}(p_{3})\rightarrow H+g_{+}(p_{2})$
Table 3: Different kinematic crossings of the $H\to q\bar{q}g$ amplitudes. The
comma is used to separate the helicities of initial and final state partons.

The three amplitudes were computed with a combination of crossings and
analytic continuation as follows (an additional helicity flip due to time
reversal is again understood after continuation to (3a)),

$\displaystyle\mathcal{M}^{LR+}_{q\bar{q}g}\,$
$\displaystyle\xrightarrow[]{\text{to (3a)}}\,\mathcal{M}^{R-,R}_{qg,q}$ (109)
$\displaystyle\mathcal{M}^{LR+}_{q\bar{q}g}\,\xrightarrow[]{p_{1}\leftrightarrow
p_{2}}\,$ $\displaystyle\mathcal{M}^{RL+}_{\bar{q}qg}\,\xrightarrow[]{\text{to
(3a)}}\,\mathcal{M}^{L-,L}_{\bar{q}g,\bar{q}}$ (110)
$\displaystyle\mathcal{M}^{LR+}_{q\bar{q}g}\,\xrightarrow[]{p_{2}\leftrightarrow
p_{3}}\,$ $\displaystyle\mathcal{M}^{L+R}_{qg\bar{q}}\,\xrightarrow[]{\text{to
(3a)}}\,\mathcal{M}^{RL,+}_{q\bar{q},g}\,.$ (111)

Finally, we stress that each crossing must also be applied to the spinor
prefactors in (30). We have cross-checked all amplitudes in all helicity
configurations against previously published results Gehrmann:2012  at the
level of the finite remainder up to $\mathcal{O}(\epsilon^{0})$ and found
perfect agreement. Conventions and normalization of helicity amplitudes have
been fixed by numerical evaluation of our result, in the $q\bar{q}g$ channel,
using OpenLoops 2 Buccioni:2019sur .

##  7 Conclusions

In this paper, we presented the calculation of the two loop corrections to the
helicity amplitudes for the $H\to ggg$ and $H\to q\bar{q}g$ up to order
$\mathcal{O}(\epsilon^{2})$ in the large top Higgs effective field theory.
These amplitudes constitute the first missing ingredient towards the
calculation of $H+$jet production to N3LO at the LHC, as they are required to
properly define the finite remainder of the corresponding three loop virtual
corrections. We follow a standard approach to compute the helicity amplitudes.
We start by decomposing the amplitude in a basis of independent tensor
structures, and we use spinor-helicity to express the helicity amplitudes in
terms of linear combinations of the corresponding scalar form factors. Next,
we derived a canonical basis for the relevant master integrals and, through
the method of differential equations, provided their solution up to weight six
in MPLs. We verified that our results for the finite remainder of the
amplitude match the literature Gehrmann:2012 , and analytically continued the
helicity amplitudes to all kinematic regions. The infrared structure of the
result was inspected in the frameworks of SCET and the Catani infrared
factorization formula. This updated result marks the first step towards
computing the N3LO corrections to $H$+jet production.

## Acknowledgments

This work was supported in part by the Excellence Cluster ORIGINS funded by
the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany’s Excellence Strategy – EXC-2094-390783311, by the Swiss National
Science Foundation (SNF) under contract 200020-204200, and by the European
Research Council (ERC) under the European Union’s research and innovation
programme grant agreements 949279 (ERC Starting Grant HighPHun) and 101019620
(ERC Advanced Grant TOPUP).

##  Appendix A Anomalous dimensions

In this appendix, we give the perturbative coefficients for the cusp anomalous
dimension $\gamma^{\text{cusp}}$, the quark anomalous dimension $\gamma_{q}$
and the gluon anomalous dimension $\gamma_{g}$ up to two loops,
$\mathcal{O}(\alpha_{s}^{2})$. The perturbative expansion reads

$\displaystyle\gamma^{\text{cusp}}=\sum_{n=0}^{\infty}\gamma_{n}^{\text{cusp}}\Big{(}\frac{\alpha_{s}}{2\pi}\Big{)}^{n+1}\,\quad\,\gamma^{i}=\sum_{n=0}^{\infty}\gamma_{n}^{i}\Big{(}\frac{\alpha_{s}}{2\pi}\Big{)}^{n+1}$
(112)

with $i=q,\overline{q},g$.
The cusp anomalous dimension was computed at two loops in korchemksy:1987

$\displaystyle\gamma_{0}^{\text{cusp}}$ $\displaystyle=2,$
$\displaystyle\gamma_{1}^{\text{cusp}}$
$\displaystyle=\Big{(}\frac{67}{9}-\frac{\pi^{2}}{3}\Big{)}C_{A}-\frac{10}{9}N_{f}\,,$
(113)

while for the (anti-)quark anomalous dimension we have vanneerven:1986 ;
matsuura:1988sm ,

$\displaystyle\gamma_{0}^{q}$ $\displaystyle=-\frac{3}{2}\,C_{F},$
$\displaystyle\gamma_{1}^{q}$
$\displaystyle=C_{F}^{2}\Big{(}-\frac{3}{8}+\frac{\pi^{2}}{2}-6\zeta_{3}\Big{)}+C_{F}C_{A}\Big{(}-\frac{961}{216}-\frac{11}{24}\pi^{2}+\frac{13}{2}\zeta_{3}\Big{)}+C_{F}N_{F}\Big{(}\frac{65}{108}+\frac{\pi^{2}}{12}\Big{)}\,,$
(114)

and finally, for the gluon, Harlander:2000 ,

$\displaystyle\gamma_{0}^{g}$ $\displaystyle=-\beta_{0},$
$\displaystyle\gamma_{1}^{g}$
$\displaystyle=C_{A}^{2}\Big{(}-\frac{173}{27}+\frac{11}{72}\pi^{2}+\frac{\zeta_{3}}{2}\Big{)}+C_{A}N_{F}\Big{(}\frac{32}{27}-\frac{\pi^{2}}{36}\Big{)}+\frac{1}{2}C_{F}N_{F}.$
(115)

##  Appendix B One-loop master integrals

In order to obtain a result for the two-loop amplitudes up to
$\mathcal{O}(\epsilon^{2})$, we ought to calculate the one-loop master
integrals and amplitudes up to $\mathcal{O}(\epsilon^{4})$. Denoting

$\displaystyle
I_{a_{1},...,a_{4}}=\int\frac{e^{\epsilon\gamma_{E}}}{i\pi^{\frac{D}{2}}}\frac{1}{(k)^{a_{i}}(k-p_{1})^{a_{i}}(k-p_{1}-p_{2})^{a_{i}}(k-p_{1}-p_{2}-p_{3})^{a_{i}}}\,d^{D}k\,,$
(116)

a canonical basis for this process is (for simplicity we set $M^{2}_{H}=-1$
again)

$\displaystyle\vec{I}=\begin{pmatrix}-(1-y-z)\epsilon\,I_{2,0,1,0}\\\
-z\epsilon\,I_{0,2,0,1}\\\ \epsilon\,I_{2,0,0,1}\\\
(1-y-z)z\epsilon^{2}I_{1,1,1,1}\end{pmatrix}\,,$ (117)

representing 3 bubble diagrams with a squared propagator in the variables
$s_{12}$, $s_{23}$ and $M_{H}^{2}$, respectively, and ultimately the box
diagram. The former three can be related to simple bubbles by IBPs:

$\displaystyle I_{2,0,1,0}$
$\displaystyle=-\frac{1-2\epsilon}{-1+y+z}\mathrm{Bub}_{s_{12}}$
$\displaystyle=\frac{(1-y-z)^{-1-\epsilon}}{\epsilon}C(\epsilon)\,,$ (118)
$\displaystyle I_{0,2,0,1}$
$\displaystyle=\frac{1-2\epsilon}{z}\mathrm{Bub}_{s_{23}}$
$\displaystyle=\frac{(z)^{-1-\epsilon}}{\epsilon}C(\epsilon)\,,$ (119)
$\displaystyle I_{2,0,0,1}$
$\displaystyle=(1-2\epsilon)\mathrm{Bub}_{m_{H}^{2}}$
$\displaystyle=\frac{1}{\epsilon}C(\epsilon)\,,$ (120)

where we used the well-known expression for the bubble,

$\displaystyle\mathrm{Bub}_{s}=s^{-\epsilon}\frac{1}{\epsilon(1-2\epsilon)}\left(-e^{\epsilon\gamma_{E}}\frac{\Gamma(1+\epsilon)\Gamma(1-\epsilon)^{2}}{\Gamma(1-2\epsilon)}\right)=s^{-\epsilon}\frac{1}{\epsilon(1-2\epsilon)}C(\epsilon)\,.$
(121)

The differential equations take the form $\partial_{y}\vec{I}=A\vec{I}$ and
$\partial_{z}\vec{I}=B\vec{I}$ where

$\displaystyle A$ $\displaystyle=\begin{pmatrix}\frac{1}{x}&0&0&0\\\
0&0&0&0\\\ 0&0&0&0\\\
\left(\frac{2}{y}-\frac{2}{y+z}\right)&\frac{2}{y}&\left(-\frac{2}{y}+\frac{2}{y+z}\right)&\left(\frac{1}{y}+\frac{1}{x}\right)\end{pmatrix}\;,$
(122) $\displaystyle B$ $\displaystyle=\begin{pmatrix}\frac{1}{x}&0&0&0\\\
0&-\frac{1}{z}&0&0\\\ 0&0&0&0\\\
-\frac{2}{y+z}&\frac{-2}{-1+z}&\left(\frac{2}{-1+z}+\frac{2}{y+z}\right)&\left(-\frac{1}{z}+\frac{1}{x}\right)\end{pmatrix}$
(123)

and $x=1-y-z$. They confirm the scaling of the bubbles given in (118, 119,
120), and allow us to compute the box. The boundary condition obtained by
multiplying the first equation by $y$ and sending it to 0 reads

$\displaystyle
0=2I_{2,0,1,0}\big{|}_{y=0}+2I_{0,2,0,1}\big{|}_{y=0}-2I_{2,0,0,1}\big{|}_{y=0}+I_{1,1,1,1}\big{|}_{y=0}$
(124)

and is the only condition needed to find a unique solution for the box
integral. We could have obtained this requirement also via considering the
limit $y\to 0$ for $\epsilon<0$. The solution to the DE in $y$ near the point
$y=0$ is

$\displaystyle e^{\epsilon\log{y}\lim_{y\to
0}[yA_{y}]}\vec{I}\big{|}_{y=0}=\begin{pmatrix}1&0&0&0\\\ 0&1&0&0\\\
0&0&1&0\\\
2(-1+y^{\epsilon})&2(-1+y^{\epsilon})&-2(-1+y^{\epsilon})&y^{\epsilon}\end{pmatrix}\vec{I}\big{|}_{y=0}\,.$
(125)

Clearly requiring the non-appearance of terms $y^{\epsilon}$ in the solution
replicates the condition (124). The strength of the second approach is that it
applies to regular as well as singular limits. Expanded solutions for all four
integrals can be found in the supplementary material.

## References

  * (1) ATLAS collaboration, _Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC_ , _Phys. Lett. B_ 716 (2012) 1 [1207.7214].
  * (2) CMS collaboration, _Observation of a New Boson at a Mass of 125 GeV with the CMS Experiment at the LHC_ , _Phys. Lett. B_ 716 (2012) 30 [1207.7235].
  * (3) A. Dainese, M. Mangano, A.B. Meyer, A. Nisati, G. Salam and M.A. Vesterinen, eds., _Report on the Physics at the HL-LHC,and Perspectives for the HE-LHC_ , vol. 7/2019 of _CERN Yellow Reports: Monographs_ , CERN, Geneva, Switzerland (2019), 10.23731/CYRM-2019-007.
  * (4) G. Heinrich, _Collider Physics at the Precision Frontier_ , _Phys. Rept._ 922 (2021) 1 [2009.00516].
  * (5) S. Dawson, _Radiative corrections to Higgs boson production_ , _Nucl. Phys. B_ 359 (1991) 283.
  * (6) M. Spira, A. Djouadi, D. Graudenz and P.M. Zerwas, _Higgs boson production at the LHC_ , _Nucl. Phys. B_ 453 (1995) 17 [hep-ph/9504378].
  * (7) M. Czakon, R.V. Harlander, J. Klappert and M. Niggetiedt, _Exact Top-Quark Mass Dependence in Hadronic Higgs Production_ , _Phys. Rev. Lett._ 127 (2021) 162002 [2105.04436].
  * (8) F. Wilczek, _Decays of heavy vector mesons into higgs particles_ , _Phys. Rev. Lett._ 39 (1977) 1304.
  * (9) M. Shifman, A. Vainshtein and V. Zakharov, _Remarks on higgs-boson interactions with nucleons_ , _Physics Letters B_ 78 (1978) 443.
  * (10) T. Inami, T. Kubota and Y. Okada, _Effective Gauge Theory and the Effect of Heavy Quarks in Higgs Boson Decays_ , _Z. Phys. C_ 18 (1983) 69.
  * (11) C. Anastasiou, C. Duhr, F. Dulat, F. Herzog and B. Mistlberger, _Higgs Boson Gluon-Fusion Production in QCD at Three Loops_ , _Phys. Rev. Lett._ 114 (2015) 212001 [1503.06056].
  * (12) B. Mistlberger, _Higgs boson production at hadron colliders at N 3LO in QCD_, _JHEP_ 05 (2018) 028 [1802.00833].
  * (13) C. Anastasiou, C. Duhr, F. Dulat, E. Furlan, T. Gehrmann, F. Herzog, A. Lazopoulos and B. Mistlberger, _High precision determination of the gluon fusion Higgs boson cross-section at the LHC_ , _JHEP_ 05 (2016) 058 [1602.00695].
  * (14) X. Chen, T. Gehrmann, E.W.N. Glover, A. Huss, B. Mistlberger and A. Pelloni, _Fully Differential Higgs Boson Production to Third Order in QCD_ , _Phys. Rev. Lett._ 127 (2021) 072002 [2102.07607].
  * (15) F. Bishara, U. Haisch, P.F. Monni and E. Re, _Constraining Light-Quark Yukawa Couplings from Higgs Distributions_ , _Phys. Rev. Lett._ 118 (2017) 121801 [1606.09253].
  * (16) R. Boughezal, F. Caola, K. Melnikov, F. Petriello and M. Schulze, _Higgs boson production in association with a jet at next-to-next-to-leading order_ , _Phys. Rev. Lett._ 115 (2015) 082003 [1504.07922].
  * (17) R. Boughezal, C. Focke, W. Giele, X. Liu and F. Petriello, _Higgs boson production in association with a jet at NNLO using jettiness subtraction_ , _Phys. Lett. B_ 748 (2015) 5 [1505.03893].
  * (18) F. Caola, K. Melnikov and M. Schulze, _Fiducial cross sections for Higgs boson production in association with a jet at next-to-next-to-leading order in QCD_ , _Phys. Rev. D_ 92 (2015) 074032 [1508.02684].
  * (19) X. Chen, J. Cruz-Martinez, T. Gehrmann, E.W.N. Glover and M. Jaquier, _NNLO QCD corrections to Higgs boson production at large transverse momentum_ , _JHEP_ 10 (2016) 066 [1607.08817].
  * (20) J.M. Campbell, R.K. Ellis and S. Seth, _H + 1 jet production revisited_ , _JHEP_ 10 (2019) 136 [1906.01020].
  * (21) J.M. Lindert, K. Kudashkin, K. Melnikov and C. Wever, _Higgs bosons with large transverse momentum at the LHC_ , _Phys. Lett. B_ 782 (2018) 210 [1801.08226].
  * (22) S.P. Jones, M. Kerner and G. Luisoni, _Next-to-Leading-Order QCD Corrections to Higgs Boson Plus Jet Production with Full Top-Quark Mass Dependence_ , _Phys. Rev. Lett._ 120 (2018) 162001 [1802.00349].
  * (23) X. Chen, A. Huss, S.P. Jones, M. Kerner, J.N. Lang, J.M. Lindert and H. Zhang, _Top-quark mass effects in H+jet and H+2 jets production_ , _JHEP_ 03 (2022) 096 [2110.06953].
  * (24) R. Bonciani, V. Del Duca, H. Frellesvig, J.M. Henn, F. Moriello and V.A. Smirnov, _Two-loop planar master integrals for Higgs $\to 3$ partons with full heavy-quark mass dependence_, _JHEP_ 12 (2016) 096 [1609.06685].
  * (25) R. Bonciani, V. Del Duca, H. Frellesvig, J.M. Henn, M. Hidding, L. Maestri, F. Moriello, G. Salvatori and V.A. Smirnov, _Evaluating a family of two-loop non-planar master integrals for Higgs + jet production with full heavy-quark mass dependence_ , _JHEP_ 01 (2020) 132 [1907.13156].
  * (26) R. Bonciani, V. Del Duca, H. Frellesvig, M. Hidding, V. Hirschi, F. Moriello, G. Salvatori, G. Somogyi and F. Tramontano, _Next-to-leading-order QCD Corrections to Higgs Production in association with a Jet_ , 2206.10490.
  * (27) F. Caola, W. Chen, C. Duhr, X. Liu, B. Mistlberger, F. Petriello, G. Vita and S. Weinzierl, _The Path forward to N 3LO_, in _2022 Snowmass Summer Study_ , 3, 2022 [2203.06730].
  * (28) C. Duhr, _Hopf algebras, coproducts and symbols: an application to Higgs boson amplitudes_ , _JHEP_ 08 (2012) 043 [1203.0454].
  * (29) T. Gehrmann, M. Jaquier, E.W.N. Glover and A. Koukoutsakis, _Two-loop QCD corrections to the helicity amplitudes for H $\rightarrow$ 3 partons_, _Journal of High Energy Physics_ 2012 (2012) .
  * (30) X. Chen, T. Gehrmann, E.W.N. Glover, A. Huss, Y. Li, D. Neill, M. Schulze, I.W. Stewart and H.X. Zhu, _Precise QCD Description of the Higgs Boson Transverse Momentum Spectrum_ , _Phys. Lett. B_ 788 (2019) 425 [1805.00736].
  * (31) T. Peraro and L. Tancredi, _Physical projectors for multi-leg helicity amplitudes_ , _JHEP_ 07 (2019) 114 [1906.03298].
  * (32) T. Peraro and L. Tancredi, _Tensor decomposition for bosonic and fermionic scattering amplitudes_ , _Phys. Rev. D_ 103 (2021) 054042 [2012.00820].
  * (33) K.G. Chetyrkin and F.V. Tkachov, _Integration by Parts: The Algorithm to Calculate beta Functions in 4 Loops_ , _Nucl. Phys. B_ 192 (1981) 159.
  * (34) S. Laporta, _High precision calculation of multiloop Feynman integrals by difference equations_ , _Int. J. Mod. Phys. A_ 15 (2000) 5087 [hep-ph/0102033].
  * (35) A.V. Kotikov, _Differential equations method: New technique for massive Feynman diagrams calculation_ , _Phys. Lett. B_ 254 (1991) 158.
  * (36) A.V. Kotikov, _Differential equations method: The Calculation of vertex type Feynman diagrams_ , _Phys. Lett. B_ 259 (1991) 314.
  * (37) A.V. Kotikov, _Differential equation method: The Calculation of N point Feynman diagrams_ , _Phys. Lett. B_ 267 (1991) 123.
  * (38) T. Gehrmann and E. Remiddi, _Differential equations for two loop four point functions_ , _Nucl. Phys. B_ 580 (2000) 485 [hep-ph/9912329].
  * (39) T. Gehrmann and E. Remiddi, _Two loop master integrals for $\gamma^{*}\to$ 3 jets: The Planar topologies_, _Nucl. Phys. B_ 601 (2001) 248 [hep-ph/0008287].
  * (40) T. Gehrmann and E. Remiddi, _Two loop master integrals for $\gamma^{*}\to$ 3 jets: The Nonplanar topologies_, _Nucl. Phys. B_ 601 (2001) 287 [hep-ph/0101124].
  * (41) N. Arkani-Hamed, J.L. Bourjaily, F. Cachazo and J. Trnka, _Local Integrals for Planar Scattering Amplitudes_ , _JHEP_ 06 (2012) 125 [1012.6032].
  * (42) A.V. Kotikov, _The Property of maximal transcendentality in the N=4 Supersymmetric Yang-Mills_ , pp. 150–174, 5, 2010 [1005.5029].
  * (43) A.V. Kotikov and L.N. Lipatov, _On the highest transcendentality in N=4 SUSY_ , _Nucl. Phys. B_ 769 (2007) 217 [hep-th/0611204].
  * (44) A.V. Kotikov and L.N. Lipatov, _DGLAP and BFKL evolution equations in the N=4 supersymmetric gauge theory_ , in _35th Annual Winter School on Nuclear and Particle Physics_ , 12, 2001 [hep-ph/0112346].
  * (45) J.M. Henn, _Multiloop integrals in dimensional regularization made simple_ , _Phys. Rev. Lett._ 110 (2013) 251601 [1304.1806].
  * (46) A.B. Goncharov, _Multiple polylogarithms, cyclotomy and modular complexes_ , _Math. Res. Lett._ 5 (1998) 497 [1105.2076].
  * (47) E. Remiddi and J.A.M. Vermaseren, _Harmonic polylogarithms_ , _Int. J. Mod. Phys. A_ 15 (2000) 725 [hep-ph/9905237].
  * (48) J. Vollinga and S. Weinzierl, _Numerical evaluation of multiple polylogarithms_ , _Comput. Phys. Commun._ 167 (2005) 177 [hep-ph/0410259].
  * (49) S. Di Vita, P. Mastrolia, U. Schubert and V. Yundin, _Three-loop master integrals for ladder-box diagrams with one massive leg_ , _JHEP_ 09 (2014) 148 [1408.3107].
  * (50) R.K. Ellis, I. Hinchliffe, M. Soldate and J.J. van der Bij, _Higgs Decay to $\tau^{+}\tau^{-}$: A Possible Signature of Intermediate Mass Higgs Bosons at the SSC_, _Nucl. Phys. B_ 297 (1988) 221.
  * (51) U. Baur and E. Glover, _Higgs boson production at large transverse momentum in hadronic collisions_ , _Nucl. Phys. B_ 339 (1990) 38.
  * (52) B.A. Kniehl and M. Spira, _Low-energy theorems in Higgs physics_ , _Z. Phys. C_ 69 (1995) 77 [hep-ph/9505225].
  * (53) K.G. Chetyrkin, B.A. Kniehl and M. Steinhauser, _Decoupling relations to O $(\alpha_{s}^{3})$ and their connection to low-energy theorems_, _Nucl. Phys. B_ 510 (1998) 61 [hep-ph/9708255].
  * (54) D. Graudenz, M. Spira and P.M. Zerwas, _QCD corrections to Higgs-boson production at proton-proton colliders_ , _Phys. Rev. Lett._ 70 (1993) 1372.
  * (55) M. Spira, _QCD effects in Higgs physics_ , _Fortsch. Phys._ 46 (1998) 203 [hep-ph/9705337].
  * (56) T. Gehrmann and E. Remiddi, _Analytic continuation of massless two loop four point functions_ , _Nucl. Phys. B_ 640 (2002) 379 [hep-ph/0207020].
  * (57) P. Nogueira, _Automatic Feynman graph generation_ , _J. Comput. Phys._ 105 (1993) 279.
  * (58) C. Studerus, _Reduze-Feynman Integral Reduction in C++_ , _Comput. Phys. Commun._ 181 (2010) 1293 [0912.2546].
  * (59) A. von Manteuffel and C. Studerus, _Reduze 2 - Distributed Feynman Integral Reduction_ , 1201.4330.
  * (60) J.A.M. Vermaseren, _New features of FORM_ , math-ph/0010025.
  * (61) P. Maierhöfer, J. Usovitsch and P. Uwer, _Kira—A Feynman integral reduction program_ , _Comput. Phys. Commun._ 230 (2018) 99 [1705.05610].
  * (62) J. Klappert, F. Lange, P. Maierhöfer and J. Usovitsch, _Integral reduction with Kira 2.0 and finite field methods_ , _Comput. Phys. Commun._ 266 (2021) 108024 [2008.06494].
  * (63) C. Meyer, _Algorithmic transformation of multi-loop master integrals to a canonical basis with CANONICA_ , _Comput. Phys. Commun._ 222 (2018) 295 [1705.06252].
  * (64) O. Gituliar and V. Magerya, _Fuchsia: a tool for reducing differential equations for Feynman master integrals to epsilon form_ , _Comput. Phys. Commun._ 219 (2017) 329 [1701.04269].
  * (65) J. Henn, B. Mistlberger, V.A. Smirnov and P. Wasser, _Constructing d-log integrands and computing master integrals for three-loop four-particle scattering_ , _JHEP_ 04 (2020) 167 [2002.09492].
  * (66) A. Primo and L. Tancredi, _On the maximal cut of Feynman integrals and the solution of their differential equations_ , _Nucl. Phys. B_ 916 (2017) 94 [1610.08397].
  * (67) P.A. Baikov, _Explicit solutions of the multiloop integral recurrence relations and its application_ , _Nucl. Instrum. Meth. A_ 389 (1997) 347 [hep-ph/9611449].
  * (68) T. Gehrmann, T. Huber and D. Maitre, _Two-loop quark and gluon form-factors in dimensional regularisation_ , _Phys. Lett. B_ 622 (2005) 295 [hep-ph/0507061].
  * (69) C. Duhr and F. Dulat, _PolyLogTools — polylogs for the masses_ , _JHEP_ 08 (2019) 135 [1904.07279].
  * (70) H. Ferguson and D. Bailey, _A polynomial time, numerically stable integer relation algorithm_ , _RNR Technical Report_ (1992) .
  * (71) P. Wasser, _Scattering Amplitudes and Logarithmic Differential Forms_ , Ph.D. thesis, Mainz U., Johannes Gutenberg-Universität Mainz, 2022.  10.25358/openscience-6801.
  * (72) B. Jantzen, A.V. Smirnov and V.A. Smirnov, _Expansion by regions: revealing potential and Glauber regions automatically_ , _Eur. Phys. J. C_ 72 (2012) 2139 [1206.0546].
  * (73) A.V. Smirnov, _FIESTA4: Optimized Feynman integral calculations with GPU support_ , _Comput. Phys. Commun._ 204 (2016) 189 [1511.03614].
  * (74) S. Borowka, G. Heinrich, S. Jahn, S.P. Jones, M. Kerner, J. Schlenk and T. Zirke, _pySecDec: a toolbox for the numerical evaluation of multi-scale integrals_ , _Comput. Phys. Commun._ 222 (2018) 313 [1703.09692].
  * (75) T. Kinoshita, _Mass singularities of Feynman amplitudes_ , _J. Math. Phys._ 3 (1962) 650.
  * (76) T.D. Lee and M. Nauenberg, _Degenerate Systems and Mass Singularities_ , _Phys. Rev._ 133 (1964) B1549.
  * (77) S. Catani and M. Seymour, _A general algorithm for calculating jet cross sections in NLO QCD_ , _Nuclear Physics B_ 485 (1997) 291.
  * (78) S. Catani, _The singular behaviour of QCD amplitudes at two-loop order_ , _Physics Letters B_ 427 (1998) 161.
  * (79) G.F. Sterman and M.E. Tejeda-Yeomans, _Multiloop amplitudes and resummation_ , _Phys. Lett. B_ 552 (2003) 48 [hep-ph/0210130].
  * (80) S.M. Aybat, L.J. Dixon and G.F. Sterman, _The Two-loop anomalous dimension matrix for soft gluon exchange_ , _Phys. Rev. Lett._ 97 (2006) 072001 [hep-ph/0606254].
  * (81) S.M. Aybat, L.J. Dixon and G.F. Sterman, _The Two-loop soft anomalous dimension matrix and resummation at next-to-next-to leading pole_ , _Phys. Rev. D_ 74 (2006) 074004 [hep-ph/0607309].
  * (82) T. Becher and M. Neubert, _Infrared singularities of scattering amplitudes in perturbative QCD_ , _Phys. Rev. Lett._ 102 (2009) 162001 [0901.0722].
  * (83) T. Becher and M. Neubert, _On the Structure of Infrared Singularities of Gauge-Theory Amplitudes_ , _JHEP_ 06 (2009) 081 [0903.1126].
  * (84) L.J. Dixon, _Matter Dependence of the Three-Loop Soft Anomalous Dimension Matrix_ , _Phys. Rev. D_ 79 (2009) 091501 [0901.3414].
  * (85) E. Gardi and L. Magnea, _Factorization constraints for soft anomalous dimensions in QCD scattering amplitudes_ , _JHEP_ 03 (2009) 079 [0901.1091].
  * (86) E. Gardi and L. Magnea, _Infrared singularities in QCD amplitudes_ , _Nuovo Cim. C_ 32N5-6 (2009) 137 [0908.3273].
  * (87) O. Almelid, C. Duhr and E. Gardi, _Three-loop corrections to the soft anomalous dimension in multileg scattering_ , _Phys. Rev. Lett._ 117 (2016) 172002 [1507.00047].
  * (88) F. Buccioni, J.-N. Lang, J.M. Lindert, P. Maierhöfer, S. Pozzorini, H. Zhang and M.F. Zoller, _OpenLoops 2_ , _Eur. Phys. J. C_ 79 (2019) 866 [1907.13071].
  * (89) G. Korchemsky and A. Radyushkin, _Renormalization of the wilson loops beyond the leading order_ , _Nuclear Physics B_ 283 (1987) 342.
  * (90) W.L. van Neerven, _Dimensional regularization of mass and infrared singularities in two’loop on-shell vertex functions_ , _Nucl. Phys. B_ 268 (1986) 453.
  * (91) T. Matsuura, S.C. van der Marck and W.L. van Neerven, _The calculation of the second order soft and virtual contributions to the Drell-Yan cross-cection_ , _Nucl. Phys. B_ 319 (1989) 570.
  * (92) R.V. Harlander, _Virtual corrections to g g $to$ H to two loops in the heavy top limit_, _Phys. Lett. B_ 492 (2000) 74 [hep-ph/0007289].