aainstitutetext:  Dipartimento di Fisica and Arnold-Regge Center, Università
di Torino, and INFN, Sezione di Torino, Via P. Giuria 1, I-10125 Torino, Italy
bbinstitutetext:  Cavendish Laboratory, University of Cambridge, Cambridge CB3
0HE, United Kingdom ccinstitutetext:  Institute for Particle Physics
Phenomenology, Department of Physics, Durham University, Durham DH1 3LE,
United Kingdom

# Two-loop leading-colour QCD helicity amplitudes for Higgs boson production
in association with a bottom-quark pair at the LHC

Simon Badger b Heribertus Bayu Hartanto a,c Jakub Kryś a Simone Zoia
simondavid.badger@unito,it<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We compute the two-loop QCD helicity amplitudes for the production of a Higgs
boson in association with a bottom quark pair at a hadron collider. We take
the approximations of leading colour and work in the five flavour scheme,
where the bottom quarks are massless while the Yukawa coupling is non-zero. We
extract analytic expressions from multiple numerical evaluations over finite
fields and present the results in terms of an independent set of special
functions that can be reliably evaluated over the full phase space.

††preprint: CAVENDISH-HEP-21/11

##  1 Introduction

Precise theoretical predictions are an essential ingredient in the search for
subtle deviations from the Standard Model (SM) at current and future collider
experiments. The enormous amount of data gathered at the Large Hadron Collider
(LHC) is already challenging the theoretical precision for many processes and
the pressure will increase as data continues to pour in. It has been clear for
some time Badger:2016bpw ; Bendavid:2018nar ; Amoroso:2020lgh  that, for a
large class of two- and three-particle final states, at least the Next-to-
Next-to-Leading Order corrections in Quantum Chromodynamics (NNLO QCD) will be
required for fully differential cross sections. The high multiplicity final
states, particularly those with many kinematic scales, pose the biggest
technical challenge owing to the enormous analytic and algebraic complexity.
Nevertheless, thanks to a huge effort across the theoretical community, new
tools and methods have been developed that have produced the first NNLO QCD
predictions for $2\to 3$ processes Chawdhry:2019bji ; Kallweit:2020gcp ;
Chawdhry:2021hkp , most notably for 3-jet production Czakon:2021mjy . This
remarkable progress has been driven both by the advancements in the scattering
amplitude computation and by efficient subtraction schemes. By now, all the
two-loop master integrals – which are one of the important ingredients in the
amplitude computation – are available for the fully massless five-particle
processes Papadopoulos:2015jft ; Gehrmann:2018yef ; Abreu:2018rcw ;
Chicherin:2018mue ; Chicherin:2018old ; Chicherin:2020oor , allowing for
several two-loop QCD amplitudes to be derived analytically Gehrmann:2015bfy ;
Badger:2018enw ; Abreu:2018zmy ; Abreu:2019odu ; Badger:2019djh ;
Abreu:2020cwb ; Chawdhry:2020for ; Agarwal:2021grm ; Abreu:2021fuk ;
Chawdhry:2021mkw ; Agarwal:2021vdh ; Badger:2021imn , improving on previous
results that were obtained numerically Badger:2017jhb ; Abreu:2017hqn ;
Badger:2018gip ; Abreu:2018jgq . These new results have been achieved thanks
to technological breakthroughs in the method of differential equations
Kotikov:1990kg ; Bern:1993kr ; Remiddi:1997ny ; Gehrmann:1999as ; Henn:2013pwa
, integral reduction algorithms Gluza:2010ws ; Schabinger:2011dz ; Ita:2015tya
; Larsen:2015ped ; Boehm:2017wjc  and the use of finite-field arithmetic to
tame the algebraic complexity of multi-leg and multi-scale problems
vonManteuffel:2014ixa ; Peraro:2016wsq ; Klappert:2019emp ; Peraro:2019svx ;
Klappert:2020aqs ; Klappert:2020nbg .

A different class of $2\to 3$ processes that is of great interest involves an
external massive leg. The two-loop planar helicity amplitudes for $W+4$-parton
scattering (contributing to the prediction for $W+2$-jet production at NNLO
QCD) have been previously studied numerically Hartanto:2019uvl . All two-loop
integrals needed for this type of process are available for the planar
contributions Papadopoulos:2015jft ; Abreu:2020jxa ; Canko:2020ylt ;
Syrrakos:2020kba , and the first analytic result was derived for $Wb\bar{b}$
production in the leading colour, massless $b$-quark and on-shell $W$-boson
approximations Badger:2021nhg . Very recently, complete analytic results for
one of the non-planar integral topologies have also become available
Papadopoulos:2019iam ; abreu2021twoloop .

In this article we consider the two-loop amplitudes relevant for the
production of a Higgs boson in association with a bottom-quark pair, i.e.
$pp\to b\bar{b}H$, in the leading colour approximation. $b\bar{b}H$ production
at the LHC has been a subject of great phenomenological interest due to its
potential in directly measuring the bottom-quark Yukawa coupling. In the
Standard Model (SM), the coupling strengths of the Higgs boson to the fermions
and vector bosons are proportional to their mass, causing the rate of the
$b\bar{b}H$ production to be suppressed with respect to, for example, Higgs
production in gluon fusion ($gg\to H$) or vector boson fusion ($pp\to Hjj$),
associated production with a vector boson ($pp\to VH$), and associated
production with a top-quark pair ($pp\to t\bar{t}H$). In addition, the
presence of $b$-tagging further suppresses the $b\bar{b}H$ production rate. In
some new physics scenarios, such as the Two Higgs Doublet Models (2HDM’s) and
the Minimal Supersymmetric Standard Model (MSSM), the bottom-quark Yukawa
coupling can be dramatically enhanced, resulting in a considerable increase of
the $b\bar{b}H$ production rate Balazs:1998nt ; Dawson:2005vi . Thus, the
study of the $b\bar{b}H$ production will allow to constrain supersymmetric
models and other extensions of the SM that modify the bottom-quark Yukawa
coupling. Recent studies on the interplay between $b\bar{b}H$ signal and
backgrounds can be found in Refs. Pagani:2020rsg ; Grojean:2020ech ;
Konar:2021nkk .

The theoretical approach to obtaining predictions for the $pp\to b\bar{b}H$
process has been subject of much discussion in the community. This is due to
the fact that, in the presence of bottom quarks, a theoretical prediction can
be computed in either the four-flavour scheme (4FS) or the five-flavour scheme
(5FS). In the 4FS computation, bottom quarks are treated as massive and they
do not contribute to the parton distributions functions (PDFs), hence only
appearing in the final state. Large logarithms of the form $\log(m_{b}/Q)$
with $Q\propto m_{H}$ arise when the integration over the bottom-quark phase
space is performed, and such contributions may spoil the convergence of the
perturbative series. These large logarithms can be resummed to all orders by
introducing the bottom-quark parton distribution functions. The 5FS approach
stems from this prescription, where the bottom-quarks are included in the
parton distribution functions, allowed to appear in the initial state, and
treated as massless. We refer the reader to Ref. Maltoni:2012pa  for further
discussion on the 5FS and 4FS approaches. In 5FS the inclusive $b\bar{b}H$
production (where the tree level process is $b\bar{b}\to H$) has been computed
up to N3LO QCD Dicus:1998hs ; Balazs:1998sb ; Maltoni:2003pn ;
Harlander:2010cz ; Buehler:2012cu ; Harlander:2012pb ; H:2019nsw ;
Duhr:2019kwi ; Mondini:2021nck , while for the case where a single bottom
quark is observed NLO QCD Campbell:2002zm , weak Dawson:2010yz  and SUSY QCD
Dawson:2011pe  corrections are available. In 4FS the $b\bar{b}H$ production
has been calculated up to NLO QCD Dittmaier:2003ej ; Dawson:2003kb ;
Dawson:2004sh ; Wiesemann:2014ioa ; Jager:2015hka ; Deutschmann:2018avk , and
the supersymmetric QCD corrections Liu:2012qu ; Dittmaier:2014sva  are also
known. There have also been efforts in matching the 5FS and 4FS calculations
to obtain accurate predictions across the entire kinematic region
Harlander:2011aa ; Bonvini:2015pxa ; Forte:2015hba ; Forte:2016sja ;
Duhr:2020kzd . A first step towards a massive version of the five-flavour
scheme (5FMS) has been devised to naturally connect the 4FS and 5FS approaches
Krauss:2017wmx ; Figueroa:2018chn .

In this work we compute the two-loop QCD corrections to the $gg\rightarrow
b\bar{b}H$, $q\bar{q}/\bar{q}q\rightarrow b\bar{b}H$,
$b\bar{b}/\bar{b}b\rightarrow b\bar{b}H$, $bb\rightarrow bbH$ and
$\bar{b}\bar{b}\rightarrow\bar{b}\bar{b}H$ reactions in the 5FS. These two-
loop amplitudes enter the computation of $pp(b\bar{b})\to H$ at N4LO, $pp\to
b(\bar{b})H$ at N3LO when one $b$-jet is tagged, and $pp\to b\bar{b}H$ at NNLO
when two $b$-jets are required in the final state. We note that beyond NLO,
for the computation with massless bottom quarks, a flavoured jet algorithm
Banfi:2006hf  would have to be employed when identifying the $b$-jets, since
the use of conventional $k_{T}$ or anti-$k_{T}$ jet algorithms would render
the fixed order computation infrared unsafe. We further remark that the two-
loop amplitudes for $b\bar{b}H$ production derived in this work can also be
used in the computation of Higgs decaying into a bottom-quark pair in 5FS, by
crossing initial partons to the final state. Specifically, they will
contribute in the N4LO $H\to b\bar{b}$, N3LO $H\to b\bar{b}j$ and NNLO $H\to
b\bar{b}jj$ computations. In addition, by crossing the $b\bar{b}$ pair to the
initial state and the $gg/q\bar{q}$ pair to the final state we obtain the
contribution of the bottom-quark initiated channel to $H+2j$ production
($b\bar{b}\to Hjj$).

We present analytic results for the finite remainders after ultraviolet (UV)
and infrared (IR) poles have been subtracted. This is possible using a basis
of independent special functions recently identified in the context of
$Wb\bar{b}$ production Badger:2021nhg . We obtain numerical results valid
across the full phase space by applying the generalised series expansion
approach Francesco:2019yqt ; Abreu:2020jxa ; Hidding:2020ytt  to the
differential equations satisfied by the special functions appearing in the
finite remainders.

This paper is organised as follows. We begin by describing the structure of
the $b\bar{b}H$ amplitudes at leading colour in Section 2, followed by a
discussion of the methodology used in deriving the analytic expression of the
amplitudes in Section 3. A description of the basis of special functions is
presented in Section 4 and a number of validations that we performed on the
results derived in this work are discussed in Section 5. We present benchmark
numerical evaluations together with evaluations on a physical phase-space
slice in Section 6. Finally we draw our conclusions in Section 7.

##  2 Structure of the $pp\to b\bar{b}H$ Amplitudes at Leading Colour

We compute the two-loop QCD corrections in the leading colour approximation
for the following subprocesses

$\displaystyle
0\rightarrow\bar{b}(p_{1})+b(p_{2})+g(p_{3})+g(p_{4})+H(p_{5})\,,$ (1)
$\displaystyle
0\rightarrow\bar{b}(p_{1})+b(p_{2})+\bar{q}(p_{3})+q(p_{4})+H(p_{5})\,,$ (2)
$\displaystyle
0\rightarrow\bar{b}(p_{1})+b(p_{2})+\bar{b}(p_{3})+b(p_{4})+H(p_{5})\,,$ (3)

where all momenta are taken as outgoing,

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

We work in the five-flavour scheme, where the bottom quark is taken as
massless while its Yukawa coupling to the Higgs boson is kept finite, so that

$\displaystyle p_{1}^{2}=p_{2}^{2}=p_{3}^{2}=p_{4}^{2}=0\,,\qquad\quad
p_{5}^{2}=m_{H}^{2}\,,$ (5)

where $m_{H}$ is the mass of the Higgs boson. The kinematics is described by
six independent scalar products, which we choose as

$(s_{12},s_{23},s_{34},s_{45},s_{15},p_{5}^{2})\,,$

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

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

The latter is connected to the Gram determinant of the four linearly
independent momenta, $\Delta=\mathrm{det}(s_{ij})$ with $1\leq i,j\leq 4$, via
$\Delta=(\text{tr}_{5})^{2}$.

The colour decomposition of the $L$-loop amplitudes in the leading colour
approximation is given by

$\displaystyle\mathcal{A}^{(L)}(1_{\bar{b}},2_{b},3_{g},4_{g},5_{H})=n^{L}g_{s}^{2}y_{b}\bigg{[}(T^{a_{3}}T^{a_{4}})_{i_{2}}^{\;\;\bar{i}_{1}}\,A^{(L)}(1_{\bar{b}},2_{b},3_{g},4_{g},5_{H})+(3\leftrightarrow
4)\bigg{]}\,,$
$\displaystyle\mathcal{A}^{(L)}(1_{\bar{b}},2_{b},3_{\bar{q}},4_{q},5_{H})=n^{L}g_{s}^{2}y_{b}\delta_{i_{4}}^{\;\;\bar{i}_{1}}\delta_{i_{2}}^{\;\;\bar{i}_{3}}\,A^{(L)}(1_{\bar{b}},2_{b},3_{\bar{q}},4_{q},5_{H})\,,$
(7)
$\displaystyle\mathcal{A}^{(L)}(1_{\bar{b}},2_{b},3_{\bar{b}},4_{b},5_{H})=n^{L}g_{s}^{2}y_{b}\bigg{[}\delta_{i_{4}}^{\;\;\bar{i}_{1}}\delta_{i_{2}}^{\;\;\bar{i}_{3}}\,\biggl{(}A^{(L)}(1_{\bar{b}},2_{b},3_{\bar{q}},4_{q},5_{H})+A^{(L)}(3_{\bar{b}},4_{b},1_{\bar{q}},2_{q},5_{H})\biggr{)}\,$
$\displaystyle\hskip
128.0374pt-\delta_{i_{2}}^{\;\;\bar{i}_{1}}\delta_{i_{4}}^{\;\;\bar{i}_{3}}\,\biggl{(}A^{(L)}(1_{\bar{b}},4_{b},3_{\bar{q}},2_{q},5_{H})+A^{(L)}(3_{\bar{b}},2_{b},1_{\bar{q}},4_{q},5_{H})\biggr{)}\bigg{]}\,,$

where $n=m_{\epsilon}\alpha_{s}/(4\pi),\ \alpha_{s}=g_{s}^{2}/(4\pi)$,
$m_{\epsilon}=i(4\pi)^{\epsilon}e^{-\epsilon\gamma_{E}}$, $T^{a}$ are the
fundamental generators of $SU(N_{c})$ normalised such that $\operatorname{\rm
tr}(T^{a}T^{b})=\delta^{ab}$, while $g_{s}$ and $y_{b}$ are the strong
coupling constant and bottom-quark Yukawa coupling, respectively. We further
decompose the partial amplitudes at one and two loops based on their closed
fermion loop contributions, as

$\displaystyle A^{(1)}=N_{c}A^{(1),1}+n_{f}A^{(1),n_{f}}\,,$ (8)
$\displaystyle
A^{(2)}=N_{c}^{2}A^{(2),1}+N_{c}n_{f}A^{(2),n_{f}}+n_{f}^{2}A^{(2),n_{f}^{2}}\,,$
(9)

where $n_{f}$ is the number of light quarks circulating in the loop. The
Feynman diagrams with the Higgs boson directly coupled to a closed bottom-
quark loop vanish since, for a massless bottom-quark, they contain a Dirac
trace with an odd number of $\gamma$ matrices. In our computation we do not
consider the closed top-quark loop contribution.

Figure 1: Sample Feynman diagrams corresponding to the various closed fermion
loop contributions at one loop as specified in Eq. (8). Black-dashed, red,
black-spiralled and black lines represent Higgs bosons, bottom quarks, gluons
and light quarks (bottom quarks included), respectively. Figure 2: Sample
Feynman diagrams corresponding to the various closed fermion loop
contributions at two loops as specified in Eq. (9). Black-dashed, red, black-
spiralled and black lines represent Higgs bosons, bottom quarks, gluons and
light quarks (bottom quarks included), respectively.

The pole structure of the unrenormalised amplitudes in the ’t Hooft-Veltman
(tHV) scheme at one and two loops is given by

