††institutetext: Department of Physics,
University at Buffalo, The State University of New York, Buffalo 14260, USA

# Pseudoscalar Higgs plus jet production at Next-to-Next-to-Leading Order in
QCD

Youngjin Kim  and Ciaran Williams<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We present a calculation of pseudoscalar Higgs production in association with
a jet at Next-to-Next-to Leading Order (NNLO) accuracy in QCD. We work in an
effective field theory in which $m_{t}\rightarrow\infty$ resulting in
effective operators which couple the pseudoscalar to gluons and (massless)
quarks. We have calculated all of the relevant amplitudes for the two-loop,
one-loop and tree-level contributions. As a cross-check of our calculation we
have re-calculated all of the scalar Higgs plus parton amplitudes and perform
a detailed comparison to the literature. In order to regulate the infra-red
singularities present at this order we employ the $N-$jettiness slicing
method. In addition to a detailed validation of our calculation at this order
we investigate LHC phenomenology for a selection of pseudoscalar Higgs masses.
Our results are implemented into the parton-level Monte Carlo code MCFM.

##  1 Introduction

The continued study of the Higgs boson, discovered over a decade ago
ATLAS:2012yve ; CMS:2012qbp , remains one of the top priorities for high
energy physics. Over the next twenty years the full data set of the HL-LHC
will allow for unprecedented glimpse at the Higgs boson’s coupling to matter,
other bosons and ultimately itself. In all instances the hope is that
deviations from the predictions of the Standard Model (SM) may be observed
suggestive of new interactions which may elucidate the origins of Electroweak
symmetry breaking. In order to achieve this aim, precision predictions for the
Higgs boson will be essential. For the inclusive Higgs situation calculations
have achieved a high level of precision with N3LO predictions completed
Anastasiou:2015vya ; Mistlberger:2018etf ; Chen:2021isd . At the LHC an
observable of great interest is the transverse momentum of the recoiling Higgs
boson. The calculation of the production of a Higgs boson at non-zero $p_{T}$
has been a successful area of research over the last ten years, with several
groups pushing the accuracy to NNLO in QCD in an effective field theory in
which the top quark mass goes to infinity Boughezal:2013uia ; Chen:2014gva ;
Boughezal:2015dra ; Boughezal:2015aha ; Caola:2015wna ; Chen:2016zka ;
Campbell:2019gmd . The Higgs effective theory with $m_{t}\to\infty$ provides,
for the most part, an excellent description of the full theory. However,
results from the EFT breakdown in certain regions of phase space (typically
high-$p_{T}$) where a sensitivity to the top mass returns. The exact
calculation the Higgs production including the top-quark mass dependence is
known up to NLO in QCD Spira:1995rr . The EFT prediction can be improved by
including sub-leading top-quark mass corrections systematically as an
expansion in the inverse top-quark mass Pak:2009dg ; Harlander:2009mq .

Although the SM is a beautiful quantum field theory, and the discovery of the
Higgs provided the last missing piece, there are well motivated reasons to
believe it is not the fundamental theory of nature. Leading concerns include,
the lack of a dark matter candidate, the ability to explain neutrino masses,
or the amount of CP violation observed in the universe. For these reasons and
more many theories Beyond the Standard Model (BSM) have been constructed to
address and solve these problems. A common feature of a vast majority of BSM
models includes an extended Higgs sector, for instance, any type of two Higgs
doublet models or the Minimal Supersymmetric Standard Model Fayet:1974pd ;
Dimopoulos:1981zb ; Sakai:1981gr ; Inoue:1982pi ; Inoue:1983pp ; Inoue:1982ej
introduce an extra scalar doublet. In BSM theories which introduce a second
Higgs doublet, after electroweak symmetry breaking, there are five physical
scalar bosons which include two CP even Higgs states $h,H$, a CP odd
pseudoscalar $A$ and two charged Higgs bosons $H^{\pm}$. Given the SM-like
nature of (one of the) CP even Higgs bosons, extensions of the SM are
naturally constrained such that the spectrum of states is in an alignment
limit with the SM.

Multiple new scalar bosons is thus a simple, and somewhat model independent,
prediction of extended Higgs sectors, and this motivated many dedicated
experimental searches and theoretical work. In this paper we are primarily
interested in the psuedoscalar $A$ boson (focusing on the situation where the
ratio between the vacuum expectation of the Higgs doublets, $\tan{\beta}$ is
small). Given the mathematical similarities between the CP even and CP odd
Higgs bosons, precision predictions for psuedoscalar bosons tend to track
their CP-even Higgs counterparts. For the effective field theory induced by
taking $m_{t}\rightarrow\infty$ scalar Higgs production at next to leading
order Dawson:1990zj ; Djouadi:1991tka  was calculated in the early 1990s, with
the corresponding studies for the pseudoscalar Higgs production case following
shortly after Kauffman:1993nv 111The exact top mass dependence for pseudo
scalar Higgs is known up to NLO Spira:1993bb .. Later, in the early 2000s, the
Next-to-Next-to-leading order (NNLO) predictions were completed
Harlander:2002wh ; Anastasiou:2002yz ; Ravindran:2003um ; Anastasiou:2005qj ,
immediately followed by the pseudoscalar Higgs case Harlander:2002vv ;
Anastasiou:2002wq . Inclusive Higgs production has recently been computed at
N3LO Anastasiou:2015vya ; Mistlberger:2018etf ; Chen:2021isd , and some
progress has been made towards the full N3LO predictions for pseudoscalar
production Ahmed:2015qda ; Ahmed:2016otz . The corresponding studies of the
(pseudo)scalar Higgs boson produced with an associated jet show a similar
pattern. NLO predictions can be found in Refs. deFlorian:1999zd ;
Ravindran:2002dc ; Field:2002pb ; Bernreuther:2010uw . However, after a series
of recent calculations of the production of the scalar Higgs plus jet at Next-
to-Next-to leading order Boughezal:2013uia ; Chen:2014gva ; Boughezal:2015dra
; Boughezal:2015aha ; Caola:2015wna ; Chen:2016zka ; Campbell:2019gmd , no
such calculation for pseudoscalar Higgs bosons is available. Addressing this
gap is the primary motivation of our paper.

The NNLO calculations of Higgs plus jet have relied on a series of
breakthroughs over the last decade. For a long time NNLO calculations at
hadron colliders involving final state jets were beyond reach due to the
technical complexities of infra-red (IR) limits. However, over the last ten
years dramatic improvements have been made in this regard. Two popular
techniques for dealing with IR singularities include N-jettiness slicing
Gaunt:2015pea ; Boughezal:2015dva  and antenna subtraction Gehrmann-
DeRidder:2005btv ; Daleo:2006xa ; Daleo:2009yj ; Gehrmann:2011wi ;
Boughezal:2010mc ; Gehrmann-DeRidder:2012too ; Currie:2013vh  and these
methods have led to dramatic improvements in the number of NNLO predictions
available. In our calculation, we have taken advantage of the well-established
MCFM Boughezal:2016wmq ; Campbell:2019dru  framework (which uses $N$-jettiness
slicing) that has been successfully deployed for many calculations at NNLO,
the most relevant for this work being the calculation of $H+j$ at NNLO
Campbell:2019gmd  and bottom induced Higgs plus jet production Mondini:2021nck
.

Our paper proceeds as follows. In Section 2 we give a general overview of the
calculation, while a detailed discussion of our two-loop hard function
calculation is presented in Section 3, and the calculation of $A+2j$ at NLO is
discussed in Section 4. We discuss the results of our Monte Carlo
implementation of MCFM in Section 5 including the phenomenology for this
process at the 13 TeV LHC. Then, lastly, we draw our conclusions in Section 6.
A series of useful formulae, and a detailed comparison with the $H+j$
calculation, as well as a description of ancillary files make up a series of
Appendices.

##  2 Overview of the calculation

In this paper we consider the production of a pseudoscalar Higgs boson $(A)$
and an additional jet up to NNLO in QCD.

(a) Virtual-Virtual

(b) Real-Virtual

(c) Real-Real

Figure 1: Representative Feynman diagrams for the production of an $A$ boson
with an additional jet at NNLO.

Representative Feynman diagrams for this process at NNLO are shown in Fig. 1
and Fig. 2.

(a) $f=ggg$

(b) $f=q\bar{q}g$

Figure 2: Representative $O_{J}$ Feynman diagrams which contribute at NNLO.

In each diagram the $A$ boson couples to partons through an effective
Lagrangian. The effective Lagrangian describing the coupling of the $A$ boson
to gluons (and massless quarks) arises from integrating out the top quark and
can be written as Chetyrkin:1998mw 222Note we factored $\frac{1}{8}$ and
$\frac{1}{2}$ in the each Wilson coefficients respectively. See eq.(2.2) in
Banerjee:2017faz .,

$\displaystyle\mathcal{L}^{A}_{\rm
eff}=-A\Big{[}{C}_{G}O_{G}(x)+{C}_{J}O_{J}(x)\Big{]},$ (1)

where the operators are defined as

$O_{G}(x)=G^{\mu\nu}_{a}\tilde{G}_{a,\mu\nu}\equiv\epsilon_{\mu\nu\rho\sigma}G^{\mu\nu}_{a}G^{\rho\sigma}_{a}\,,\qquad
O_{J}(x)=\partial_{\mu}\left(\bar{\psi}\gamma^{\mu}\gamma_{5}\psi\right)\,.$
(2)

The Wilson coefficients $C_{G}$ and $C_{J}$ are obtained by integrating out
the top quark loops and are defined as follows:

$\displaystyle
C_{G}=-\frac{\alpha_{s}}{2\pi}\frac{1}{v}\left(\frac{1}{8}\right),\
\text{when}\ \tan{\beta}\sim 1$ (3) $\displaystyle
C_{J}=-\left[\left(\frac{\alpha_{s}}{2\pi}\right)\frac{C_{F}}{4}\left(\frac{3}{2}-3\ln\frac{\mu_{R}^{2}}{m_{t}^{2}}\right)+\left(\frac{\alpha_{s}}{2\pi}\right)^{2}C_{J}^{(2)}+\cdots\right]C_{G}\,.$
(4)

The Adler-Bardeen theorem Adler:1969gk  constraints the higher order
corrections to $C_{G}$, in our calculation we need $C_{J}$ expanded to
$O(\alpha_{s}^{2})$. The $C_{J}$ operator is one order higher in $\alpha_{s}$
than $C_{G}$ and as such the contributions from this operator are simpler and
correspond to one-loop boxes. The full complexities of the NNLO calculation
are therefore felt in the $C_{G}$ piece of the calculation. It is interesting
to note that the tree-level contribution from the $C_{J}$ pieces (i.e. the
right hand diagram in Fig. 2 without the loop gluon) enter at
$\mathcal{O}(\epsilon)$. Therefore as $\epsilon\rightarrow 0$ there is no
contribution from $C_{J}$ at $\mathcal{O}(\alpha_{s}^{4})$, but these pieces
do play a role in correctly performing the renormalization at NNLO.

The topologies presented in Fig. 1 are examples of double-virtual (2-loop),
real-virtual (1-loop) and real-real (tree-level) corrections which appear at
$\mathcal{O}(\alpha_{s}^{2})$ in perturbation theory. Each topology
individually contains copious infra-red (IR) singularities which cancel only
upon combinations of the sub-topologies into IR-safe observables. Over the
last decade significant progress has been achieved on regulating IR
divergences. In our calculation we will employ the $N$-jettiness slicing
method Gaunt:2015pea ; Boughezal:2015dva . This proceeds by defining the
1-jettiness observable as follows:

$\displaystyle\tau_{1}=\sum_{j=1,m}\min_{i}\left\\{\frac{2q_{i}\cdot
p_{j}}{Q_{i}}\right\\}\,,$ (5)

where the index $j$ runs over the momenta of all partons $p_{j}$ while $q_{i}$
represents a momenta of the jet-clustered phase space. For small values of
$\tau_{1}$ (defined by $\tau^{\rm{cut}}_{1}$) a factorization theorem exists
from Soft-Collinear-Effective-Theory (SCET) Stewart:2010tn ; Stewart:2009yx :

$\displaystyle\sigma(\tau_{1}\leq\tau^{\rm{cut}}_{1})=\int_{0}^{\tau^{\rm{cut}}_{1}}d\tau_{1}\;\left(\mathcal{S}\otimes\mathcal{J}\otimes\prod_{a=1,2}\mathcal{B}_{a}\otimes\mathcal{H}\right)+\mathcal{F}(\tau^{\rm{cut}}_{1}),$
(6)

The above equation is valid up to power corrections (denoted by the
$\mathcal{F}(\tau^{\rm{cut}}_{1})$ term), which vanish in the limit
$\tau^{\rm{cut}}_{1}\rightarrow 0$. At NLO the leading power corrections are
well described by the form $\tau^{\rm{cut}}_{1}\log(\tau^{\rm{cut}}_{1}/Q)$,
and at NNLO the leading power corrections have the form
$\tau^{\rm{cut}}_{1}\log^{3}(\tau^{\rm{cut}}_{1}/Q_{i})$ (where in both cases
$Q_{i}$ is a hard scale associated with the process). The factorization
theorem contains the universal soft ($\mathcal{S}$), jet ($\mathcal{J}$), and
beam ($\mathcal{B}$) functions, for which calculations accurate to
$\mathcal{O}(\alpha_{s}^{2})$ needed for our calculation can be found in Refs.
Campbell:2017hsw ; Boughezal:2015eha ; Becher:2010pd ; Becher:2006qw ;
Gaunt:2014xga ; Gaunt:2014cfa . Specifically we use the implementation in MCFM
described in Refs. Boughezal:2016wmq ; Campbell:2019dru . The process-specific
hard function $\mathcal{H}$ must be calculated on a case-by-case basis and is
obtained from the double-virtual type topologies in Fig. 1.In our calculation
we can therefore recycle much of the results presented for $H+j$ at NNLO
Campbell:2019gmd , which given the similarity to our final state of interest,
differs only by the hard function. We describe the calculation of the hard
function in detail at $\mathcal{O}(\alpha_{s}^{2})$ in section 3. Finally, we
note that there are various options regarding the choice of the hard scale
$Q_{i}$, appearing in the definitions above. The most common choice of $Q_{i}$
is $Q_{i}=2E_{i}$ which is the so-called geometric measure Jouttenus:2011wh ;
Jouttenus:2013hs . However in Ref. Campbell:2019gmd , for higher order
corrections to Higgs plus jet it was shown that a boosted frame, in which the
final-state system is at rest, resulted in a smaller impact of power
corrections. We also have confirmed this boosted choice shows better
asymptotic behavior compared to the geometric measure, and therefore will
present results using the boosted definition in this paper.

The SCET-factorization theorem describes the cross section in region below the
cut, leaving the above-cut region $\tau_{1}>\tau^{\rm{cut}}_{1}$ to be
determined. The primary advantage of the slicing method resides in the fact
that in the above-cut region there is at least one resolved emission, and
therefore the phase space contains only NLO IR-singularities. Therefore in
order to describe this region, a NLO calculation of $A+2j$ is all that is
required. We discuss this calculation in detail in section 4.

##  3 Hard function for $A\rightarrow ggg$ and $A\rightarrow q\overline{q}g$
at NNLO

###  3.1 Overview

In this section we describe the calculation of the hard function for the
process $A\rightarrow ggg$ and $A\rightarrow q\overline{q}g$ at NNLO accuracy
($\mathcal{O}(\alpha_{s}^{5})$) in the $m_{t}\rightarrow\infty$ limit. We
begin by considering the decay of an $A$ boson to partons. At LO these
reactions can be denoted as follows:

$\displaystyle A\to g(p_{1})\,g(p_{2})\,g(p_{3})\quad{\rm{and}}\quad\ A\to
q(p_{1})\,\bar{q}(p_{2})\,g(p_{3})\,.$

We define Mandelstam invariants for this process as

$\displaystyle s=(p_{1}+p_{2})^{2}>0\,,\hskip
42.67912ptt=(p_{1}+p_{3})^{2}>0\,,\hskip 42.67912ptu=(p_{2}+p_{3})^{2}>0\,,$

and from momentum conservation these satisfy $s+t+u=m_{A}^{2}$, where $m_{A}$
denotes the mass of the pseudoscalar boson. We also introduce the
dimensionless quantities

$\displaystyle x=\frac{s}{m_{A}^{2}}\,,\hskip
28.45274pty=\frac{t}{m_{A}^{2}}\,,\hskip 28.45274ptz=\frac{u}{m_{A}^{2}}\,,$
(7)

which satisfy $0<x<1$, $0<y<1$, $0<z<1$, and $x+y+z=1$ (for decay kinematics).
In order to obtain predictions for the collider production of an $A$ boson
with an additional jet, the decay amplitudes must be crossed to generate the
relevant partons in the initial state, this alters the allowed regions in
$(x,y,z)$ space, which we will discuss shortly. We note that the calculation
for the decay kinematics $A\rightarrow$ 3 partons has been considered in the
literature Banerjee:2017faz . Although to allow for an efficient mechanism of
generating the crossed process we have re-calculated the 2-loop amplitudes
independently. However, where possible we try to maintain the notation from
Banerjee:2017faz  to facilitate an easy comparison, whenever we deviate from
the notation in the literature we will mention it333For instance, we note that
$t\leftrightarrow u$ and $y\leftrightarrow z$ in Banerjee:2017faz ..

###  3.2 Calculation

We now discuss the calculation of the amplitudes up to two-loop order. The
effective Lagrangian contains both $\epsilon_{\mu\nu\rho\sigma}$ and
$\gamma_{5}$ which are inherently four-dimensional objects. Care must
therefore be taken when calculating using the dimensional regulator
$d=4-2\epsilon$. We proceed as follows, $\gamma_{5}$ is defined as

$\displaystyle\gamma_{5}=\frac{i}{4!}\varepsilon_{\mu\nu\rho\sigma}\gamma^{\mu}\gamma^{\nu}\gamma^{\rho}\gamma^{\sigma}\,,$
(8)
$\displaystyle\varepsilon_{\mu_{1}\nu_{1}\rho_{1}\sigma_{1}}\,\varepsilon^{\mu_{2}\nu_{2}\rho_{2}\sigma_{2}}=\large{\left|\begin{array}[]{cccc}\delta_{\mu_{1}}^{\mu_{2}}&\delta_{\mu_{1}}^{\nu_{2}}&\delta_{\mu_{1}}^{\rho_{2}}&\delta_{\mu_{1}}^{\sigma_{2}}\\\
\delta_{\nu_{1}}^{\mu_{2}}&\delta_{\nu_{1}}^{\nu_{2}}&\delta_{\nu_{1}}^{\rho_{2}}&\delta_{\nu_{1}}^{\sigma_{2}}\\\
\delta_{\rho_{1}}^{\mu_{2}}&\delta_{\rho_{1}}^{\nu_{2}}&\delta_{\rho_{1}}^{\rho_{2}}&\delta_{\rho_{1}}^{\sigma_{2}}\\\
\delta_{\sigma_{1}}^{\mu_{2}}&\delta_{\sigma_{1}}^{\nu_{2}}&\delta_{\sigma_{1}}^{\rho_{2}}&\delta_{\sigma_{1}}^{\sigma_{2}}\end{array}\right|}$
(13)

where all the Lorentz indices are $d$-dimensional Larin:1993tq . This is
commonly referred to as the Larin prescription. While being straightforward to
implement in $d$ dimensions, the above definition fails to satisfy the Ward
identities, and therefore an additional, finite, renormalization of
$\gamma_{5}$ is required in order to restore the Ward Identity.

Aside from the issues with $\gamma_{5}$ there are additional complications
arising from operator mixing. The essential details of the UV renormalization
procedure are well summarized in Refs. Ahmed:2015qpa ; Banerjee:2017faz , here
for brevity we provide a brief overview and we refer the reader to the
literature for further details. We renormalize the bare strong coupling
constant by performing the replacement

$\hat{\alpha}_{s}\to\alpha_{s}\,S_{\epsilon}\,Z_{\alpha}\,,$ (14)

with $S_{\epsilon}=\frac{\exp{(\epsilon\gamma_{E})}}{(4\pi)^{\epsilon}}$,
$\alpha_{s}\equiv\alpha_{s}(\mu_{R})$ at the renormalization scale $\mu_{R}$
which keeps the coupling constant dimensionless in $d=4-2\epsilon$. The
renormalization factors are given by

$\displaystyle Z_{\alpha}$
$\displaystyle=1+\left(\frac{\alpha_{s}}{2\pi}\right)r_{1}+\left(\frac{\alpha_{s}}{2\pi}\right)^{2}r_{2}+\mathcal{O}(\alpha_{s}^{3})\,,$
(15)

where explicit formula for the terms $r_{1}$, $r_{2}$ are defined in Appendix
A. After coupling renormalization the sub-amplitudes for the two operators
$O_{\Lambda=G}$ and $O_{\Lambda=J}$ for each process $f=ggg,\ q\bar{q}g$ are
written as,

$\displaystyle\mathcal{M}_{f}^{\Sigma,(0)}$
$\displaystyle=\hat{\mathcal{M}}_{f}^{\Sigma,(0)}\,,$ (16)
$\displaystyle\mathcal{M}_{f}^{\Sigma,(1)}$
$\displaystyle=S_{\epsilon}\hat{\mathcal{M}}_{f}^{\Sigma,(1)}+\frac{r_{1}}{2}\hat{\mathcal{M}}_{f}^{\Sigma,(0)}\,,$
(17) $\displaystyle\mathcal{M}_{f}^{\Sigma,(2)}$
$\displaystyle=S_{\epsilon}^{2}\hat{\mathcal{M}}_{f}^{\Sigma,(2)}+\frac{3r_{1}}{2}S_{\epsilon}\hat{\mathcal{M}}_{f}^{\Sigma,(1)}$
$\displaystyle\quad+\left(\frac{r_{2}}{2}-\frac{r_{1}^{2}}{8}\right)\hat{\mathcal{M}}_{f}^{\Sigma,(0)}\,,$
(18)

where $\Sigma=G,J$. Next, applying the renormalization of $\gamma_{5}$ and
operator mixing described previously we define

$\displaystyle\mathcal{M}_{f}^{\Lambda}=\sum_{\Sigma=G,J}Z_{\Lambda\Sigma}\left(\mathcal{M}_{f}^{\Sigma,(0)}+\left(\frac{\alpha_{s}}{2\pi}\right)\mathcal{M}_{f}^{\Sigma,(1)}+\left(\frac{\alpha_{s}}{2\pi}\right)^{2}\mathcal{M}_{f}^{\Sigma,(2)}+\mathcal{O}(\alpha_{s}^{3})\right),$
(19)

all the necessary renormalization coefficients can be found in Larin:1993tq ;
Zoller:2013ixa  and we just simply present them with our notation here,

$\displaystyle Z_{GG}$
$\displaystyle=1+\left(\frac{\alpha_{s}}{2\pi}\right)z_{GG,1}+\left(\frac{\alpha_{s}}{2\pi}\right)^{2}z_{GG,2}+\mathcal{O}(\alpha_{s}^{3})$
(20) $\displaystyle Z_{GJ}$
$\displaystyle=\left(\frac{\alpha_{s}}{2\pi}\right)z_{GJ,1}+\left(\frac{\alpha_{s}}{2\pi}\right)^{2}z_{GJ,2}+\mathcal{O}(\alpha_{s}^{3})$
(21) $\displaystyle Z_{JG}$ $\displaystyle=0$ (22) $\displaystyle Z_{JJ}$
$\displaystyle=1+\left(\frac{\alpha_{s}}{2\pi}\right)z_{JJ,1}+\left(\frac{\alpha_{s}}{2\pi}\right)^{2}z_{JJ,2}+\mathcal{O}(\alpha_{s}^{3}),$
(23)

where $z_{\Lambda,\Sigma}$ can be found in Appendix A. After applying the
prescription defined above we finally obtain UV-finite amplitudes. These
amplitudes can be written as

$\displaystyle\mathcal{A}_{f}=\left(4\pi\alpha_{s}\right)^{1/2}\sum_{\Lambda=G,J}C_{\Lambda}(\alpha_{s})\mathcal{M}_{f}^{\Lambda}.$
(24)

Such that, with (19) and (24), we obtain all the necessary ingredients to
define the squared amplitudes. For $f=ggg$,

$\displaystyle\mathcal{A}_{ggg}\mathcal{A}_{ggg}^{*}$
$\displaystyle=(4\pi\alpha_{s})\left(\frac{\alpha_{s}}{2\pi}\right)^{2}\left(C_{G}^{(1)}\right)^{2}\Bigg{[}\mathcal{M}_{ggg}^{G,(0)}\mathcal{M}_{ggg}^{G,(0)*}+\left(\frac{\alpha_{s}}{2\pi}\right)2\,\text{Re}\left(\mathcal{M}_{ggg}^{G,(1)}\mathcal{M}_{ggg}^{G,(0)*}\right)$
$\displaystyle+\left(\frac{\alpha_{s}}{2\pi}\right)^{2}\bigg{[}2\,\text{Re}\left(\mathcal{M}_{ggg}^{G,(2)}\mathcal{M}_{ggg}^{G,(0)*}\right)+\mathcal{M}_{ggg}^{G,(1)}\mathcal{M}_{ggg}^{G,(1)*}+2C_{J}^{(1)}\,\text{Re}\left(\mathcal{M}_{ggg}^{J,(1)}\mathcal{M}_{ggg}^{G,(0)*}\right)\bigg{]}\Bigg{]}$
(25)

and for $f=q\bar{q}g$,

$\displaystyle\mathcal{A}_{q\bar{q}g}\mathcal{A}_{q\bar{q}g}^{*}$
$\displaystyle=(4\pi\alpha_{s})\left(\frac{\alpha_{s}}{2\pi}\right)^{2}\left(C_{G}^{(1)}\right)^{2}\Bigg{[}\mathcal{M}_{q\bar{q}g}^{G,(0)}\mathcal{M}_{q\bar{q}g}^{G,(0)*}$
$\displaystyle+\left(\frac{\alpha_{s}}{2\pi}\right)\bigg{[}2\,\text{Re}\left(\mathcal{M}_{q\bar{q}g}^{G,(1)}\mathcal{M}_{q\bar{q}g}^{G,(0)*}\right)+2C_{J}^{(1)}\mathcal{M}_{q\bar{q}g}^{G,(0)}\mathcal{M}_{q\bar{q}g}^{J,(0)*}\bigg{]}$
$\displaystyle+\left(\frac{\alpha_{s}}{2\pi}\right)^{2}\bigg{[}2\,\text{Re}\left(\mathcal{M}_{q\bar{q}g}^{G,(2)}\mathcal{M}_{q\bar{q}g}^{G,(0)*}\right)+\mathcal{M}_{q\bar{q}g}^{G,(1)}\mathcal{M}_{q\bar{q}g}^{G,(1)*}$
$\displaystyle+2C_{J}^{(1)}\left(\mathcal{M}_{q\bar{q}g}^{G,(1)}\mathcal{M}_{q\bar{q}g}^{J,(0)*}+\mathcal{M}_{q\bar{q}g}^{G,(0)}\mathcal{M}_{q\bar{q}g}^{J,(1)*}\right)$
$\displaystyle+\left(C_{J}^{(1)}\right)^{2}\mathcal{M}_{q\bar{q}g}^{J,(0)}\mathcal{M}_{q\bar{q}g}^{J,(0)*}+2C_{J}^{(2)}\mathcal{M}_{q\bar{q}g}^{G,(0)}\mathcal{M}_{q\bar{q}g}^{J,(0)*}\bigg{]}\Bigg{]}\
,$ (26)

where

$\displaystyle C_{G}^{(1)}=-\frac{1}{v}\left(\frac{1}{8}\right)\,,\quad\quad
C_{J}^{(1)}=-\frac{C_{F}}{4}\left(\frac{3}{2}-3\ln\frac{\mu_{R}^{2}}{m_{t}^{2}}\right).$
(27)

The sum over the polarization states of the external gluons is defined by the
completeness relation

$\displaystyle\sum_{\text{pol}}\epsilon^{\mu}(p_{i})\epsilon^{\nu*}(p_{i})=-g^{\mu\nu}+\frac{p_{i}^{\mu}q^{\nu}+q^{\mu}p_{i}^{\nu}}{q\cdot
p_{i}}$ (28)

where $q$ is an auxiliary vector. In our calculation we choose
$q=p_{2},p_{3},p_{1}$ for $p_{i}=p_{1},p_{2},p_{3}$ respectively.

The workflow of our calculation proceeds in the following manner. Initially we
generated the tree-level, one-loop, and two-loop Feynman diagrams using
FeynArts Hahn:2000kx . We confirmed the results using an independent
implementation using QGRAF Nogueira:1991ex . For the FeynArts model file
generation we have used FeynRules Christensen:2008py ; Alloul:2013bka . The
Feynman rules have been implemented using FeynCalc Mertig:1990an ;
Shtabovenko:2020gxv . The resulting one-loop and two-loop amplitudes can be
written in terms of scalar loop integrals. These integrals are then further
reduced to a smaller set of Master Integrals (MIs). The reduction to MI’s was
performed using LiteRed Lee:2013mka  and KIRA Maierhofer:2017gsa . At the one-
loop level all the integrals are reduced to the bubble and box integral which
can be found in Appendix A of Ref. Garland:2001tf . At the two-loop level all
integrals can be expressed in terms of MIs which occur in the similar
$H\rightarrow$ 3 partons process, which are presented in Ref. Gehrmann:2000zt
; Gehrmann:2001ck . All of the MIs are expressed in terms of HPLs
Remiddi:1999ew  and two-dimensional HPLs (2dHPLs) Gehrmann:2000zt ;
Gehrmann:2001ck . We also use GiNaC Bauer:2000cp  and HandyG Naterop:2019xaf
to numerically evaluate the HPLs and 2dHPLs for checking purposes.

As a cross check of our calculation we calculated both the desired
$A\rightarrow 3$ parton process and $H\rightarrow 3$ partons, the later of
which has been studied extensively in the literature. We also checked the
various crossings of our amplitudes to the physical kinematic regions in the
same way. While the (re-)calculation of $H+j$ is not the primary focus of this
paper or discussion herein we provide a detailed summary of our Higgs
calculation (and comparison to the literature) in Appendix D.

###  3.3 IR subtraction and conversion to $\overline{\text{MS}}$ scheme

In order to define the hard function utilized in the SCET factorization
theorem, the amplitudes must be converted to a suitable form. In order to
obtain the hard function we remove the explicit soft and collinear divergences
from the UV-renormalized coefficients. The IR structure of one-loop and two-
loop QCD amplitudes is universally known Catani:1998bh  and can be written
using Catani’s subtraction operators $\bm{I}^{(\ell)}(\epsilon)$. The finite
coefficients $\mathcal{M}_{f}^{\Lambda,\text{fin}}$ are defined as

$\displaystyle\mathcal{M}_{f}^{\Lambda,(0),\text{fin}}$
$\displaystyle=\mathcal{M}_{f}^{\Lambda,(0)}\,,$ (29)
$\displaystyle\mathcal{M}_{f}^{\Lambda,(1),\text{fin}}$
$\displaystyle=\mathcal{M}_{f}^{\Lambda,(1)}-\bm{I}_{f}^{(1)}(\epsilon)\mathcal{M}_{f}^{\Lambda,(0)}\,,$
(30) $\displaystyle\mathcal{M}_{f}^{\Lambda,(2),\text{fin}}$
$\displaystyle=\mathcal{M}_{f}^{\Lambda,(2)}-\bm{I}_{f}^{(1)}(\epsilon)\mathcal{M}_{f}^{\Lambda,(1)}-\bm{I}_{f}^{(2)}(\epsilon)\mathcal{M}_{f}^{\Lambda,(0)}\,.$
(31)

We confirmed the IR-singularities in our amplitude cancel after this operation
as desired. The explicit expressions of the subtraction operators for
$A\rightarrow ggg,q\bar{g}g$ can be found in Appendix A. In Appendix C, we
provide full results for the tree and one-loop expressions for
$\mathcal{M}_{f}^{\Lambda,\text{fin}}$ for a quick comparison for readers. The
full results for the two-loop matrix elements can be obtained from the
attached ancillary files with the arXiv submission. The description of the
files can be found in Appendix E.

Finally, following the discussion in Ref. Becher:2013vva , we obtain the
$\overline{\text{MS}}$-renormalized coefficients
$\mathcal{M}_{f}^{\Lambda,\text{ren}}$ in the following way:

$\displaystyle\mathcal{M}_{f}^{\Lambda,(0),\text{ren}}$
$\displaystyle=\mathcal{M}_{f}^{\Lambda,(0),\text{fin}}\,,$ (32)
$\displaystyle\mathcal{M}_{f}^{\Lambda,(1),\text{ren}}$
$\displaystyle=\mathcal{M}_{f}^{\Lambda,(1),\text{fin}}+\left(\bm{I}_{f}^{(1)}(\epsilon)+\bm{Z}_{f}^{(1)}(\epsilon)\right)\mathcal{M}_{f}^{\Lambda,(0),\text{fin}}\,,$
(33) $\displaystyle\mathcal{M}_{f}^{\Lambda,(2),\text{ren}}$
$\displaystyle=\mathcal{M}_{f}^{\Lambda,(2),\text{fin}}+\left(\bm{I}_{f}^{(1)}(\epsilon)+\bm{Z}_{f}^{(1)}(\epsilon)\right)\mathcal{M}_{f}^{\Lambda,(1),\text{fin}}$
$\displaystyle+\Bigg{[}\bm{I}_{f}^{(2)}(\epsilon)+\left(\bm{I}_{f}^{(1)}(\epsilon)+\bm{Z}_{f}^{(1)}(\epsilon)\right)\bm{I}_{f}^{(1)}(\epsilon)+\bm{Z}_{f}^{(2)}(\epsilon)\Bigg{]}\mathcal{M}_{f}^{\Lambda,(0),\text{fin}}\,,$
(34)

where $\bm{Z}_{f}^{(1)}$ and $\bm{Z}_{f}^{(2)}$ for $f=ggg,q\bar{q}g$ are
presented in Appendix A.

As mentioned previously, in order to evaluate the amplitude in a kinematic
region relevant for scattering at the LHC the amplitudes must be crossed to
ensure the correct partons are in the initial state. The detailed procedure on
how to achieve this can be found in Refs. Gehrmann:2011ab ; Gehrmann:2002zr .
In this article, on the other hand, we followed coproduct method Duhr:2014woa
; Duhr:2012fh  as outlined in Ref. Mondini:2019vub . We have confirmed the
consistency of our crossed MI’s using a numerical evaluation of the MI’s in
the scattering region using the package AMFlow Liu:2022chg .

In the following we present numerical results of hard functions at each
channel using the parameters defined in Eq.(131) with $m_{A}=0.1\,\text{TeV}$:

  * •

$gg$-channel:

$\displaystyle\text{Re}\left(\mathcal{M}_{gg\to
Ag}^{G,\text{ren}}\,\mathcal{M}_{gg\to Ag}^{G*,\text{ren}}\right)$
$\displaystyle=4.85921+(38.5429)\,\alpha_{s}$
$\displaystyle+(187.3091)\,\alpha_{s}^{2}\,,$
$\displaystyle\text{Re}\left(\mathcal{M}_{gg\to
Ag}^{J,\text{ren}}\,\mathcal{M}_{gg\to Ag}^{G*,\text{ren}}+\mathcal{M}_{gg\to
Ag}^{G,\text{ren}}\,\mathcal{M}_{gg\to Ag}^{J*,\text{ren}}\right)$
$\displaystyle=0+(0)\,\alpha_{s}+(1.23085)\,\alpha_{s}^{2}\,,$

  * •

$q\bar{q}$-channel:

$\displaystyle\text{Re}\left(\mathcal{M}_{q\bar{q}\to
Ag}^{G,\text{ren}}\,\mathcal{M}_{q\bar{q}\to Ag}^{G*,\text{ren}}\right)$
$\displaystyle=0.50810+(1.38529)\,\alpha_{s}$
$\displaystyle+(3.67559)\,\alpha_{s}^{2}\,,$
$\displaystyle\text{Re}\left(\mathcal{M}_{q\bar{q}\to
Ag}^{J,\text{ren}}\,\mathcal{M}_{q\bar{q}\to
Ag}^{G*,\text{ren}}+\mathcal{M}_{q\bar{q}\to
Ag}^{G,\text{ren}}\,\mathcal{M}_{q\bar{q}\to Ag}^{J*,\text{ren}}\right)$
$\displaystyle=0+(0)\,\alpha_{s}+(0.12870)\,\alpha_{s}^{2}\,,$
$\displaystyle\text{Re}\left(\mathcal{M}_{q\bar{q}\to
Ag}^{J,\text{ren}}\,\mathcal{M}_{q\bar{q}\to Ag}^{J*,\text{ren}}\right)$
$\displaystyle=(0)+(0)\,\alpha_{s}+(0)\,\alpha_{s}^{2}\,,$

  * •

$qg$-channel:

$\displaystyle\text{Re}\left(\mathcal{M}_{qg\to
Ag}^{G,\text{ren}}\,\mathcal{M}_{qg\to Ag}^{G*,\text{ren}}\right)$
$\displaystyle=-1.96610+(-16.4889)\,\alpha_{s}$
$\displaystyle+(-77.2028)\,\alpha_{s}^{2}\,,$
$\displaystyle\text{Re}\left(\mathcal{M}_{qg\to
Ag}^{J,\text{ren}}\,\mathcal{M}_{qg\to Ag}^{G*,\text{ren}}+\mathcal{M}_{qg\to
Ag}^{G,\text{ren}}\,\mathcal{M}_{qg\to Ag}^{J*,\text{ren}}\right)$
$\displaystyle=0+(0)\,\alpha_{s}+(-0.49802)\,\alpha_{s}^{2}\,,$
$\displaystyle\text{Re}\left(\mathcal{M}_{qg\to
Ag}^{J,\text{ren}}\,\mathcal{M}_{qg\to Ag}^{J*,\text{ren}}\right)$
$\displaystyle=(0)+(0)\,\alpha_{s}+(0)\,\alpha_{s}^{2}\,.$ (37)

We note that the above results have not been normalized by the leading order
results.

###  3.4 Factorization properties of the two-loop amplitude

As this section has detailed, the calculation of the second order hard
function is quite intricate and involves several distinct stages. Through
diagram generation, reduction to MI’s, the definition of a UV finite
amplitude, extraction of IR poles, evaluation and analytically continuing the
decay amplitudes to the scattering region. Therefore it is natural to search
for methods of validating the correctness of our approach. In this section we
outline the checks which we have performed on our result, which correspond to
testing the IR singular limits of our calculation against the know
factorization properties of QCD Badger:2004uk ; Gehrmann:2011ab .

A particularly good test is the collinear limit Badger:2004uk , we consider
the two cases for $f=ggg,q\bar{q}g$ separately. For the $f=ggg$ case, we
consider a phase space point in which one of the gluons ($g=g(p_{3})$) becomes
collinear to an another gluon ($g=g(p_{1})$), as a result of which the
invariant $t$ vanishes which means $y\to 0$ while $x,z\neq 0$. For the
$f=q\bar{q}g$ case, we consider phase space point where the $q\overline{q}$
pair becomes collinear, and as such the invariant $s$ vanishes which means
$x\to 0$ while $y,z\neq 0$. The collinear limit at two loops reads:

$\displaystyle\hat{\mathcal{M}}^{G,(2)}_{f}\hat{\mathcal{M}}^{G,(0)*}_{f}$
$\displaystyle\to C^{(2)}_{f}=P_{f}^{(0)}\cdot\hat{\mathcal{M}}^{G,(2)}_{A\to
gg}\hat{\mathcal{M}}^{G,(0)*}_{A\to gg}$
$\displaystyle\qquad+P_{f}^{(1)}\cdot\hat{\mathcal{M}}^{G,(1)}_{A\to
gg}\hat{\mathcal{M}}^{G,(0)*}_{A\to gg}$
$\displaystyle\qquad+P_{f}^{(2)}\cdot\hat{\mathcal{M}}^{G,(0)}_{A\to
gg}\hat{\mathcal{M}}^{G,(0)*}_{A\to gg}\,.$ (38)

The splitting functions $P_{f}^{(\ell)}(y,z)$ and required amplitudes
$\mathcal{M}^{G,(\ell)}_{A\to gg}\mathcal{M}^{G,(0)*}_{A\to gg}$ for
$\ell=0,1,2$ are given in Appendix B. We compared our result for
$\hat{\mathcal{M}}^{G,(2)}_{f}\hat{\mathcal{M}}^{G,(0)*}_{f}$ as a series in
$\epsilon$ with $C^{(2)}_{f}$. We multiply both expressions by a factor of $x$
or $y$ to remove the leading divergence. The numerical results are displayed
in table 1. We observe excellent agreement between our result and the known
collinear limit. Additionally we also investigated the soft limit
$p_{3}\rightarrow 0$ for $f=ggg$ by following the procedure outlined in
section 8. of Ref. Badger:2004uk  finding excellent agreement, we note this
check is much more intricate for $f=q\bar{q}g$ case due to the effective
operator structure, so we did not pursue it here.

Coefficient | $y\,C^{(2)}_{ggg}$ | $y\,\hat{\mathcal{M}}^{G,(2)}_{ggg}\hat{\mathcal{M}}^{G,(0)*}_{ggg}$ | $x\,C^{(2)}_{q\bar{q}g}$ | $x\,\hat{\mathcal{M}}^{G,(2)}_{q\bar{q}g}\hat{\mathcal{M}}^{G,(0)*}_{q\bar{q}g}$
---|---|---|---|---
$\epsilon^{-4}$ | $1.20981960\cdot 10^{6}$ | $1.20981960\cdot 10^{6}$ | $4.05026555\cdot 10^{2}$ | $4.05026555\cdot 10^{2}$
$\epsilon^{-3}$ | $1.58228295\cdot 10^{7}$ | $1.58228295\cdot 10^{7}$ | $-2.59019027\cdot 10^{3}$ | $-2.59019027\cdot 10^{3}$
$\epsilon^{-2}$ | $2.36283980\cdot 10^{8}$ | $2.36283980\cdot 10^{8}$ | $-1.20976857\cdot 10^{4}$ | $-1.20976857\cdot 10^{4}$
$\epsilon^{-1}$ | $2.58965014\cdot 10^{9}$ | $2.58966527\cdot 10^{9}$ | $5.16726263\cdot 10^{4}$ | $5.16726262\cdot 10^{4}$
$\epsilon^{0}$ | $2.19247701\cdot 10^{10}$ | $2.19253448\cdot 10^{10}$ | $2.38532152\cdot 10^{5}$ | $2.38475465\cdot 10^{5}$
Table 1: Numerical comparison between our two-loop results and the known
collinear limit at phase space points $(y,z)=(10^{-10},0.23)$ and
$(x,z)=(10^{-10},0.23)$ for $f=ggg,q\bar{q}g$ respectively with
$\mu_{R}=m_{A}=125\,\text{GeV}$. For readability here we have only presented
the real part, although we note that we have found excellent agreement for the
imaginary part also.

###  3.5 Comparison with existing results

As a further check of our results, we can compare our calculation for the
squared amplitudes for $\mathcal{M}_{f}^{\Lambda,(\ell),\text{fin}}$ up to
$\ell=2$ with the existing results in the literature Banerjee:2017faz , for
decay kinematics. We have performed this check and we have found that at tree-
and one-loop level our two calculations are perfectly in agreement. The two-
loop level comparisons are non-trivial. We begin by taking account of the
following differences in notation;

$\displaystyle S_{q}^{G,(2)}\quad\text{Eq.(2.32)
in~{}\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{Banerjee:2017faz}{\@@citephrase{(}}{\@@citephrase{)}}}}$
$\displaystyle\leftrightarrow\left(\frac{1}{8}\right)^{-2}\left(\frac{1}{2}\right)^{-2}\Bigg{[}2\,\text{Re}\left(\mathcal{M}_{q\bar{q}g}^{G,(2)}\mathcal{M}_{q\bar{q}g}^{G,(0)*}\right)+\mathcal{M}_{q\bar{q}g}^{G,(1)}\mathcal{M}_{q\bar{q}g}^{G,(1)*}\Bigg{]}$
(39) $\displaystyle S_{g}^{G,(2)}\quad\text{Eq.(2.30)
in~{}\cite[cite]{\@@bibref{Authors
Phrase1YearPhrase2}{Banerjee:2017faz}{\@@citephrase{(}}{\@@citephrase{)}}}}$
$\displaystyle\leftrightarrow\left(\frac{1}{8}\right)^{-2}\left(\frac{1}{2}\right)^{-2}\Bigg{[}2\,\text{Re}\left(\mathcal{M}_{ggg}^{G,(2)}\mathcal{M}_{ggg}^{G,(0)*}\right)+\mathcal{M}_{ggg}^{G,(1)}\mathcal{M}_{ggg}^{G,(1)*}\Bigg{]}$
(40)

where the factors $(1/8)^{-2}$ and $(1/2)^{-2}$ are due to the different
definitions of the Lagrangian (1) and expansion in units of $\alpha_{s}/2\pi$
versus of $\alpha_{s}/4\pi$ as in Banerjee:2017faz . Our results should be
expanded with a factor $(-1)^{-2\epsilon}$ before comparison due to our
convention $\mu_{R}^{2}=m_{A}^{2}$ as opposed to the choice
$\mu^{2}=-m_{A}^{2}$ in the literature result. Also, note that we normalized
the amplitudes by Born color factor $C_{A}C_{F}$, $2C_{A}^{2}C_{F}$ for
$f=q\bar{q}g,ggg$ respectively as in the literature. After making these
adjustments we have found perfect agreement with the literature results for
the case $f=q\bar{q}g$, and for the case $f=ggg$ we found perfect agreement
after a couple of minor typographical errors were fixed in Ref.
Banerjee:2017faz 444We thank the authors of Banerjee:2017faz  for assistance
with this comparison..

###  3.6 Summary

In this section we have presented the calculation of the hard function
required to construct the $\tau_{1}<\tau_{1}^{\text{cut}}$ part of our NNLO
calculation. We have verified that our expressions reproduce the known IR
limits at this order, and were able to reproduce a known result in the
literature for decay kinematics. Further we used our methodology to recompute
the known processes for $H\rightarrow$ 3 partons. As a result of these checks
we are therefore confident in using our results for the phenomenology
presented in the subsequent sections of this paper.

##  4 The $\tau_{1}>\tau^{\text{cut}}_{1}$ contribution

The remaining part of our calculation corresponds to the above cut piece where
$\tau_{1}>\tau_{1}^{\text{cut}}$. This piece corresponds to the NLO
calculation of the process $A+2j$, allowing for a very soft jet requirement on
the second jet to probe the regions with small $\tau_{1}$. The one-loop
amplitude for the virtual part of the $A+2j$ process can be readily obtained
from the literature by virtue of how the Higgs plus 4 parton one-loop
amplitudes have been calculated. The Higgs plus 4 parton one-loop amplitudes
have been presented in a compact form several years ago Dixon:2004za ;
Badger:2006us ; Dixon:2009uk ; Badger:2009hw ; Badger:2009vh . One of the
reasons that the results are written so compactly is due to the amplitudes
being calculated in terms of the following complex fields;

$\displaystyle\phi=\frac{H+iA}{2},\qquad\phi^{\dagger}=\frac{H-iA}{2},$ (41)

such that the amplitude for a Higgs plus parton amplitude is the sum of $\phi$
and $\phi^{\dagger}$ amplitudes

$\displaystyle\mathcal{M}^{(1)}(H;\\{p_{k}\\})=\mathcal{M}^{(1)}(\phi;\\{p_{k}\\})+\mathcal{M}^{(1)}(\phi^{\dagger};\\{p_{k}\\}).$
(42)

On the other hand we can generate pseudoscalar amplitudes from the difference
of $\phi$ and $\phi^{\dagger}$ components,

$\displaystyle\mathcal{M}^{(1)}(A;\\{p_{k}\\})=\frac{1}{i}\left(\mathcal{M}^{(1)}(\phi;\\{p_{k}\\})-\mathcal{M}^{(1)}(\phi^{\dagger};\\{p_{k}\\})\right).$
(43)

The simplicity arises from the observation that the $\phi$ field couples only
to the self-dual part of the (trace of) the Gluon field strengths and the
$\phi^{\dagger}$ couples only to the anti-self dual components. Due to this
the helicity amplitudes for $\phi$ ($\phi^{\dagger}$) plus partons are
considerably simpler than the corresponding $H$ or $A$ equivalents
Dixon:2004za . By taking advantage of this fact, we have implemented the
virtual part of the $A+2j$ process in MCFM with slight modifications of the
scalar Higgs case. Specifically, for the NMHV helicity configurations in the
literature only the $H$ amplitude has been defined and not the individual
$\phi$ and $\phi^{\dagger}$ components. However (as a result of the helicity
arguments discussed above) for an NMHV configuration only one of the
amplitudes is complicated, the other being a simple rational function. We were
therefore able to modify the results of the literature to extract the $\phi$
and $\phi^{\dagger}$ pieces separately and construct the $A+4$ parton
amplitudes accordingly.

The double-real part of our NNLO calculation requires the tree-level
amplitudes for $A\to 5p$, which we have calculated using the BCFW recursion
relations Britto:2005fq , obtaining compact analytic expressions. All the
amplitudes present in the calculation have been checked against Madgraph
Alwall:2014hca  and GOSAM GoSam:2014iqq  at random phase space points, finding
excellent agreement.

The regularization of IR singularities present in the above-cut region has
been performed using the dipole subtraction method Catani:1996vz . In the
dipole method there are user-defined “$\alpha$ parameters” Nagy:1998bb  which
determine the amount of (non-singular) phase space integrated over by the
subtraction counter-terms. An advantage of the method is that each type of
dipole has a unique unphysical parameter which can be varied individually. If
the cancellation has proceeded correctly the individual real and virtual phase
spaces depend on this choice of $\alpha$ but the combined physical predictions
do not. Therefore following the procedure described in Ref. Campbell:2019gmd
we define the following quantities

$\epsilon^{ab}=\frac{\sigma(\alpha_{ab}=1)-\sigma(\alpha_{ab}=0.01)}{\sigma(\alpha_{ab}=1)}\,,$
(44)

where $a$ and $b$ are drawn from the set of initial and final $\\{I,F\\}$ type
dipoles. We then proceed to calculate the cross sections using the following
phase space requirements

$\displaystyle{\rm LHC},~{}\sqrt{s}=13~{}{\rm TeV},\qquad$
$\displaystyle\mu_{R}=\mu_{F}=m_{A}=125~{}{\rm GeV},$ $\displaystyle
p_{T}^{\rm jet}>20~{}{\rm GeV},\qquad$ $\displaystyle\Delta R=0.4$
$\displaystyle{\rm anti-}k_{T},\qquad$ $\displaystyle\text{no explicit cut on
rapidities}\,.$ (45)

Figure 3: The dependence of the $A+2j$ cross-section on the $\alpha$
parameters, for each partonic channel. The points represent the deviation from
$\alpha_{ab}=1$ to $10^{-2}$. The cross-sections obtained at default
parameters ($\alpha_{ab}=1$), are indicated in the plots. The horizontal lines
represent the uncertainty on a fit of the results to a constant. The plots
show excellent $\alpha$-independence within the MC uncertainty range.

Our results for the $\alpha$-(in)dependence are presented in Fig. 3, where we
have broken down cross sections into individual partonic sub-channels and
investigated each type of dipole separately. In all cases excellent agreement
is obtained with $\epsilon^{ab}=0$ within the (sub)per-mille level MC
uncertainty.

In summary we have calculated all helicity amplitudes at one-loop level for
$A+4$ partons and tree-level for $A+5$ partons and implemented the results
into MCFM. We have checked the validity of the individual amplitudes against
Madgraph, and performed intricate tests of the dipole cancellation. We are
therefore in position to put the above cut pieces together with the below cut
pieces discussed in the previous section and test our full NNLO calculation.
We do this in the next section.

##  5 Results

The cross-section expanded to NNLO accuracy can be defined as follows,

$\displaystyle\sigma_{NLO}$ $\displaystyle=\sigma_{LO}+\delta\sigma_{NLO}\,,$
(46) $\displaystyle\sigma_{NNLO}$
$\displaystyle=\sigma_{LO}+\delta\sigma_{NLO}+\delta\sigma_{NNLO}\,.$ (47)