$\displaystyle P^{(1)}$
$\displaystyle=2I_{1}(\epsilon)+\frac{\beta_{0}}{\epsilon}-2s_{1},$ (10)
$\displaystyle P^{(2)}$
$\displaystyle=2I_{1}(\epsilon)\bigg{(}\hat{A}^{(1)}-\frac{\beta_{0}}{\epsilon}+2s_{1}\bigg{)}+4I_{2}(\epsilon)+\bigg{(}\frac{2\beta_{0}}{\epsilon}-2s_{1}\bigg{)}\hat{A}^{(1)}$
$\displaystyle\hskip
11.38092pt-\frac{\beta_{0}^{2}}{\epsilon^{2}}+\frac{\beta_{1}}{2\epsilon}-4s_{2}+\frac{2s_{1}\beta_{0}}{\epsilon}\,,$
(11)

where $\hat{A}^{(1)}$ is the unrenormalised one-loop amplitude normalised to
the tree-level amplitude. $s_{1}$ and $s_{2}$ are the bottom-quark Yukawa
renormalisation constants, and their expressions can be found in Appendix A.
We used a mixed renormalisation scheme where the strong coupling $\alpha_{s}$
and the bottom-Yukawa coupling $y_{b}$ are renormalised in the
$\overline{\text{MS}}$ scheme, while the bottom-quark mass and wave function
are renormalised in the on-shell ($\mathrm{OS}$) scheme. This allows to keep
$y_{b}$ finite while taking the bottom-quark mass smoothly to zero
($m_{b}^{\mathrm{OS}}\rightarrow 0$) Mondini:2021nck . Such a mixed
renormalisation scheme can be used so long as pure QCD corrections are
considered. In fact, using the $\overline{\text{MS}}$ scheme to renormalise
$y_{b}$ allows us to better control the convergence of the perturbative
corrections by resumming the large logarithms that appear in the OS scheme by
running $y_{b}$ to a scale close to the Higgs mass. In the presence of
electroweak (EW) corrections, however, the relationship between $y_{b}$ and
$m_{b}$ must be imposed to guarantee the cancellation of UV singularities
Pagani:2020rsg .

The $I_{2}(\epsilon)$ operator is defined by

$I_{2}(\epsilon)=-\frac{1}{2}I_{1}(\epsilon)\left[I_{1}(\epsilon)+\frac{\beta_{0}}{\epsilon}\right]+\frac{N(\epsilon)}{N(2\epsilon)}\left[\frac{\beta_{0}}{2\epsilon}+\frac{\gamma_{1}^{{\mathrm{cusp}}}}{8}\right]I_{1}(2\epsilon)+H^{(2)}(\epsilon)\,,$
(12)

while the $I_{1}(\epsilon)$ operators for $b\bar{b}H$ production in both the
$gg$ and the $q\bar{q}$ channels are given at leading colour by

$\displaystyle I^{\bar{b}b\bar{q}qH}_{1}(\epsilon)$
$\displaystyle=-N_{c}\frac{N(\epsilon)}{2}\bigg{(}\frac{1}{\epsilon^{2}}+\frac{3}{2\epsilon}\bigg{)}\big{[}\left(-s_{23}\right)^{-\epsilon}+\left(-s_{14}\right)^{-\epsilon}\big{]},$
(13) $\displaystyle I^{\bar{b}bggH}_{1}(\epsilon)$
$\displaystyle=-N_{c}\frac{N(\epsilon)}{2}\bigg{\\{}\bigg{(}\frac{1}{\epsilon^{2}}+\frac{3}{4\epsilon}+\frac{\beta_{0}}{4\epsilon}\bigg{)}\big{[}\left(-s_{23}\right)^{-\epsilon}+\left(-s_{14}\right)^{-\epsilon}\big{]}+\bigg{(}\frac{1}{\epsilon^{2}}+\frac{\beta_{0}}{2\epsilon}\bigg{)}\left(-s_{34}\right)^{-\epsilon}\bigg{\\}}\,,$
(14)

where $N(\epsilon)={e^{\epsilon\gamma_{E}}}/{\Gamma(1-\epsilon)}$,
$s_{14}=s_{23}-s_{15}-s_{45}+p_{5}^{2}$ and

$\displaystyle H_{\bar{b}b\bar{q}qH}^{(2)}(\epsilon)$
$\displaystyle=\frac{1}{16\epsilon}\bigg{\\{}4\gamma_{1}^{q}-\gamma_{1}^{\mathrm{cusp}}\gamma_{0}^{q}+\frac{\pi^{2}}{4}\beta_{0}\gamma_{0}^{\mathrm{cusp}}C_{F}\bigg{\\}}\,,$
(15) $\displaystyle H_{\bar{b}bggH}^{(2)}(\epsilon)$
$\displaystyle=\frac{1}{16\epsilon}\bigg{\\{}2\left(\gamma_{1}^{q}+\gamma_{1}^{g}\right)-\frac{1}{2}\gamma_{1}^{\mathrm{cusp}}\left(\gamma_{0}^{q}+\gamma_{0}^{g}\right)+\frac{\pi^{2}}{8}\beta_{0}\gamma_{0}^{\mathrm{cusp}}\left(C_{F}+C_{A}\right)\bigg{\\}}\,.$
(16)

The $\beta$ function coefficients and anomalous dimensions are given in
Appendix A. The finite remainder of the $L$-loop partial amplitude is then
obtained by subtracting the poles $P^{(L)}$ (which include both the
ultraviolet and infrared singularities) from the unrenormalised partial
amplitude $A^{(L)}$ and setting $\epsilon$ to $0$,

$F^{(L)}=\lim_{\epsilon\to 0}\left[A^{(L)}-P^{(L)}A^{(0)}\right]\,.$ (17)

The partial finite remainders $F^{(L)}$ inherit from the partial amplitudes
the decomposition in powers of $n_{f}$,

$\displaystyle F^{(1)}=N_{c}F^{(1),1}+n_{f}F^{(1),n_{f}}\,,$ (18)
$\displaystyle
F^{(2)}=N_{c}^{2}F^{(2),1}+N_{c}n_{f}F^{(2),n_{f}}+n_{f}^{2}F^{(2),n_{f}^{2}}\,.$
(19)

The full finite remainders $\mathcal{F}^{(L)}$ are obtained from the partial
ones $F^{(L)}$ through a colour decomposition analogous to that given in Eq.
(7) for the bare amplitudes,

$\displaystyle\mathcal{F}^{(L)}(1_{\bar{b}},2_{b},3_{g},4_{g},5_{H})=n^{L}g_{s}^{2}y_{b}\bigg{[}(T^{a_{3}}T^{a_{4}})_{i_{2}}^{\;\;\bar{i}_{1}}\,F^{(L)}(1_{\bar{b}},2_{b},3_{g},4_{g},5_{H})+(3\leftrightarrow
4)\bigg{]}\,,$
$\displaystyle\mathcal{F}^{(L)}(1_{\bar{b}},2_{b},3_{\bar{q}},4_{q},5_{H})=n^{L}g_{s}^{2}y_{b}\delta_{i_{4}}^{\;\;\bar{i}_{1}}\delta_{i_{2}}^{\;\;\bar{i}_{3}}\,F^{(L)}(1_{\bar{b}},2_{b},3_{\bar{q}},4_{q},5_{H})\,,$
(20)
$\displaystyle\mathcal{F}^{(L)}(1_{\bar{b}},2_{b},3_{\bar{b}},4_{b},5_{H})=n^{L}g_{s}^{2}y_{b}\bigg{[}\delta_{i_{4}}^{\;\;\bar{i}_{1}}\delta_{i_{2}}^{\;\;\bar{i}_{3}}\,\biggl{(}F^{(L)}(1_{\bar{b}},2_{b},3_{\bar{q}},4_{q},5_{H})+F^{(L)}(3_{\bar{b}},4_{b},1_{\bar{q}},2_{q},5_{H})\biggr{)}\,$
$\displaystyle\hskip
128.0374pt-\delta_{i_{2}}^{\;\;\bar{i}_{1}}\delta_{i_{4}}^{\;\;\bar{i}_{3}}\,\biggl{(}F^{(L)}(1_{\bar{b}},4_{b},3_{\bar{q}},2_{q},5_{H})+F^{(L)}(3_{\bar{b}},2_{b},1_{\bar{q}},4_{q},5_{H})\biggr{)}\bigg{]}\,.$

###  2.1 Tree-Level Amplitudes

The tree-level amplitudes can be obtained using the BCFW recursion relations
Britto:2004ap ; Britto:2005fq  within the spinor helicity formalism. In the
$\bar{b}bggH$ case, we choose to shift the momenta of gluons $3$ and $4$,
while in the $\bar{b}b\bar{q}qH$ case we choose particles $1$ and $4$ to avoid
shifting the momenta of adjacent quarks of the same flavour. Moreover, we
ensure that the shifted brackets $[\hat{i}\rangle,|\hat{j}]$ do not belong to
particles of helicities $i^{-},j^{+}$. These choices are necessary for the
validity of the recursion relations as they prevent the shifted amplitude from
having poles at infinity.

For the $\bar{b}bggH$ channel we obtain the following non-vanishing tree-level
partial amplitudes,

$\displaystyle A^{(0)}(1^{+}_{\bar{b}},2^{+}_{b},3^{+}_{g},4^{+}_{g},5_{H})$
$\displaystyle=\frac{s_{5}}{\langle 23\rangle\langle 34\rangle\langle
41\rangle}\,,$ (21) $\displaystyle
A^{(0)}(1^{+}_{\bar{b}},2^{+}_{b},3^{-}_{g},4^{-}_{g},5_{H})$
$\displaystyle=-\frac{[12]^{2}}{[23][34][41]}\,,$ (22) $\displaystyle
A^{(0)}(1^{+}_{\bar{b}},2^{+}_{b},3^{+}_{g},4^{-}_{g},5_{H})$
$\displaystyle=\frac{\langle 24\rangle\langle 4|5|1]^{2}}{s_{234}\langle
23\rangle\langle 34\rangle\langle
2|5|1]}-\frac{s_{5}[13]^{3}}{s_{134}[14][34]\langle 2|5|1]}\,.$ (23)

The $\bar{b},b$ quarks need to have the same helicity as that is the only way
they can couple to the Higgs boson. For the $\bar{b}bggH$ channel the “all-
plus” and MHV configurations vanish, and we are left with

$\displaystyle
A^{(0)}(1^{+}_{\bar{b}},2^{+}_{b},3^{+}_{\bar{q}},4^{-}_{q},5_{H})$
$\displaystyle=\frac{\langle 4|5|1]^{2}}{s_{234}\langle 34\rangle\langle
2|5|1]}+\frac{s_{5}[31]^{2}}{s_{134}[34]\langle 2|5|1]}\,.$ (24)

In both cases, due to the colour decomposition of the full amplitudes given by
Eq. (7), the $A^{(0)}(1^{+},2^{+},3^{-},4^{+},5_{H})$ partial amplitude is
related to $A^{(0)}(1^{+},2^{+},3^{+},4^{-},5_{H})$ by swapping the particles
$1\leftrightarrow 2,3\leftrightarrow 4$ , and flipping the overall sign for
the subprocess $\bar{b}bggH$. The remaining non-vanishing helicity
configurations can be obtained by parity transformations, that is by swapping
the brackets $\langle\phantom{1}\rangle\leftrightarrow[\phantom{1}]$.

##  3 Amplitude Reduction and Reconstruction

We derive the analytic form of the two-loop QCD helicity amplitudes using the
framework discussed in Refs. Badger:2021owl ; Badger:2021imn , which combines
the Feynman diagram approach with numerical sampling over finite fields and
functional reconstruction techniques. We generate a set of Feynman diagrams
contributing to both the $\bar{b}bggH$ and $\bar{b}b\bar{q}qH$ processes using
QGRAF Nogueira:1991ex , and perform diagram filtering, topology identification
and colour decomposition using a collection of in-house Mathematica and Form
Kuipers:2012rf ; Ruijl:2017dtg  scripts. In the leading colour approximation
there are 749 (264) Feynman diagrams contributing to the two-loop
$\bar{b}bggH$ ($\bar{b}b\bar{q}qH$) QCD amplitudes, including all closed
fermion loop contributions. The numerators of the loop amplitudes are then
computed for each independent helicity configuration, and the spinor-helicity
manipulations are carried out using the library Spinney Cullen:2010jv . The
helicity dependent loop numerators are written in terms of spinor products
($\langle ij\rangle$, $[ij]$), scalar products ($p_{i}\cdot p_{j}$,
$k_{i}\cdot k_{j}$, $k_{i}\cdot p_{j}$) and spinor strings ($\langle
i|k_{i}|j]$, $\langle i|p_{5}|j]$, $\langle i|k_{i}p_{5}|j\rangle$,
$[i|k_{i}p_{5}|j]$), where $p_{i}$’s and $k_{i}$’s are the external and loop
momenta, respectively. The angle $|i\rangle$ and square $|i]$ brackets denote
the usual holomorphic and anti-holomorphic bi-spinors built from the external
massless momenta. The momenta appearing between them in the spinor strings are
intended as the matrices obtained by contracting the corresponding four-
momenta with the vector of the Pauli matrices, e.g. $\langle
i|p_{5}|j]=\langle i|\sigma_{\mu}|j]\,p^{\mu}_{5}$. The loop numerators are
therefore expressed as linear combinations of monomials of loop-momentum
dependent objects that multiply coefficients which depend only on the external
kinematics. These coefficients contain spinor products and loop momentum
independent spinor strings, and such objects not only are not independent (due
to momentum conservation and Schouten identity), but they are also
incompatible with the finite field arithmetic framework.

In order to allow for the use of finite field sampling in the computation of
the helicity amplitudes, we make use of an explicit rational parametrisation
of the external kinematics. This is constructed with help of the momentum
twistor Hodges:2009hk  and spinor-helicity formalisms. To obtain a
configuration for the off-shell five-particle configuration we begin with a
massless configuration for six particles. While the exact parametrisation is
not important, the form presented in Ref. Badger:2016uuq  is a useful starting
point. We can think of the massless process as the result of the off-shell leg
decaying into a pair of massless particles. There are 8 independent variables
in the six-particle process, but we can reduce the degrees of freedom by
choosing one of the decay direction to be collinear with one of the other
massless legs in the five-point kinematics. Since the momentum twistors are
associated with complex momenta we consider the angle and square bracket
spinors products to be independent. We can write this explicitly as follows:
the six massless momenta $q_{i}^{\mu}$ are related to the off-shell five
particle momenta $p_{i}^{\mu}$ by

$\displaystyle p_{1}^{\mu}=q_{1}^{\mu}\,,\ \qquad p_{2}^{\mu}=q_{2}^{\mu}\,,\
\qquad p_{3}^{\mu}=q_{3}^{\mu}\,,\ \qquad p_{4}^{\mu}=q_{4}^{\mu}\,,\ \qquad
p_{5}^{\mu}=q_{5}^{\mu}+q_{6}^{\mu}\,.$ (25)

We impose the constraints $\langle q_{2}q_{6}\rangle=0$ and $[q_{2}q_{6}]=0$
to reduce the independent variables to 6. We then apply a transformation onto
the following choice of variables:

$\displaystyle x_{1}$ $\displaystyle=s_{12},$ (26) $\displaystyle x_{2}$
$\displaystyle=-\frac{\operatorname{\rm tr}_{+}(1234)}{s_{12}s_{34}},$
$\displaystyle x_{3}$ $\displaystyle=\frac{\operatorname{\rm
tr}_{+}(134152)}{s_{13}\operatorname{\rm tr}_{+}{(1452)}},$ $\displaystyle
x_{4}$ $\displaystyle=\frac{s_{23}}{s_{12}},$ $\displaystyle x_{5}$
$\displaystyle=-\frac{\operatorname{\rm
tr}_{-}(1(2+3)(1+5)523)}{s_{23}\operatorname{\rm tr}_{-}(1523)},$
$\displaystyle x_{6}$ $\displaystyle=\frac{s_{45}}{s_{12}},$

where $\operatorname{\rm tr}_{\pm}(ij\cdots kl)=\frac{1}{2}\operatorname{\rm
tr}((1\pm\gamma_{5})\not{p}_{i}\not{p}_{j}\cdots\not{p}_{k}\not{p}_{l})$. Some
features of this parametrisation are particularly convenient. For example, the
only dimensionful quantity is $x_{1}$, and $\operatorname{\rm tr}_{5}$ is
rational. However, such a choice does break some symmetries in the problem and
so different helicity configurations will have different complexities. For the
processes considered here the polynomial complexity was manageable using this
form. Applications to other off-shell five-particle processes may require
further thought.