In the following section we will present and discuss results at NLO and NNLO.
We have taken the calculations described in Sections 3 and 4, and implemented
them into the parton level Monte Carlo code MCFM, making use of the existing
implementation of the $H+j$ process Campbell:2019gmd  where possible. In order
to facilitate a comparison with the production of a scalar Higgs we introduce
the following phase space selection requirements;

$\displaystyle{\rm LHC},~{}\sqrt{s}=13~{}{\rm TeV},$
$\displaystyle\quad\mu_{R}=\mu_{F}=m_{A}=125~{}{\rm GeV}\,,$ $\displaystyle
p_{T}^{\rm jet}>30~{}{\rm GeV},$ $\displaystyle\quad{\rm anti-}k_{T}~{}{\rm
algorithm},~{}\Delta R=0.4\,,$ (48) $\displaystyle{\rm PDF~{}set:}$
$\displaystyle\quad{\tt PDF4LHC15\\_nnlo\\_30}\,.$ (49)

We note that these cuts match those used to study various NNLO calculations of
$H+j$ at NNLO in Refs. Campbell:2019gmd ; Chen:2014gva ; Chen:2016zka ;
Bizon:2018foh . We will use these cuts throughout this section.

###  5.1 Validation

We begin by validating our calculation. Since the SCET-based factorization
theorem for the below cut pieces neglects power suppressed terms, a natural
check is to ensure that the Monte Carlo code can be run in a manner within the
on-set of asymptotic behavior. We use the same definition for
$\tau_{1}^{\text{cut}}$ as Ref. Campbell:2019gmd  namely:

$\displaystyle\tau^{\text{cut}}_{1}=\epsilon\times\sqrt{m_{A}^{2}+\left(p_{T}^{j_{1}}\right)^{2}}\,.$
(50)

Typically we will draw values of $\epsilon$ from the range, $2.5\times
10^{-5}\leq\epsilon\leq 5\times 10^{-4}$. It has been shown Campbell:2019gmd
that if the $N$-jettiness variable is evaluated in the so-called boosted frame
(corresponding to the rest frame of the Higgs-jet final state system) the
resulting dependence on the unphysical parameter $\tau^{\text{cut}}_{1}$ is
softened, and asymptotic behavior is reached sooner. Therefore in this paper
we evaluate the 1-jettiness in the boosted frame.
The parametric form of the leading power corrections are well known and have
the following structure

$\displaystyle\sigma_{NLO}(\epsilon)$
$\displaystyle=\sigma_{NLO}^{0}+c_{0}\epsilon\log(\epsilon)+\cdots\,,$ (51)
$\displaystyle\delta\sigma_{NNLO}(\epsilon)$
$\displaystyle=\delta\sigma_{NNLO}^{0}+c_{0}\epsilon\log^{3}(\epsilon)+\cdots\,,$
(52)

where $\epsilon$ is defined through Eq. (50).

Figure 4: $\tau$-dependence of the total NLO cross section, $\sigma_{NLO}$.
The plot is made in the boosted frame. The blue solid line corresponds to the
fitted curve from in Eq. (51), with the blue zone representing the errors from
the fitted result. The dipole subtraction is shown as the black dashed line.

In Fig. 4 we present the results for the $A+j$ cross section at NLO using the
default set of cuts. Shown in the figure is the result obtained using dipole
subtraction (as the dashed line), the results from our calculation, and a fit
to the results as described above. Using our parametric fit one can extract
the following results

$\displaystyle\sigma^{0}_{NLO}$ $\displaystyle=31.674\pm 0.022\ \text{pb}\,,$
(53) $\displaystyle\sigma^{dipole}_{NLO}$ $\displaystyle=31.675\pm 0.031\
\text{pb}\,.$ (54)

which show the excellent agreement between the two methodologies at NLO
accuracy.

Figure 5:  $\tau$-dependence of $\delta\sigma_{NNLO}$ using our default set of
cuts. The blue solid line corresponds to the fit form in Eq. (52), with the
shaded blue color band representing the fitting errors.

We turn our attention to the NNLO coefficient in Fig. 5. Again the results
from our calculation are well modelled by the parametric fit (Eq. (52)).
Performing the fit described above allows to extract the coefficient in the
$\tau^{\text{cut}}_{1}\rightarrow 0$ limit, obtaining:

$\displaystyle\delta\sigma^{0}_{NNLO}$ $\displaystyle=6.435\pm 0.083\
\text{pb}\,.$ (55)

The general $\tau_{1}^{\rm{cut}}$ dependence of the NNLO coefficient is
unsurprisingly almost identical to that reported in the calculation of $H+j$
Campbell:2019gmd . While the asymptotic limit is harder to reach due to the
more intricate NNLO phase space, as $\epsilon\sim 2-3\times 10^{-5}$ the
results agree with the limit within the reported uncertainties. We therefore
take these values as within the asymptotic region and used $\epsilon=2.5\times
10^{-5}$ as our default value for the phenomenological studies in the next
section.

###  5.2 Phenomenology

In this section we explore the phenomenology of psuedoscalar Higgs plus jet
production at NNLO. We keep the same fiducial selection criteria described in
the previous section. NNLO predictions shown in this section are calculated
using $\epsilon=2.5\times 10^{-5}$ with the boosted $\tau$ definition. In our
predictions we take $\mu_{F}=\mu_{R}=m_{A}$ as a central scale choice. We then
perform a six-point scale variation i.e. we compute distributions taking the
extremeums from the set $(\mu_{R}/m_{A},\mu_{F}/m_{A})=(\alpha,\beta)$ where
$(\alpha,\beta)\in\\{(1,2),(1,1/2),(1/2,1),(2,1),(2,2),(1/2,1/2)\\}$ and we
exclude choices in which the scales are increased and decreased in opposite
directions.

Figure 6: Total cross sections and ratios of the cross sections across various
pseudoscalar Higgs masses, using the our default cuts. The error bands are
obtained via variations of the renormalization and the factorization factor as
described in the text.

In Fig. 6 we present the total cross section for $A+j$ as a function of the
mass of the pseudoscalar Higgs with $m_{A}$ taken in the range 0-1 TeV. As is
expected from the analogous case in which the scalar Higgs is considered, the
NNLO $K$-factor (the ratio of the NNLO to NLO cross sections) is sizable,
around a factor of 1.2. The $K$-factor is reasonably flat across the mass
range (with a gentle decrease as the mass increases). Going from NLO to NNLO
significantly improves the scale variation, for instance at $m_{A}=125$ GeV
the scale variation at NLO is 38 %, which drops to 22 % at NNLO. We also
remind the reader that the PDF choice is the same across all predictions,
larger differences are to be expected if lower order predictions matched to
the relevant PDF set are used. It should be noted that in the region
$m_{A}>2m_{t}$ the effective field theory prescription becomes unreliable due
to the resolved top quark. Significant work has been undertaken in the case of
the scalar Higgs, and it has been shown that, if the cross section is weighted
by the ratio of exact LO results including the mass, then the combined
prediction does a good job of incorporating mass effects and higher order QCD
corrections. Since a prediction resolving the top quark mass requires a (model
dependent) definition of the $At\overline{t}$ vertex, we do not include such a
re-weighting in this work (where we are primarily interested in the technical
applications of the NNLO calculation). When performing a more detailed study
of LHC phenomenology this is worth bearing in mind, but we leave this
interesting study to future work.

Figure 7: $p_{T}^{A}$ distribution of the pseudoscalar Higgs ($m_{A}=125$ GeV)
produced in association with an additional jet up to NNLO using our standard
selection criteria. The green, blue, and red colored histograms represent the
LO, NLO, and NNLO predictions respectively. The middle subplot shows
$K$-factors which are the ratios of NLO over LO and NNLO over NLO, scale
uncertainties are computed using a 7 point variation described in the text.
For each ratio the scale choice in the denominator is fixed at
$\mu_{R}=\mu_{F}=m_{A}$. The bottom panel shows the relative impact of the
pieces proportional to the $C_{J}$ operator.  Figure 8: $p_{T}^{A}$
distribution of the pseudoscalar Higgs evaluated with a mass of $m_{A}=700$
GeV. The remaining parameters are same as in Fig. 7.

Figs. 7-8 show distribution of the transverse momentum of the pseudoscalar
Higgs up to NNLO for two different mass choices $m_{A}=125$ and 700 GeV.
Including the NNLO corrections reduce the scale dependence compared the NLO
prediction as expected. For instance, in the bulk phase space ($p_{T}^{A}>30$
GeV) the typical scale variation at NLO is around 44 %, reduced to 13 % at
NNLO. We also note that the Sudakov shoulder effect Catani:1997xc  is observed
around $p_{T}^{A}$ = 30 GeV from the distribution as in the scalar Higgs case
Boughezal:2015aha . In the region below $p_{T}^{A}<p_{T}^{j}$ only the real-
virtual and real-real parts of the NNLO calculation contribute, so here the
prediction is effectively NLO with the corresponding large $K$-factor and
scale variation associated with predictions at this level. Finally in the
lowest sub-plot we isolate the contributions sensitive to the $C_{J}$
operator. We recall that these pieces arise first at NNLO and contribute via
an interference effect with the LO amplitude. As such, they reside in the
$2\rightarrow 2$ phase space. These pieces are a small part of the total NNLO
cross-section (around 0.5% across the differential distribution), and they
soften slightly at higher transverse momentum, which can be understood as an
increasing relevance of the higher multiplicity phase spaces.

Figure 9: The rapidity distribution of the pseudoscalar Higgs ($m_{A}=125$
GeV) produced in association with an additional jet up to NNLO. The green,
blue, and red colored histograms represent the LO, NLO, and NNLO respectively,
scales are varied as described in the text. The middle subplot shows
$K$-factors which are the ratios of NLO over LO and NNLO over NLO. The bottom
subplot shows the relative impact of the $C_{J}$ operator.  Figure 10: The
rapidity distribution of the pseudoscalar Higgs with $m_{A}$=700 GeV, all
other parameters are the same as Fig. 10.

Finally in Figs. 9 and 10 we present the rapidity distribution of the
pseudoscalar Higgs at the same mass choices ($125$ and 700 GeV). The NNLO
distribution is obtained in the same way as the $p^{A}_{T}$ distribution.
Whilst the effect the NLO and the NNLO corrections are approximately constant
in rapidity as in the standard Higgs case Boughezal:2015aha ; Campbell:2019gmd
, the NNLO correction significantly reduces the scale dependence as desired.
Again the $C_{J}$ pieces represent a reasonably flat correction across the
distribution, at around 0.5% of the total NNLO prediction.

##  6 Conclusions

This paper presents the NNLO calculation of a psuedoscalar Higgs boson in
association with an additional jet. The calculation of a psuedoscalar (or
scalar) Higgs boson in association with a jet at this order is rather
intricate, requiring a delicate combination of loop calculations and IR
regulation. We have performed the calculation by independently re-calculating
all components, including a recalculation of Higgs plus jet at this order to
validate our psuedoscalar Higgs calculation at every stage. We calculated the
two-loop amplitude for $A(H)\rightarrow$ 3 partons and the relevant crossings
for LHC physics. We performed several checks of our calculation against both
the literature and known IR limits of QCD amplitudes.

In order to regulate the IR divergences present in the NNLO calculation we
used the $N$-jettiness slicing method. Here one imposes a cut on the
$1$-jettiness variable to separate the calculation into two regions. When the
$1$-jettiness variable is sufficiently large the calculation has maximally one
unresolved emission and resembles a traditional NLO calculation. Below the cut
a factorization theorem from SCET allows the amplitude to be written in a
convenient form (up to suppressed power corrections).

We validated our calculation in several ways including; comparing to the
literature results where appropriate, checking the dipole cancellation in the
NLO $A+2j$ calculation and finally investigating the slicing dependence at NLO
and NNLO to ensure asymptotic behaviour. We implemented the results of
calculation into MCFM and used the code to produce phenomenological results at
13 TeV LHC. In particular we studied the transverse momentum and rapidity
distribution of the pseudoscalar focusing on two representative masses. A
natural extension of the work presented here is to include decays of the boson
and investigate the impact on constraints on current searches at the LHC.

###### Acknowledgements.

We thank John Campbell, Prasanna Dhani, Thomas Gehrmann, Roberto Mondini and
Ulrich Schubert for useful discussions. The authors are supported by the
National Science Foundation through awards NSF-PHY-2014021 and NSF-
PHY-2310363. Support provided by the Center for Computational Research at the
University at Buffalo.

##  Appendix A Formulae for renormalization and IR subtraction

The renormalization coefficients of Eqs. (15) and (20)-(23) are defined as

$\displaystyle r_{1}$ $\displaystyle=-\frac{\beta_{0}}{\epsilon}\,,$ (56)
$\displaystyle r_{2}$
$\displaystyle=\frac{\beta_{0}^{2}}{\epsilon^{2}}-\frac{\beta_{1}}{2\epsilon}\,,$
(57) $\displaystyle z_{GG,1}$
$\displaystyle=\frac{-11C_{A}}{6\epsilon}+\frac{N_{F}}{3\epsilon}\,,$ (58)
$\displaystyle z_{GG,2}$
$\displaystyle=\frac{1}{4\epsilon^{2}}\left(\frac{121C_{A}^{2}}{9}-\frac{44C_{A}N_{F}}{9}+\frac{4N_{F}^{2}}{9}\right)+\frac{1}{4\epsilon}\left(\frac{-17C_{A}^{2}}{3}+\frac{5C_{A}N_{F}}{3}+C_{F}N_{F}\right)\,,$
(59) $\displaystyle z_{GJ,1}$ $\displaystyle=\frac{6C_{F}}{\epsilon}\,,$ (60)
$\displaystyle z_{GJ,2}$
$\displaystyle=\frac{1}{2\epsilon^{2}}\left(-22C_{A}C_{F}+4C_{F}N_{F}\right)+\frac{1}{2\epsilon}\left(\frac{71C_{A}C_{F}}{3}-21C_{F}^{2}-\frac{2}{3}C_{F}N_{F}\right)\,,$
(61) $\displaystyle z_{JJ,1}$ $\displaystyle=-2C_{F}\,,$ (62) $\displaystyle
z_{JJ,2}$
$\displaystyle=\frac{1}{4\epsilon}\left(\frac{22C_{A}C_{F}}{3}+\frac{5C_{F}N_{F}}{3}\right)+\frac{1}{4}\left(22C_{F}^{2}-\frac{107C_{A}C_{F}}{9}+\frac{31C_{F}N_{F}}{18}\right)\,,$
(63)

with

$\displaystyle\beta_{0}$ $\displaystyle=\frac{11C_{A}-2N_{F}}{6}\,,$ (64)
$\displaystyle\beta_{1}$
$\displaystyle=\frac{17C_{A}^{2}-10C_{A}T_{R}N_{F}-6C_{F}T_{R}N_{F}}{6}\,,$
(65)

and $T_{R}=\frac{1}{2}$, $C_{A}=N_{c}$, $C_{F}=\frac{N_{c}^{2}-1}{2N_{c}}$.

The subtraction operators $\bm{I}_{f}^{(\ell)}(\epsilon)$ for generic QCD
processes can be found in Ref. Catani:1998bh . For completeness, we show here
the explicit expressions for the subtraction operators in CDR.

For the process $A\rightarrow ggg,q\bar{q}g$ at the one-loop order:

$\displaystyle\bm{I}_{ggg}^{(1)}(\epsilon)$
$\displaystyle=\frac{-e^{\epsilon\gamma_{E}}}{2\Gamma(1-\epsilon)}\left(-\frac{m_{H}^{2}}{\mu^{2}}\right)^{-\epsilon}\bigg{[}N_{C}\left(\frac{1}{\epsilon^{2}}+\frac{\beta_{0}}{\epsilon
N_{C}}\left(x^{-\epsilon}+y^{-\epsilon}+z^{-\epsilon}\right)\right)\bigg{]}\,,$
(66) $\displaystyle\bm{I}_{q\bar{q}g}^{(1)}(\epsilon)$
$\displaystyle=\frac{-e^{\epsilon\gamma_{E}}}{2\Gamma(1-\epsilon)}\left(-\frac{m_{H}^{2}}{\mu^{2}}\right)^{-\epsilon}\bigg{[}N_{C}\left(\frac{1}{\epsilon^{2}}+\frac{3}{4\epsilon}+\frac{\beta_{0}}{2\epsilon
N_{C}}\right)\left(y^{-\epsilon}+z^{-\epsilon}\right)-\frac{1}{N_{C}}\left(\frac{1}{\epsilon^{2}}+\frac{3}{2\epsilon}\right)x^{-\epsilon}\bigg{]}\,.$
(67)

At the two-loop order for each $f=ggg,q\bar{q}g$:

$\displaystyle\bm{I}_{f}^{(2)}(\epsilon)$
$\displaystyle=-\frac{1}{2}\bm{I}_{f}^{(1)}(\epsilon)\bm{I}_{f}^{(1)}(\epsilon)-\frac{\beta_{0}}{\epsilon}\bm{I}_{f}^{(1)}(\epsilon)+e^{-\epsilon\gamma_{E}}\frac{\Gamma(1-2\epsilon)}{\Gamma(1-\epsilon)}\left(\frac{\beta_{0}}{\epsilon}+K\right)\bm{I}_{f}^{(1)}(2\epsilon)+\bm{H}_{f}^{(2)}(\epsilon)\,,$
(68)

where