Having set up a rational parametrisation of the external kinematics, the
helicity dependent loop numerators are constructed analytically and ready to
be further processed. This is the starting point of our numerical algorithm in
the finite field setup. In order to write the loop amplitude in terms of
Feynman scalar integrals, we first define the integral families for the
maximal topologies, which are the topologies with the maximum number of loop
propagators (8 in this case). There are 15 maximal topologies for both the
$\bar{b}bggH$ and $\bar{b}b\bar{q}qH$ processes, and they are shown in Fig. 3
of Ref. Badger:2021nhg . Each diagram topology appearing in the amplitude can
be mapped onto at least one of the maximal topologies. After that is done, the
objects which depend on the loop momenta in the helicity dependent loop
numerators are written in terms of the 11 propagators associated with the
chosen maximal topology. These mapping procedures are performed numerically
within the FiniteFlow framework Peraro:2019svx . External kinematics that
involve in the mapping are already expressed in terms of momentum twistor
variables. At this stage, the loop amplitude is expressed as a linear
combination of scalar integrals that will further be reduced to a set of
master integrals. The coefficients associated with these scalar integrals are
functions only of the external kinematics in the form of momentum twistor
variables.

The integrand is now ready to be reduced to a set of master integrals using
the Integration-by-Parts (IBP) reduction method Chetyrkin:1981qh  within the
finite field setup. The IBP relations are generated using LiteRed Lee:2012cn
and solved numerically over finite fields using the Laporta algorithm
Laporta:2001dd . The IBP reduction is performed directly to a set of master
integrals with uniform transcendental (UT) weight identified in Ref.
Abreu:2020jxa . The UT master integrals are further decomposed into a basis of
special functions $f$ that is built out of the UT master integral components
as proposed in Ref. Badger:2021nhg . We refer to this basis of special
functions as the master integral function basis. Subtracting the UV and IR
poles from the bare helicity amplitudes and performing Laurent expansion in
the dimensional regulator $\epsilon=(4-d)/2$, we arrive at the following
expression for the finite remainders,

$F^{(L)}=\sum_{i}r_{i}(x)m_{i}(f)\,,$ (27)

where $m_{i}(f)$ are monomials of the special functions and $r_{i}$ are
rational functions of the momentum twistor variables $x$. It is important to
note that the definition of the UT master integrals involve three square
roots. One is related to the pseudo-scalar invariant $\operatorname{\rm
tr}_{5}$, which captures the parity-odd part of the spinor expressions and is
therefore present already in the rational coefficients. This square root is
rationalised by the momentum-twistor inspired parameterisation (26). The other
two square roots are associated with the three-mass triangle Gram
determinants,

$\displaystyle\Delta_{1}$ $\displaystyle=2\sqrt{(p_{23}\cdot
p_{5})^{2}-p_{5}^{2}s_{23}}\,,$ $\displaystyle\Delta_{2}$
$\displaystyle=2\sqrt{(p_{12}\cdot p_{5})^{2}-p_{5}^{2}s_{12}}\,,$ (28)

where $p_{ij}=p_{i}+p_{j}$. They are not rationalised by our parameterisation
and may therefore be problematic for the finite field setup. We overcome this
issue by absorbing the three square roots in the definition of the UT master
integrals, which is possible because they appear only as overall
normalisations of the latter. As a result, they are contained in the monomials
$m_{i}(f)$ in Eq. (27) and do not appear in the amplitude reconstruction. We
note that in our computation, while the $\operatorname{\rm tr}_{5}$
originating from the UT master integral definition is absorbed in $m_{i}(f)$,
the $\operatorname{\rm tr}_{5}$ already present in the rational coefficients
is rationalised and present in the amplitude reconstruction.

At this stage we have an algorithm which computes the coefficients of the
special function monomials, $r_{i}(x)$ in Eq. (28), numerically over finite
fields for each of the independent helicity configurations of the two
processes $\bar{b}bggH$ and $\bar{b}b\bar{q}qH$. We emphasise that from the
start of our numerical algorithm until the evaluation of $r_{i}(x)$, the
computation is done within one system of FiniteFlow graphs. The last step
which remains to be done is the functional reconstruction of the rational
coefficients. This task is made challenging by the high polynomial degrees,
which are shown in the third column of Tables 1 and 2. We tackle the
complexity of the reconstruction through the strategy already used in Refs.
Badger:2021nhg ; Badger:2021imn . We refer to Ref. Badger:2021imn  for a
thorough discussion, and give here only a brief outline.

The first step of the strategy consists in fitting the linear relations among
the rational coefficients. These linear relations are then used to express the
rational coefficients in terms of a set of linearly independent coefficients,
which are chosen to have the lowest possible degrees. The degrees of the
independent rational coefficients are given in the fourth column of Tables 1
and 2.

The second step of our reconstruction strategy consists in determining the
factors in the denominators of the rational coefficients. The analytic
structure of the special functions is determined by a set of algebraic
functions of the kinematics called letters. In other words, the singularities
and branch cuts of the special functions are located on the hypersurfaces
where any of the letters vanishes (or goes to infinity). It is therefore
natural to expect that the rational coefficients which multiply the special
functions should feature poles which are similarly linked to the letters,
given for this problem in Ref. Abreu:2020jxa . Inspired from the letters we
can therefore guess the factors in the denominators of the rational
coefficients. Our guess is given by

$\displaystyle\biggl{\\{}$ $\displaystyle\langle ij\rangle\,,[ij]\,,\langle
i|p_{5}|j]\,,s_{ij}\,,s_{ij}-s_{kl}\,,s_{i5}-p_{5}^{2}\,,p_{5}^{2}\,,\operatorname{\rm
tr}_{5}\,,\Delta_{1}\,,\Delta_{2}\,,$ (29) $\displaystyle
s_{15}(s_{13}+s_{34})-p_{5}^{2}s_{34}\,,s_{25}(s_{24}+s_{34})-p_{5}^{2}s_{34}\biggr{\\}}\,,$

where the indices $i,j$ in the spinor expressions vary from $1$ to $4$, while
the indices $i,j,k,l$ in the Mandelstam invariants vary from $1$ to $5$. The
various free indices in each of the factors are understood to be different
from each other. We then match an ansatz made of the factors in Eq. (29)
against the rational coefficients reconstructed on a random univariate slice
of the phase space. This allows to determine entirely the denominators, as
well as some factors in the numerators.

Having determined the denominators of the rational coefficients we can proceed
to the third and last step of our reconstruction strategy: a univariate
partial fraction decomposition on the fly. We find that, with the
parameterisation of the external kinematics given by Eq. (26), partial
fractioning with respect to $x_{5}$ is the most convenient choice. The partial
fraction decomposition of multivariate functions has recently drawn a lot
interest as a powerful tool to simplify the analytic expressions of scattering
amplitudes. A number of new approaches have been proposed, which make use of
algebraic geometry techniques to handle the multivariate case efficiently
Leinartas:1978 ; Raichev:2012 ; Abreu:2019odu ; Boehm:2020ijp ; Heller:2021qkz
. These algorithms are however applied to known analytic expressions, namely
after the rational reconstruction has been performed. Conversely, the
univariate partial fraction decomposition can be performed within the finite
field setup to simplify the reconstruction of the rational coefficients. The
algorithm makes use of the knowledge of the denominators to construct an
ansatz for the decomposition into partial fractions with respect to the chosen
variable. The coefficients of this ansatz can then be solved for and
reconstructed within the finite field workflow. In order to reconstruct the
coefficients of the partial fraction decomposition we perform a further
matching of the factors from the ansatz (29). The determined factors are then
removed, and the remaining functions are reconstructed with FiniteFlow’s
reconstruction algorithm. The degrees of the rational coefficients which
remain to be reconstructed are shown in the fifth column of Tables 1 and 2.
Note that after partial fractioning the denominators are not entirely
determined from the ansatz (29). The univariate partial fraction decomposition
in fact introduces spurious factors in the denominators. The latter could be
determined as well, but we refrain from doing so as it does not reduce the
complexity of the reconstruction.

$\bar{b}bggH$ |  helicity configurations |  $r_{i}(x)$ |  independent $r_{i}(x)$ |  partial fraction in $x_{5}$ |  number of points
---|---|---|---|---|---
$F^{(2),1}$ |  $++++$ |  63/57 |  52/46 |  20/6 |  3361
|  $+++-$ |  135/134 |  119/120 |  28/22 |  24901
|  $++--$ |  105/111 |  105/111 |  22/12 |  4797
$F^{(2),n_{f}}$ |  $++++$ |  45/41 |  45/41 |  16/6 |  1381
|  $+++-$ |  94/95 |  94/95 |  17/6 |  1853
|  $++--$ |  89/95 |  62/69 |  18/3 |  2492
$F^{(2),n_{f}^{2}}$ |  $++++$ |  12/8 |  9/7 |  0/0 |  3
|  $+++-$ |  11/16 |  11/16 |  3/0 |  22
|  $++--$ |  12/20 |  8/16 |  8/0 |  242
Table 1:  Maximum numerator/denominator polynomial degrees of the finite remainder coefficients $r_{i}(x)$ in Eq. (27) at each stage of our reconstruction steps, together with the number of sample points needed for the analytic reconstruction in the $\bar{b}bggH$ subprocess, for the various closed fermion loop contributions. $\bar{b}b\bar{q}qH$ |  helicity configurations |  $r_{i}(x)$ |  independent $r_{i}(x)$ |  partial fraction in $x_{5}$ |  number of points
---|---|---|---|---|---
$F^{(2),1}$ |  $+++-$ |  82/81 |  69/70 |  24/16 |  10326
$F^{(2),n_{f}}$ |  $+++-$ |  28/30 |  25/24 |  8/6 |  379
$F^{(2),n_{f}^{2}}$ |  $+++-$ |  6/11 |  6/11 |  3/0 |  22
Table 2:  Maximum numerator/denominator polynomial degrees of the finite
remainder coefficients $r_{i}(x)$ in Eq. (27) at each stage of our
reconstruction steps, together with the number of sample points needed for the
analytic reconstruction in the $\bar{b}b\bar{q}qH$ subprocess, for the various
closed fermion loop contributions.

Following the strategy outlined above we reconstructed the partial finite
remainders for the independent mostly-plus helicity configurations of the
subprocesses $\bar{b}bggH$ and $\bar{b}b\bar{q}qH$. Reconstruction data are
shown in Tables 1 and 2. The mostly-minus helicity finite remainders can be
obtained by parity conjugation. Moreover, the finite remainders for the
helicity configuration $++-+$ can be obtained by swapping $1\leftrightarrow
2$, $3\leftrightarrow 4$ in the $+++-$ finite remainders, as discussed in
Section 2.1 for the tree-level amplitudes. This symmetry follows from the
colour structure and thus holds at any loop order. We could therefore evaluate
the $++-+$ finite remainders by permuting numerically the $+++-$ ones, as we
do in order to get all the other helicity configurations from the independent
ones. By permuting numerically we mean that we obtain the values of the
permuted finite remainders by evaluating the un-permuted ones at permuted
points. This is possible because our approach to the evaluation of the special
functions, which we discuss in Section 4, handles the analytic continuation to
any region automatically. Each numerical permutation therefore amounts to one
more evaluation for each point. The permutation that takes from $+++-$ to
$++-+$ is however peculiar, as it is covered by the basis of special functions
defined in Ref. Badger:2021nhg . We can thus reconstruct the finite remainders
for both helicity configurations directly in terms of the same basis of
special functions. This is much more convenient, as it reduces the number of
permutations which need to be carried out numerically, this way decreasing the
global evaluation time of the finite remainders. For this reason we
reconstructed the analytic expression of the $++-+$ finite remainders as well,
and include the $++-+$ configuration in the independent helicity set in the
following sections. The relation with the $+++-$ configuration constitutes a
non-trivial check of our results, which we discuss in Section 5.

##  4 A Custom Basis of Special Functions for the Finite Remainders

The one and two-loop finite remainders are expressed as combinations of
rational coefficients – functions of the momentum twistor variables (26) – and
monomials of square roots and elements of the master integral function basis
$\\{f^{(w)}_{i}\\}$. The latter were classified in Ref. Badger:2021nhg  so as
to span the cyclic permutations of the planar five-particle integrals with one
massive off-shell leg up to two loops. The function space of the finite
remainders is however simpler than that of the integrals and of the
amplitudes. This becomes particularly clear when we express the special
functions in terms of Chen’s iterated integrals Chen:1977oja  (see the notes
Brown:2013qva  for a thorough discussion of their properties). Given a set of
logarithmic integration kernels $\left\\{d\log W_{i}\right\\}$, where the
arguments $W_{i}$ are algebraic functions of the external kinematics, Chen’s
iterated integrals can be defined iteratively as

$\displaystyle\left[W_{i_{1}},\ldots,W_{i_{n}}\right]_{s_{0}}(s)=\int_{0}^{1}dt\,\frac{d\log
W_{i_{n}}\left(\gamma(t)\right)}{dt}\left[W_{i_{1}},\ldots,W_{i_{n-1}}\right]_{s_{0}}\left(\gamma(t)\right)\,,$
(30)

where we denote cumulatively by $s$ the dependence on the external kinematics,
$s_{0}$ is an arbitrary reference point, and $\gamma$ is an arbitrary contour
in the phase space going from $\gamma(0)=s_{0}$ to $\gamma(1)=s$. The
iteration starts with $[]_{s_{0}}(s):=1$, and the number of iterated
integrations – $n$ in Eq. (30) – is called transcendental weight. At two loops
up to order $\epsilon^{0}$ iterated integrals with transcendental weight up to
four are required. In our case, the integration kernels are given by the
letters $\\{W_{i}\\}$ of the alphabet defined in Ref. Abreu:2020jxa . Although
they may look rather abstract, Chen’s iterated integrals offer several
important advantages. Most importantly, they implement automatically all the
functional relations which would otherwise be hidden in a representation in
terms of other types of functions, such as Goncharov polylogarithms. This
property was exploited in Ref. Badger:2021nhg  to construct the master
integral function basis $\\{f^{(w)}_{i}\\}$. The expression in terms of Chen’s
iterated integrals of the master integral basis functions can be obtained
through their definition in terms of master integrals components given in Ref.
Badger:2021nhg , and the canonical differential equations for the master
integrals given in Ref. Abreu:2020jxa  (see Ref. Badger:2021nhg  for a
thorough discussion). Another important benefit of using iterated integrals is
that their analytic structure is beautifully manifest: an iterated integral
may have singularities or branch points only where one of its letters vanishes
or diverges. Expressing the special functions in the finite remainders in
terms of Chen’s iterated integrals therefore highlights their analytic
structure. Indeed this unveils important simplifications: certain letters,
which are present in the expressions of the master integrals, are absent in
the finite remainders. In other words, certain branch cuts of the integrals
drop out of the finite remainders. We observe the same pattern noticed for the
$Wb\bar{b}$ amplitudes in Ref. Badger:2021nhg . One letter,
$W_{49}=\text{tr}_{5}$, is present in the bare amplitudes but absent in the
finite remainders. This pattern of the pseudo-scalar invariant $\text{tr}_{5}$
has already been observed explicitly in many massless two-loop five-particle
finite remainders Badger:2018enw ; Abreu:2018zmy ; Abreu:2018aqd ;
Chicherin:2018yne ; Chicherin:2019xeg ; Abreu:2019rpt ; Abreu:2019odu ;
Badger:2019djh ; Caron-Huot:2020vlo ; Abreu:2020cwb ; Chawdhry:2020for , and
is linked to an underlying cluster algebra structure Chicherin:2020umh . In
addition, six letters, $W_{i}$ with $i=16,17,27,28,29,30$, are present in the
two-loop integrals but drop out already at the level of the bare amplitudes
(truncated at order $\epsilon^{0}$ at two loops). It would be of great
interest to uncover the physical principle underlying these simplifications.

The simplifications of the analytic structure mentioned above require the
interplay among various elements of the master integral function basis
$\\{f^{(w)}_{i}\\}$ of Ref. Badger:2021nhg . In other words, the separate
terms of the finite remainders may have spurious branch cuts which cancel out
in the sum. It is therefore convenient to construct a new, ad hoc basis of
special functions where the properties of the finite remainders are manifest.
In addition to being more elegant from the theoretical point of view, such a
basis is also much more convenient from the practical point of view.
Evaluating the master integral function basis $\\{f^{(w)}_{i}\\}$ in fact
requires handling integration kernels – letters – which eventually cancel out
in the objects we are interested in evaluating. It is desirable that these
cancellations take place analytically rather than numerically, so that the
spurious kernels are avoided altogether. For this reason we need a new basis
of special functions where these properties are manifest, and an approach to
evaluate it which by-passes the master integral function basis
$\\{f^{(w)}_{i}\\}$ and its unnecessary complexity.

We thus construct a new basis of special functions, which we label by
$\\{h^{(w)}_{i}\\}$, where $w=1,2,3,4$ denotes the transcendental weight. This
new basis contains only those linearly independent combinations of functions
$f^{(w)}_{i}$ which are actually present in the finite remainders, and whose
iterated integral expression is free of the letters $W_{i}$ with
$i=16,17,27,28,29,30,49$. We dub this basis finite remainder function basis,
as opposed to the master integral function basis $\\{f^{(w)}_{i}\\}$. This
leads to an important simplification in the expressions of the finite
remainders, which are then expressed as combinations of rational coefficients
and monomials in the $h^{(w)}_{i}$’s and the square roots. In particular, only
$23$ weight-4 functions $\\{h^{(4)}_{i}\\}_{i=1}^{23}$ are required, to be
compared with the $113$ weight-4 functions in the master integral function
basis. Since the evaluation of the weight-4 functions is the most expensive
step in the evaluation of the finite remainders, this reduction has a strong
positive impact on the total evaluation time. Note that we have also improved
the master integral basis $\\{f^{(w)}_{i}\\}$ of Ref. Badger:2021nhg  by
identifying relations among the higher weight functions and products of lower
weight ones which were originally missed. We achieved this following the
approach of Ref. Badger:2021nhg , but evaluating the boundary values with
higher accuracy.

In order to evaluate numerically the finite remainder function basis
$\\{h^{(w)}_{i}\\}$ we apply the method of the generalised power series
expansion Francesco:2019yqt . This approach has already found several
successful applications to the evaluation of master integrals from the
differential equations they satisfy Francesco:2019yqt ; Bonciani:2019jyb ;
Frellesvig:2019byn ; Abreu:2020jxa ; Becchetti:2020wof ; Bonciani:2021zzf ;
abreu2021twoloop . Here, following Ref. Badger:2021nhg , we apply it directly
to the basis of special functions. The main advantage is that having a basis
of special functions allows to subtract from the bare amplitudes the IR and UV
poles analytically, which makes it possible to reconstruct the finite
remainders directly over finite fields.

The method of the generalised power series expansion can be applied to the
finite remainder and master integral function basis because they too, like the
master integral bases they stem from, satisfy systems of differential
equations in the canonical form Henn:2013pwa . This follows from the fact that
the functions in the bases $\\{h^{(w)}_{i}\\}$ and $\\{f^{(w)}_{i}\\}$ are by
construction pure, that is they have the following two properties. First, they
have uniform transcendental weight, i.e. they are expressed by linear
combinations of Chen’s iterated integrals having the same transcendental
weight with constant rational coefficients. Explicitly, any of the functions
$h^{(w)}_{i}$, for $w>0$, is given in terms of Chen’s iterated integrals by

$\displaystyle
h^{(w)}_{i}(s)=\sum_{I=(i_{1},\ldots,i_{w})}c_{I}^{(i)}\,\left[W_{i_{1}},\ldots,W_{i_{w}}\right]_{s_{0}}(s)\,,$
(31)

for some rational constant coefficients $c_{I}^{(i)}$111In general the
iterated integrals in Eq. (31) include also transcendental constants, such as
$\pi$ and $\zeta(3)$. We neglect them here to simplify the presentation..
Second, the derivative of a pure weight-$w$ function is a uniform
transcendental function with weight $w-1$. This follows straightforwardly from
the differential of Chen’s iterated integrals,

$\displaystyle d\left[W_{i_{1}},\ldots,W_{i_{n}}\right]_{s_{0}}(s)=d\log
W_{i_{n}}(s)\,\left[W_{i_{1}},\ldots,W_{i_{n-1}}\right]_{s_{0}}(s)\,,$ (32)

where we note that the derivatives of $\log W_{i}$ are algebraic functions and
hence have transcendental weight $0$. The transcendental purity of the
functions in the basis therefore implies that the vector of all weight-$w$
functions in the finite remainder basis,
$\vec{h}^{(w)}=\bigl{(}h^{(w)}_{1},h^{(w)}_{2},\ldots\bigr{)}^{T}$, satisfies
a differential equation of the form

$\displaystyle d\vec{h}^{(w)}=\left(\sum_{j}b_{j}^{(w)}\,d\log
W_{j}\right)\,\vec{h}^{(w-1)}\,,$ (33)

for $w>1$, where $b_{j}$ are constant rational matrices. The matrices
$b_{j}^{(w)}$ can be determined easily from the iterated integral expression
of the functions through Eq. (32). In general the derivatives of the
weight-$w$ functions may involve weight-$(w-1)$ functions which are not needed
to express the finite remainders. These must be included in the basis as well,
in order to be able to write down the differential equations. So, strictly
speaking, the finite remainder function basis $\\{h^{(w)}_{i}\\}$ contains all
linearly independent functions which appear in the finite remainders and in
the derivatives of the higher-weight functions of the basis itself. The
differential equations (33) for all weights can be put together in one system
by defining the vector of all the functions in the finite remainder basis as

$\displaystyle\vec{h}=\left(\epsilon^{4}\vec{h}^{(4)}\,,\epsilon^{3}\vec{h}^{(3)}\,,\epsilon^{2}\vec{h}^{(2)}\,,\epsilon\vec{h}^{(1)}\,,1\right)^{T}\,,$
(34)

where $\epsilon$ is an auxiliary parameter with transcendental weight $-1$, so
that the vector $\vec{h}$ is pure with transcendental weight $0$. The full
basis $\vec{h}$ then satisfies a system of differential equations in the
canonical form Henn:2013pwa ,

$\displaystyle d\vec{h}=\epsilon\,\left(\sum_{i=1}^{49}a_{i}\,d\log
W_{i}\right)\,\vec{h}\,,$ (35)

where the constant matrices $a_{i}$ are strictly upper block triangular, with
the blocks given by the matrices $b_{i}^{(w)}$ in the differential equations
(33) for the various weights. This system of differential equations is however
much simpler than that for the master integrals. First, it contains only the
letters which do appear in the finite remainders, i.e. $a_{i}=0$ for
$i=16,17,27,28,29,30,49$. Second, while the differential equations for the
master integrals contain information about all the orders in $\epsilon$, the
system (35) encodes only those orders which are relevant for the finite
remainders, i.e. up to $\epsilon^{4}$. The constant matrices $a_{i}$ are in
fact nilpotent with degree $5$, i.e. $(a_{i})^{5}=0$ for any $i$, which
follows from their strictly upper triangular form.

We evaluate the finite remainder function basis $\\{h^{(w)}_{i}\\}$ by solving
the system of differential equations (35) numerically with the method of the
generalised power series expansions Francesco:2019yqt . For this purpose we
make use of the public Mathematica package DiffExp Hidding:2020ytt . In order
to compute the boundary values, i.e. the values of the functions at some base
point, we relate the finite remainder functions to master integral components
using the definition of the master integral function basis Badger:2021nhg ,
and evaluate the latter through their analytic expression in terms of
Goncharov polylogarithms 2011arXiv1105.2076G ; Remiddi:1999ew ;
2001math……3059G  provided in Refs. Papadopoulos:2015jft ; Syrrakos:2020kba ;
Canko:2020ylt . We evaluate the Goncharov polylogarithms numerically with the
C++ library GiNaC Vollinga:2004sn . We provide in ancillary files the
differential equations (35) and the values of the finite remainder functions
at six points, one for each $2\to 3$ scattering region with massless incoming
particles, with (at least) $200$-digit accuracy. Using this information, the
generalised power series expansion method allows to evaluate the finite
remainder function basis reliably anywhere in the phase space. This technique
in fact handles the analytic continuation automatically, so that also the
permutations required to evaluate the complete finite remainders starting from
the partial ones – as shown in Eqs. (7) for the amplitudes – as well as the
other helicity configurations can be performed numerically straightforwardly.
We assessed the reliability of this evaluation approach by checking the
convergence of the finite remainders close to the spurious poles of the
rational coefficients. We discuss this and other checks in the next section.

##  5 Further Validation

The finite remainders are defined by subtracting the expected UV and IR poles
from the bare amplitudes. Therefore, the fact that all the poles in $\epsilon$
cancel out, so that our expressions for the finite remainders are indeed
finite, is already a strong check of the validity of our results. We have also
checked that our amplitudes are correctly normalised by comparing the spin and
colour averaged squared tree-level amplitudes with full colour dependence
against MadGraph Alwall:2014hca . On top of that, we performed a number of
additional checks, which we discuss in the next subsections.

###  5.1 Direct computation of the squared finite remainders

We checked the helicity amplitudes derived in this paper against a squared
matrix element computation, carried out independently following the approach
used in the previous work on $u\bar{d}\to Wb\bar{b}$ Badger:2021nhg . In the
squared matrix element computation the two-loop amplitude is interfered with
the tree-level amplitude analytically. After manipulating the Dirac traces,
the loop numerators contain only scalar products ($k_{i}\cdot k_{j}$,
$k_{i}\cdot p_{j}$, $p_{i}\cdot p_{j}$), that can be mapped easily onto
propagator denominators. This way we derive an analytic form of the two-loop
squared matrix element written as a linear combination of scalar Feynman
integrals. The squared matrix element is then processed further through IBP
reduction to obtain the master integral representation, followed by mapping of
the master integrals onto the basis of special functions, subtraction of UV
and IR singularities, and finally Laurent expansion in $\epsilon$. All these
steps are carried out numerically over finite fields within the FiniteFlow
framework. The squared matrix element computation uses the Conventional
Dimensional Regularisation (CDR) scheme, where internal and external live in
$d=4-2\epsilon$ dimensions. We find complete agreement with this approach and
the helicity amplitudes in the tHV scheme at the level of the squared finite
remainders.

###  5.2 Renormalisation scale dependence

The analytic expressions of the one- and two-loop finite remainders have been
derived with the renormalisation scale $\mu$ set to $\mu=1$. The dependence of
the finite remainders on the renormalisation scale can be determined by
restoring the $\mu$ dependence of the strong and Yukawa couplings
($\alpha_{s}\rightarrow\alpha_{s}\mu^{2\epsilon}$, $y_{b}\rightarrow
y_{b}\mu^{\epsilon}$), and by replacing
$(-s_{ij})^{-\epsilon}\rightarrow(-\mu^{2}/s_{ij})^{\epsilon}$ in the $I_{1}$
operators defined in Eqs. (13) – (14), which enter the pole terms in Eqs. (10)
– (11). In order to capture the scale dependence of the finite remainders we
define the difference

$\delta F^{(L),i}(\mu^{2})=F^{(L),i}(\mu^{2})-F^{(L),i}(\mu^{2}=1)\,,$ (36)

where the dependence on the kinematic variables is understood. The difference
$\delta F^{(L),i}(\mu^{2})$ is entirely determined by the finite remainders
evaluated at $\mu^{2}=1$ – which we reconstructed analytically – and
logarithms of $\mu^{2}$. For the $\bar{b}bggH$ finite remainders it is given
by

$\displaystyle\delta F^{(1),1}(\mu^{2})$
$\displaystyle=\frac{31}{6}A^{(0)}\log(\mu^{2})\,,$ (37) $\displaystyle\delta
F^{(1),n_{f}}(\mu^{2})$ $\displaystyle=-\frac{2}{3}A^{(0)}\log(\mu^{2})\,,$
(38) $\displaystyle\delta F^{(2),1}(\mu^{2})$
$\displaystyle=\log(\mu_{R}^{2})\bigg{\\{}\bigg{(}\frac{1261}{54}-\frac{11}{72}\pi^{2}+9\zeta_{3}\bigg{)}A^{(0)}+\frac{53}{6}F^{(1),1}(1)\bigg{\\}}+\frac{1643}{72}A^{(0)}\log^{2}(\mu^{2})\,,$
(39) $\displaystyle\delta F^{(2),n_{f}}(\mu^{2})$
$\displaystyle=\log(\mu_{R}^{2})\bigg{\\{}\bigg{(}-\frac{241}{27}+\frac{\pi^{2}}{36}\bigg{)}A^{(0)}-\frac{4}{3}F^{(1),1}(1)+\frac{53}{6}F^{(1),n_{f}}(1)\bigg{\\}}$
$\displaystyle\quad-\frac{115}{18}A^{(0)}\log^{2}(\mu^{2})\,,$ (40)
$\displaystyle\delta F^{(2),n_{f}^{2}}(\mu^{2})$
$\displaystyle=\log(\mu^{2})\bigg{(}\frac{20}{27}A^{(0)}-\frac{4}{3}F^{(1),n_{f}}(1)\bigg{)}+\frac{4}{9}A^{(0)}\log^{2}(\mu^{2})\,,$
(41)

while for the $\bar{b}b\bar{q}qH$ finite remainders it is given by

$\displaystyle\delta F^{(1),1}(\mu^{2})$
$\displaystyle=\frac{31}{6}A^{(0)}\log(\mu^{2})\,,$ (42) $\displaystyle\delta
F^{(1),n_{f}}(\mu^{2})$ $\displaystyle=-\frac{2}{3}A^{(0)}\log(\mu^{2})\,,$
(43) $\displaystyle\delta F^{(2),1}(\mu^{2})$
$\displaystyle=\log(\mu^{2})\bigg{\\{}\bigg{(}\frac{5093}{216}-\frac{11}{12}\pi^{2}+14\zeta_{3}\bigg{)}A^{(0)}+\frac{53}{6}F^{(1),1}(1)\bigg{\\}}+\frac{1643}{72}A^{(0)}\log^{2}(\mu^{2})\,,$
(44) $\displaystyle\delta F^{(2),n_{f}}(\mu^{2})$
$\displaystyle=\log(\mu^{2})\bigg{\\{}\bigg{(}-\frac{329}{54}+\frac{\pi^{2}}{6}\bigg{)}A^{(0)}-\frac{4}{3}F^{(1),1}(1)+\frac{53}{6}F^{(1),n_{f}}(1)\bigg{\\}}$
$\displaystyle\quad-\frac{115}{18}A^{(0)}\log^{2}(\mu^{2})\,,$ (45)
$\displaystyle\delta F^{(2),n_{f}^{2}}(\mu^{2})$
$\displaystyle=-\frac{4}{3}F^{(1),n_{f}}(1)\log(\mu^{2})+\frac{4}{9}A^{(0)}\log^{2}(\mu^{2})\,.$
(46)

To check that our results for the finite remainders have the correct scale
dependence, we evaluate them at two kinematic points related by a rescaling by
some positive factor $a$,

$\displaystyle\vec{s}$
$\displaystyle=(s_{12},s_{23},s_{34},s_{45},s_{15},p_{5}^{2})\,,$ (47)
$\displaystyle\vec{s}^{\,{}^{\prime}}$
$\displaystyle=a\,\vec{s}=(as_{12},as_{23},as_{34},as_{45},as_{15},ap_{5}^{2})\,.$

We then verify numerically that the finite remainders satisfy the following
relation,

$\frac{F^{(L),i}(1,a\,\vec{s})+\delta
F^{(L),i}(a,a\,\vec{s})}{A^{(0)}(a\,\vec{s})}=\frac{F^{(L),i}(1,\vec{s})}{A^{(0)}(\vec{s})}\,,$
(48)

where we extended the notation of the finite remainders and tree-level
amplitudes to take into account their dependence on the phase space point
$\vec{s}$.

###  5.3 Relation between $+++-$ and $++-+$

The partial finite remainders for the single-minus helicity configurations,
$+++-$ and $++-+$, are related by a permutation of the external particles,

$\displaystyle
F^{(L)}(1^{+}_{\bar{b}},2^{+}_{b},3^{-}_{g},4^{+}_{g},5_{H})=-\,F^{(L)}(2^{+}_{\bar{b}},1^{+}_{b},4^{+}_{g},3^{-}_{g},5_{H})\,,$
(49) $\displaystyle
F^{(L)}(1^{+}_{\bar{b}},2^{+}_{b},3^{-}_{\bar{q}},4^{+}_{q},5_{H})=F^{(L)}(2^{+}_{\bar{b}},1^{+}_{b},4^{+}_{\bar{q}},3^{-}_{q},5_{H})\,.$