$\displaystyle K$
$\displaystyle=\left(\frac{67}{18}-\frac{\pi^{2}}{6}\right)C_{A}-\frac{10}{9}T_{R}N_{F}\,,$
(69) $\displaystyle\bm{H}_{ggg}^{(2)}(\epsilon)$
$\displaystyle=\frac{e^{\epsilon\gamma_{E}}}{4\epsilon\Gamma(1-\epsilon)}\left[3H^{(2)}_{g}\right]\,,$
(70) $\displaystyle\bm{H}_{q\bar{q}g}^{(2)}(\epsilon)$
$\displaystyle=\frac{e^{\epsilon\gamma_{E}}}{4\epsilon\Gamma(1-\epsilon)}\left[2H^{(2)}_{q}+H^{(2)}_{g}\right]\,,$
(71)

with

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

Lastly, we present the expressions for $\bm{Z}_{f}^{(1)}$ and
$\bm{Z}_{f}^{(2)}$ in Eqs. (33) and (34).
For $f=q\bar{q}g$:

$\displaystyle\bm{Z}_{q\bar{q}g}^{(1)}$
$\displaystyle=\frac{1}{\epsilon^{2}}\left(N_{C}-\frac{1}{2N_{C}}\right)+\frac{1}{\epsilon}\left[\frac{N_{C}}{6}\left(-3L(y)-3L(z)+10\right)+\frac{1}{4N_{C}}\left(2L(x)-3\right)-\frac{N_{F}}{6}\right]\,,$
(74) $\displaystyle\bm{Z}_{q\bar{q}g}^{(2)}$
$\displaystyle=\frac{1}{\epsilon^{4}}\Bigg{[}\frac{N_{C}^{2}}{2}+\frac{1}{8N_{C}^{2}}-\frac{1}{2}\Bigg{]}$
$\displaystyle+\frac{1}{\epsilon^{3}}\Bigg{[}\frac{N_{C}^{2}}{24}\bigg{(}-12L(y)-12L(z)+7\bigg{)}+\frac{1}{N_{C}^{2}}\bigg{(}3-2L(x)\bigg{)}$
$\displaystyle+\frac{1}{48}\bigg{(}24L(x)+12L(y)+12L(z)-43\bigg{)}+\frac{N_{C}N_{F}}{12}-\frac{N_{F}}{24N_{C}}\Bigg{]}$
$\displaystyle+\frac{1}{\epsilon^{2}}\Bigg{[}\frac{N_{C}^{2}}{24}\bigg{(}L(y)(6L(z)-9)+3L(y)^{2}+3L(z)^{2}-9L(z)-\pi^{2}+19\bigg{)}$
$\displaystyle+\frac{1}{32N_{C}^{2}}\bigg{(}3-2L(x)\bigg{)}^{2}+\frac{N_{C}N_{F}}{72}+\frac{5N_{F}}{72N_{C}}-\frac{N_{F}^{2}}{72}$
$\displaystyle+\frac{1}{144}\bigg{(}-18L(x)(2L(y)+2L(z)-3)+54L(y)+54L(z)+3\pi^{2}-148\bigg{)}\Bigg{]}$
$\displaystyle+\frac{1}{\epsilon}\Bigg{[}\frac{N_{C}^{2}}{1728}\bigg{(}-1728\zeta_{3}+24\left(3\pi^{2}-67\right)L(y)+24\left(3\pi^{2}-67\right)L(z)+24\pi^{2}+4771\bigg{)}$
$\displaystyle+\frac{1}{64N_{C}^{2}}\bigg{(}48\zeta_{3}-4\pi^{2}+3\bigg{)}+\frac{1}{864}\bigg{(}108\zeta_{3}+\left(804-36\pi^{2}\right)L(x)+9\pi^{2}-1042\bigg{)}$
$\displaystyle+\frac{N_{C}N_{F}}{216}\bigg{(}30L(y)+30L(z)-3\pi^{2}-110\bigg{)}+\frac{N_{F}}{432N_{C}}\bigg{(}-60L(x)+9\pi^{2}+92\bigg{)}\Bigg{]}\,,$
(75)

and for $f=ggg$:

$\displaystyle\bm{Z}_{ggg}^{(1)}$
$\displaystyle=\frac{1}{\epsilon^{2}}\left(\frac{3N_{C}}{2}\right)+\frac{1}{\epsilon}\left(\frac{N_{C}}{4}(-2L(x)-2L(y)-2L(z)+11)-\frac{N_{F}}{2}\right)\,,$
(76) $\displaystyle\bm{Z}_{ggg}^{(2)}$
$\displaystyle=\frac{1}{\epsilon^{4}}\Bigg{[}\frac{9N_{C}^{2}}{8}\Bigg{]}$
$\displaystyle+\frac{1}{\epsilon^{3}}\Bigg{[}\frac{-3N_{C}^{2}}{16}\bigg{(}4L(x)+4L(y)+4L(z)-11\bigg{)}-\frac{3}{8}N_{C}N_{F}\Bigg{]}$
$\displaystyle+\frac{1}{\epsilon^{2}}\Bigg{[}\frac{N_{C}^{2}}{96}\bigg{(}8L(x)(3L(y)+3L(z)-11)+12L(x)^{2}+8L(y)(3L(z)-11)$
$\displaystyle+12L(y)^{2}+12L(z)^{2}-88L(z)-6\pi^{2}+255\bigg{)}$
$\displaystyle+\frac{N_{C}N_{F}}{6}(L(x)+L(y)+L(z)-4)+\frac{N_{F}^{2}}{24}\Bigg{]}$
$\displaystyle+\frac{1}{\epsilon}\Bigg{[}\frac{N_{C}^{2}}{288}\bigg{(}-108\zeta_{3}+4\left(3\pi^{2}-67\right)L(x)+4\left(3\pi^{2}-67\right)L(y)$
$\displaystyle+12\pi^{2}L(z)-268L(z)-33\pi^{2}+1384\bigg{)}$
$\displaystyle+\frac{N_{C}N_{F}}{144}\bigg{(}20L(x)+20L(y)+20L(z)+3\pi^{2}-155\bigg{)}+\frac{3N_{F}}{16N_{C}}\Bigg{]}\,,$
(77)

where for brevity
$L(a)=\text{ln}\\!\left(-\frac{m_{A}^{2}}{\mu^{2}}\right)+\text{ln}\,a$.

##  Appendix B Formulae for soft and collinear limits

We list the unrenormalized $A\to gg$ ($\Lambda=G$) matrix elements that are
needed for the collinear limit checks of the two-loop $A\to ggg,q\bar{q}g$
($\Lambda=G$) amplitudes. The matrix elements in CDR in which $D=4-2\epsilon$
and $\mu_{R}^{2}=m_{A}^{2}$ read:

$\displaystyle\hat{\mathcal{M}}^{G}_{A\to gg}\hat{\mathcal{M}}^{G*}_{A\to gg}$
$\displaystyle=C_{A}C_{F}\left(C_{G}\right)^{2}m_{A}^{2}\Bigg{[}\hat{\mathcal{M}}^{G,(0)}_{A\to
gg}\hat{\mathcal{M}}^{G,(0)*}_{A\to
gg}+\left(\frac{\alpha_{s}}{2\pi}\right)(-1)^{\epsilon}S_{\epsilon}\hat{\mathcal{M}}^{G,(0)}_{A\to
gg}\hat{\mathcal{M}}^{G,(1)*}_{A\to gg}$
$\displaystyle+\left(\frac{\alpha_{s}}{2\pi}\right)^{2}(-1)^{2\epsilon}S_{\epsilon}^{2}\hat{\mathcal{M}}^{G,(0)}_{A\to
gg}\hat{\mathcal{M}}^{G,(2)*}_{A\to gg}\Bigg{]}\,,$ (78)

where

$\displaystyle\hat{\mathcal{M}}^{G,(0)}_{A\to
gg}\hat{\mathcal{M}}^{G,(0)*}_{A\to gg}$
$\displaystyle=1-3\epsilon+2\epsilon^{2}\,,$ (79)
$\displaystyle\hat{\mathcal{M}}^{G,(1)}_{A\to
gg}\hat{\mathcal{M}}^{G,(0)*}_{A\to gg}$
$\displaystyle=N_{C}\Bigg{[}\frac{-1}{\epsilon^{2}}+\frac{3}{\epsilon}+\frac{\pi^{2}}{12}+\epsilon\left(\frac{7\zeta_{3}}{3}-\frac{\pi^{2}}{4}\right)+\epsilon^{2}\left(\frac{47\pi^{4}}{1440}-7\zeta_{3}\right)\Bigg{]}+\mathcal{O}(\epsilon^{3})\,,$
(80) $\displaystyle\hat{\mathcal{M}}^{G,(2)}_{A\to
gg}\hat{\mathcal{M}}^{G,(0)*}_{A\to gg}$
$\displaystyle=\frac{N_{C}^{2}}{2\epsilon^{4}}+\frac{1}{\epsilon^{3}}\Bigg{[}\frac{N_{C}N_{F}}{12}-\frac{47N_{C}^{2}}{24}\Bigg{]}$
$\displaystyle+\frac{1}{\epsilon^{2}}\Bigg{[}-\frac{N_{C}N_{F}}{9}+N_{C}^{2}\left(-\frac{5}{9}-\frac{\pi^{2}}{24}\right)\Bigg{]}$
$\displaystyle+\frac{1}{\epsilon}\Bigg{[}N_{C}^{2}\left(-\frac{25\zeta_{3}}{12}+\frac{17\pi^{2}}{48}+\frac{811}{216}\right)+N_{C}N_{F}\left(-\frac{161}{108}-\frac{\pi^{2}}{24}\right)+\frac{3N_{F}}{4N_{C}}\Bigg{]}$
$\displaystyle+\Bigg{[}N_{C}^{2}\bigg{(}\frac{59\zeta_{3}}{9}+\frac{\pi^{2}}{36}-\frac{7\pi^{4}}{240}+\frac{5941}{324}\bigg{)}+N_{C}N_{F}\bigg{(}-\frac{19\zeta_{3}}{18}+\frac{\pi^{2}}{18}-\frac{3965}{648}\bigg{)}$
$\displaystyle+\frac{N_{F}}{N_{C}}\left(\frac{71}{24}-\zeta_{3}\right)\Bigg{]}+\mathcal{O}(\epsilon)\,,$
(81)

which also can be found from Ref. Ahmed:2015qpa  and
$S_{\epsilon}=\frac{\exp{(\epsilon\gamma_{E})}}{(4\pi)^{\epsilon}}$. The
collinear functions in Eq. (38) are for $f=ggg$,

$\displaystyle P_{ggg}^{(0)}(y,z)$
$\displaystyle=\frac{\left(z^{2}-z+1\right)^{2}}{(1-z)z}\,,$ (82)
$\displaystyle P_{ggg}^{(1)}(y,z)$
$\displaystyle=\frac{1}{(1-z)z}\text{Sp}_{ggg,1}^{(1)}(\epsilon)+\left(\frac{z^{3}}{(1-z)}+\frac{(1-z)^{3}}{z}\right)\text{Sp}_{ggg,2}^{(1)}(\epsilon)\,,$
(83) $\displaystyle P_{ggg}^{(2)}(y,z)$
$\displaystyle=\frac{1}{(1-z)z}\text{Sp}_{ggg,1}^{(2)}(\epsilon)+\left(\frac{z^{3}}{(1-z)}+\frac{(1-z)^{3}}{z}\right)\text{Sp}_{ggg,2}^{(2)}(\epsilon)\,,$
(84)

where $\text{Sp}_{ggg,n}^{(1)}(\epsilon)$ at $\mu_{R}^{2}=m_{A}^{2}$ are
defined as(see section 5. in Badger:2004uk .)

$\displaystyle\text{Sp}_{ggg,1}^{(1)}(\epsilon)$
$\displaystyle=(4\pi)^{\epsilon}c_{\Gamma}(\epsilon)(-y)^{-\epsilon}\Bigg{[}\frac{N_{C}}{2\epsilon^{2}}\Bigg{\\{}-\Gamma(1-\epsilon)\Gamma(1+\epsilon)\left(\frac{z}{1-z}\right)^{\epsilon}$
$\displaystyle+\sum_{m=1}^{\infty}2\epsilon^{2m-1}\text{Li}_{2m-1}\left(\frac{1-z}{-z}\right)\Bigg{\\}}+\frac{z(1-z)}{(1-2\epsilon)(2-2\epsilon)(3-2\epsilon)}(N_{C}-N_{C}\epsilon-
N_{F})\Bigg{]}\,,$ (85) $\displaystyle\text{Sp}_{ggg,2}^{(1)}(\epsilon)$
$\displaystyle=(4\pi)^{\epsilon}c_{\Gamma}(\epsilon)(-y)^{-\epsilon}\Bigg{[}\frac{N_{C}}{2\epsilon^{2}}\Bigg{\\{}-\Gamma(1-\epsilon)\Gamma(1+\epsilon)\left(\frac{z}{1-z}\right)^{\epsilon}$
$\displaystyle+\sum_{m=1}^{\infty}2\epsilon^{2m-1}\text{Li}_{2m-1}\left(\frac{1-z}{-z}\right)\Bigg{\\}}\Bigg{]}\,.$
(86)

For $f=q\bar{q}g$,

$\displaystyle P_{q\bar{q}g}^{(0)}(x,z)$
$\displaystyle=1-\frac{2(1-z)z}{1-\epsilon}\,,$ (87) $\displaystyle
P_{q\bar{q}g}^{(1)}(x,z)$
$\displaystyle=\left(1-\frac{2(1-z)z}{1-\epsilon}\right)\text{Sp}_{q\bar{q}g,1}^{(1)}(\epsilon)\,,$
(88) $\displaystyle P_{q\bar{q}g}^{(2)}(x,z)$
$\displaystyle=\left(1-\frac{2(1-z)z}{1-\epsilon}\right)\text{Sp}_{q\bar{q}g,1}^{(2)}(\epsilon)\,.$
(89)

where $\text{Sp}_{q\bar{q}g,1}^{(1)}(\epsilon)$ at $\mu_{R}^{2}=m_{A}^{2}$ is
defined as(see section 6. in Badger:2004uk .)

$\displaystyle\text{Sp}_{q\bar{q}g,1}^{(1)}(\epsilon)$
$\displaystyle=(4\pi)^{\epsilon}c_{\Gamma}(\epsilon)(-x)^{-\epsilon}\Bigg{[}\frac{N_{C}}{2\epsilon^{2}}\Bigg{\\{}1-\Gamma(1-\epsilon)\Gamma(1+\epsilon)\left(\frac{z}{1-z}\right)^{\epsilon}$
$\displaystyle+\sum_{m=1}^{\infty}2\epsilon^{2m-1}\text{Li}_{2m-1}\left(\frac{1-z}{-z}\right)\Bigg{\\}}$
$\displaystyle+N_{C}\left(\frac{13}{12\epsilon(1-2\epsilon)}+\frac{1}{6(1-2\epsilon)(3-2\epsilon)}\right)+\frac{1}{N_{C}}\left(\frac{1}{2\epsilon^{2}}+\frac{3}{4\epsilon(1-2\epsilon)}+\frac{1}{2(1-2\epsilon)}\right)$
$\displaystyle+N_{F}\left(\frac{-1}{3\epsilon(1-2\epsilon)}+\frac{1}{3(1-2\epsilon)(3-2\epsilon)}\right)\Bigg{]}\,.$
(90)

The two-loop splitting functions $\text{Sp}^{(2)}_{f,n}(\epsilon)$ in CDR are
defined as:

$\displaystyle\text{Sp}^{(2)}_{f,n}(\epsilon)$
$\displaystyle=\frac{1}{2}\left(\text{Sp}^{(1)}_{f,n}(\epsilon)\right)^{2}+\frac{e^{-\epsilon\gamma_{E}}\,c_{\Gamma}(\epsilon)}{c_{\Gamma}(2\epsilon)}\left(\frac{\beta_{0}}{\epsilon}+K\right)\text{Sp}^{(1)}_{f,n}(2\epsilon)+H_{f}(\epsilon)+\text{Sp}^{(2),\text{fin}}_{f,n}+\mathcal{O}(\epsilon)\,,$
(91)

where

$\displaystyle H_{ggg}(\epsilon)$
$\displaystyle=\frac{e^{-\epsilon\gamma_{E}}\,c_{\Gamma}(\epsilon)}{4\epsilon}(4\pi)^{2\epsilon}\,\left(-y\right)^{-2\epsilon}\left[z(1-z)\right]^{-2\epsilon}\left(H^{(2)}_{g}-\beta_{0}K+\beta_{1}\right)\,,$
(92) $\displaystyle H_{q\bar{q}g}(\epsilon)$
$\displaystyle=\frac{e^{-\epsilon\gamma_{E}}\,c_{\Gamma}(\epsilon)}{4\epsilon}(4\pi)^{2\epsilon}\,\left(-x\right)^{-2\epsilon}\left[z(1-z)\right]^{-2\epsilon}\left(2H^{(2)}_{q}-H^{(2)}_{g}-\beta_{0}K+\beta_{1}\right)\,,$
(93) $\displaystyle c_{\Gamma}(\epsilon)$
$\displaystyle=\frac{\Gamma(1-\epsilon)^{2}\,\Gamma(1+\epsilon)}{\Gamma(1-2\epsilon)}\,.$
(94)

with $H^{(2)}_{g},H^{(2)}_{q}$ and $K$ as defined in Appendix A. Finally, the
functions $\text{Sp}^{(2),\text{fin}}_{f,n}$ correspond to Eqs. (5.18–5.20)
and (6.8) of Ref. Badger:2004uk  respectively with the replacement $w\to z$.