Permuting the special functions is however non-trivial, and may in general
require analytic continuation. The permuted special functions are defined by

$\displaystyle\left(\sigma\circ\vec{h}\right)(\\{p_{i}\\})=\vec{h}(\\{p_{\sigma(i)}\\})\,,$
(50)

where $\sigma$ denotes the permutation $(12345)\to(21435)$ of the external
particles, and $\\{p_{i}\\}$ denotes the dependence on the external momenta.
In order to check the relations (49) analytically, we need to express the
permuted functions $(\sigma\circ\vec{h})$ in terms of the ones in the standard
orientation, $\vec{h}$. To perform this operation in an algorithmic and robust
way we resort to the system of differential equations (35) satisfied by the
finite remainder function basis. The permuted functions in fact satisfy the
permuted differential equations

$\displaystyle
d(\sigma\circ\vec{h})=\epsilon\left[\sum_{i=1}^{49}a_{i}\,d\log\left(\sigma\circ
W_{i}\right)\right]\,(\sigma\circ\vec{h})\,.$ (51)

Permuting the letters $W_{i}$ is however trivial, as they only involve
algebraic functions. Since the alphabet is closed by construction under this
permutation, we obtain an explicit system of differential equations for the
permuted special functions,

$\displaystyle
d(\sigma\circ\vec{h})=\epsilon\left(\sum_{i=1}^{49}a_{i}^{\prime}\,d\log
W_{i}\right)\,(\sigma\circ\vec{h})\,.$ (52)

The latter can straightforwardly be solved in terms of Chen’s iterated
integrals (see Ref. Badger:2021nhg  for a thorough discussion). In order to be
comparable with the solution of the system (35), which defines the finite
remainder function basis $\vec{h}$, we must make sure that the same boundary
point is used when solving both systems of differential equations in terms of
iterated integrals. The boundary values can be obtained numerically with
arbitrary precision using the expressions in terms of Goncharov polylogarithms
of Refs. Papadopoulos:2015jft ; Syrrakos:2020kba ; Canko:2020ylt , as
discussed at the end of Section 4. Using the differential equations (52) we
can then express the permuted finite remainder special functions in terms of
the ones in the original orientation through Chen’s iterated integrals. Once
that is done, the right-hand sides of Eqs. (49) are written in terms of the
same special function basis as the left-hand sides. Since the rational
coefficients can be permuted analytically trivially, it is immediate to verify
that our results for the one- and two-loop finite remainders satisfy the
relations given by Eqs. (49).

###  5.4 Convergence near spurious poles

The rational coefficients in the finite remainders contain spurious poles,
namely poles which are not related to the physical singularities of the
amplitudes. When evaluating at a phase space point infinitesimally close to a
spurious pole – but at a finite distance from all physical poles – the
rational coefficients diverge, whereas the finite remainders must stay finite.
This can only occur through large numerical cancellations involving both the
rational coefficients and the special functions. Verifying this behaviour
constitutes both a check on our analytic expressions for the finite remainders
and a stress test of our evaluation approach for the special functions. We do
it as follows. First, we use the factor matching strategy discussed in Section
3 with the ansatz given by Eq. (29). This allows us to determine that the
spurious poles in the rational coefficients are associated with the following
factors in their denominators,

$\displaystyle\left\\{P_{i}\right\\}=\biggl{\\{}$ $\displaystyle\langle
i|p_{5}|j]\,,\text{tr}_{5}\,,\Delta_{1}\,,\Delta_{2}\,,s_{15}-s_{23}\,,s_{15}-s_{24}\,,s_{13}-s_{25}\,,s_{14}-s_{25}\,,s_{15}-s_{34}\,,$
(53) $\displaystyle
s_{25}-s_{34}\,,s_{14}-s_{35}\,,s_{23}-s_{45}\,,s_{15}(s_{13}+s_{34})-p_{5}^{2}s_{34}\,,s_{25}(s_{24}+s_{34})-p_{5}^{2}s_{34}\biggr{\\}}\,,$

where $i,j\in\\{1,2,3,4\\}$ with $i\neq j$. Next, for each factor in Eq. (53),
we construct a one-dimensional slice of the phase space, parametrised by a
parameter $\delta$, such that $P_{i}=\delta$. As we probe the small-$\delta$
region all the other factors in Eq. (53) and the factors associated with the
physical singularities – $(p_{i}+p_{j})^{2}$ and $p_{i}\cdot p_{j}$ – must
stay finite, i.e. they must neither diverge nor vanish. This ensures that we
target a specific spurious pole, rather than multiple at once, and that we
stay away from the physical singularities. Finally, we evaluate the finite
remainders on these slices for increasingly small values of $\delta$. We
evaluate the special functions with $64$-digit accuracy. Our analytic
expressions for the finite remainders exhibit the expected behaviour: as
$\delta$ approaches zero the rational coefficients diverge, while the finite
remainders converge to finite values.

##  6 Results

We provide the analytic expressions of the independent one- and two-loop
finite remainders in the ancillary files. We present them as combinations of
linearly independent rational coefficients and monomials of the square roots
and of the finite remainder basis functions. The rational coefficients are
expressed in terms of the momentum twistor variables defined in Eq. (26). We
evaluate the finite remainder function basis numerically by integrating the
differential equations (35) with the method of the generalised series
expansions. We also provide Mathematica scripts which illustrate how to
evaluate the finite remainders interfered with the tree-level amplitudes for
all the partonic channels contributing to the process $pp\to b\bar{b}H$, which
we label as

$\displaystyle\mathbf{gg}:\quad
g(-p_{3})+g(-p_{4})\rightarrow\bar{b}(p_{1})+b(p_{2})+H(p_{5})\,,$ (54)
$\displaystyle\mathbf{q\bar{q}}:\quad
q(-p_{3})+\bar{q}(-p_{4})\rightarrow\bar{b}(p_{1})+b(p_{2})+H(p_{5})\,,$
$\displaystyle\mathbf{\bar{q}q}:\quad\bar{q}(-p_{3})+q(-p_{4})\rightarrow\bar{b}(p_{1})+b(p_{2})+H(p_{5})\,,$
$\displaystyle\mathbf{b\bar{b}}:\quad
b(-p_{3})+\bar{b}(-p_{4})\rightarrow\bar{b}(p_{1})+b(p_{2})+H(p_{5})\,,$
$\displaystyle\mathbf{\bar{b}b}:\quad\bar{b}(-p_{3})+b(-p_{4})\rightarrow\bar{b}(p_{1})+b(p_{2})+H(p_{5})\,,$
$\displaystyle\mathbf{bb}:\quad b(-p_{3})+b(-p_{4})\rightarrow
b(p_{1})+b(p_{2})+H(p_{5})\,,$
$\displaystyle\mathbf{\bar{b}\bar{b}}:\quad\bar{b}(-p_{3})+\bar{b}(-p_{4})\rightarrow\bar{b}(p_{1})+\bar{b}(p_{2})+H(p_{5})\,.$

The interference between the finite remainders and the tree-level amplitudes
summed over colour and helicity is given at leading colour by

$\displaystyle\sum_{\text{colour}}\sum_{\text{helicity}}\mathcal{A}^{(0)*}\mathcal{F}^{(L)}=:g_{s}^{4}y_{b}^{2}n^{L}N_{c}^{\alpha}\,\mathcal{H}^{(L)}\,,$
(55)

where $\alpha=3$ for $\mathbf{gg}$ and $\alpha=2$ for all the other channels,
and the reduced squared finite remainders, $\mathcal{H}^{(L)}$, are given
explicitly for each channel by

$\displaystyle\mathcal{H}^{(L)}_{\mathbf{gg}}$
$\displaystyle=\sum_{\mathrm{helicity}}\big{[}A^{(0)}(1_{\bar{b}},2_{b},3_{g},4_{g},5_{H})\big{]}^{*}F^{(L)}(1_{\bar{b}},2_{b},3_{g},4_{g},5_{H})$
$\displaystyle\quad+\sum_{\mathrm{helicity}}\big{[}A^{(0)}(1_{\bar{b}},2_{b},4_{g},3_{g},5_{H})\big{]}^{*}F^{(L)}(1_{\bar{b}},2_{b},4_{g},3_{g},5_{H})\,,$
(56) $\displaystyle\mathcal{H}^{(L)}_{\mathbf{q\bar{q}}}$
$\displaystyle=\sum_{\mathrm{helicity}}\big{[}A^{(0)}(1_{\bar{b}},2_{b},3_{\bar{q}},4_{q},5_{H})\big{]}^{*}F^{(L)}(1_{\bar{b}},2_{b},3_{\bar{q}},4_{q},5_{H})\,,$
(57) $\displaystyle\mathcal{H}^{(L)}_{\mathbf{\bar{q}q}}$
$\displaystyle=\sum_{\mathrm{helicity}}\big{[}A^{(0)}(1_{\bar{b}},2_{b},4_{\bar{q}},3_{q},5_{H})\big{]}^{*}F^{(L)}(1_{\bar{b}},2_{b},4_{\bar{q}},3_{q},5_{H})\,,$
(58) $\displaystyle\mathcal{H}^{(L)}_{\mathbf{b\bar{b}}}$
$\displaystyle=\sum_{\mathrm{helicity}}\big{[}A^{(0)}(1_{\bar{b}},2_{b},3_{\bar{q}},4_{q},5_{H})+A^{(0)}(3_{\bar{b}},4_{b},1_{\bar{q}},2_{q},5_{H})\big{]}^{*}$
$\displaystyle\hskip
42.67912pt\times\big{[}F^{(L)}(1_{\bar{b}},2_{b},3_{\bar{q}},4_{q},5_{H})+F^{(L)}(3_{\bar{b}},4_{b},1_{\bar{q}},2_{q},5_{H})\big{]}$
$\displaystyle+\sum_{\mathrm{helicity}}\big{[}A^{(0)}(1_{\bar{b}},4_{b},3_{\bar{q}},2_{q},5_{H})+A^{(0)}(3_{\bar{b}},2_{b},1_{\bar{q}},4_{q},5_{H})\big{]}^{*}$
$\displaystyle\hskip
42.67912pt\times\big{[}F^{(L)}(1_{\bar{b}},4_{b},3_{\bar{q}},2_{q},5_{H})+F^{(L)}(3_{\bar{b}},2_{b},1_{\bar{q}},4_{q},5_{H})\big{]}\,,$
(59) $\displaystyle\mathcal{H}^{(L)}_{\mathbf{\bar{b}b}}$
$\displaystyle=\sum_{\mathrm{helicity}}\big{[}A^{(0)}(1_{\bar{b}},2_{b},4_{\bar{q}},3_{q},5_{H})+A^{(0)}(4_{\bar{b}},3_{b},1_{\bar{q}},2_{q},5_{H})\big{]}^{*}$
$\displaystyle\hskip
42.67912pt\times\big{[}F^{(L)}(1_{\bar{b}},2_{b},4_{\bar{q}},3_{q},5_{H})+F^{(L)}(4_{\bar{b}},3_{b},1_{\bar{q}},2_{q},5_{H})\big{]}$
$\displaystyle+\sum_{\mathrm{helicity}}\big{[}A^{(0)}(1_{\bar{b}},3_{b},4_{\bar{q}},2_{q},5_{H})+A^{(0)}(4_{\bar{b}},2_{b},1_{\bar{q}},3_{q},5_{H})\big{]}^{*}$
$\displaystyle\hskip
42.67912pt\times\big{[}F^{(L)}(1_{\bar{b}},3_{b},4_{\bar{q}},2_{q},5_{H})+F^{(L)}(4_{\bar{b}},2_{b},1_{\bar{q}},3_{q},5_{H})\big{]}\,,$
(60) $\displaystyle\mathcal{H}^{(L)}_{\mathbf{bb}}$
$\displaystyle=\sum_{\mathrm{helicity}}\big{[}A^{(0)}(3_{\bar{b}},1_{b},4_{\bar{q}},2_{q},5_{H})+A^{(0)}(4_{\bar{b}},2_{b},3_{\bar{q}},1_{q},5_{H})\big{]}^{*}$
$\displaystyle\hskip
42.67912pt\times\big{[}F^{(L)}(3_{\bar{b}},1_{b},4_{\bar{q}},2_{q},5_{H})+F^{(L)}(4_{\bar{b}},2_{b},3_{\bar{q}},1_{q},5_{H})\big{]}$
$\displaystyle+\sum_{\mathrm{helicity}}\big{[}A^{(0)}(3_{\bar{b}},2_{b},4_{\bar{q}},1_{q},5_{H})+A^{(0)}(4_{\bar{b}},1_{b},3_{\bar{q}},2_{q},5_{H})\big{]}^{*}$
$\displaystyle\hskip
42.67912pt\times\big{[}F^{(L)}(3_{\bar{b}},2_{b},4_{\bar{q}},1_{q},5_{H})+F^{(L)}(4_{\bar{b}},1_{b},3_{\bar{q}},2_{q},5_{H})\big{]}\,,$
(61) $\displaystyle\mathcal{H}^{(L)}_{\mathbf{\bar{b}\bar{b}}}$
$\displaystyle=\sum_{\mathrm{helicity}}\big{[}A^{(0)}(1_{\bar{b}},3_{b},2_{\bar{q}},4_{q},5_{H})+A^{(0)}(2_{\bar{b}},4_{b},1_{\bar{q}},3_{q},5_{H})\big{]}^{*}$
$\displaystyle\hskip
42.67912pt\times\big{[}F^{(L)}(1_{\bar{b}},3_{b},2_{\bar{q}},4_{q},5_{H})+F^{(L)}(2_{\bar{b}},4_{b},1_{\bar{q}},3_{q},5_{H})\big{]}$
$\displaystyle+\sum_{\mathrm{helicity}}\big{[}A^{(0)}(1_{\bar{b}},4_{b},3_{\bar{q}},2_{q},5_{H})+A^{(0)}(2_{\bar{b}},3_{b},1_{\bar{q}},4_{q},5_{H})\big{]}^{*}$
$\displaystyle\hskip
42.67912pt\times\big{[}F^{(L)}(1_{\bar{b}},4_{b},3_{\bar{q}},2_{q},5_{H})+F^{(L)}(2_{\bar{b}},3_{b},1_{\bar{q}},4_{q},5_{H})\big{]}\,.$
(62)

We evaluate the permutations of the finite remainders numerically. The
analytic continuation is performed by adding a small positive imaginary part
to each $s_{ij}$ and to $p_{5}^{2}$, which is done automatically by DiffExp.

To facilitate future checks, we provide here the benchmark evaluation at a
physical point corresponding to the $pp\to b\bar{b}H$ scattering process,

$\displaystyle s_{12}=\frac{49}{576}\,,\quad
s_{23}=-\frac{15337}{2048}\,,\quad s_{34}=\frac{63}{4}\,,\quad
s_{45}=-\frac{288491}{38912}\,,$ $\displaystyle\qquad\quad
s_{15}=\frac{455}{64}\,,\quad
p_{5}^{2}=7\,,\quad\mathrm{tr}_{5}=i\frac{49\sqrt{50998583}}{1400832}\,,$ (63)

which corresponds to the following values of the momentum twistor variables,

$\displaystyle x_{1}=\frac{49}{576}\,,\quad
x_{2}=-\frac{77405}{76608}-i\frac{\sqrt{50998583}}{76608}\,,\quad
x_{3}=-\frac{557251}{411874}+i\frac{95\sqrt{50998583}}{411874}\,,\quad$
$\displaystyle\qquad\quad x_{4}=-\frac{2817}{32}\,,\quad
x_{5}=-\frac{11290}{41629}-i\frac{2\sqrt{50998583}}{41629}\,,\quad
x_{6}=-\frac{370917}{4256}\,.$ (64)

The values of the bare two-loop amplitudes normalised by the tree-level
amplitudes, $\hat{A}^{(L)}=A^{(L)}/A^{(0)}$, for the independent mostly-plus
helicity configurations of the subprocesses $0\to\bar{b}bggH$ and
$0\to\bar{b}b\bar{q}qH$ are given in Tables 3 and 4. Table 5 shows the values
of the two-loop reduced squared finite remainders $\mathcal{H}^{(2)}$.
Analogous tables for the one-loop amplitudes are given in Appendix B.

$\bar{b}bggH$ |  helicity |  $\epsilon^{-4}$ |  $\epsilon^{-3}$ |  $\epsilon^{-2}$ |  $\epsilon^{-1}$ |  $\epsilon^{0}$
---|---|---|---|---|---|---
$\hat{A}^{(2),1}$ |  $++++$ |  $4.5$ |  $-11.9857+9.42478i$ |  $1.38005-40.6951i$ |  $37.5629+74.9878i$ |  $-160.364-16.4633i$
|  $+++-$ |  $4.5$ |  $-11.9857+9.42478i$ |  $11.3257-12.3672i$ |  $-26.8161+82.1522i$ |  $-142.327-160.925i$
|  $++-+$ |  $4.5$ |  $-11.9857+9.42478i$ |  $2.69154-41.4561i$ |  $35.9446+68.3748i$ |  $-132.233-11.7912i$
|  $++--$ |  $4.5$ |  $-11.9857+9.42478i$ |  $21.8803-71.2779i$ |  $85.0932+67.5004i$ |  $-293.742+11.2118i$
$\hat{A}^{(2),n_{f}}$ |  $++++$ |  $0$ |  $0.5$ |  $-0.177826+1.04486i$ |  $-0.769158-3.80277i$ |  $-5.39544+7.05528i$
|  $+++-$ |  $0$ |  $0.5$ |  $-0.192856+1.0472i$ |  $1.4513+2.42621i$ |  $-3.57357+44.5555i$
|  $++-+$ |  $0$ |  $0.5$ |  $-0.192856+1.0472i$ |  $-0.467396-4.03798i$ |  $-3.83854+2.69906i$
|  $++--$ |  $0$ |  $0.5$ |  $0.987798+0.631652i$ |  $3.957-5.16329i$ |  $33.7155-38.6759i$
$\hat{A}^{(2),n_{f}^{2}}$ |  $++++$ |  $0$ |  $0$ |  $0$ |  $0.00334-0.000519914i$ |  $0.00266436+0.0210796i$
|  $+++-$ |  $0$ |  $0$ |  $0$ |  $0$ |  $0$
|  $++-+$ |  $0$ |  $0$ |  $0$ |  $0$ |  $0$
|  $++--$ |  $0$ |  $0$ |  $0$ |  $0.262368-0.0923434i$ |  $0.532893+1.66516i$
Table 3:  Numerical values of the bare $\bar{b}bggH$ partial amplitudes at two loops (normalised to the tree-level amplitude) at the kinematic point in Eq. (64) for the four independent helicity configurations and the various closed fermion loops contributions. $\bar{b}b\bar{q}qH$ |  helicity |  $\epsilon^{-4}$ |  $\epsilon^{-3}$ |  $\epsilon^{-2}$ |  $\epsilon^{-1}$ |  $\epsilon^{0}$
---|---|---|---|---|---|---
$\hat{A}^{(2),1}$ |  $+++-$ |  $2$ |  $-6.81012$ |  $25.5694+17.9036i$ |  $-60.3404-6.4188i$ |  $48.2991-125.381i$
|  $++-+$ |  $2$ |  $-6.81012$ |  $22.4573+14.9001i$ |  $-60.7797+3.42105i$ |  $96.4449-180.941i$
$\hat{A}^{(2),n_{f}}$ |  $+++-$ |  $0$ |  $1.66667$ |  $-4.60863+4.18879i$ |  $13.2979-5.52188i$ |  $4.96804+95.7191i$
|  $++-+$ |  $0$ |  $1.66667$ |  $-4.60863+4.18879i$ |  $11.2232-7.52422i$ |  $-1.06892+93.2862i$
$\hat{A}^{(2),n_{f}^{2}}$ |  $+++-$ |  $0$ |  $0$ |  $0.444444$ |  $-0.969043+2.79253i$ |  $-6.91677-6.08868i$
|  $++-+$ |  $0$ |  $0$ |  $0.444444$ |  $-0.969043+2.79253i$ |  $-6.91677-6.08868i$
Table 4:  Numerical values of the bare $\bar{b}b\bar{q}qH$ partial amplitudes at two loops (normalised to the tree level amplitude) at the kinematic point in Eq. (64) for the four independent helicity configurations and the various closed fermion loops contributions.  channel |  $\mathrm{Re}\;\mathcal{H}^{(2),1}$ |  $\mathrm{Re}\;\mathcal{H}^{(2),n_{f}}$ |  $\mathrm{Re}\;\mathcal{H}^{(2),n_{f}^{2}}$
---|---|---|---
$\mathbf{gg}$ |  $156680.6267$ |  $-41215.80337$ |  $405.9379563$
$\mathbf{q\bar{q}}$ |  $0.09391314268$ |  $-0.02045942258$ |  $-0.004225713438$
$\mathbf{\bar{q}q}$ |  $0.3494872243$ |  $-0.08069122736$ |  $-0.004225713438$
$\mathbf{b\bar{b}}$ |  $48640.80398$ |  $-26530.01855$ |  $2458.442153$
$\mathbf{\bar{b}b}$ |  $-141130.5373$ |  $42183.03094$ |  $3711.445449$
$\mathbf{bb}/\mathbf{\bar{b}\bar{b}}$ |  $-53679.25708$ |  $1988.662899$ |  $894.7895467$
Table 5:  Numerical values of the two-loop reduced squared finite remainders
$\mathcal{H}^{(2)}$ defined in Eqs. (6)-(6) at the kinematic point in Eq. (64)
for the various closed fermion loops contributions and the scattering channels
specified in Eq. (54).

To demonstrate the suitability of our results for phenomenological
applications we present the evaluation of the finite remainders interfered
with the tree-level amplitudes on a univariate slice of the phase space. For
this purpose we parametrise the momenta for the scattering processes (1) in
terms of angles and energy fractions of the final state as

$\displaystyle
p_{1}=\frac{y_{1}\sqrt{s}}{2}\left(1\,,1\,,0\,,0\right)\,,\quad$ (65)
$\displaystyle
p_{2}=\frac{y_{2}\sqrt{s}}{2}\left(1\,,\cos\theta\,,-\sin\theta\sin\phi\,,-\sin\theta\cos\phi\right)\,,$
$\displaystyle p_{3}=\frac{\sqrt{s}}{2}\left(-1\,,0\,,0\,,-1\right)\,,$
$\displaystyle p_{4}=\frac{\sqrt{s}}{2}\left(-1\,,0\,,0\,,1\right)\,,$

while $p_{5}$ is fixed by momentum conservation,
$p_{5}=-p_{1}-p_{2}-p_{3}-p_{4}$. We have chosen the particle with momentum
$p_{1}$ to be produced at an angle of $\pi/2$ with respect to the beam axis.
Requiring that the Higgs boson is on-shell, $p_{5}^{2}=m_{H}^{2}$, constrains
the angle $\theta$ as

$\displaystyle\cos\theta=1+\frac{2}{y_{1}y_{2}}\left(1-y_{1}-y_{2}-\frac{m_{H}^{2}}{s}\right)\,.$
(66)

To restrict the kinematics to a one-dimensional slice we choose

$\displaystyle s=1\,,\quad\qquad\phi=\frac{1}{10}\,,\quad\qquad
y_{1}=\frac{3}{5}\,,\quad\qquad m_{H}=\frac{1}{10}\,.$ (67)

The reality of the angle $\theta$ then restricts the free parameter $y_{2}$ to
the interval $y_{2}\in[\frac{39}{100},\frac{39}{40}]$. In order to evaluate
the finite remainders for all the processes shown in Eq. (54) we need to
compute the finite remainder special functions at $16$ permutations of each
phase space point, as can be seen explicitly in Eqs.(6) - (6). We do this by
integrating the system of canonical differential equations (35) with DiffExp
using rationalised values of the invariants. For each permutation of each
point we compute the “distance” (technically speaking, the number of segments
into which the integration contour is divided by DiffExp) from the six
reference points provided in the ancillary files, and choose the closest one
as initial point for the integration. In Figure 3 we plot the values of the
one- and two-loop reduced squared finite remainders for all the processes
defined in Eq. (54), as functions of the parameter $y_{2}$.

(a) $\mathbf{gg}$

(b) $\mathbf{q\bar{q}}$

(c) $\mathbf{\bar{q}q}$

(d) $\mathbf{bb}$

(e) $\mathbf{b\bar{b}}$

(f) $\mathbf{\bar{b}b}$

Figure 3: Reduced squared finite remainders $\mathcal{H}^{(L)}$ at tree level,
one and two loops evaluated on the one-dimensional phase space slice defined
in Eq. (65), as functions of the variable $y_{2}$, for the channels defined in
Eq. (54).

We stress that the purpose of the plots in Figure 3 is merely to demonstrate
that our results for the finite remainders can be evaluated reliably in the
physical scattering region. Nothing can be inferred about the convergence of
the perturbative expansion at the cross section level. One interesting feature
which can be appreciated from Figure 3 is the appearance of a loop-induced
peak in the finite remainders for the channel $\mathbf{\bar{q}q}$. The peak is
absent at tree level for the same channel and up to two loops for
$\mathbf{q\bar{q}}$. The latter channel is related to $\mathbf{\bar{q}q}$ by
the exchange $3\leftrightarrow 4$ of the external particles. We observe that
the peak stems from the values of the finite remainder function basis, while
the rational coefficients are not enhanced in that region. In order to
pinpoint more precisely the origin of this phenomenon, we construct the
explicit analytic expressions for some of the special functions which exhibit
the peak, starting from their iterated integral expression obtained by solving
the system of canonical differential equations (35). For instance,

$\displaystyle h^{(2)}_{4}$
$\displaystyle=2\,\text{Li}_{2}\left(1-\frac{s_{15}}{p_{5}^{2}}\right)+2\,\text{Li}_{2}\left(1-\frac{s_{23}}{s_{15}}\right)-\frac{\pi^{2}}{4}+\frac{1}{2}\log^{2}\left(p_{5}^{2}\right)+\frac{1}{2}\log^{2}\left(-s_{45}\right)+2\log^{2}\left(s_{15}\right)$
(68)
$\displaystyle-\frac{1}{2}\log^{2}\left(-s_{23}\right)-2\log\left(s_{15}\right)\log\left(s_{15}-s_{23}\right)+\log^{2}\left(s_{15}-s_{23}\right)+\log\left(p_{5}^{2}\right)\log\left(-s_{23}\right)$
$\displaystyle-2\log\left(p_{5}^{2}\right)\log\left(s_{15}\right)-\log\left(p_{5}^{2}\right)\log\left(-s_{14}\right)+2\log\left(s_{15}\right)\log\left(-s_{14}\right)-\log\left(-s_{23}\right)\log\left(-s_{14}\right)$
$\displaystyle+i\pi\biggl{[}\log\left(p_{5}^{2}\right)-2\log\left(s_{15}\right)+2\log\left(s_{15}-s_{23}\right)-\log\left(-s_{23}\right)-\log\left(-s_{14}\right)\biggr{]}\,,$

which is well defined in the $s_{34}$ physical scattering region. The analytic
continuation to any other region is obtained by adding a small positive
imaginary part to each $s_{ij}$ and to $p_{5}^{2}$. We checked that the values
of $h^{(2)}_{4}$ as given by Eq. (68) (and of its permutations) agree with the
evaluation through the generalised series expansion. The function
$h^{(2)}_{4}$ exhibits no peak on the phase space sliced defined by Eqs. (65)
and (67), and indeed the finite remainder for the channel $\mathbf{q\bar{q}}$
does not exhibit such feature. The permutation $3\leftrightarrow 4$ of
$h^{(2)}_{4}$, which contributes to $\mathbf{\bar{q}q}$, is instead peaked
around $y_{2}\approx 0.5566$. Thanks to the analytic expression (68) we can
identify the source of the peak in the logarithms of $s_{24}$, which originate
from the $\log(-s_{23})$’s in Eq. (68) upon swapping $3\leftrightarrow 4$.
Indeed the tree-level amplitude for the subprocess $0\to\bar{b}b\bar{q}qH$,
given by Eq. (24), is manifestly free of $1/\langle 23\rangle$ poles. The
$0\to\bar{b}b\bar{q}qH$ diagrams with a $1/s_{23}$ pole would come with a loop
and so end up scaling as $\log(-s_{23})$. While $\log(-s_{23})$ is not
enhanced on the one-dimensional slice under consideration, its
$3\leftrightarrow 4$ permutation $\log(-s_{24})$, which contributes to
$\mathcal{H}^{(L)}_{\mathbf{\bar{q}q}}$ for $L=1,2$, is peaked at
$y_{2}\approx 0.5566$, where $s_{24}$ reaches its minimum absolute value on
the slice. The tree-level amplitudes for the channel $\mathbf{gg}$ instead do
have poles at $s_{23}=0$, which can be seen explicitly in Eqs. (21). Since
$\mathcal{H}^{(L)}_{\mathbf{gg}}$ with $L=1,2$ receive contributions from the
partial finite remainders in both the standard orientation and with the swap
$3\leftrightarrow 4$, as shown in Eq. (6), their plot in Figure 3 (a) exhibits
this peak already at tree level. The same holds for the $\mathbf{bb}$ and
$\mathbf{\bar{b}b}$ channels, as can be seen in Figures 3 (d) and (f).

Also in Figure 3, we observe divergences at $y_{2}=39/100$ for the processes
$\mathbf{b\bar{b}}$ and $\mathbf{\bar{b}b}$. This divergence is associated
with the propagator $1/s_{12}$, which can only appear in processes with two
pairs of bottom quarks. In Eqs. (6) and (6) we can see the $q\bar{q}$ fermion
pairs can appear with momenta $p_{1}$ and $p_{2}$, which is not the case for
other processes. All the other features of the plots in Figure 3 can be
similarly understood in terms of tree-level propagators going on shell.

##  7 Conclusions

In this article we have presented an analytic form for the two-loop QCD
corrections to the process $pp\to b\bar{b}H$. To the best of our knowledge
this is the first complete set of helicity amplitudes provided for a $2\to 3$
scattering process with an off-shell leg. A special function basis for the
finite remainder was identified obeying a canonical form differential
equation. The method for constructing such a basis was presented in the recent
$pp\to Wb\bar{b}$ computation Badger:2021nhg  by three of the present authors.
The finite remainders were then extracted from multiple evaluations over
finite fields using a rational phase-space and IBP reduction. We obtained
relatively compact results after determining the linear relations between the
rational coefficients and performing a univariate partial fraction
decomposition on the fly. The final expressions were evaluated using the
method of generalised series expansions as implemented in the DiffExp code
Hidding:2020ytt .

The expressions have been validated in a number of ways and we observe that
they exhibit a smooth behaviour in all scattering regions. Evaluation times
appear to be suitable for phenomenological applications, especially since we
have not tried to optimise the route through the phase-space evaluations as
has been done in other applications of the method Frellesvig:2019byn ;
Abreu:2020jxa ; Becchetti:2020wof ; Bonciani:2021zzf ; abreu2021twoloop .

The techniques presented here show promise for applications to other important
scattering processes such as $pp\to V+2j$ and $pp\to H+2j$, although the non-
planar sector is not entirely known and remains a high priority.

The analytic expressions together with Mathematica scripts to evaluate them
numerically are provided in the ancillary files accompanying this article.

###### Acknowledgements.

SZ wishes to thank Dmitry Chicherin and Vasily Sotnikov for useful
discussions. This project received funding from the European Union’s Horizon
2020 research and innovation programmes New level of theoretical precision for
LHC Run 2 and beyond (grant agreement No 683211), High precision multi-jet
dynamics at the LHC (grant agreement No 772009), and Novel structures in
scattering amplitudes (grant agreement No 725110). HBH was partially supported
by STFC consolidated HEP theory grant ST/T000694/1. SZ gratefully acknowledges
the computing resources provided by the Max Planck Institute for Physics and
by the Max Planck Computing & Data Facility.

##  Appendix A Renormalisation Constants

The one- and two-loop bottom-quark Yukawa renormalisation constants entering
the UV counterterm in Eqs. (10) and (11) are

$\displaystyle s_{1}$ $\displaystyle=-\frac{3C_{F}}{2\epsilon}\,,$ (69)
$\displaystyle s_{2}$
$\displaystyle=\frac{3}{8\epsilon^{2}}\big{(}3C_{F}^{2}+\beta_{0}C_{F}\big{)}-\frac{1}{8\epsilon}\bigg{(}\frac{3C_{F}^{2}}{2}+\frac{97}{6}C_{F}C_{A}-\frac{10}{3}C_{F}T_{F}n_{f}\bigg{)}\,,$
(70)