##  Appendix C IR-finite squared amplitudes and hard functions for $A\to
ggg,q\bar{q}g$

In this appendix we present the IR-finite squared amplitudes for
$f=ggg,q\bar{q}g$ with $\mu_{R}^{2}=m_{A}^{2}$555Note the results presented in
Ref. Banerjee:2017faz  are with the set $\mu_{R}^{2}=-m_{A}^{2}$. . We
normalized the amplitudes by Born color factor $2C_{A}^{2}C_{F}$ for $f=ggg$
and $C_{A}C_{F}$ for $q\bar{q}g$. The results for the one-loop and two-loop
hard functions in each kinematic region are also available in the attached
auxiliary Mathematica files. To save space, we only present up to the one-loop
order for quick comparisons in this appendix. The higher order that the two-
loop and the one-loop self interference amplitudes are provided in the
attached auxiliary files.

###  C.1 $f=ggg$ case

The tree and one-loop results are given by

$\displaystyle\mathcal{M}_{ggg}^{G,(0),\text{fin}}\mathcal{M}_{ggg}^{G,(0),\text{fin}*}$
$\displaystyle=-2C_{1,g}\,,$ (95)
$\displaystyle\text{Re}\left(\mathcal{M}_{ggg}^{G,(1),\text{fin}}\mathcal{M}_{ggg}^{G,(0),\text{fin}*}\right)$
$\displaystyle=N_{C}\Bigg{[}C_{1,g}\bigg{\\{}H(0,y)H(0,z)-H(0,y)H(1,z)-\frac{11}{6}H(1-z,y)$
$\displaystyle-H(0,z)H(1-z,y)+2H(1,z)H(z,y)-H(0,1-z,y)$
$\displaystyle-H(1-z,0,y)+2H(z,1-z,y)+\frac{11}{6}H(0,y)+2H(1,0,y)$
$\displaystyle+\frac{11}{6}H(0,z)-\frac{11}{6}H(1,z)+H(0,1,z)+H(1,0,z)+\frac{\pi^{2}}{6}-4\bigg{\\}}$
$\displaystyle+\frac{C_{2,g}}{6}\Bigg{]}+N_{F}\Bigg{[}\frac{C_{1,g}}{6}\bigg{\\{}2H(1-z,y)-2H(0,y)-2H(0,z)$
$\displaystyle+2H(1,z)\bigg{\\}}-\frac{C_{2,g}}{6}\Bigg{]}\,,$ (96)
$\displaystyle\text{Im}\left(\mathcal{M}_{ggg}^{G,(1),\text{fin}}\mathcal{M}_{ggg}^{G,(0),\text{fin}*}\right)$
$\displaystyle=C_{1,g}\left(-\frac{11\pi}{2}N_{C}\right)+C_{1,g}N_{F}\pi\,,$
(97)
$\displaystyle\text{Re}\left(\mathcal{M}_{ggg}^{J,(1),\text{fin}}\mathcal{M}_{ggg}^{G,(0),\text{fin}*}\right)$
$\displaystyle=-2C_{1,g}N_{F}\,,$ (98)
$\displaystyle\text{Im}\left(\mathcal{M}_{ggg}^{J,(1),\text{fin}}\mathcal{M}_{ggg}^{G,(0),\text{fin}*}\right)$
$\displaystyle=0\,,$ (99)

and the hard functions are given by

$\displaystyle\mathcal{M}_{ggg}^{G,(0),\text{ren}}\mathcal{M}_{ggg}^{G,(0),\text{ren}*}$
$\displaystyle=-2C_{1,g}\,,$ (100)
$\displaystyle\text{Re}\left(\mathcal{M}_{ggg}^{G,(1),\text{ren}}\mathcal{M}_{ggg}^{G,(0),\text{ren}*}\right)$
$\displaystyle=N_{C}\Bigg{[}C_{1,g}\bigg{\\{}H(0,y)H(0,z)-H(0,z)H(1-z,y)$
$\displaystyle-H(0,y)H(1,z)+H(1,z)H(1-z,y)+2H(1,z)H(z,y)$
$\displaystyle-H(0,1-z,y)-H(1-z,0,y)+H(1-z,1-z,y)$
$\displaystyle+2H(z,1-z,y)+H(0,0,y)+2H(1,0,y)+H(0,0,z)$
$\displaystyle+H(0,1,z)+H(1,0,z)+H(1,1,z)\bigg{\\}}-\frac{19\pi^{2}}{12}C_{1,g}$
$\displaystyle-\frac{C_{3,g}}{6}\bigg{\\{}24y^{4}+48y^{3}z-48y^{3}+72y^{2}z^{2}-143y^{2}z+71y^{2}$
$\displaystyle+48yz^{3}-143yz^{2}+142yz-47y+24z^{4}-48z^{3}+71z^{2}$
$\displaystyle-47z+24\bigg{\\}}\Bigg{]}-\frac{N_{F}}{6}C_{2,g}\,,$ (101)
$\displaystyle\text{Im}\left(\mathcal{M}_{ggg}^{G,(1),\text{ren}}\mathcal{M}_{ggg}^{G,(0),\text{ren}*}\right)$
$\displaystyle=C_{1,g}N_{C}\bigg{[}H(1-z,y)-H(0,y)-H(0,z)+H(1,z)\bigg{]}\pi\,,$
(102)
$\displaystyle\text{Re}\left(\mathcal{M}_{ggg}^{J,(1),\text{ren}}\mathcal{M}_{ggg}^{G,(0),\text{ren}*}\right)$
$\displaystyle=-2C_{1,g}N_{F}\,,$ (103)
$\displaystyle\text{Im}\left(\mathcal{M}_{ggg}^{J,(1),\text{ren}}\mathcal{M}_{ggg}^{G,(0),\text{ren}*}\right)$
$\displaystyle=0\,,$ (104)

where the constants $C_{i,g}$, $i=1,2,3$ are defined as

$\displaystyle C_{1,g}$
$\displaystyle=m_{A}^{2}\frac{y^{4}+2y^{3}(z-1)+3y^{2}(z-1)^{2}+2y(z-1)^{3}+\left(z^{2}-z+1\right)^{2}}{yz(y+z-1)}\,,$
(105) $\displaystyle C_{2,g}$
$\displaystyle=m_{A}^{2}\frac{(1-y)(1-z)(y+z)}{yz(-y-z+1)}\,,$ (106)
$\displaystyle C_{3,g}$ $\displaystyle=\frac{m_{A}^{2}}{yz(y+z-1)}\,.$ (107)

The amplitudes at NNLO (two-loop and one-loop2) from the operator $O_{G}$ are
decompose as following color terms:

$\displaystyle\mathcal{M}_{ggg}^{G,(2(1))}\mathcal{M}_{ggg}^{G,(0(1))*}$
$\displaystyle=N_{C}^{2}\left(\mathcal{M}_{ggg}^{G,(2(1))}\mathcal{M}_{ggg}^{G,(0(1))*}\right)^{N_{C}^{2}}+\frac{N_{F}}{N_{C}}\left(\mathcal{M}_{ggg}^{G,(2(1))}\mathcal{M}_{ggg}^{G,(0(1))*}\right)^{\frac{N_{F}}{N_{C}}}$
$\displaystyle+N_{F}N_{C}\left(\mathcal{M}_{ggg}^{G,(2(1))}\mathcal{M}_{ggg}^{G,(0(1))*}\right)^{N_{F}N_{C}}+N_{F}^{2}\left(\mathcal{M}_{ggg}^{G,(2(1))}\mathcal{M}_{ggg}^{G,(0(1))*}\right)^{N_{F}^{2}}\,.$
(108)

Explicit expressions of the IR-finite amplitudes and the hard functions for
each term can be found in the attached Mathematica files.

###  C.2 $f=q\bar{q}g$ case

The tree and one-loop results are given by

$\displaystyle\mathcal{M}_{q\bar{q}g}^{G,(0),\text{fin}}\mathcal{M}_{q\bar{q}g}^{G,(0),\text{fin}*}$
$\displaystyle=-C_{1,q}\,,$ (109)
$\displaystyle\mathcal{M}_{q\bar{q}g}^{G,(0),\text{fin}}\mathcal{M}_{q\bar{q}g}^{J,(0),\text{fin}*}$
$\displaystyle=0\,,$ (110)
$\displaystyle\mathcal{M}_{q\bar{q}g}^{J,(0),\text{fin}}\mathcal{M}_{q\bar{q}g}^{J,(0),\text{fin}*}$
$\displaystyle=0\,,$ (111)
$\displaystyle\text{Re}\left(\mathcal{M}_{q\bar{q}g}^{G,(1),\text{fin}}\mathcal{M}_{q\bar{q}g}^{G,(0),\text{fin}*}\right)$
$\displaystyle=N_{C}\Bigg{[}\frac{C_{1,q}}{36}\bigg{\\{}-18H(0,y)H(1,z)-18H(0,z)H(1-z,y)$
$\displaystyle-39H(1-z,y)+36H(1,z)H(z,y)-18H(0,1-z,y)$
$\displaystyle-18H(1-z,0,y)+36H(z,1-z,y)+30H(0,y)$
$\displaystyle+18H(1,0,y)+30H(0,z)-39H(1,z)+18H(0,1,z)\bigg{\\}}$
$\displaystyle-\frac{C_{2,q}}{36}\bigg{\\{}143y^{2}-18yz+9y+143z^{2}+9z\bigg{\\}}\Bigg{]}$
$\displaystyle+\frac{1}{N_{C}}\Bigg{[}\frac{C_{1,q}}{2}\bigg{\\{}-H(0,y)H(0,z)-H(1,0,y)-H(1,0,z)\bigg{\\}}$
$\displaystyle-\frac{C_{2,q}}{12}\bigg{\\{}\pi^{2}\left(y^{2}+z^{2}\right)+21y^{2}-6yz+3y+21z^{2}+3z\bigg{\\}}\Bigg{]}$
$\displaystyle+\frac{C_{1,q}}{36}N_{F}\Bigg{[}12H(1-z,y)-3H(0,y)-3H(0,z)$
$\displaystyle+12H(1,z)+20\Bigg{]}\,,$ (112)
$\displaystyle\text{Im}\left(\mathcal{M}_{q\bar{q}g}^{G,(1),\text{fin}}\mathcal{M}_{q\bar{q}g}^{G,(0),\text{fin}*}\right)$
$\displaystyle=-\frac{11\pi}{4}C_{1,q}N_{C}+\frac{\pi}{2}C_{1,q}N_{F}\,,$
(113)
$\displaystyle\text{Re}\left(\mathcal{M}_{q\bar{q}g}^{G,(1),\text{fin}}\mathcal{M}_{q\bar{q}g}^{J,(0),\text{fin}*}+\mathcal{M}_{q\bar{q}g}^{J,(1),\text{fin}}\mathcal{M}_{q\bar{q}g}^{G,(0),\text{fin}*}\right)$
$\displaystyle=-C_{1,q}N_{F}\,,$ (114)
$\displaystyle\text{Im}\left(\mathcal{M}_{q\bar{q}g}^{G,(1),\text{fin}}\mathcal{M}_{q\bar{q}g}^{J,(0),\text{fin}*}+\mathcal{M}_{q\bar{q}g}^{J,(1),\text{fin}}\mathcal{M}_{q\bar{q}g}^{G,(0),\text{fin}*}\right)$
$\displaystyle=0\,,$ (115)

and the hard functions are given by

$\displaystyle\mathcal{M}_{q\bar{q}g}^{G,(0),\text{ren}}\mathcal{M}_{q\bar{q}g}^{G,(0),\text{ren}*}$
$\displaystyle=-C_{1,q}\,,$ (116)
$\displaystyle\mathcal{M}_{q\bar{q}g}^{G,(0),\text{ren}}\mathcal{M}_{q\bar{q}g}^{J,(0),\text{ren}*}$
$\displaystyle=0\,,$ (117)
$\displaystyle\mathcal{M}_{q\bar{q}g}^{J,(0),\text{ren}}\mathcal{M}_{q\bar{q}g}^{J,(0),\text{ren}*}$
$\displaystyle=0\,,$ (118)
$\displaystyle\text{Re}\left(\mathcal{M}_{q\bar{q}g}^{G,(1),\text{ren}}\mathcal{M}_{q\bar{q}g}^{G,(0),\text{ren}*}\right)$
$\displaystyle=N_{C}\frac{C_{1,q}}{12}\Bigg{[}-6H(0,y)H(1,z)+12H(1,z)H(z,y)$
$\displaystyle-6H(0,z)H(1-z,y)-13H(1-z,y)-6H(0,1-z,y)$
$\displaystyle-6H(1-z,0,y)+12H(z,1-z,y)+6H(0,0,y)+6H(1,0,y)$
$\displaystyle-13H(1,z)+6H(0,0,z)+6H(0,1,z)-7\pi^{2}-\frac{143}{3}+3C_{3,q}\Bigg{]}$
$\displaystyle+\frac{1}{N_{C}}\frac{C_{1,q}}{4}\Bigg{[}-2H(0,y)H(0,z)-2H(1,z)H(1-z,y)$
$\displaystyle-3H(1-z,y)-2H(1-z,1-z,y)-2H(1,0,y)-3H(1,z)$
$\displaystyle-2H(1,0,z)-2H(1,1,z)+\frac{5\pi^{2}}{6}-7+C_{3,q}\Bigg{]}$
$\displaystyle+N_{F}\frac{C_{1,q}}{9}\Bigg{[}3H(1-z,y)+3H(1,z)+5\Bigg{]}\,,$
(119)
$\displaystyle\text{Im}\left(\mathcal{M}_{q\bar{q}g}^{G,(1),\text{ren}}\mathcal{M}_{q\bar{q}g}^{G,(0),\text{ren}*}\right)$
$\displaystyle=-\frac{C_{1,q}}{12}\pi N_{C}\bigg{(}6H(0,y)+6H(0,z)+13\bigg{)}$
$\displaystyle-\frac{C_{1,q}}{4}\pi\frac{1}{N_{C}}\bigg{(}2H(1-z,y)+2H(1,z)+3\bigg{)}+\frac{C_{1,q}}{3}\pi
N_{F}\,,$ (120)
$\displaystyle\text{Re}\left(\mathcal{M}_{q\bar{q}g}^{G,(1),\text{ren}}\mathcal{M}_{q\bar{q}g}^{J,(0),\text{ren}*}+\mathcal{M}_{q\bar{q}g}^{J,(1),\text{ren}}\mathcal{M}_{q\bar{q}g}^{G,(0),\text{ren}*}\right)$
$\displaystyle=-C_{1,q}N_{F}\,,$ (121)
$\displaystyle\text{Im}\left(\mathcal{M}_{q\bar{q}g}^{G,(1),\text{ren}}\mathcal{M}_{q\bar{q}g}^{J,(0),\text{ren}*}+\mathcal{M}_{q\bar{q}g}^{J,(1),\text{ren}}\mathcal{M}_{q\bar{q}g}^{G,(0),\text{ren}*}\right)$
$\displaystyle=0\,,$ (122)

where the constants $C_{i,q}$, $i=1,2,3$ are defined as

$\displaystyle C_{1,q}$ $\displaystyle=m_{A}^{2}\frac{y^{2}+z^{2}}{y+z-1}\,,$
(123) $\displaystyle C_{2,q}$ $\displaystyle=\frac{m_{A}^{2}}{y+z-1}\,,$ (124)
$\displaystyle C_{3,q}$ $\displaystyle=\frac{2yz-y-z}{y^{2}+z^{2}}\,.$ (125)

The amplitudes at NNLO (two-loop and one-loop2) from the operator $O_{G}$ are
decompose as following color terms:

$\displaystyle\mathcal{M}_{q\bar{q}g}^{G,(2(1))}\mathcal{M}_{q\bar{q}g}^{G,(0(1))*}$
$\displaystyle=N_{C}^{2}\left(\mathcal{M}_{q\bar{q}g}^{G,(2(1))}\mathcal{M}_{q\bar{q}g}^{G,(0(1))*}\right)^{N_{C}^{2}}+\frac{1}{N_{C}^{2}}\left(\mathcal{M}_{q\bar{q}g}^{G,(2(1))}\mathcal{M}_{q\bar{q}g}^{G,(0(1))*}\right)^{N_{C}^{-2}}$
$\displaystyle+N_{C}^{0}\left(\mathcal{M}_{q\bar{q}g}^{G,(2(1))}\mathcal{M}_{q\bar{q}g}^{G,(0(1))*}\right)^{N_{C}^{0}}+\frac{N_{F}}{N_{C}}\left(\mathcal{M}_{q\bar{q}g}^{G,(2(1))}\mathcal{M}_{q\bar{q}g}^{G,(0(1))*}\right)^{\frac{N_{F}}{N_{C}}}$
$\displaystyle+N_{F}N_{C}\left(\mathcal{M}_{q\bar{q}g}^{G,(2(1))}\mathcal{M}_{q\bar{q}g}^{G,(0(1))*}\right)^{N_{F}N_{C}}+N_{F}^{2}\left(\mathcal{M}_{q\bar{q}g}^{G,(2(1))}\mathcal{M}_{q\bar{q}g}^{G,(0(1))*}\right)^{N_{F}^{2}}\,.$
(126)

Explicit expressions of the IR-finite amplitudes and the hard functions for
each term can be found in the attached Mathematica files.

##  Appendix D Recalculation of $H\to ggg,q\bar{q}g$ Amplitudes