while the $\beta$ function coefficients and anomalous dimensions are

$\displaystyle\beta_{0}=$
$\displaystyle\;\frac{11}{3}C_{A}-\frac{4}{3}T_{F}n_{f}\,,$ (71)
$\displaystyle\beta_{1}=$
$\displaystyle\;\frac{34}{3}C_{A}^{2}-\frac{20}{3}C_{A}T_{F}n_{f}-4C_{F}T_{F}n_{f}\,,$
(72) $\displaystyle\gamma_{0}^{g}=$
$\displaystyle\;-\frac{11}{3}C_{A}+\frac{4}{3}T_{F}n_{f}\,,$ (73)
$\displaystyle\gamma_{1}^{g}=$
$\displaystyle\;C_{A}^{2}\left(-\frac{692}{27}+\frac{11\pi^{2}}{18}+2\zeta_{3}\right)+4C_{F}T_{F}n_{f}+C_{A}T_{F}n_{f}\left(\frac{256}{27}-\frac{2\pi^{2}}{9}\right)\,,$
(74) $\displaystyle\gamma_{0}^{q}=$ $\displaystyle\;-3C_{F}\,,$ (75)
$\displaystyle\gamma_{1}^{q}=$
$\displaystyle\;C_{F}^{2}\left(-\frac{3}{2}+2\pi^{2}-24\zeta_{3}\right)+C_{F}C_{A}\left(-\frac{961}{54}-\frac{11\pi^{2}}{6}+26\zeta_{3}\right)$
$\displaystyle\;+C_{F}T_{F}n_{f}\left(\frac{130}{27}+\frac{2\pi^{2}}{3}\right)\,,$
(76) $\displaystyle\gamma_{0}^{\mathrm{cusp}}=$ $\displaystyle\;4\,,$ (77)
$\displaystyle\gamma_{1}^{\mathrm{cusp}}=$
$\displaystyle\;\left(\frac{268}{9}-\frac{4\pi^{2}}{3}\right)C_{A}-\frac{80}{9}T_{F}n_{f}\,,$
(78)

where $C_{A}=N_{c}$ and $C_{F}=(N_{c}^{2}-1)/(2N_{c})$.

##  Appendix B One-Loop Results

$\bar{b}bggH$ |  helicity |  $\epsilon^{-2}$ |  $\epsilon^{-1}$ |  $\epsilon^{0}$ |  $\epsilon^{1}$ |  $\epsilon^{2}$
---|---|---|---|---|---|---
$\hat{A}^{(1),1}$ |  $++++$ |  $-3$ |  $3.07857-3.14159i$ |  $0.317351+8.42128i$ |  $-1.25257-8.56907i$ |  $25.8294-4.35648i$
|  $+++-$ |  $-3$ |  $3.07857-3.14159i$ |  $-2.99786-1.02133i$ |  $2.86487-28.7164i$ |  $30.3093-26.3373i$
|  $++-+$ |  $-3$ |  $3.07857-3.14159i$ |  $-0.119814+8.67497i$ |  $-1.43041-5.33656i$ |  $19.6373-0.110475i$
|  $++--$ |  $-3$ |  $3.07857-3.14159i$ |  $-6.51606+18.6156i$ |  $-21.7849+24.0036i$ |  $-24.6605+55.6878i$
$\hat{A}^{(1),n_{f}}$ |  $++++$ |  $0$ |  $0$ |  $-0.005010+0.000779871i$ |  $-0.00700827-0.0150298i$ |  $0.0109029-0.0163643i$
|  $+++-$ |  $0$ |  $0$ |  $0$ |  $0$ |  $0$
|  $++-+$ |  $0$ |  $0$ |  $0$ |  $0$ |  $0$
|  $++--$ |  $0$ |  $0$ |  $-0.393552+0.138515i$ |  $-0.793221-1.11035i$ |  $0.635641-1.48796i$
$\bar{b}b\bar{q}qH$ |  helicity |  $\epsilon^{-2}$ |  $\epsilon^{-1}$ |  $\epsilon^{0}$ |  $\epsilon^{1}$ |  $\epsilon^{2}$
$\hat{A}^{(1),1}$ |  $+++-$ |  $-2$ |  $2.48840$ |  $-9.99430-8.95182i$ |  $2.20899-24.3401i$ |  $4.76962-27.6604i$
|  $++-+$ |  $-2$ |  $2.48840$ |  $-8.43825-7.45006i$ |  $7.21741-24.6383i$ |  $13.6369-20.4876i$
$\hat{A}^{(1),n_{f}}$ |  $+++-$ |  $0$ |  $-0.666667$ |  $0.726782-2.09440i$ |  $2.29387+2.28325i$ |  $-2.54017+0.316127i$
|  $++-+$ |  $0$ |  $-0.666667$ |  $0.726782-2.09440i$ |  $2.29387+2.28325i$ |  $-2.54017+0.316127i$
Table 6:  Numerical values of the bare $\bar{b}bggH$ and $\bar{b}b\bar{q}qH$ partial amplitudes at one loop (normalised to the tree-level amplitude) at the kinematic point in Eq. (64) for the four independent helicity configurations and the various closed fermion loops contributions.  channel |  $\mathcal{H}^{(0)}$ |  $\mathrm{Re}\;\mathcal{H}^{(1),1}$ |  $\mathrm{Re}\;\mathcal{H}^{(1),n_{f}}$
---|---|---|---
$\mathbf{gg}$ |  $1121.375369$ |  $4905.689964$ |  $204.1069797$
$\mathbf{q\bar{q}}$ |  $0.001095232986$ |  $-0.008958148524$ |  $0.0007959961305$
$\mathbf{\bar{q}q}$ |  $0.001095232986$ |  $0.01182947634$ |  $0.0007959961305$
$\mathbf{b\bar{b}}$ |  $738.4111805$ |  $5948.275150$ |  $-2005.976183$
$\mathbf{\bar{b}b}$ |  $774.9861507$ |  $-8346.007933$ |  $-2253.325645$
$\mathbf{bb}/\mathbf{\bar{b}\bar{b}}$ |  $71.81424881$ |  $-678.1382010$ |  $-243.5040325$
Table 7:  Numerical values of the tree-level reduced squared amplitudes
$\mathcal{H}^{(0)}$ and one-loop reduced squared finite remainders
$\mathcal{H}^{(1)}$ defined in Eqs. (6)-(6) at the kinematic point in Eq. (64)
for the closed fermion loops contributions and the scattering channels
specified in Eq. (54).

We present in Table 6 the numerical values of the one-loop bare helicity
amplitudes at the kinematic point given in Eq. (64), evaluated through
$O(\epsilon^{2})$ for the different closed fermion loop contributions defined
in Eq. (8). These one-loop amplitudes are required in the computation of the
two-loop pole terms in Eqs. (10) and (11). Furthermore, Table 7 shows the
tree-level reduced squared amplitudes and one-loop reduced squared finite
remainders for the scattering channels listed in Eq. (54), at the same
kinematic point.

## References

  * (1) J. R. Andersen et al., _Les Houches 2015: Physics at TeV Colliders Standard Model Working Group Report_ , in _9th Les Houches Workshop on Physics at TeV Colliders_ , 5, 2016.  1605.04692.
  * (2) _Les Houches 2017: Physics at TeV Colliders Standard Model Working Group Report_ , 3, 2018.
  * (3) S. Amoroso et al., _Les Houches 2019: Physics at TeV Colliders: Standard Model Working Group Report_ , in _11th Les Houches Workshop on Physics at TeV Colliders: PhysTeV Les Houches_ , 3, 2020.  2003.01700.
  * (4) H. A. Chawdhry, M. Czakon, A. Mitov and R. Poncelet, _NNLO QCD corrections to three-photon production at the LHC_ , _JHEP_ 02 (2020) 057, [1911.00479].
  * (5) S. Kallweit, V. Sotnikov and M. Wiesemann, _Triphoton production at hadron colliders in NNLO QCD_ , 2010.04681.
  * (6) H. A. Chawdhry, M. Czakon, A. Mitov and R. Poncelet, _NNLO QCD corrections to diphoton production with an additional jet at the LHC_ , 2105.06940.
  * (7) M. Czakon, A. Mitov and R. Poncelet, _Tour de force in Quantum Chromodynamics: A first next-to-next-to-leading order study of three-jet production at the LHC_ , 2106.05331.
  * (8) C. G. Papadopoulos, D. Tommasini and C. Wever, _The Pentabox Master Integrals with the Simplified Differential Equations approach_ , _JHEP_ 04 (2016) 078, [1511.09404].
  * (9) T. Gehrmann, J. Henn and N. Lo Presti, _Pentagon functions for massless planar scattering amplitudes_ , _JHEP_ 10 (2018) 103, [1807.09812].
  * (10) S. Abreu, B. Page and M. Zeng, _Differential equations from unitarity cuts: nonplanar hexa-box integrals_ , _JHEP_ 01 (2019) 006, [1807.11522].
  * (11) D. Chicherin, T. Gehrmann, J. Henn, N. Lo Presti, V. Mitev and P. Wasser, _Analytic result for the nonplanar hexa-box integrals_ , _JHEP_ 03 (2019) 042, [1809.06240].
  * (12) D. Chicherin, T. Gehrmann, J. Henn, P. Wasser, Y. Zhang and S. Zoia, _All Master Integrals for Three-Jet Production at Next-to-Next-to-Leading Order_ , _Phys. Rev. Lett._ 123 (2019) 041603, [1812.11160].
  * (13) D. Chicherin and V. Sotnikov, _Pentagon Functions for Scattering of Five Massless Particles_ , _JHEP_ 12 (2020) 167, [2009.07803].
  * (14) T. Gehrmann, J. Henn and N. Lo Presti, _Analytic form of the two-loop planar five-gluon all-plus-helicity amplitude in QCD_ , _Phys. Rev. Lett._ 116 (2016) 062001, [1511.05409].
  * (15) S. Badger, C. Brønnum-Hansen, H. B. Hartanto and T. Peraro, _Analytic helicity amplitudes for two-loop five-gluon scattering: the single-minus case_ , _JHEP_ 01 (2019) 186, [1811.11699].
  * (16) S. Abreu, J. Dormans, F. Febres Cordero, H. Ita and B. Page, _Analytic Form of Planar Two-Loop Five-Gluon Scattering Amplitudes in QCD_ , _Phys. Rev. Lett._ 122 (2019) 082002, [1812.04586].
  * (17) S. Abreu, J. Dormans, F. Febres Cordero, H. Ita, B. Page and V. Sotnikov, _Analytic Form of the Planar Two-Loop Five-Parton Scattering Amplitudes in QCD_ , _JHEP_ 05 (2019) 084, [1904.00945].
  * (18) S. Badger, D. Chicherin, T. Gehrmann, G. Heinrich, J. Henn, T. Peraro et al., _Analytic form of the full two-loop five-gluon all-plus helicity amplitude_ , _Phys. Rev. Lett._ 123 (2019) 071601, [1905.03733].
  * (19) S. Abreu, B. Page, E. Pascual and V. Sotnikov, _Leading-Color Two-Loop QCD Corrections for Three-Photon Production at Hadron Colliders_ , 2010.15834.
  * (20) H. A. Chawdhry, M. Czakon, A. Mitov and R. Poncelet, _Two-loop leading-color helicity amplitudes for three-photon production at the LHC_ , 2012.13553.
  * (21) B. Agarwal, F. Buccioni, A. von Manteuffel and L. Tancredi, _Two-loop leading colour QCD corrections to $q\bar{q}\to\gamma\gamma g$ and $qg\to\gamma\gamma q$_, _JHEP_ 04 (2021) 201, [2102.01820].
  * (22) S. Abreu, F. Febres Cordero, H. Ita, B. Page and V. Sotnikov, _Leading-Color Two-Loop QCD Corrections for Three-Jet Production at Hadron Colliders_ , 2102.13609.
  * (23) H. A. Chawdhry, M. Czakon, A. Mitov and R. Poncelet, _Two-loop leading-colour QCD helicity amplitudes for two-photon plus jet production at the LHC_ , 2103.04319.
  * (24) B. Agarwal, F. Buccioni, A. von Manteuffel and L. Tancredi, _Two-loop helicity amplitudes for diphoton plus jet production in full color_ , 2105.04585.
  * (25) S. Badger, C. Brønnum-Hansen, D. Chicherin, T. Gehrmann, H. B. Hartanto, J. Henn et al., _Virtual QCD corrections to gluon-initiated diphoton plus jet production at hadron colliders_ , 2106.08664.
  * (26) S. Badger, C. Brønnum-Hansen, H. B. Hartanto and T. Peraro, _First look at two-loop five-gluon scattering in QCD_ , _Phys. Rev. Lett._ 120 (2018) 092001, [1712.02229].
  * (27) S. Abreu, F. Febres Cordero, H. Ita, B. Page and M. Zeng, _Planar Two-Loop Five-Gluon Amplitudes from Numerical Unitarity_ , _Phys. Rev. D_ 97 (2018) 116014, [1712.03946].
  * (28) S. Badger, C. Brønnum-Hansen, T. Gehrmann, H. B. Hartanto, J. Henn, N. A. Lo Presti et al., _Applications of integrand reduction to two-loop five-point scattering amplitudes in QCD_ , _PoS_ LL2018 (2018) 006, [1807.09709].
  * (29) S. Abreu, F. Febres Cordero, H. Ita, B. Page and V. Sotnikov, _Planar Two-Loop Five-Parton Amplitudes from Numerical Unitarity_ , _JHEP_ 11 (2018) 116, [1809.09067].
  * (30) A. V. Kotikov, _Differential equations method: New technique for massive Feynman diagrams calculation_ , _Phys. Lett. B_ 254 (1991) 158–164.
  * (31) Z. Bern, L. J. Dixon and D. A. Kosower, _Dimensionally regulated pentagon integrals_ , _Nucl. Phys. B_ 412 (1994) 751–816, [hep-ph/9306240].
  * (32) E. Remiddi, _Differential equations for Feynman graph amplitudes_ , _Nuovo Cim. A_ 110 (1997) 1435–1452, [hep-th/9711188].
  * (33) T. Gehrmann and E. Remiddi, _Differential equations for two loop four point functions_ , _Nucl. Phys. B_ 580 (2000) 485–518, [hep-ph/9912329].
  * (34) J. M. Henn, _Multiloop integrals in dimensional regularization made simple_ , _Phys. Rev. Lett._ 110 (2013) 251601, [1304.1806].
  * (35) J. Gluza, K. Kajda and D. A. Kosower, _Towards a Basis for Planar Two-Loop Integrals_ , _Phys. Rev. D_ 83 (2011) 045012, [1009.0472].
  * (36) R. M. Schabinger, _A New Algorithm For The Generation Of Unitarity-Compatible Integration By Parts Relations_ , _JHEP_ 01 (2012) 077, [1111.4220].
  * (37) H. Ita, _Two-loop Integrand Decomposition into Master Integrals and Surface Terms_ , _Phys. Rev. D_ 94 (2016) 116015, [1510.05626].
  * (38) K. J. Larsen and Y. Zhang, _Integration-by-parts reductions from unitarity cuts and algebraic geometry_ , _Phys. Rev. D_ 93 (2016) 041701, [1511.01071].
  * (39) J. Böhm, A. Georgoudis, K. J. Larsen, M. Schulze and Y. Zhang, _Complete sets of logarithmic vector fields for integration-by-parts identities of Feynman integrals_ , _Phys. Rev. D_ 98 (2018) 025023, [1712.09737].
  * (40) A. von Manteuffel and R. M. Schabinger, _A novel approach to integration by parts reduction_ , _Phys. Lett. B_ 744 (2015) 101–104, [1406.4513].
  * (41) T. Peraro, _Scattering amplitudes over finite fields and multivariate functional reconstruction_ , _JHEP_ 12 (2016) 030, [1608.01902].
  * (42) J. Klappert and F. Lange, _Reconstructing rational functions with FireFly_ , _Comput. Phys. Commun._ 247 (2020) 106951, [1904.00009].
  * (43) T. Peraro, _FiniteFlow: multivariate functional reconstruction using finite fields and dataflow graphs_ , _JHEP_ 07 (2019) 031, [1905.08019].
  * (44) J. Klappert, S. Y. Klein and F. Lange, _Interpolation of dense and sparse rational functions and other improvements in FireFly_ , _Comput. Phys. Commun._ 264 (2021) 107968, [2004.01463].
  * (45) J. Klappert, F. Lange, P. Maierhöfer and J. Usovitsch, _Integral reduction with Kira 2.0 and finite field methods_ , _Comput. Phys. Commun._ 266 (2021) 108024, [2008.06494].
  * (46) H. B. Hartanto, S. Badger, C. Brønnum-Hansen and T. Peraro, _A numerical evaluation of planar two-loop helicity amplitudes for a W-boson plus four partons_ , _JHEP_ 09 (2019) 119, [1906.11862].
  * (47) S. Abreu, H. Ita, F. Moriello, B. Page, W. Tschernow and M. Zeng, _Two-Loop Integrals for Planar Five-Point One-Mass Processes_ , _JHEP_ 11 (2020) 117, [2005.04195].
  * (48) D. D. Canko, C. G. Papadopoulos and N. Syrrakos, _Analytic representation of all planar two-loop five-point Master Integrals with one off-shell leg_ , _JHEP_ 01 (2021) 199, [2009.13917].
  * (49) N. Syrrakos, _Pentagon integrals to arbitrary order in the dimensional regulator_ , _JHEP_ 06 (2021) 037, [2012.10635].
  * (50) S. Badger, H. B. Hartanto and S. Zoia, _Two-Loop QCD Corrections to $Wb\bar{b}$ Production at Hadron Colliders_, _Phys. Rev. Lett._ 127 (2021) 012001, [2102.02516].
  * (51) C. G. Papadopoulos and C. Wever, _Internal Reduction method for computing Feynman Integrals_ , _JHEP_ 02 (2020) 112, [1910.06275].
  * (52) S. Abreu, H. Ita, B. Page and W. Tschernow, _Two-Loop Hexa-Box Integrals for Non-Planar Five-Point One-Mass Processes_ , 2107.14180.
  * (53) C. Balazs, J. L. Diaz-Cruz, H. J. He, T. M. P. Tait and C. P. Yuan, _Probing Higgs bosons with large bottom Yukawa coupling at hadron colliders_ , _Phys. Rev. D_ 59 (1999) 055016, [hep-ph/9807349].
  * (54) S. Dawson, C. B. Jackson, L. Reina and D. Wackeroth, _Higgs production in association with bottom quarks at hadron colliders_ , _Mod. Phys. Lett. A_ 21 (2006) 89–110, [hep-ph/0508293].
  * (55) D. Pagani, H.-S. Shao and M. Zaro, _RIP $Hb\overline{b}$: how other Higgs production modes conspire to kill a rare signal at the LHC_, _JHEP_ 11 (2020) 036, [2005.10277].
  * (56) C. Grojean, A. Paul and Z. Qian, _Resurrecting $b\overline{b}h$ with kinematic shapes_, _JHEP_ 04 (2021) 139, [2011.13945].
  * (57) P. Konar, B. Mukhopadhyaya, R. Rahaman and R. K. Singh, _Probing non-standard $b\bar{b}h$ interaction at the LHC at $\sqrt{s}=13$ TeV_, _Phys. Lett. B_ 818 (2021) 136358, [2101.10683].
  * (58) F. Maltoni, G. Ridolfi and M. Ubiali, _b-initiated processes at the LHC: a reappraisal_ , _JHEP_ 07 (2012) 022, [1203.6393].
  * (59) D. Dicus, T. Stelzer, Z. Sullivan and S. Willenbrock, _Higgs boson production in association with bottom quarks at next-to-leading order_ , _Phys. Rev. D_ 59 (1999) 094016, [hep-ph/9811492].
  * (60) C. Balazs, H.-J. He and C. P. Yuan, _QCD corrections to scalar production via heavy quark fusion at hadron colliders_ , _Phys. Rev. D_ 60 (1999) 114001, [hep-ph/9812263].
  * (61) F. Maltoni, Z. Sullivan and S. Willenbrock, _Higgs-Boson Production via Bottom-Quark Fusion_ , _Phys. Rev. D_ 67 (2003) 093005, [hep-ph/0301033].
  * (62) R. V. Harlander, K. J. Ozeren and M. Wiesemann, _Higgs plus jet production in bottom quark annihilation at next-to-leading order_ , _Phys. Lett. B_ 693 (2010) 269–273, [1007.5411].
  * (63) S. Bühler, F. Herzog, A. Lazopoulos and R. Müller, _The fully differential hadronic production of a Higgs boson via bottom quark fusion at NNLO_ , _JHEP_ 07 (2012) 115, [1204.4415].
  * (64) R. V. Harlander, S. Liebler and H. Mantler, _SusHi: A program for the calculation of Higgs production in gluon fusion and bottom-quark annihilation in the Standard Model and the MSSM_ , _Comput. Phys. Commun._ 184 (2013) 1605–1617, [1212.3249].
  * (65) A. H. Ajjath, P. Banerjee, A. Chakraborty, P. K. Dhani, P. Mukherjee, N. Rana et al., _NNLO QCD $\oplus$QED corrections to Higgs production in bottom quark annihilation_, _Phys. Rev. D_ 100 (2019) 114016, [1906.09028].
  * (66) C. Duhr, F. Dulat and B. Mistlberger, _Higgs Boson Production in Bottom-Quark Fusion to Third Order in the Strong Coupling_ , _Phys. Rev. Lett._ 125 (2020) 051804, [1904.09990].
  * (67) 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].
  * (68) J. M. Campbell, R. K. Ellis, F. Maltoni and S. Willenbrock, _Higgs-Boson Production in Association with a Single Bottom Quark_ , _Phys. Rev. D_ 67 (2003) 095002, [hep-ph/0204093].
  * (69) S. Dawson and P. Jaiswal, _Weak Corrections to Associated Higgs-Bottom Quark Production_ , _Phys. Rev. D_ 81 (2010) 073008, [1002.2672].
  * (70) S. Dawson, C. B. Jackson and P. Jaiswal, _SUSY QCD Corrections to Higgs-b Production : Is the $\Delta_{b}$ Approximation Accurate?_, _Phys. Rev. D_ 83 (2011) 115007, [1104.1631].
  * (71) S. Dittmaier, M. Krämer and M. Spira, _Higgs radiation off bottom quarks at the Tevatron and the CERN LHC_ , _Phys. Rev. D_ 70 (2004) 074010, [hep-ph/0309204].
  * (72) S. Dawson, C. B. Jackson, L. Reina and D. Wackeroth, _Exclusive Higgs boson production with bottom quarks at hadron colliders_ , _Phys. Rev. D_ 69 (2004) 074027, [hep-ph/0311067].
  * (73) S. Dawson, C. B. Jackson, L. Reina and D. Wackeroth, _Higgs boson production with one bottom quark jet at hadron colliders_ , _Phys. Rev. Lett._ 94 (2005) 031802, [hep-ph/0408077].
  * (74) M. Wiesemann, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni and P. Torrielli, _Higgs production in association with bottom quarks_ , _JHEP_ 02 (2015) 132, [1409.5301].
  * (75) B. Jager, L. Reina and D. Wackeroth, _Higgs boson production in association with b jets in the POWHEG BOX_ , _Phys. Rev. D_ 93 (2016) 014030, [1509.05843].
  * (76) N. Deutschmann, F. Maltoni, M. Wiesemann and M. Zaro, _Top-Yukawa contributions to bbH production at the LHC_ , _JHEP_ 07 (2019) 054, [1808.01660].
  * (77) N. Liu, L. Wu, P. Wu and J. M. Yang, _Complete one-loop effects of SUSY QCD in $b\bar{b}h$ production at the LHC under current experimental constraints_, _JHEP_ 01 (2013) 161, [1208.3413].
  * (78) S. Dittmaier, P. Häfliger, M. Krämer, M. Spira and M. Walser, _Neutral MSSM Higgs-boson production with heavy quarks: NLO supersymmetric QCD corrections_ , _Phys. Rev. D_ 90 (2014) 035010, [1406.5307].
  * (79) R. Harlander, M. Kramer and M. Schumacher, _Bottom-quark associated Higgs-boson production: reconciling the four- and five-flavour scheme approach_ , 1112.3478.
  * (80) M. Bonvini, A. S. Papanastasiou and F. J. Tackmann, _Resummation and matching of b-quark mass effects in $b\overline{b}H$ production_, _JHEP_ 11 (2015) 196, [1508.03288].
  * (81) S. Forte, D. Napoletano and M. Ubiali, _Higgs production in bottom-quark fusion in a matched scheme_ , _Phys. Lett. B_ 751 (2015) 331–337, [1508.01529].
  * (82) S. Forte, D. Napoletano and M. Ubiali, _Higgs production in bottom-quark fusion: matching beyond leading order_ , _Phys. Lett. B_ 763 (2016) 190–196, [1607.00389].
  * (83) C. Duhr, F. Dulat, V. Hirschi and B. Mistlberger, _Higgs production in bottom quark fusion: matching the 4- and 5-flavour schemes to third order in the strong coupling_ , _JHEP_ 08 (2020) 017, [2004.04752].
  * (84) F. Krauss and D. Napoletano, _Towards a fully massive five-flavor scheme_ , _Phys. Rev. D_ 98 (2018) 096002, [1712.06832].
  * (85) D. Figueroa, S. Honeywell, S. Quackenbush, L. Reina, C. Reuschle and D. Wackeroth, _Electroweak and QCD corrections to $Z$-boson production with one $b$ jet in a massive five-flavor scheme_, _Phys. Rev. D_ 98 (2018) 093002, [1805.01353].
  * (86) A. Banfi, G. P. Salam and G. Zanderighi, _Infrared safe definition of jet flavor_ , _Eur. Phys. J. C_ 47 (2006) 113–124, [hep-ph/0601139].
  * (87) F. Moriello, _Generalised power series expansions for the elliptic planar families of Higgs + jet production at two loops_ , _JHEP_ 01 (2020) 150, [1907.13234].
  * (88) M. Hidding, _DiffExp, a Mathematica package for computing Feynman integrals in terms of one-dimensional series expansions_ , 2006.05510.
  * (89) R. Britto, F. Cachazo and B. Feng, _New recursion relations for tree amplitudes of gluons_ , _Nucl. Phys. B_ 715 (2005) 499–522, [hep-th/0412308].
  * (90) 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].
  * (91) S. Badger, E. Chaubey, H. B. Hartanto and R. Marzucca, _Two-loop leading colour QCD helicity amplitudes for top quark pair production in the gluon fusion channel_ , _JHEP_ 06 (2021) 163, [2102.13450].
  * (92) P. Nogueira, _Automatic Feynman graph generation_ , _J. Comput. Phys._ 105 (1993) 279–289.
  * (93) J. Kuipers, T. Ueda, J. A. M. Vermaseren and J. Vollinga, _FORM version 4.0_ , _Comput. Phys. Commun._ 184 (2013) 1453–1467, [1203.6543].
  * (94) B. Ruijl, T. Ueda and J. Vermaseren, _FORM version 4.2_ , 1707.06453.
  * (95) G. Cullen, M. Koch-Janusz and T. Reiter, _Spinney: A Form Library for Helicity Spinors_ , _Comput. Phys. Commun._ 182 (2011) 2368–2387, [1008.0803].
  * (96) A. Hodges, _Eliminating spurious poles from gauge-theoretic amplitudes_ , _JHEP_ 05 (2013) 135, [0905.1473].
  * (97) S. Badger, _Automating QCD amplitudes with on-shell methods_ , _J. Phys. Conf. Ser._ 762 (2016) 012057, [1605.02172].
  * (98) K. G. Chetyrkin and F. V. Tkachov, _Integration by Parts: The Algorithm to Calculate beta Functions in 4 Loops_ , _Nucl. Phys._ B192 (1981) 159–204.
  * (99) R. N. Lee, _Presenting LiteRed: a tool for the Loop InTEgrals REDuction_ , 1212.2685.
  * (100) S. Laporta, _High precision calculation of multiloop Feynman integrals by difference equations_ , _Int. J. Mod. Phys._ A15 (2000) 5087–5159, [hep-ph/0102033].
  * (101) L. E. K., _Factorization of rational functions of several variables into partial fractions_ , _Izvestiya Vysshikh Uchebnykh Zavedenii. Matematika_ 47 (1978) .
  * (102) A. Raichev, _Leinartas’s partial fraction decomposition_ , _arXiv e-prints_ (June, 2012) arXiv:1206.4740, [1206.4740].
  * (103) J. Boehm, M. Wittmann, Z. Wu, Y. Xu and Y. Zhang, _IBP reduction coefficients made simple_ , _JHEP_ 12 (2020) 054, [2008.13194].
  * (104) M. Heller and A. von Manteuffel, _MultivariateApart: Generalized Partial Fractions_ , 2101.08283.
  * (105) K.-T. Chen, _Iterated path integrals_ , _Bull. Am. Math. Soc._ 83 (1977) 831–879.
  * (106) F. Brown, _Iterated integrals in quantum field theory_ , in _6th Summer School on Geometric and Topological Methods for Quantum Field Theory_ , 2013.  DOI.
  * (107) S. Abreu, L. J. Dixon, E. Herrmann, B. Page and M. Zeng, _The two-loop five-point amplitude in $\mathcal{N}=4$ super-Yang-Mills theory_, _Phys. Rev. Lett._ 122 (2019) 121603, [1812.08941].
  * (108) D. Chicherin, T. Gehrmann, J. Henn, P. Wasser, Y. Zhang and S. Zoia, _Analytic result for a two-loop five-particle amplitude_ , _Phys. Rev. Lett._ 122 (2019) 121602, [1812.11057].
  * (109) D. Chicherin, T. Gehrmann, J. M. Henn, P. Wasser, Y. Zhang and S. Zoia, _The two-loop five-particle amplitude in $\mathcal{N}$ = 8 supergravity_, _JHEP_ 03 (2019) 115, [1901.05932].
  * (110) S. Abreu, L. J. Dixon, E. Herrmann, B. Page and M. Zeng, _The two-loop five-point amplitude in $\mathcal{N}$ = 8 supergravity_, _JHEP_ 03 (2019) 123, [1901.08563].
  * (111) S. Caron-Huot, D. Chicherin, J. Henn, Y. Zhang and S. Zoia, _Multi-Regge Limit of the Two-Loop Five-Point Amplitudes in $\mathcal{N}=4$ Super Yang-Mills and $\mathcal{N}=8$ Supergravity_, _JHEP_ 10 (2020) 188, [2003.03120].
  * (112) D. Chicherin, J. M. Henn and G. Papathanasiou, _Cluster algebras for Feynman integrals_ , _Phys. Rev. Lett._ 126 (2021) 091603, [2012.12285].
  * (113) R. Bonciani, V. Del Duca, H. Frellesvig, J. M. Henn, M. Hidding, L. Maestri et al., _Evaluating a family of two-loop non-planar master integrals for Higgs + jet production with full heavy-quark mass dependence_ , _JHEP_ 01 (2020) 132, [1907.13156].
  * (114) H. Frellesvig, M. Hidding, L. Maestri, F. Moriello and G. Salvatori, _The complete set of two-loop master integrals for Higgs + jet production in QCD_ , _JHEP_ 06 (2020) 093, [1911.06308].
  * (115) M. Becchetti, R. Bonciani, V. Del Duca, V. Hirschi, F. Moriello and A. Schweitzer, _Next-to-leading order corrections to light-quark mixed QCD-EW contributions to Higgs boson production_ , _Phys. Rev. D_ 103 (2021) 054037, [2010.09451].
  * (116) R. Bonciani, L. Buonocore, M. Grazzini, S. Kallweit, N. Rana, F. Tramontano et al., _Mixed strong $-$electroweak corrections to the Drell$-$Yan process_, 2106.11953.
  * (117) A. B. Goncharov, _Multiple polylogarithms, cyclotomy and modular complexes_ , _arXiv e-prints_ (May, 2011) arXiv:1105.2076, [1105.2076].
  * (118) E. Remiddi and J. A. M. Vermaseren, _Harmonic polylogarithms_ , _Int. J. Mod. Phys. A_ 15 (2000) 725–754, [hep-ph/9905237].
  * (119) A. B. Goncharov, _Multiple polylogarithms and mixed Tate motives_ , _arXiv Mathematics e-prints_ (Mar., 2001) math/0103059, [math/0103059].
  * (120) J. Vollinga and S. Weinzierl, _Numerical evaluation of multiple polylogarithms_ , _Comput. Phys. Commun._ 167 (2005) 177, [hep-ph/0410259].
  * (121) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer et al., _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].