The NNLO QCD corrections to $H+j$ have been well studied in the literature,
and serve as a useful cross check of our calculation. Therefore, we have
checked our results for $A+j$ by reproducing the relevant results for $H+j$ at
every stage of our calculation. The most intricate cross check regards the
reproduction of the two-loop amplitudes for $H\to ggg,q\bar{q}g$, which were
originally presented in Ref. Gehrmann:2011aa  and more recently expanded to
higher orders in $\epsilon$, in Ref. Gehrmann:2023etk . In Refs.
Gehrmann:2011aa ; Gehrmann:2023etk , the helicity amplitudes $\alpha,\,\beta,$
and $\gamma$ are calculated for the $H\to ggg,q\bar{q}g$ process up to two-
loop order. These are expanded to second order in $\alpha_{s}$ in terms of
$\Omega=\alpha,\,\beta,\,\gamma$ as follows

$\displaystyle\Omega=C_{H}\sqrt{4\pi\alpha_{s}}\,T_{\Omega}\left[\Omega^{(0)}+\left(\frac{\alpha_{s}}{2\pi}\right)\Omega^{(1)}+\left(\frac{\alpha_{s}}{2\pi}\right)^{2}\Omega^{(2)}+\mathcal{O}(\alpha_{s}^{3})\right]\,,$
(127)

where $C_{H}$ is the effective $Hgg$ coupling Chetyrkin:1997iv ;
Chetyrkin:1997un  in the $m_{t}\to\infty$ limit and the color factor is
$T_{\alpha}=T_{\beta}=f^{a_{1}a_{2}a_{3}}$ and
$T_{\gamma}=T^{a_{3}}_{i_{1}j_{2}}$. We have been able to reproduce the same
IR-finite results for $\alpha,\,\beta,$ and $\gamma$ as presented in appendix
A and B in ref. Gehrmann:2011aa . Since we do not use a projector-based
approach for our $A+j$ calculation a second useful check it so calculate the
squared amplitudes directly. The squared amplitudes are expanded in order of
$\alpha_{s}$ as

$\displaystyle\mathcal{A}_{H\to f}\mathcal{A}_{H\to f}^{*}$
$\displaystyle=(4\pi\alpha_{s})\left(\frac{\alpha_{s}}{2\pi}\right)^{2}\left(C_{H}\right)^{2}\mathcal{C}_{f}\Bigg{[}\mathcal{M}_{H\to
f}^{(0)}\mathcal{M}_{H\to
f}^{(0)*}+\left(\frac{\alpha_{s}}{2\pi}\right)2\,\text{Re}\left(\mathcal{M}_{H\to
f}^{(1)}\mathcal{M}_{H\to f}^{(0)*}\right)$
$\displaystyle+\left(\frac{\alpha_{s}}{2\pi}\right)^{2}\bigg{\\{}2\,\text{Re}\left(\mathcal{M}_{H\to
f}^{(2)}\mathcal{M}_{H\to f}^{(0)*}\right)+\mathcal{M}_{H\to
f}^{(1)}\mathcal{M}_{H\to f}^{(1)*}\bigg{\\}}\Bigg{]}\,,$ (128)

where the Born color factor $\mathcal{C}_{f}$ is $2C_{A}^{2}C_{F}$ and
$C_{A}C_{F}$ for $f=ggg,\,q\bar{q}g$. We have made use of the same program as
described in section 3 for calculating those and could reproduce the same IR-
finite results as the squared amplitudes constructed from the
$\alpha,\,\beta,$ and $\gamma$, which gives us confidence in the related
calculation for the pseudo Higgs amplitudes. Furthermore, we have performed a
test for soft and collinear limit behavior of the helicity amplitudes and they
all are consistent with the results from Ref. Badger:2004uk  in which the
splitting functions have been extracted from Ref. Gehrmann:2011aa .
In order to utilize the amplitudes described above for LHC scattering the
partons must be crossed to the initial state via analytic continuation. We
confirmed our crossed integrals numerically with AMFlow Liu:2022chg  finding
perfect agreement. We have also checked our crossed results against the
literature and, found perfect agreement with the results from the recent paper
Gehrmann:2023etk . We note that there is a small typo in Ref. Gehrmann:2011aa
for $\gamma^{(2)}$ of the $q\bar{q}\to Hg$ kinematic region (the
region$(2a_{+})$ as denoted in Gehrmann:2002zr ). Since the results of Ref.
Gehrmann:2011aa  have been widely used, this discrepancy can give small
modifications to some other published articles. For example, in Ref.
Becher:2013vva , the authors assembled the two-loop hard functions for $H+j$.
The numerical results for the hard functions are provided in Eq.(73) in
section.5, which we present again here for the reader:

$\displaystyle\hat{H}_{gg\to Hg}(\hat{u},\hat{t},\mu)$
$\displaystyle=1+(6.02164)\,\alpha_{s}+(24.2724)\,\alpha_{s}^{2}\,,$
$\displaystyle\hat{H}_{q\bar{q}\to Hg}(\hat{u},\hat{t},\mu)$
$\displaystyle=1+(1.85023)\,\alpha_{s}+(8.15565)\,\alpha_{s}^{2}\,,$
$\displaystyle\hat{H}_{qg\to Hq}(\hat{u},\hat{t},\mu)$
$\displaystyle=1+(6.63865)\,\alpha_{s}+(24.9851)\,\alpha_{s}^{2}\,,$ (129)

at the following inputs:

$\displaystyle(\hat{s},\hat{t})$
$\displaystyle=(1\,\text{TeV}^{2},\,-0.4\,\text{TeV}^{2})\,,$ $\displaystyle
m_{H}$ $\displaystyle=0.1\,\text{TeV}\,,$ $\displaystyle\mu_{R}$
$\displaystyle=0.6\,\text{TeV}\,.$ (130)

With our calculation we have reproduced the same numerical results for
$\hat{H}_{gg\to Hg}$ and $\hat{H}_{qg\to Hq}\,$, but have obtained different
results for $\hat{H}_{q\bar{q}\to Hg}$, namely, :

$\displaystyle\hat{H}_{q\bar{q}\to Hg}(\hat{u},\hat{t},\mu)$
$\displaystyle=1+(1.85023)\,\alpha_{s}+(7.54245)\,\alpha_{s}^{2}\,.$ (131)

A natural concern is that this difference may result in a phenomenological
difference for the various calculations of $H+j$@NNLO, which used this hard
function. However, due to the relative smallness of the $q\overline{q}$
initial state (as can be seen from the breakdown of the partonic channels for
instance in Ref. Campbell:2019gmd ), the impact of the different hard function
is negligible at the level of the total cross section. We have assembled all
the IR-finite squared amplitudes and the hard functions for each kinematic
region up to NNLO for an arbitrary renormalization scale. These are attached
as auxiliary files to the arXiv submission, which may be useful in further
comparisons of $H@$N3LO and $H+j$@NNLO.

##  Appendix E Descriptions of the ancillary files

In this appendix we describe the structure and the notation of the attached
Mathematica ancillary files. The IR-finite squared amplitudes for the SM Higgs
are written in the files with a name which begins with Hggg_IR and Hqqg_IR for
the process $f=ggg$ and $f=q\bar{q}g$ respectively. Likewise, the hard
functions for the SM Higgs are presented in the files with a name which begins
Hggg_Hard and Hqqg_Hard. For the pseudoscalar Higgs case, the files are named
in a similar manner which begins with A instead of H. We have normalized the
each amplitude for $f=ggg,q\bar{q}g$ processes by $2C_{A}^{2}C_{F}$ and
$C_{A}C_{F}$ respectively. Then any one-loop and two-loop(one-loop self
interference) amplitudes can be decomposed as

$\displaystyle\mathcal{M}_{f}^{(1)}\mathcal{M}_{f}^{(0)*}$
$\displaystyle=N_{C}\left(\mathcal{M}_{f}^{(1)}\mathcal{M}_{f}^{(0)*}\right)^{N_{C}}+\frac{1}{N_{C}}\left(\mathcal{M}_{f}^{(1)}\mathcal{M}_{f}^{(0)*}\right)^{N_{C}^{-1}}+N_{F}\left(\mathcal{M}_{f}^{(1)}\mathcal{M}_{f}^{(0)*}\right)^{N_{F}}\,,$
(132) $\displaystyle\mathcal{M}_{f}^{(2(1))}\mathcal{M}_{f}^{(0(1))*}$
$\displaystyle=N_{C}^{2}\left(\mathcal{M}_{f}^{(2(1))}\mathcal{M}_{f}^{(0(1))*}\right)^{N_{C}^{2}}+N_{C}^{0}\left(\mathcal{M}_{f}^{(2(1))}\mathcal{M}_{f}^{(0(1))*}\right)^{N_{C}^{0}}$
$\displaystyle+\frac{1}{N_{C}^{2}}\left(\mathcal{M}_{f}^{(2(1))}\mathcal{M}_{f}^{(0(1))*}\right)^{N_{C}^{-2}}+\frac{N_{F}}{N_{C}}\left(\mathcal{M}_{f}^{(2(1))}\mathcal{M}_{f}^{(0(1))*}\right)^{\frac{N_{F}}{N_{C}}}$
$\displaystyle+N_{F}N_{C}\left(\mathcal{M}_{f}^{(2(1))}\mathcal{M}_{f}^{(0(1))*}\right)^{N_{F}N_{C}}+N_{F}^{2}\left(\mathcal{M}_{f}^{(2(1))}\mathcal{M}_{f}^{(0(1))*}\right)^{N_{F}^{2}}\,.$
(133)

The evolution of the renormalization scale $\mu_{R}$ has been carried for the
all orders in the attached files. For the LHC application the squared
amplitudes need to be crossed to the appropriate kinematic regions. To do so,
the required analytic continuations are from the region(1) to the region(2)
and (4)(according to the Ref. Becher:2013vva ; Gehrmann:2002zr ) such as

$\displaystyle g(p_{1})+g(p_{2})\to A(H)(p_{4})+g(-p_{3})$ (134)

from $f=ggg$ and

$\displaystyle q(p_{1})+\bar{q}(p_{2})\to
A(H)(p_{4})+g(-p_{3})\,,\quad\bar{q}(p_{2})+g(p_{3})\to
A(H)(p_{4})+\bar{q}(-p_{1})$ (135)

from $f=q\bar{q}g$. The dimensionless variables for the region(2) and the
region(4) are defined as Becher:2013vva ; Gehrmann:2002zr :

$\displaystyle u_{2}=-\frac{y}{x}\,,\qquad v_{2}=\frac{1}{x}\,,$ (136)
$\displaystyle u_{4}=-\frac{y}{z}\,,\qquad v_{4}=\frac{1}{z}\,.$ (137)

The IR-finite squared amplitudes and the hard functions for the each region is
attached to the file name ends with R1, R2, and R4 respectively. The notation
of the attached Mathematica files are as follows:

$\displaystyle N_{C}$ $\displaystyle\to\texttt{Nc}\,,\qquad
N_{F}\to\texttt{Nf},$ $\displaystyle u_{2,4}$
$\displaystyle\to\texttt{u}\,,\qquad v_{2,4}\to\texttt{v},$
$\displaystyle\zeta_{2},\zeta_{3},\zeta_{4}$ $\displaystyle\to\texttt{zeta2,
zeta3, zeta4}\,,$ $\displaystyle\frac{\mu_{R}^{2}}{m_{H}^{2}}$
$\displaystyle\to\texttt{muSqOMH2}\,,\quad\frac{\mu_{R}^{2}}{m_{A}^{2}}\to\texttt{muSqOMA2}\,.$
(138)

For the convention of the $1D$ and $2D$ HPLs we followed the notation
introduced in the Ref. Gehrmann:2000zt . For instance, we denote

$\displaystyle H(0,1,1,0;x)\to\texttt{H[0,1,1,0,x]}\,,\qquad
H(1,z,1-z;y)\to\texttt{H[1,z,1-z,y]}\,.$ (139)

The required terms at $\text{LO}(\mathcal{O}(\alpha_{s}^{3}))$,
$\text{NLO}(\mathcal{O}(\alpha_{s}^{4}))$, and
$\text{NNLO}(\mathcal{O}(\alpha_{s}^{5}))$ expanded as (D), (25), and (26) for
the Higgs and the pseudoscalar Higgs are written in the files for the each
kinematic region. The notation for the terms for the Higgs case are written
as:

$\displaystyle\mathcal{M}_{H\to f}^{(0)}\mathcal{M}_{H\to f}^{(0)*}$
$\displaystyle\to\texttt{M0M0}\,,\qquad\mathcal{M}_{H\to
f}^{(1)}\mathcal{M}_{H\to f}^{(0)*}\to\texttt{M0M1}\,,$
$\displaystyle\mathcal{M}_{H\to f}^{(1)}\mathcal{M}_{H\to f}^{(1)*}$
$\displaystyle\to\texttt{M1M1}\,,\qquad\mathcal{M}_{H\to
f}^{(2)}\mathcal{M}_{H\to f}^{(0)*}\to\texttt{M0M2}\,,$ (140)

for the both $f=ggg,q\bar{q}g$. In a similar manner, the pseudo scalar Higgs
amplitudes are written as,

$\displaystyle\mathcal{M}_{q\bar{q}g}^{G,(0)}\mathcal{M}_{q\bar{q}g}^{G,(0)*}$
$\displaystyle\to\texttt{MG0MG0}\,,\qquad\mathcal{M}_{q\bar{q}g}^{G,(1)}\mathcal{M}_{q\bar{q}g}^{G,(0)*}\to\texttt{MG0MG1}\,,$
$\displaystyle\mathcal{M}_{q\bar{q}g}^{G,(1)}\mathcal{M}_{q\bar{q}g}^{G,(1)*}$
$\displaystyle\to\texttt{MG1MG1}\,,\qquad\mathcal{M}_{q\bar{q}g}^{G,(2)}\mathcal{M}_{q\bar{q}g}^{G,(0)*}\to\texttt{MG0MG2}\,,$
$\displaystyle\mathcal{M}_{q\bar{q}g}^{G,(0)}\mathcal{M}_{q\bar{q}g}^{J,(0)*}$
$\displaystyle\to\texttt{MG0MJ0}\,,\qquad\mathcal{M}_{q\bar{q}g}^{J,(0)}\mathcal{M}_{q\bar{q}g}^{J,(0)*}\to\texttt{MJ0MJ0}\,,$
$\displaystyle\mathcal{M}_{q\bar{q}g}^{G,(1)}\mathcal{M}_{q\bar{q}g}^{J,(0)*}$
$\displaystyle+\mathcal{M}_{q\bar{q}g}^{G,(0)}\mathcal{M}_{q\bar{q}g}^{J,(1)*}\to\texttt{MG0MJ1plusMG1MJ0}\,,$
(141)

for $f=q\bar{q}g$. The $f=ggg$ amplitudes are written in the same fashion.
Readers can download the files and import the results into a Mathematica
environment using Get or <<.

## References

  * (1) ATLAS Collaboration, G. Aad et. al., 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–29 [1207.7214].
  * (2) CMS Collaboration, S. Chatrchyan et. al., Observation of a New Boson at a Mass of 125 GeV with the CMS Experiment at the LHC, Phys. Lett. B 716 (2012) 30–61 [1207.7235].
  * (3) 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].
  * (4) B. Mistlberger, Higgs boson production at hadron colliders at N3LO in QCD, JHEP 05 (2018) 028 [1802.00833].
  * (5) 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), no. 7 072002 [2102.07607].
  * (6) 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 in perturbative QCD, JHEP 06 (2013) 072 [1302.6216].
  * (7) X. Chen, T. Gehrmann, E. W. N. Glover and M. Jaquier, Precise QCD predictions for the production of Higgs + jet final states, Phys. Lett. B 740 (2015) 147–150 [1408.5325].
  * (8) 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), no. 8 082003 [1504.07922].
  * (9) 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–8 [1505.03893].
  * (10) 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), no. 7 074032 [1508.02684].
  * (11) 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].
  * (12) J. M. Campbell, R. K. Ellis and S. Seth, H + 1 jet production revisited, JHEP 10 (2019) 136 [1906.01020].
  * (13) M. Spira, A. Djouadi, D. Graudenz and P. M. Zerwas, Higgs boson production at the LHC, Nucl. Phys. B 453 (1995) 17–82 [hep-ph/9504378].
  * (14) A. Pak, M. Rogal and M. Steinhauser, Finite top quark mass effects in NNLO Higgs boson production at LHC, JHEP 02 (2010) 025 [0911.4662].
  * (15) R. V. Harlander and K. J. Ozeren, Finite top mass effects for hadronic Higgs production at next-to-next-to-leading order, JHEP 11 (2009) 088 [0909.3420].
  * (16) P. Fayet, Supergauge Invariant Extension of the Higgs Mechanism and a Model for the electron and Its Neutrino, Nucl. Phys. B 90 (1975) 104–124.
  * (17) S. Dimopoulos and H. Georgi, Softly Broken Supersymmetry and SU(5), Nucl. Phys. B 193 (1981) 150–162.
  * (18) N. Sakai, Naturalness in Supersymmetric Guts, Z. Phys. C 11 (1981) 153.
  * (19) K. Inoue, A. Kakuto, H. Komatsu and S. Takeshita, Aspects of Grand Unified Models with Softly Broken Supersymmetry, Prog. Theor. Phys. 68 (1982) 927. [Erratum: Prog.Theor.Phys. 70, 330 (1983)].
  * (20) K. Inoue, A. Kakuto, H. Komatsu and S. Takeshita, Renormalization of Supersymmetry Breaking Parameters Revisited, Prog. Theor. Phys. 71 (1984) 413.
  * (21) K. Inoue, A. Kakuto, H. Komatsu and S. Takeshita, Low-Energy Parameters and Particle Masses in a Supersymmetric Grand Unified Model, Prog. Theor. Phys. 67 (1982) 1889.
  * (22) S. Dawson, Radiative corrections to Higgs boson production, Nucl. Phys. B 359 (1991) 283–300.
  * (23) A. Djouadi, M. Spira and P. M. Zerwas, Production of Higgs bosons in proton colliders: QCD corrections, Phys. Lett. B 264 (1991) 440–446.
  * (24) R. P. Kauffman and W. Schaffer, QCD corrections to production of Higgs pseudoscalars, Phys. Rev. D 49 (1994) 551–554 [hep-ph/9305279].
  * (25) M. Spira, A. Djouadi, D. Graudenz and P. M. Zerwas, SUSY Higgs production at proton colliders, Phys. Lett. B 318 (1993) 347–353.
  * (26) R. V. Harlander and W. B. Kilgore, Next-to-next-to-leading order Higgs production at hadron colliders, Phys. Rev. Lett. 88 (2002) 201801 [hep-ph/0201206].
  * (27) C. Anastasiou and K. Melnikov, Higgs boson production at hadron colliders in NNLO QCD, Nucl. Phys. B 646 (2002) 220–256 [hep-ph/0207004].
  * (28) V. Ravindran, J. Smith and W. L. van Neerven, NNLO corrections to the total cross-section for Higgs boson production in hadron hadron collisions, Nucl. Phys. B 665 (2003) 325–366 [hep-ph/0302135].
  * (29) C. Anastasiou, K. Melnikov and F. Petriello, Fully differential Higgs boson production and the di-photon signal through next-to-next-to-leading order, Nucl. Phys. B 724 (2005) 197–246 [hep-ph/0501130].
  * (30) R. V. Harlander and W. B. Kilgore, Production of a pseudoscalar Higgs boson at hadron colliders at next-to-next-to leading order, JHEP 10 (2002) 017 [hep-ph/0208096].
  * (31) C. Anastasiou and K. Melnikov, Pseudoscalar Higgs boson production at hadron colliders in NNLO QCD, Phys. Rev. D 67 (2003) 037501 [hep-ph/0208115].
  * (32) T. Ahmed, M. C. Kumar, P. Mathews, N. Rana and V. Ravindran, Pseudo-scalar Higgs boson production at threshold N3 LO and N3 LL QCD, Eur. Phys. J. C 76 (2016), no. 6 355 [1510.02235].
  * (33) T. Ahmed, M. Bonvini, M. C. Kumar, P. Mathews, N. Rana, V. Ravindran and L. Rottoli, Pseudo-scalar Higgs boson production at N3 LO${}_{\text{A}}$ +N3 LL ′, Eur. Phys. J. C 76 (2016), no. 12 663 [1606.00837].
  * (34) D. de Florian, M. Grazzini and Z. Kunszt, Higgs production with large transverse momentum in hadronic collisions at next-to-leading order, Phys. Rev. Lett. 82 (1999) 5209–5212 [hep-ph/9902483].
  * (35) V. Ravindran, J. Smith and W. L. Van Neerven, Next-to-leading order QCD corrections to differential distributions of Higgs boson production in hadron hadron collisions, Nucl. Phys. B 634 (2002) 247–290 [hep-ph/0201114].
  * (36) B. Field, J. Smith, M. E. Tejeda-Yeomans and W. L. van Neerven, NLO corrections to differential cross-sections for pseudoscalar Higgs boson production, Phys. Lett. B 551 (2003) 137–145 [hep-ph/0210369].
  * (37) W. Bernreuther, P. Gonzalez and M. Wiebusch, Pseudoscalar Higgs Bosons at the LHC: Production and Decays into Electroweak Gauge Bosons Revisited, Eur. Phys. J. C 69 (2010) 31–43 [1003.5585].
  * (38) J. Gaunt, M. Stahlhofen, F. J. Tackmann and J. R. Walsh, N-jettiness Subtractions for NNLO QCD Calculations, JHEP 09 (2015) 058 [1505.04794].
  * (39) R. Boughezal, C. Focke, X. Liu and F. Petriello, $W$-boson production in association with a jet at next-to-next-to-leading order in perturbative QCD, Phys. Rev. Lett. 115 (2015), no. 6 062002 [1504.02131].
  * (40) A. Gehrmann-De Ridder, T. Gehrmann and E. W. N. Glover, Antenna subtraction at NNLO, JHEP 09 (2005) 056 [hep-ph/0505111].
  * (41) A. Daleo, T. Gehrmann and D. Maitre, Antenna subtraction with hadronic initial states, JHEP 04 (2007) 016 [hep-ph/0612257].
  * (42) A. Daleo, A. Gehrmann-De Ridder, T. Gehrmann and G. Luisoni, Antenna subtraction at NNLO with hadronic initial states: initial-final configurations, JHEP 01 (2010) 118 [0912.0374].
  * (43) T. Gehrmann and P. F. Monni, Antenna subtraction at NNLO with hadronic initial states: real-virtual initial-initial configurations, JHEP 12 (2011) 049 [1107.4037].
  * (44) R. Boughezal, A. Gehrmann-De Ridder and M. Ritzmann, Antenna subtraction at NNLO with hadronic initial states: double real radiation for initial-initial configurations with two quark flavours, JHEP 02 (2011) 098 [1011.6631].
  * (45) A. Gehrmann-De Ridder, T. Gehrmann and M. Ritzmann, Antenna subtraction at NNLO with hadronic initial states: double real initial-initial configurations, JHEP 10 (2012) 047 [1207.5779].
  * (46) J. Currie, E. W. N. Glover and S. Wells, Infrared Structure at NNLO Using Antenna Subtraction, JHEP 04 (2013) 066 [1301.4693].
  * (47) R. Boughezal, J. M. Campbell, R. K. Ellis, C. Focke, W. Giele, X. Liu, F. Petriello and C. Williams, Color singlet production at NNLO in MCFM, Eur. Phys. J. C 77 (2017), no. 1 7 [1605.08011].
  * (48) J. Campbell and T. Neumann, Precision Phenomenology with MCFM, JHEP 12 (2019) 034 [1909.09117].
  * (49) R. Mondini and C. Williams, Bottom-induced contributions to Higgs plus jet at next-to-next-to-leading order, JHEP 05 (2021) 045 [2102.05487].
  * (50) K. G. Chetyrkin, B. A. Kniehl, M. Steinhauser and W. A. Bardeen, Effective QCD interactions of CP odd Higgs bosons at three loops, Nucl. Phys. B 535 (1998) 3–18 [hep-ph/9807241].
  * (51) P. Banerjee, P. K. Dhani and V. Ravindran, Two loop QCD corrections for the process Pseudo-scalar Higgs $\rightarrow 3$ partons, JHEP 10 (2017) 067 [1708.02387].
  * (52) S. L. Adler, Axial vector vertex in spinor electrodynamics, Phys. Rev. 177 (1969) 2426–2438.
  * (53) I. W. Stewart, F. J. Tackmann and W. J. Waalewijn, N-Jettiness: An Inclusive Event Shape to Veto Jets, Phys. Rev. Lett. 105 (2010) 092002 [1004.2489].
  * (54) I. W. Stewart, F. J. Tackmann and W. J. Waalewijn, Factorization at the LHC: From PDFs to Initial State Jets, Phys. Rev. D 81 (2010) 094035 [0910.0467].
  * (55) J. M. Campbell, R. K. Ellis, R. Mondini and C. Williams, The NNLO QCD soft function for 1-jettiness, Eur. Phys. J. C 78 (2018), no. 3 234 [1711.09984].
  * (56) R. Boughezal, X. Liu and F. Petriello, $N$-jettiness soft function at next-to-next-to-leading order, Phys. Rev. D 91 (2015), no. 9 094035 [1504.02540].
  * (57) T. Becher and G. Bell, The gluon jet function at two-loop order, Phys. Lett. B 695 (2011) 252–258 [1008.1936].
  * (58) T. Becher and M. Neubert, Toward a NNLO calculation of the anti-B —$>$ X(s) gamma decay rate with a cut on photon energy. II. Two-loop result for the jet function, Phys. Lett. B 637 (2006) 251–259 [hep-ph/0603140].
  * (59) J. R. Gaunt, M. Stahlhofen and F. J. Tackmann, The Quark Beam Function at Two Loops, JHEP 04 (2014) 113 [1401.5478].
  * (60) J. Gaunt, M. Stahlhofen and F. J. Tackmann, The Gluon Beam Function at Two Loops, JHEP 08 (2014) 020 [1405.1044].
  * (61) T. T. Jouttenus, I. W. Stewart, F. J. Tackmann and W. J. Waalewijn, The Soft Function for Exclusive N-Jet Production at Hadron Colliders, Phys. Rev. D 83 (2011) 114030 [1102.4344].
  * (62) T. T. Jouttenus, I. W. Stewart, F. J. Tackmann and W. J. Waalewijn, Jet mass spectra in Higgs boson plus one jet at next-to-next-to-leading logarithmic order, Phys. Rev. D 88 (2013), no. 5 054031 [1302.0846].
  * (63) S. A. Larin, The Renormalization of the axial anomaly in dimensional regularization, Phys. Lett. B 303 (1993) 113–118 [hep-ph/9302240].
  * (64) T. Ahmed, T. Gehrmann, P. Mathews, N. Rana and V. Ravindran, Pseudo-scalar Form Factors at Three Loops in QCD, JHEP 11 (2015) 169 [1510.01715].
  * (65) M. F. Zoller, OPE of the pseudoscalar gluonium correlator in massless QCD to three-loop order, JHEP 07 (2013) 040 [1304.2232].
  * (66) T. Hahn, Generating Feynman diagrams and amplitudes with FeynArts 3, Comput. Phys. Commun. 140 (2001) 418–431 [hep-ph/0012260].
  * (67) P. Nogueira, Automatic Feynman graph generation, J. Comput. Phys. 105 (1993) 279–289.
  * (68) N. D. Christensen and C. Duhr, FeynRules - Feynman rules made easy, Comput. Phys. Commun. 180 (2009) 1614–1641 [0806.4194].
  * (69) A. Alloul, N. D. Christensen, C. Degrande, C. Duhr and B. Fuks, FeynRules 2.0 - A complete toolbox for tree-level phenomenology, Comput. Phys. Commun. 185 (2014) 2250–2300 [1310.1921].
  * (70) R. Mertig, M. Bohm and A. Denner, FEYN CALC: Computer algebraic calculation of Feynman amplitudes, Comput. Phys. Commun. 64 (1991) 345–359.
  * (71) V. Shtabovenko, R. Mertig and F. Orellana, FeynCalc 9.3: New features and improvements, Comput. Phys. Commun. 256 (2020) 107478 [2001.04407].
  * (72) R. N. Lee, LiteRed 1.4: a powerful tool for reduction of multiloop integrals, J. Phys. Conf. Ser. 523 (2014) 012059 [1310.1145].
  * (73) P. Maierhöfer, J. Usovitsch and P. Uwer, Kira—A Feynman integral reduction program, Comput. Phys. Commun. 230 (2018) 99–112 [1705.05610].
  * (74) L. W. Garland, T. Gehrmann, E. W. N. Glover, A. Koukoutsakis and E. Remiddi, The Two loop QCD matrix element for e+ e- —$>$ 3 jets, Nucl. Phys. B 627 (2002) 107–188 [hep-ph/0112081].
  * (75) T. Gehrmann and E. Remiddi, Two loop master integrals for gamma* —$>$ 3 jets: The Planar topologies, Nucl. Phys. B 601 (2001) 248–286 [hep-ph/0008287].
  * (76) T. Gehrmann and E. Remiddi, Two loop master integrals for gamma* –$>$ 3 jets: The Nonplanar topologies, Nucl. Phys. B 601 (2001) 287–317 [hep-ph/0101124].
  * (77) E. Remiddi and J. A. M. Vermaseren, Harmonic polylogarithms, Int. J. Mod. Phys. A 15 (2000) 725–754 [hep-ph/9905237].
  * (78) C. W. Bauer, A. Frink and R. Kreckel, Introduction to the GiNaC framework for symbolic computation within the C++ programming language, J. Symb. Comput. 33 (2002) 1–12 [cs/0004015].
  * (79) L. Naterop, A. Signer and Y. Ulrich, handyG —Rapid numerical evaluation of generalised polylogarithms in Fortran, Comput. Phys. Commun. 253 (2020) 107165 [1909.01656].
  * (80) S. Catani, The Singular behavior of QCD amplitudes at two loop order, Phys. Lett. B 427 (1998) 161–171 [hep-ph/9802439].
  * (81) T. Becher, G. Bell, C. Lorentzen and S. Marti, Transverse-momentum spectra of electroweak bosons near threshold at NNLO, JHEP 02 (2014) 004 [1309.3245].
  * (82) T. Gehrmann and L. Tancredi, Two-loop QCD helicity amplitudes for $q\bar{q}\to W^{\pm}\gamma$ and $q\bar{q}\to Z^{0}\gamma$, JHEP 02 (2012) 004 [1112.1531].
  * (83) T. Gehrmann and E. Remiddi, Analytic continuation of massless two loop four point functions, Nucl. Phys. B 640 (2002) 379–411 [hep-ph/0207020].
  * (84) C. Duhr, Mathematical aspects of scattering amplitudes, in Theoretical Advanced Study Institute in Elementary Particle Physics: Journeys Through the Precision Frontier: Amplitudes for Colliders, pp. 419–476, 2015.  1411.7538.
  * (85) C. Duhr, Hopf algebras, coproducts and symbols: an application to Higgs boson amplitudes, JHEP 08 (2012) 043 [1203.0454].
  * (86) R. Mondini and C. Williams, $H\to b\overline{b}j$ at next-to-next-to-leading order accuracy, JHEP 06 (2019) 120 [1904.08961].
  * (87) X. Liu and Y.-Q. Ma, AMFlow: A Mathematica package for Feynman integrals computation via auxiliary mass flow, Comput. Phys. Commun. 283 (2023) 108565 [2201.11669].
  * (88) S. D. Badger and E. W. N. Glover, Two loop splitting functions in QCD, JHEP 07 (2004) 040 [hep-ph/0405236].
  * (89) L. J. Dixon, E. W. N. Glover and V. V. Khoze, MHV rules for Higgs plus multi-gluon amplitudes, JHEP 12 (2004) 015 [hep-th/0411092].
  * (90) S. D. Badger and E. W. N. Glover, One-loop helicity amplitudes for H —$>$ gluons: The All-minus configuration, Nucl. Phys. B Proc. Suppl. 160 (2006) 71–75 [hep-ph/0607139].
  * (91) L. J. Dixon and Y. Sofianatos, Analytic one-loop amplitudes for a Higgs boson plus four partons, JHEP 08 (2009) 058 [0906.0008].
  * (92) S. Badger, E. W. Nigel Glover, P. Mastrolia and C. Williams, One-loop Higgs plus four gluon amplitudes: Full analytic results, JHEP 01 (2010) 036 [0909.4475].
  * (93) S. Badger, J. M. Campbell, R. K. Ellis and C. Williams, Analytic results for the one-loop NMHV Hqqgg amplitude, JHEP 12 (2009) 035 [0910.4481].
  * (94) R. Britto, F. Cachazo, B. Feng and E. Witten, Direct proof of tree-level recursion relation in Yang-Mills theory, Phys. Rev. Lett. 94 (2005) 181602 [hep-th/0501052].
  * (95) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli and M. Zaro, The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations, JHEP 07 (2014) 079 [1405.0301].
  * (96) GoSam Collaboration, G. Cullen et. al., G$\scriptsize{O}$S$\scriptsize{AM}$-2.0: a tool for automated one-loop calculations within the Standard Model and beyond, Eur. Phys. J. C 74 (2014), no. 8 3001 [1404.7096].
  * (97) S. Catani and M. H. Seymour, A General algorithm for calculating jet cross-sections in NLO QCD, Nucl. Phys. B 485 (1997) 291–419 [hep-ph/9605323]. [Erratum: Nucl.Phys.B 510, 503–504 (1998)].
  * (98) Z. Nagy and Z. Trocsanyi, Next-to-leading order calculation of four jet observables in electron positron annihilation, Phys. Rev. D 59 (1999) 014020 [hep-ph/9806317]. [Erratum: Phys.Rev.D 62, 099902 (2000)].
  * (99) W. Bizoń, X. Chen, A. Gehrmann-De Ridder, T. Gehrmann, N. Glover, A. Huss, P. F. Monni, E. Re, L. Rottoli and P. Torrielli, Fiducial distributions in Higgs and Drell-Yan production at N3LL+NNLO, JHEP 12 (2018) 132 [1805.05916].
  * (100) S. Catani and B. R. Webber, Infrared safe but infinite: Soft gluon divergences inside the physical region, JHEP 10 (1997) 005 [hep-ph/9710333].
  * (101) T. Gehrmann, M. Jaquier, E. W. N. Glover and A. Koukoutsakis, Two-Loop QCD Corrections to the Helicity Amplitudes for $H\to$ 3 partons, JHEP 02 (2012) 056 [1112.3554].
  * (102) T. Gehrmann, P. Jakubčík, C. C. Mella, N. Syrrakos and L. Tancredi, Two-loop helicity amplitudes for $H+$jet production to higher orders in the dimensional regulator, JHEP 04 (2023) 016 [2301.10849].
  * (103) K. G. Chetyrkin, B. A. Kniehl and M. Steinhauser, Hadronic Higgs decay to order alpha-s**4, Phys. Rev. Lett. 79 (1997) 353–356 [hep-ph/9705240].
  * (104) 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–87 [hep-ph/9708255].