P3H-19-047, TTP19-041

Higgs decay into massive $b$-quarks at NNLO QCD
in the nested soft-collinear subtraction scheme

Arnd Behringa, Wojciech Bizońa,b.

a Institut für Theoretische Teilchenphysik (TTP), KIT, 76128 Karlsruhe,
Germany
b Institut für Kernphysik (IKP), KIT, 76344 Eggenstein-Leopoldshafen, Germany

Abstract

We present a fully differential description of a decay of a scalar Higgs boson
into massive $b$-quarks valid at next-to-next-to-leading order (NNLO) in
perturbative quantum chromodynamics (QCD). We work within the nested soft-
collinear subtraction scheme extended to accommodate massive partons. We
include the loop-induced contribution involving a Higgs coupling to a top
quark. We test our calculation against results existing in the literature,
comparing the predictions for the total decay width and jet rates.

###### Contents

  1. 1 Introduction
  2. 2 General considerations
    1. 2.1 Notation
    2. 2.2 Outline of the subtraction scheme
    3. 2.3 IR poles of virtual amplitudes
    4. 2.4 Phase-space parametrisation
    5. 2.5 Pole vs. $\overline{\textrm{MS}}$ Yukawa coupling
  3. 3 $H\to{b\bar{b}{}}$ decay at NLO
    1. 3.1 Real contribution
    2. 3.2 Virtual contribution
    3. 3.3 Pole cancellation
  4. 4 $H\to{b\bar{b}{}}$ decay at NNLO
    1. 4.1 Double-real contribution
    2. 4.2 Real-virtual contribution
    3. 4.3 Double-virtual contribution
    4. 4.4 Pole cancellation
    5. 4.5 Top-quark contribution to the $H\to{b\bar{b}{}}$ decay
  5. 5 Results
    1. 5.1 Overview of the calculation
    2. 5.2 Total width of the $H\to{b\bar{b}{}}$ decay at NNLO
    3. 5.3 Jet rates for $H\to{b\bar{b}{}}$ at NNLO
  6. 6 Conclusions
  7. A Renormalisation
  8. B Useful formulae
    1. B.1 Anomalous dimensions for IR factorisation
    2. B.2 Coefficients for on-shell to $\overline{\textrm{MS}}$-scheme conversion
  9. C Factorisation formulae
    1. C.1 Single-collinear factorisation
    2. C.2 Single-soft factorisation (tree-level)
    3. C.3 Single-soft factorisation (one-loop)
    4. C.4 Double-soft factorisation
  10. D Integrated subtraction terms
    1. D.1 Single-collinear subtraction terms
    2. D.2 Single-soft subtraction terms (tree-level)
    3. D.3 Single-soft subtraction terms (one-loop)
    4. D.4 Double-soft subtraction terms

##  1 Introduction

After the Higgs boson discovery by the ATLAS and CMS collaborations in 2012,
the study of Higgs boson properties has become one of the major research
avenues in particle physics. Since the mass of the Higgs boson has already
been precisely measured [1], all couplings between the Higgs boson and other
Standard Model (SM) particles can be accurately predicted. Nevertheless, these
couplings can be modified by New Physics that lies beyond the SM. Hence,
actual measurements of those couplings can provide important constraints on
many extensions of the SM.

The Higgs boson decay into a pair of $b$-quarks is the most common decay
channel of the Higgs boson and it is essential to many New Physics searches.
Indeed, it plays a particularly important role when considering rare Higgs
boson production modes, which benefit from the large $H\to{b\bar{b}{}}$
branching fraction. Although such measurements are often very challenging, due
to overwhelming QCD backgrounds, the $H\to{b\bar{b}{}}$ decay has already been
observed by both ATLAS and CMS [2, 3]. The upcoming years of data taking at
the Large Hadron Collider (LHC) will allow for further exploration and use of
this Higgs boson decay channel.

In order to fully utilise the data collected at the LHC, a good theoretical
understanding of the $H\to{b\bar{b}{}}$ process is required. The next-to-
leading order (NLO) QCD corrections have been available for a long time [4, 5,
6, 7, 8]. Currently, corrections to the total decay width are known up to
$\mathcal{O}(\alpha_{s}^{4})$ in the limit of massless $b$-quarks [9]. Mass
effects at $\mathcal{O}(\alpha_{s}^{2})$ have been estimated using a large
momentum expansion [10]. Furthermore, the impact of a separate class of
corrections from diagrams arising at $\mathcal{O}(\alpha_{s}^{2})$ which
involve the top-quark Yukawa coupling has been calculated in the limit of a
large top-quark mass [11, 12] as well as for general values of the masses
[13]. Recently, a set of two-loop master integrals required for mixed QCD-
electroweak corrections has been computed [14]. In the limit of massless
$b$-quarks a number of fully differential next-to-next-to-leading order (NNLO)
calculations of the $H\to{b\bar{b}{}}$ decay have been presented [15, 16, 17,
18] with first N3LO QCD results appearing recently [19]. The $b$-quark mass
effects for differential observables have been studied at NNLO QCD in Ref.
[20].

In this paper, we present an independent calculation of the $b$-quark mass
effects in the $H\to{b\bar{b}{}}$ decay at NNLO QCD. Although we believe such
a calculation is interesting in its own right and serves as a useful check of
the results presented in Ref. [20], it is also an essential step towards
studying mass effects in associated Higgs production with a vector boson,
$pp\to HV\to{b\bar{b}{}}V$. NNLO QCD corrections to this process have already
been studied in the limit of massless $b$-quarks [21, 17, 18], and large
effects, related to radiative corrections, have been reported. The impact of
the $b$-quark mass may be sizeable in certain regions of the phase space.
Higher-order effects in that process have also been investigated using parton
showers [22, 23, 24]. Another contribution to the Higgs decay width at NNLO is
mediated by top-quark loops. In the context of differential distributions, it
has been discussed in Ref. [17] and subsequently investigated in Ref. [13]. A
consistent treatment of these contributions requires keeping the $b$-quark
mass finite [17].

We work within the nested soft-collinear subtraction scheme [25, 17, 26, 27],
which is an extension of the original sector-improved residue subtraction
scheme [28, 29, 30, 31]. To incorporate the $b$-quark masses into the
calculation, we rely on the treatment of massive particles outlined in Ref.
[30].

The paper is organised as follows. In Section 2 we introduce the notation and
discuss the main steps of the subtraction scheme. We also describe the
infrared (IR) poles appearing in virtual amplitudes and touch upon the
relation between the pole and $\overline{\textrm{MS}}$ Yukawa couplings. In
Section 3 and Section 4 we review the NLO QCD calculation of the
$H\to{b\bar{b}{}}$ decay and present our calculation of the NNLO corrections,
including the treatment of the top-quark induced corrections. Finally, in
Section 5, we thoroughly test our predictions against results available in the
literature. We summarise our findings in Section 6. The appendices contain a
number of expressions used throughout our calculation.

##  2 General considerations

We are interested in decays of a scalar Higgs boson into a pair of $b$-quarks.
Our goal is to treat $b$-quarks as massive particles throughout the
calculation and achieve NNLO accuracy in perturbative QCD while working within
the nested soft-collinear subtraction scheme [25, 17, 26, 27].

###  2.1 Notation

We start with a short introduction that will set the stage for our
calculation. We consider the Higgs boson decaying into a pair of $b$-quarks

$\displaystyle H(q_{1})\longrightarrow b(q_{2})+\bar{b}(q_{3})\,,$ (2.1)

with $q_{1}^{\,2}=M_{H}^{\,2}$ and $q_{2}^{\,2}=q_{3}^{\,2}=m_{b}^{\,2}$. The
leading order (LO) partial decay width of the Higgs boson into a pair of
$b$-quarks is

$\displaystyle\Gamma_{{\mathrm{LO}}}={}$
$\displaystyle\frac{1}{2M_{H}}\int\mathrm{d}\Phi_{{b\bar{b}{}}}(q_{1})|\mathcal{M}^{(0)}_{{b\bar{b}{}}}|^{2}\,,$
(2.2)

where $\mathrm{d}\Phi_{{b\bar{b}{}}}(q_{1})$ is the two-particle phase-space
volume element for the production of two $b$-quarks with total momentum
$q_{1}$, and $|\mathcal{M}^{(0)}_{{b\bar{b}{}}}|^{2}$ is the squared tree-
level matrix element. For brevity, we list only final-state particles in the
sub- and superscripts, since the initial state is always an on-shell Higgs
boson at rest. Upon integration we obtain

$\displaystyle\Gamma_{{\mathrm{LO}}}={}$
$\displaystyle\frac{N_{c}}{16\pi}y_{b}^{2}M_{H}\beta^{3}\,,$ (2.3)

where $N_{c}=3$ is the number of colours,
$\beta=\sqrt{1-4m_{b}^{2}/M_{H}^{2}}$ and $y_{b}$ stands for the $b$-quark
Yukawa coupling, $y_{b}=m_{b}(2\sqrt{2}G_{F})^{1/2}$.

The Higgs boson decay width into $b$-quarks receives radiative corrections
that can be systematically calculated in perturbative QCD. We use the
following notation for the perturbative expansion of the width

$\displaystyle\Gamma^{{b\bar{b}{}}}={}$
$\displaystyle\Gamma_{{\mathrm{LO}}}^{{b\bar{b}{}}}\left[{1+\left({\frac{\alpha_{s}}{\pi}}\right)\gamma_{1}^{{b\bar{b}{}}}+\left({\frac{\alpha_{s}}{\pi}}\right)^{2}\gamma_{2}^{{b\bar{b}{}}}+\mathcal{O}(\alpha_{s}^{3})}\right]\,.$
(2.4)

In Eq. 2.4, $\alpha_{s}$ is the $\overline{\textrm{MS}}$ QCD strong coupling
constant defined in a theory with $n_{f}=n_{l}+1$ quark flavours, where
$n_{l}$ is the number of massless flavours.

It is useful to define a shorthand notation that denotes an integral over the
Lorentz-invariant phase space of particles involved in a particular
(sub)process. Similar to the notation in Ref. [25], we define

$\displaystyle{}F_{LM}({b\bar{b}{X}})=\mathrm{d}\Phi_{{b\bar{b}{X}}}(q_{1})|\mathcal{M}^{(0)}_{{b\bar{b}{X}}}|^{2}\mathcal{F}_{\textrm{kin}}({b\bar{b}{X}})\,,$
(2.5)

where ${b\bar{b}{X}}$ denotes the constituents of the final state of a
considered subprocess, $\mathrm{d}\Phi_{{b\bar{b}{X}}}(q_{1})$ and
$|\mathcal{M}^{(0)}_{{b\bar{b}{X}}}|^{2}$ stand for the Lorentz-invariant
phase-space measure and the squared amplitude of the $H\to{b\bar{b}{X}}$
process, respectively. The momentum $q_{1}$ refers to the initial state Higgs
boson, while $\mathcal{F}_{\textrm{kin}}$ is an infrared-safe observable that
depends on the kinematical configuration of the particles involved in the
process. We will use the notation

$\displaystyle\langle A\rangle$ $\displaystyle=\int\mathrm{d}\Phi A$ (2.6)

to denote the integration of some quantity $A$ over the phase space,
$\mathrm{d}\Phi$.

###  2.2 Outline of the subtraction scheme

One of the challenges in higher-order QCD calculations is the appearance of
infrared singularities when massless particles become soft or collinear.
Dimensional regularisation [32, 33, 34, 35, 36] can be used to regulate these
singularities which show up as poles in the dimensional regularisation
parameter $\epsilon=(4-d)/2$ in both real and virtual amplitudes. For an
infrared-safe observable these poles cancel between real and virtual
corrections and collinear factorisation contributions once the loop and phase-
space integrals over the singular regions are performed [37, 38]. The
observables depend on the momenta of the real emission partons so that
numerical integration of the phase-space integrals is desirable from the
standpoint of flexibility and often also required due to the complexity of the
observable, which may involve, e.g., complicated kinematic constraints and jet
algorithms.

It follows from the factorisation theorems of QCD that the integrand of the
cross-section, i.e., the combination of squared matrix elements and phase-
space measure, scales as $E_{i}^{-1+a\epsilon}\mathrm{d}{E_{i}}$ in soft
limits and as $\theta_{ij}^{-1+b\epsilon}\mathrm{d}\theta_{ij}$ in collinear
limits of massless partons, where $E_{i}$ is the energy of the soft parton,
$\theta_{ij}$ is the angle between the collinear partons and
$a,b\in\mathbb{R}$. Thus, $\epsilon$ acts as a regulator for logarithmic
divergences in the real-emission phase-space integrals. For the numerical
integration it is necessary to explicitly extract and cancel the poles in
$\epsilon^{-1}$ and to regulate the integrals in such a way that the expansion
in $\epsilon$ can be performed at the integrand level.

A number of methods have been developed to accomplish that. Here, we follow
the nested soft-collinear subtraction scheme [25, 17, 26, 27] which is closely
related to the sector-improved residue subtraction scheme [28, 29, 30, 31]. As
with all subtraction schemes, the general idea is to introduce subtraction
terms for each singular limit which regulate the integrand and to add back
these subtraction terms integrated over the unresolved phase space.
Schematically, for an integral of a function $F$ over the phase space, we
write

$\displaystyle\langle F\rangle$ $\displaystyle=\langle F-OF\rangle+\langle
OF\rangle\,,$ (2.7)

where $O$ is an operator which extracts the asymptotic behaviour of $F$ and
the phase space in a singular limit. The term $\langle OF\rangle$, which is
integrated over the unresolved phase space in $d$ dimensions, then carries
explicit poles in $\epsilon^{-1}$, while the regulated term $\langle
F-OF\rangle$ is expanded in $\epsilon$ at the integrand level. We apply Eq.
2.7 recursively to regulate all singular limits.

The nested soft-collinear subtraction scheme consists of the following steps.

  1. 1.

Introduce subtraction terms for the soft limits. At NNLO this involves up to
two single-soft limits and the double-soft limit.

  2. 2.

Introduce a partition of unity for the phase space that isolates the collinear
limits, $1=\sum_{\\{p\\}}w_{\\{p\\}}$, where the sum runs over the sets of
partons that can produce collinear singularities and the functions
$w_{\\{p\\}}$ go to zero whenever two partons which are not in ${p}$ become
collinear, thereby regulating integrand in that limit.

  3. 3.

In each collinear partition, choose a suitable phase-space parametrisation in
terms of angles and energies of the partons that can become unresolved.

  4. 4.

Use sector decomposition [39, 40, 41] to map all singularities to the
boundaries of the region of integration so that the singularities can be
easily extracted upon using Eq. 2.7. In order to generate the limits of the
matrix elements, we use the standard QCD factorisation formulae for the soft
and collinear limits. All necessary expressions up to NNLO can be found, e.g.,
in the appendices of Ref. [30].

The $H\to{b\bar{b}{}}$ process with massive $b$-quarks is particularly simple
in this context since there are no triple- or double-collinear limits that
involve $b$-quarks so that step 2 can be avoided. Moreover, by choosing an
appropriate phase-space parametrisation in step 3, the sector decomposition of
step 4 only yields a single collinear subsector.

The calculation will be subdivided into pieces so that the cancellation of
$\epsilon^{-1}$ poles can be shown without making reference to the explicit
form of the matrix elements. In particular, we will organise the contributions
into sets according to the multiplicity of the resolved final state and the
loop order of the matrix elements. By combining this with a suitable phase-
space parametrisation which decouples the integrations over the unresolved and
resolved parts of phase space, we demonstrate pole cancellation for each
phase-space point of the resolved configuration separately. In Sections 3 and
4 we explain the application of this scheme to the $H\to{b\bar{b}{}}$ process
in greater detail.

###  2.3 IR poles of virtual amplitudes

The one- and two-loop virtual amplitudes that we encounter in an NNLO
calculation feature ultraviolet (UV) divergences that can be removed using a
suitable renormalisation procedure. We employ a hybrid scheme in which we
renormalise quark and gluon fields, the quark masses and the Yukawa coupling
in an on-shell scheme, while we use the $\overline{\textrm{MS}}$ scheme for
the strong coupling. The details of our renormalisation choice are described
in Appendix A.

At this point, the renormalised amplitudes are free of UV divergences.
Nevertheless, they still contain poles in $\epsilon^{-1}$ which are of IR
origin. These poles can be predicted from general considerations [42, 43, 44,
45, 46, 47, 48, 49, 50]. They factorise in the form

$\displaystyle|\mathcal{M}\rangle$
$\displaystyle=\mathbf{Z}|\mathcal{F}\rangle\,,$ (2.8)

where $|\mathcal{M}\rangle$ is the UV-renormalised amplitude, $\mathbf{Z}$ is
an operator in colour space which contains poles in $\epsilon^{-1}$ and
$|\mathcal{F}\rangle$ is a finite remainder which does not contain any poles.
Expanding all pieces in the strong coupling, we find

$\displaystyle|\mathcal{M}\rangle={}$
$\displaystyle|\mathcal{M}^{(0)}\rangle+\frac{\alpha_{s}}{4\pi}\left(\mathbf{Z}^{(1)}|\mathcal{M}^{(0)}\rangle+|\mathcal{F}^{(1)}\rangle\right)$
$\displaystyle+\left(\frac{\alpha_{s}}{4\pi}\right)^{2}\left(\mathbf{Z}^{(2)}|\mathcal{M}^{(0)}\rangle+\mathbf{Z}^{(1)}|\mathcal{F}^{(1)}\rangle+|\mathcal{F}^{(2)}\rangle\right)\,.$
(2.9)

Note that this notation leaves all powers of the strong coupling related to
real emissions implicit inside the amplitudes. On the one hand, we can use
Section 2.3 as a prediction in order to check the $\epsilon^{-1}$ poles of the
UV-renormalised amplitudes. On the other hand, we can also use Section 2.3 to
define the finite remainders, i.e.

$\displaystyle|\mathcal{F}^{(0)}\rangle$
$\displaystyle=|\mathcal{M}^{(0)}\rangle\,,$ (2.10)
$\displaystyle|\mathcal{F}^{(1)}\rangle$
$\displaystyle=|\mathcal{M}^{(1)}\rangle-\mathbf{Z}^{(1)}|\mathcal{M}^{(0)}\rangle\,,$
(2.11) $\displaystyle|\mathcal{F}^{(2)}\rangle$
$\displaystyle=|\mathcal{M}^{(2)}\rangle-\mathbf{Z}^{(1)}|\mathcal{F}^{(1)}\rangle-\mathbf{Z}^{(2)}|\mathcal{M}^{(0)}\rangle\,,$
(2.12)

and express all formulae in terms of these. This is useful for showing pole
cancellation since it allows us to make the pole terms explicit without
specifying the matrix elements that they multiply. The pole terms are
multiplied by lower order quantities, as expected. Note that we only include
the $\epsilon^{-1}$ poles in the definition of the $\mathbf{Z}$ operator, cf.
Refs. [45, 50].

In general, $\mathbf{Z}$ is an operator acting on vectors in colour space,
which expresses non-trivial correlations between different colour
configurations [51, 50]. However, in our case, the coefficients can be
expressed in terms of simple colour factors since we only require virtual
amplitudes with up to three coloured particles ($H\to{b\bar{b}{}}$ and
$H\to{b\bar{b}{g}}$). The expansion coefficients
$\mathbf{Z}^{(k)}_{{b\bar{b}{}}}$ for the $H\to{b\bar{b}{}}$ process are given
by

$\displaystyle\mathbf{Z}^{(1)}_{{b\bar{b}{}}}={}$
$\displaystyle\frac{1}{2\epsilon}(C_{F}\gamma_{\mathrm{cusp,Q}}^{(0)}(v_{23})+2\gamma_{Q}^{(0)})\,,$
(2.13) $\displaystyle\mathbf{Z}^{(2)}_{{b\bar{b}{}}}={}$
$\displaystyle\frac{1}{8\epsilon^{2}}(C_{F}\gamma_{\mathrm{cusp,Q}}^{(0)}(v_{23})+2\gamma_{Q}^{(0)})(C_{F}\gamma_{\mathrm{cusp,Q}}^{(0)}(v_{23})+2\gamma_{Q}^{(0)}-2\beta_{0}(n_{l}))$
$\displaystyle+\frac{1}{4\epsilon}(C_{F}\gamma_{\mathrm{cusp,Q}}^{(1)}(v_{23})+2\gamma_{Q}^{(1)})+\frac{1}{2\epsilon}(C_{F}\gamma_{\mathrm{cusp,Q}}^{(0)}(v_{23})+2\gamma_{Q}^{(0)})\beta_{0,Q}\ln\left({\frac{\mu_{R}^{2}}{m_{b}^{2}}}\right)\,,$
(2.14)

where $\gamma_{\mathrm{cusp}}^{(i)}$ and
$\gamma_{\mathrm{cusp,Q}}^{(0)}(v_{23})$ are the massless and massive cusp
anomalous dimensions, $v_{23}=\sqrt{1-m_{b}^{4}/({q_{2}}\cdot{q_{3}})^{2}}$,
$\beta_{0}(n_{l})$ is the zeroth-order coefficient of the QCD $\beta$-function
with $n_{l}$ massless flavours,
$\beta_{0}(n_{l})=\frac{11}{3}C_{A}-\frac{4}{3}T_{F}n_{l}$ and
$\beta_{0,Q}=-\frac{4}{3}T_{F}$. The $\gamma_{Q}^{(i)}$ denote the expansion
coefficients of the anomalous dimensions of the massive quark and $\mu_{R}$ is
the renormalisation scale. We collect the necessary formulae in Section B.1.
For the $H\to{b\bar{b}{g}}$ process, where we only need the one-loop
amplitude, we find

$\displaystyle\mathbf{Z}^{(1)}_{{b\bar{b}{g}}}={}$
$\displaystyle\frac{1}{4\epsilon^{2}}(-C_{A}\gamma_{\mathrm{cusp}}^{(0)})+\frac{1}{2\epsilon}\left[\gamma_{g}^{(0)}+2\gamma_{Q}^{(0)}+\left(C_{F}-\frac{C_{A}}{2}\right)\gamma_{\mathrm{cusp,Q}}^{(0)}(v_{23})\right.$
$\displaystyle\left.-\frac{C_{A}}{2}\gamma_{\mathrm{cusp}}^{(0)}\left(\ln\left(\frac{m_{b}\mu_{R}}{2({q_{2}}\cdot{q_{4}})}\right)+\ln\left(\frac{m_{b}\mu_{R}}{2({q_{3}}\cdot{q_{4}})}\right)+2i\pi\right)\right]\,.$
(2.15)

The gluon anomalous dimension $\gamma_{g}^{(0)}$ is also given in Section B.1.

###  2.4 Phase-space parametrisation

In this section we outline the parametrisation of the real-emission phase
space that we employ throughout the calculation. NNLO corrections to Higgs
decays involve contributions with up to two real emissions accompanying the
Born process. In our case, the Born process consists of the Higgs boson
decaying into massive $b$-quarks ($H\to{b\bar{b}{}}$) and we include final
states with one additional gluon ($H\to{b\bar{b}{g}}$) as well as two
additional massless partons ($H\to{b\bar{b}{gg}}$ or
$H\to{b\bar{b}{q\bar{q}{}}}$)111The $H\to{b\bar{b}{{b\bar{b}{}}}}$
contribution is finite. Thus, it can be integrated using standard techniques
and we do not discuss it here..

The guiding principle behind the construction outlined in this section is,
first, to explicitly parametrise the energies and angles that are responsible
for the soft and collinear singularities and, second, to decouple the real-
emission phase space from the phase space of the reduced process once a parton
becomes unresolved. We note that we work in the Higgs boson rest frame
throughout the paper.

The phase-space measure for an emission of a single massless parton in
$d=4-2\epsilon$ space-time dimensions reads

$\displaystyle[\mathrm{d}{}q_{i}]={}$
$\displaystyle(\mu_{R}^{2})^{\epsilon}S_{\epsilon}\frac{\mathrm{d}^{d-1}q_{i}}{(2\pi)^{d-1}(2E_{i})}\,,$
(2.16)

where we denote the parton momentum by $q_{i}$ and its energy by $E_{i}$. Note
that we also include a global factor $(\mu_{R}^{2})^{\epsilon}S_{\epsilon}$
that originates from the strong coupling renormalisation, see discussion below
Appendix A. We do not introduce an upper bound on the energy of the emitted
gluon since it naturally appears due to the energy-momentum conserving
$\delta$-function once the measure in Eq. 2.16 is considered as a part of a
specific process.

#### Single-emission phase space

We start with the process

$\displaystyle H(q_{1})\longrightarrow b(q_{2})+\bar{b}(q_{3})+g(q_{4})\,,$
(2.17)

and discuss its phase-space parametrisation.

The phase-space measure reads

$\displaystyle\int\mathrm{d}\Phi_{{b\bar{b}{g}}}(q_{1})$
$\displaystyle={}\int[\mathrm{d}{}q_{4}]\int\mathrm{d}\Phi_{{b\bar{b}{}}}(q_{1}-q_{4})$
$\displaystyle={}(\mu_{R}^{2})^{\epsilon}S_{\epsilon}\int\frac{\mathrm{d}^{d-2}\hat{q}_{4}}{2(2\pi)^{d-1}}\int\mathrm{d}E_{4}(E_{4})^{d-3}\int\mathrm{d}\Phi_{{b\bar{b}{}}}(q_{1}-q_{4}),$
(2.18)

where $\mathrm{d}\Phi_{{b\bar{b}{}}}(Q)$ stands for the Born phase space of
the two $b$-quarks with total momentum $Q$, and $\hat{q}_{4}$ determines the
direction of the gluon momentum,

$\displaystyle\hat{q}_{4}^{\mu}={}$ $\displaystyle q_{4}^{\mu}/E_{4}\,.$
(2.19)

We further parametrise the gluon energy as

$\displaystyle E_{4}$ $\displaystyle=E_{\mathrm{max}}\,\xi_{1},$ (2.20)

with $E_{\mathrm{max}}=\tfrac{1}{2}\beta^{2}M_{H}$ and
$\beta=\sqrt{1-4m_{b}^{2}/M_{H}^{2}}$. This finally leads us to

$\displaystyle\int\mathrm{d}\Phi_{{b\bar{b}{g}}}(q_{1})={}$ $\displaystyle
2^{1-2\epsilon}E_{\mathrm{max}}^{2}\left({\frac{\mu_{R}^{2}}{E_{\mathrm{max}}^{2}}}\right)^{\epsilon}S_{\epsilon}\int\frac{\mathrm{d}\Omega_{4}^{(2-2\epsilon)}}{2(2\pi)^{3-2\epsilon}}$
$\displaystyle\times\int_{0}^{1}\frac{\mathrm{d}\eta}{(\eta(1-\eta))^{\epsilon}}\int_{0}^{1}\mathrm{d}\xi_{1}\xi_{1}^{1-2\epsilon}\int\mathrm{d}\Phi_{{b\bar{b}{}}}(q_{1}-q_{4})\,.$
(2.21)

The $\mathrm{d}\Omega_{4}$ element denotes the angular integral over the
direction of $\hat{q}_{4}$. The remaining angular integrals can be performed
using

$\displaystyle\int\mathrm{d}\Omega^{(a)}$
$\displaystyle={}\frac{2\pi^{a/2}}{\Gamma(a/2)}\,.$ (2.22)

For a single gluon emission, the only unresolved limit is the single-soft one,
i.e. $\xi_{1}\to 0$. Obviously, this removes $q_{4}$ from the overall momentum
conservation and the integration over the unresolved phase space of $q_{4}$
decouples from the Born phase space. Thus, in that limit we just replace
$\int\mathrm{d}\Phi_{{b\bar{b}{}}}(q_{1}-q_{4})$ by
$\int\mathrm{d}\Phi_{{b\bar{b}{}}}(q_{1})$ in Section 2.4.

#### Double-emission phase space

We now focus on the parametrisation of the $H\to{b\bar{b}{gg}}$ phase space.
Note that, since we consider $b$-quarks to be massive, there are no
singularities associated with kinematic configurations where gluons become
collinear to $b$-quarks. Therefore, we do not need to partition the phase
space into subsectors, which are usually necessary to disentangle the
collinear singularities. Instead, we work with a global parametrisation. In
this section we focus on the two-gluon emission case since the parametrisation
of the $q\bar{q}$ emission phase space is nearly identical. We comment on the
differences where necessary.

We consider the process

$\displaystyle H(q_{1})\longrightarrow
b(q_{2})+\bar{b}(q_{3})+g(q_{4})+g(q_{5})\,.$ (2.23)

We denote the sum of the gluon momenta by $q_{45}=q_{4}+q_{5}$. The phase
space measure then reads

$\displaystyle\int\mathrm{d}\Phi_{{b\bar{b}{gg}}}(q_{1})={}$
$\displaystyle\int[\mathrm{d}{}q_{4}]\int[\mathrm{d}{}q_{5}]\int\mathrm{d}\Phi_{{b\bar{b}{}}}(q_{1}-q_{45})$
$\displaystyle={}$
$\displaystyle(\mu_{R}^{2})^{2\epsilon}S_{\epsilon}^{2}\int\frac{\mathrm{d}^{d-2}\hat{q}_{4}}{2(2\pi)^{d-1}}\int\frac{\mathrm{d}^{d-2}\hat{q}_{5}}{2(2\pi)^{d-1}}$
$\displaystyle\hphantom{(\mu_{R}^{2})^{2\epsilon}}\times\int\mathrm{d}E_{4}(E_{4})^{d-3}\int\mathrm{d}E_{5}(E_{5})^{d-3}\int\mathrm{d}\Phi_{{b\bar{b}{}}}(q_{1}-q_{45}),$
(2.24)

where the vectors $\hat{q}_{4}$ and $\hat{q}_{5}$ determine the directions of
the two gluons and the limits of the energy integrals are so that the whole
phase space is covered.

It is convenient to introduce an energy ordering among the gluons by
partitioning the phase space via

$\displaystyle 1$ $\displaystyle={}\Theta(E_{4}-E_{5})+\Theta(E_{5}-E_{4})\,$
(2.25)

which leads to the split

$\displaystyle\int\mathrm{d}\Phi_{{b\bar{b}{gg}}}(q_{1})={}$
$\displaystyle\int\mathrm{d}\Phi_{{b\bar{b}{gg}}}^{E_{4}>E_{5}}(q_{1})+\int\mathrm{d}\Phi_{{b\bar{b}{gg}}}^{E_{5}>E_{4}}(q_{1})\,.$
(2.26)

Throughout the article, we describe calculations only for the region with
$E_{4}>E_{5}$; the other region can easily be covered by performing the same
steps with the gluon momenta swapped, $q_{4}\leftrightarrow q_{5}$.222Note
that thanks to the symmetry between the gluons we have
$\int\mathrm{d}\Phi_{{b\bar{b}{gg}}}(q_{1})=2\int\mathrm{d}\Phi_{{b\bar{b}{gg}}}^{E_{4}>E_{5}}(q_{1})$.
However when considering a $q\bar{q}$ emission such a simplification may only
be used if the symmetry $q\leftrightarrow{\bar{q}}$ also holds for the
observable under consideration. Otherwise the two parts of the phase space,
introduced in Eq. 2.26, have to be considered separately.  Hence, we
parametrise the gluon energies as [31]

$\displaystyle E_{4}$
$\displaystyle={}E_{45,\mathrm{max}}\,\xi_{1}\left(1-\frac{\xi_{2}}{2}\right)\,,$
$\displaystyle E_{5}$
$\displaystyle={}E_{45,\mathrm{max}}\,\xi_{1}\frac{\xi_{2}}{2}\,,$ (2.27)

where $E_{45,\mathrm{max}}$ is to be chosen such that the integration ranges
$\xi_{1}\in[0,1]$ and $\xi_{2}\in[0,1]$ span the whole phase space. Using
momentum conservation and considering a configuration where the two $b$-quarks
are produced at threshold, we obtain

$\displaystyle E_{45,\mathrm{max}}$
$\displaystyle={}\frac{M_{H}\beta^{2}}{1+\sqrt{1-\beta^{2}\bar{q}_{45}^{2}}}\,,$
(2.28)

where $\bar{q}_{45}=q_{45}/q^{0}_{45}$, which is only a light-like momentum
when $q_{5}$ is soft or $q_{4}$ and $q_{5}$ are collinear. In this way, we
effectively parametrise the sum of the gluon energies
($\xi_{1}=E_{45}/E_{45,\mathrm{max}}$) and their ratio
($\xi_{2}=2E_{5}/E_{45}$).

We also explicitly parametrise the angle $\theta_{45}$ between the two gluons
as follows

$\displaystyle\eta$ $\displaystyle=\frac{1}{2}(1-\cos\theta_{45})\,.$ (2.29)

The first step in the phase-space construction is to choose a direction for
$\hat{q}_{4}$. Here, we explicitly parametrise angle $\theta_{4}$ between the
emission and the $\hat{z}$-axis. Next, we fix the direction of $\hat{q}_{5}$
relative to $\hat{q}_{4}$ using the angle $\theta_{45}$ between them and the
angle $\phi$ which is the azimuthal angle of $\hat{q}_{5}$ around the
direction of $\hat{q}_{4}$. Given the directions $\hat{q}_{4}$ and
$\hat{q}_{5}$ as well as $\xi_{2}$, we calculate the vector $\bar{q}_{45}$ via

$\displaystyle\bar{q}_{45}$
$\displaystyle=\left(1-\frac{\xi_{2}}{2}\right)\hat{q}_{4}+\frac{\xi_{2}}{2}\hat{q}_{5}\,.$
(2.30)

This is sufficient to use Eq. 2.28 in order to calculate the upper bound on
the energy and from that also the individual energies $E_{4}$ and $E_{5}$,
which fully determines $q_{4}$ and $q_{5}$. Then we generate a Born phase-
space configuration with invariant mass $Q^{2}=(q_{1}-q_{45})^{2}$ in its rest
frame; this is a back-to-back configuration of the two $b$-quarks. Finally, we
boost the Born configuration to have total momentum $Q=q_{1}-q_{45}$ in the
Higgs rest frame in order to restore momentum conservation. The corresponding
phase-space measure reads

$\displaystyle{\int\mathrm{d}\Phi_{{b\bar{b}{gg}}}^{E_{4}>E_{5}}(q_{1})={}2^{2\epsilon-2}(\mu_{R}^{2})^{2\epsilon}S_{\epsilon}^{2}}$
$\displaystyle\times\int\frac{\mathrm{d}\Omega_{4}^{(2-2\epsilon)}}{2(2\pi)^{3-2\epsilon}}\int_{0}^{\pi}\mathrm{d}\theta_{4}(\sin\theta_{4})^{1-2\epsilon}\int\frac{\mathrm{d}\Omega_{5}^{(1-2\epsilon)}}{2(2\pi)^{3-2\epsilon}}\int_{0}^{1}\frac{\mathrm{d}\eta}{\eta^{\epsilon}(1-\eta)^{\epsilon}}\int_{0}^{\pi}\mathrm{d}\phi(\sin\phi)^{-2\epsilon}$
$\displaystyle\times\int_{0}^{1}\mathrm{d}\xi_{1}\int_{0}^{1}\mathrm{d}\xi_{2}E_{45,\mathrm{max}}^{4-4\epsilon}\xi_{1}^{3-4\epsilon}\xi_{2}^{1-2\epsilon}(2-\xi_{2})^{1-2\epsilon}\int\mathrm{d}\Phi_{{b\bar{b}{}}}(q_{1}-q_{45})\,.$
(2.31)

The unparametrised angles in $\mathrm{d}\Omega_{4}$ and $\mathrm{d}\Omega_{5}$
can be integrated in the end using Eq. 2.22.

The parametrisation shown in Section 2.4 achieves the desired decoupling of
the emission phase space in unresolved limits. In the collinear ($\eta\to 0$)
and single-soft ($\xi_{2}\to 0$) limits the energy bound $E_{45,\mathrm{max}}$
simplifies to $E_{\mathrm{max}}=\frac{1}{2}\beta^{2}M_{H}$ and becomes
independent of $\eta$ and $\xi_{2}$. In the single-soft limit the momentum
$q_{5}$ decouples from the energy-momentum conserving $\delta$-function and
the integrations over $\eta$ and $\xi_{2}$ decouple from the resolved phase
space. In the collinear limit the momentum conservation depends only on the
sum of momenta $q_{45}$, which is the on-shell momentum of the massless parent
parton of the splitting and is independent of $\xi_{2}$. Again, the
integrations over $\eta$ and $\xi_{2}$ decouple from the resolved phase space.
Note that the collinear and the single-soft limits both yield the same
resolved configuration. Furthermore, in the double-soft limit, both gluons
decouple from the momentum conservation and the integrals over $\xi_{1}$,
$\xi_{2}$ and $\eta$ can be carried out for a fixed Born configuration.

###  2.5 Pole vs. $\overline{\textrm{MS}}$ Yukawa coupling

Already in the first calculation of the radiative corrections to the Higgs
boson decay rate to fermions discussed in Ref. [4], it has been recognised
that the result expressed in terms of the on-shell Yukawa coupling contains
large logarithms in the limit $m_{b}/M_{H}\ll 1$. It has also been shown there
that these large logarithms can be avoided by reexpressing the result in terms
of the $\overline{\textrm{MS}}$ Yukawa coupling evaluated at the
renormalisation scale $\mu=M_{H}$.

The Yukawa coupling in the $\overline{\textrm{MS}}$ scheme
$\overline{y}_{b}(\mu)$ is related to the $\overline{\textrm{MS}}$ mass
$\overline{m}_{b}(\mu)$ via

$\displaystyle\overline{y}_{b}(\mu)={}$
$\displaystyle(2\sqrt{2}G_{F})^{1/2}\,\overline{m}_{b}(\mu)\,,$ (2.32)

where $\mu$ is the renormalisation scale. Thus, the relation between the on-
shell and $\overline{\textrm{MS}}$ Yukawa couplings can be deduced from the
corresponding relation between the masses, which we need up to
$\operatorname{\mathcal{O}}\left(\alpha_{s}^{2}\right)$ [52, 53]. It reads

$\displaystyle y_{b}^{2}={}$
$\displaystyle\overline{y}_{b}^{2}(\mu)\left[{1+\left({\frac{\alpha_{s}(\mu)}{\pi}}\right)r_{1}(m_{b},\mu)+\left({\frac{\alpha_{s}(\mu)}{\pi}}\right)^{2}r_{2}(m_{b},\mu)+\mathcal{O}(\alpha_{s}^{3})}\right]\,,$
(2.33)

where the coefficients $r_{i}(m_{b},\mu)$ are presented in Section B.2.

Note that we only reexpress the overall Yukawa coupling in this way, but we
keep the mass dependence of the matrix elements and kinematical invariants in
terms of the pole mass, similar to Ref. [20].

The total decay width and its expansion coefficients computed with the
$\overline{\textrm{MS}}$ Yukawa coupling are denoted with a bar, i.e.

$\displaystyle\overline{\Gamma}^{{b\bar{b}{}}}={}$
$\displaystyle\overline{\Gamma}_{{\mathrm{LO}}}^{{b\bar{b}{}}}\left[{1+\left({\frac{\alpha_{s}}{\pi}}\right)\overline{\gamma}_{1}^{{b\bar{b}{}}}+\left({\frac{\alpha_{s}}{\pi}}\right)^{2}\overline{\gamma}_{2}^{{b\bar{b}{}}}+\mathcal{O}(\alpha_{s}^{3})}\right]\,,$
(2.34)

where the expansion coefficients in the two schemes are related by

$\displaystyle\overline{\gamma}_{1}^{{b\bar{b}{}}}={}$
$\displaystyle\gamma_{1}^{{b\bar{b}{}}}+r_{1}\,,$ (2.35)
$\displaystyle\overline{\gamma}_{2}^{{b\bar{b}{}}}={}$
$\displaystyle\gamma_{2}^{{b\bar{b}{}}}+r_{1}\gamma_{1}^{{b\bar{b}{}}}+r_{2}\,.$
(2.36)

As discussed before, the large logarithmic corrections to the total decay
width can be mitigated by reexpressing the result in terms of the running
$\overline{\textrm{MS}}$ mass of the $b$-quark. However, in a fully
differential calculation these logarithms partially enter through corrections
related to real emissions. In this case, they arise during phase-space
integration of the emission. A priori, since the mass of the $b$-quark is
small compared to the Higgs mass, one could be worried about possible
numerical instabilities when working with a fully differential calculation.
However, it turns out that in our implementation they do not pose serious
numerical problems.

##  3 $H\to{b\bar{b}{}}$ decay at NLO

In this section we briefly describe the calculation of the NLO QCD corrections
to the $H\to{b\bar{b}{}}$ decay with massive $b$-quarks. Although such a
calculation is straightforward, we find it useful to review it in order to
clarify our notation and conventions. At this order of perturbation theory we
need to consider real (R) and virtual corrections (V).

Nowadays the fully differential treatment of this decay mode can easily be
obtained using the FKS [54, 55] or Catani-Seymour [51, 56, 57] subtraction
schemes. At NLO, our approach is essentially equivalent to the FKS subtraction
method.

###  3.1 Real contribution

At NLO we consider one real emission in addition to the Born process, which
means that we need to integrate the function

$\displaystyle{}F_{LM}({b\bar{b}{g}})={}$
$\displaystyle\mathrm{d}\Phi_{{b\bar{b}{g}}}(q_{1})|\mathcal{M}^{(0)}_{{b\bar{b}{g}}}|^{2}\mathcal{F}_{\textrm{kin}}({b\bar{b}{g}})\,,$
(3.1)

over the ${b\bar{b}{g}}$ phase space. This integral is divergent in four
dimensions, due to the soft singularity of the gluon. However, there are no
collinear singularities since they are regulated by the $b$-quark mass.

We define a projection operator that allows us to extract the soft divergence
and to regulate the limit. Given a quantity $A$ that depends on the momenta,
we define a projection operator for the soft limit of momentum $q_{4}$ as

$\displaystyle S_{4}A={}$ $\displaystyle\lim_{\xi_{1}\to 0}A\,,$ (3.2)

where $\xi_{1}$ refers to the parametrisation of Section 2.4. We define the
operator to act on all quantities to the right of the $S_{4}$ symbol,
extracting the leading asymptotic behaviour in $\xi_{1}$ of the quantity $A$
if the actual limit does not exist.

Denoting the identity operation by $I$, we can immediately write

$\displaystyle{}F_{LM}({b\bar{b}{g}})={}$
$\displaystyle{(I-S_{4})}F_{LM}({b\bar{b}{g}})+{S_{4}}F_{LM}({b\bar{b}{g}})\,,$
(3.3)

where the first term is now regularised in the soft limit and can be
integrated in four dimensions. The soft singularity is exposed in the second
term in Eq. 3.3, which therefore needs to be evaluated in $d$ dimensions. Once
the soft limit is taken, the only remaining dependence on $\xi_{1}$ in
${S_{4}}F_{LM}({b\bar{b}{g}})$ is the leading behaviour
$\xi_{1}^{-1-2\epsilon}$. Thus, the $\xi_{1}$ integration becomes trivial and
the soft singularity manifests itself as an explicit $\epsilon^{-1}$ pole.

In order to show pole cancellation pointwise in the Born phase space, we split
the real emission contribution into two parts. The finite contribution is
given by

$\displaystyle
2M_{H}\,\langle{\mathrm{d}\Gamma}_{\mathrm{R}}^{\mathrm{F}}{({b\bar{b}{g}})}\rangle={}$
$\displaystyle\langle{(I-S_{4})}{}F_{LM}({b\bar{b}{g}})\rangle$
$\displaystyle={}$
$\displaystyle\left\langle{}F_{LM}({b\bar{b}{g}})+g_{s}^{2}C_{F}\left({\mathcal{S}^{(0)}_{22,4}-2\mathcal{S}^{(0)}_{23,4}+\mathcal{S}^{(0)}_{33,4}}\right)S_{4}[\mathrm{d}{}q_{4}]{}F_{LM}({b\bar{b}{}})\right\rangle\,.$
(3.4)

Here, we use the factorisation formula for the soft limit as discussed in
Section C.2, where also the eikonal factors $\mathcal{S}^{(0)}_{ij,k}$ are
defined. Moreover, we have the unresolved contribution, which contains the
integrated subtraction term and reads

$\displaystyle
2M_{H}\,\langle{\mathrm{d}\Gamma}_{\mathrm{R}}^{\mathrm{U}}{({b\bar{b}{g}})}\rangle={}$
$\displaystyle\langle{S_{4}}{}F_{LM}({b\bar{b}{g}})\rangle$ $\displaystyle={}$
$\displaystyle-\left\langle
g_{s}^{2}C_{F}\left({\mathcal{S}^{(0)}_{22,\mathrm{int}}-2\mathcal{S}^{(0)}_{23,\mathrm{int}}+\mathcal{S}^{(0)}_{33,\mathrm{int}}}\right){}F_{LM}({b\bar{b}{}})\right\rangle\,,$
(3.5)

with the integrated eikonal factors $\mathcal{S}^{(0)}_{ij,\mathrm{int}}$
given in Section D.2. Here, the integral over the unresolved phase space of
the gluon was performed and we are only left with the phase-space integral
over the underlying Born process.

###  3.2 Virtual contribution

For the virtual contribution, the phase-space integration is the same as that
for the Born process, but we need to consider a one-loop virtual amplitude.
Although this amplitude has an $\epsilon^{-1}$ pole, the singular part can be
written as a product of a tree-level matrix element and a kinematics-dependent
coefficient, as indicated in Section 2.3. We have

$\displaystyle
2\operatorname{Re}\langle\mathcal{M}^{(0)}_{{b\bar{b}{}}}|\mathcal{M}^{(1)}_{{b\bar{b}{}}}\rangle={}$
$\displaystyle
2\operatorname{Re}\langle\mathcal{M}^{(0)}_{{b\bar{b}{}}}|\mathcal{F}^{(1)}_{{b\bar{b}{}}}\rangle+2\operatorname{Re}(\mathbf{Z}^{(1)}_{{b\bar{b}{}}})\langle\mathcal{M}^{(0)}_{{b\bar{b}{}}}|\mathcal{M}^{(0)}_{{b\bar{b}{}}}\rangle\,,$
(3.6)

where the term with the $\mathbf{Z}$ operator contains an explicit
$\epsilon^{-1}$ pole, while the second term is finite. As a shorthand we
introduce

$\displaystyle{}F_{LV}^{\mathrm{fin}}({b\bar{b}{}})={}$
$\displaystyle\mathrm{d}\Phi_{{b\bar{b}{}}}(q_{1})\left({\frac{\alpha_{s}}{4\pi}}\right)2\operatorname{Re}\langle\mathcal{M}^{(0)}_{{b\bar{b}{}}}|\mathcal{F}^{(1)}_{{b\bar{b}{}}}\rangle\mathcal{F}_{\textrm{kin}}({b\bar{b}{}})\,.$
(3.7)

Accordingly, we define two contributions to the virtual correction: the
virtual finite contribution

$\displaystyle
2M_{H}\,\langle{\mathrm{d}\Gamma}_{\mathrm{V}}^{\mathrm{F}}{({b\bar{b}{}})}\rangle={}$
$\displaystyle\langle{}F_{LV}^{\mathrm{fin}}({b\bar{b}{}})\rangle$ (3.8)

and the virtual unresolved contribution

$\displaystyle
2M_{H}\,\langle{\mathrm{d}\Gamma}_{\mathrm{V}}^{\mathrm{U}}{({b\bar{b}{}})}\rangle={}$
$\displaystyle\left\langle\left({\frac{\alpha_{s}}{4\pi}}\right)2\operatorname{Re}(\mathbf{Z}^{(1)}_{{b\bar{b}{}}}){}F_{LM}({b\bar{b}{}})\right\rangle\,.$
(3.9)

The expansion coefficients of the $\mathbf{Z}$ operator are given in Eq. 2.13.

###  3.3 Pole cancellation

At this point we can combine all contributions that enter the NLO calculation.
We remind the reader that the
${\mathrm{d}\Gamma}_{\mathrm{R}}^{\mathrm{F}}{({b\bar{b}{g}})}$ and
${\mathrm{d}\Gamma}_{\mathrm{V}}^{\mathrm{F}}{({b\bar{b}{}})}$ terms are free
of $\epsilon^{-1}$ poles, while the
${\mathrm{d}\Gamma}_{\mathrm{R}}^{\mathrm{U}}{({b\bar{b}{g}})}$ and
${\mathrm{d}\Gamma}_{\mathrm{V}}^{\mathrm{U}}{({b\bar{b}{}})}$ terms feature
$\epsilon^{-1}$ poles.

Expanding the explicit results for the real unresolved contribution from Eq.
3.5 up to $\operatorname{\mathcal{O}}\left(\epsilon^{-1}\right)$ we find

$\displaystyle
2M_{H}\,\langle{\mathrm{d}\Gamma}_{\mathrm{R}}^{\mathrm{U}}{({b\bar{b}{g}})}\rangle={}$
$\displaystyle\frac{1}{\epsilon}\left[{\left({\frac{\alpha_{s}}{4\pi}}\right)4C_{F}\left[{1+\frac{1+\beta^{2}}{2\beta}\log\left({\frac{1-\beta}{1+\beta}}\right)}\right]\langle{}{}F_{LM}({b\bar{b}{}})\rangle}\right]+\operatorname{\mathcal{O}}\left(\epsilon^{0}\right)\,.$
(3.10)

Analogously, the virtual unresolved contribution, defined in Eq. 3.9, yields

$\displaystyle
2M_{H}\,\langle{\mathrm{d}\Gamma}_{\mathrm{V}}^{\mathrm{U}}{({b\bar{b}{}})}\rangle={}$
$\displaystyle-\frac{1}{\epsilon}\left[{\left({\frac{\alpha_{s}}{4\pi}}\right)4C_{F}\left[{1+\frac{1+\beta^{2}}{2\beta}\log\left({\frac{1-\beta}{1+\beta}}\right)}\right]\langle{}{}F_{LM}({b\bar{b}{}})\rangle}\right]\,.$
(3.11)

Note that the constant term in $\epsilon$ is absent since our definition of
the $\mathbf{Z}$ operator of Eq. 2.8 includes only $\epsilon^{-1}$ poles.

Obviously, the poles of the two contributions cancel. A crucial part of the
argument is to notice that the
${\mathrm{d}\Gamma}_{\mathrm{R}}^{\mathrm{U}}{({b\bar{b}{g}})}$ part, after
integrating out the real emission, is a function of a Born-like phase-space
configuration, as is the
${\mathrm{d}\Gamma}_{\mathrm{V}}^{\mathrm{U}}{({b\bar{b}{}})}$ contribution.
This allows us to demonstrate the pole cancellation for any point of the
$\mathrm{d}\Phi_{{b\bar{b}{}}}$ phase space and without specifying the
explicit form of the $\mathrm{LO}$ matrix element.

##  4 $H\to{b\bar{b}{}}$ decay at NNLO

We now consider the NNLO QCD corrections to the $H\to{b\bar{b}{}}$ process
keeping the full dependence on the $b$-quark mass. To this end, we need to
consider several contributions including

  * •

the double-real contribution (RR) – where the leading-order decay is
accompanied by an emission of a pair of massless partons ($H\to{b\bar{b}{gg}}$
and $H\to{b\bar{b}{q\bar{q}{}}}$) or an additional pair of $b$-quarks
($H\to{b\bar{b}{{b\bar{b}{}}}}$);

  * •

the real-virtual contribution (RV) – where we consider one-loop virtual
corrections to the process $H\to{b\bar{b}{g}}$;

  * •

the double-virtual contribution (VV) – where we consider two-loop virtual
corrections to the Born process $H\to{b\bar{b}{}}$.

Except for the $H\to{b\bar{b}{{b\bar{b}{}}}}$ subprocess, all these
contributions are divergent in four dimensions, due to soft and collinear
singularities. For that reason, we follow the general method recapitulated in
Section 2.2 and adopt a subtraction scheme that allows us to regulate all
singular limits and treat the divergent integrals analytically in
$d=4-2\epsilon$ dimensions.

Apart from these divergent contributions we distinguish the subprocess
$H\to{b\bar{b}{{b\bar{b}{}}}}$ which enters the calculation at
$\operatorname{\mathcal{O}}\left(\alpha_{s}^{2}\right)$. Indeed, this
subprocess is finite in four dimensions because of the non-zero $b$-quark
mass. Hence, it does not require any regularisation. It is calculated by
directly integrating the squared tree-level amplitude over the phase space
$\mathrm{d}\Phi_{{b\bar{b}b\bar{b}}}$. In our implementation we use the
sequential algorithm [58] to generate kinematic configurations and the phase-
space measure.

Finally, a distinct class of corrections to the $H\to{b\bar{b}{}}$ decay that
appears at second order of perturbation theory is related to Feynman diagrams
where the $H\to{b\bar{b}{}}$ or $H\to{b\bar{b}{g}}$ transition is induced by
the Higgs boson coupling to gluons via a top-quark loop. This contribution is
finite on its own and, hence, can be studied separately – we defer the
discussion of these top-quark mediated corrections to Section 4.5.

All tree-level amplitudes that we use in this paper are calculated using the
spinor-helicity formalism.333For a review, see e.g. Ref. [59]. The treatment
of massive external particles follows along the lines of Appendix A of Ref.
[60]. The one-loop amplitudes are calculated using a combination of Passarino-
Veltman reduction [61], to express them through one-loop scalar integrals, and
spinor-helicity techniques, to treat spinor structures appearing in the
amplitudes. The one-loop scalar integrals are evaluated using the library
QCDLoop [62, 63]. We assemble the two-loop using the two-loop scalar heavy-
quark form factor from Ref. [64]; equivalent results can be obtained using the
expressions from Ref. [65]. The form factor is expressed in terms of harmonic
polylogarithms [66], which we evaluate using HPLOG [67].

It is useful to stress that we show cancellation of all $\epsilon^{-1}$ poles
without referring to the specific form of the matrix elements. The expressions
for all amplitudes are only needed to calculate finite corrections to the
considered process and, hence, we restrict ourselves to the construction of
only four-dimensional matrix elements.

###  4.1 Double-real contribution

In this section we only focus on the $H\to{b\bar{b}{gg}}$ and
$H\to{b\bar{b}{q\bar{q}{}}}$ subprocesses. The $H\to{b\bar{b}{{b\bar{b}{}}}}$
process is completely finite on its own and does not require any
regularisation procedure. For the record we write

$\displaystyle
2M_{H}\,\langle{\mathrm{d}\Gamma}_{{b\bar{b}b\bar{b}}}{}\rangle={}$
$\displaystyle\langle{}{}F_{LM}({b\bar{b}b\bar{b}})\rangle\,,$ (4.1)

where the shorthand ${}F_{LM}({b\bar{b}b\bar{b}})$ is introduced in Eq. 2.5.

Thanks to the non-zero $b$-quark mass, the singularity structure of the
double-real contribution is simple. Indeed, we only need to take into account
three possible limits:

  * •

the soft limit ($S_{5}$) – where the energy of one of the partons vanishes,
i.e. $\xi_{2}\to 0$;

  * •

the double-soft limit($S_{45}$) – where the energies of both additional
partons vanish at a similar rate, i.e. $\xi_{1}\to 0$;

  * •

the collinear limit ($C_{45}$) – where the momenta of the two additional
partons become collinear to each other, i.e. $\eta\to 0$.

The variables $\xi_{1}$, $\xi_{2}$ and $\eta$ refer to the parametrisation
introduced in Section 2.4.

We define projection operators that allow us to extract divergences in each of
the singular regions. Given a quantity $A$ which depends on the momenta of the
$b$-quarks and gluons, we define the action of the projection operators as
follows

$\displaystyle S_{45}A={}$ $\displaystyle\lim_{\xi_{1}\to 0}A\,,$
$\displaystyle S_{5}A={}$ $\displaystyle\lim_{\xi_{2}\to 0}A\,,$
$\displaystyle C_{45}A={}$ $\displaystyle\lim_{\eta\to 0}A\,.$ (4.2)

Again, we note that taking limits in Eq. 4.2 should be understood as
extracting the most singular part of the quantity $A$ in a particular limit
whenever the limit in the conventional sense does not exist.

With these operators, we construct a nested subtraction formula which extracts
all singularities of the double-real contribution,

$\displaystyle{}F_{LM}({b\bar{b}{gg}})={}$
$\displaystyle{(I-S_{45})}F_{LM}({b\bar{b}{gg}})+{S_{45}}F_{LM}({b\bar{b}{gg}})$
$\displaystyle={}$
$\displaystyle{(I-S_{5})(I-S_{45})}F_{LM}({b\bar{b}{gg}})+{S_{5}(I-S_{45})}F_{LM}({b\bar{b}{gg}})+{S_{45}}F_{LM}({b\bar{b}{gg}})$
$\displaystyle={}$
$\displaystyle{(I-S_{5})(I-S_{45})(I-C_{45})}F_{LM}({b\bar{b}{gg}})+{S_{5}(I-S_{45})(I-C_{45})}F_{LM}({b\bar{b}{gg}})$
$\displaystyle+{(I-S_{5})(I-S_{45})C_{45}}F_{LM}({b\bar{b}{gg}})+{S_{5}(I-S_{45})C_{45}}F_{LM}({b\bar{b}{gg}})$
$\displaystyle+{S_{45}}F_{LM}({b\bar{b}{gg}})\,.$ (4.3)

To derive Section 4.1, we start by regularising the double-soft singularity
($S_{45}$), followed by further regularisation of the single-soft limit
($S_{5}$). In the last step we introduce a subtraction term for the collinear
singularity ($C_{45}$). Moreover, as these operators commute with each other
[25], we have a freedom to choose which limit to take first. Note that in
Section 4.1 we start with the collinear limit, where appropriate, having in
mind the simplicity of the corresponding factorisation formulae.

We subdivide Section 4.1 into separate contributions according to the final
state multiplicity. This allows us to discuss pole cancellation for each of
these contributions separately.

The first term on the right-hand side of Section 4.1 represents the fully
regulated double-real contribution. Therefore, it can be evaluated in four
dimensions using standard numerical techniques. We write

$\displaystyle
2M_{H}\,\langle{\mathrm{d}\Gamma}_{\mathrm{RR}}^{\mathrm{F}}{({b\bar{b}{gg}})}\rangle={}$
$\displaystyle\langle{(I-S_{5})(I-S_{45})(I-C_{45})}{}F_{LM}({b\bar{b}{gg}})\rangle\,,$
(4.4)

and similarly for the $q\bar{q}$ emission

$\displaystyle
2M_{H}\,\langle{\mathrm{d}\Gamma}_{\mathrm{RR}}^{\mathrm{F}}{({b\bar{b}{q\bar{q}}})}\rangle={}$
$\displaystyle
n_{l}\,\langle{(I-S_{45})(I-C_{45})}{}F_{LM}({b\bar{b}{q\bar{q}}})\rangle+(q\leftrightarrow\bar{q})\,,$
(4.5)

where $n_{l}$ is the number of massless quark flavours and the last term
corresponds to the phase-space region with $E_{5}>E_{4}$. Note that in case of
the $q\bar{q}$ emission we do not subtract the single-soft limit ($S_{5}$)
since this limit is not singular in case of a $g\to q\bar{q}$ splitting. The
explicit form of the subtraction terms generated by the limit operators can be
constructed using the factorisation formulae collected in Appendix C.

The next three terms on the right-hand side of Section 4.1 are regulated in
the double-soft limit, but they are evaluated in at least one of the other two
limits ($C_{45}$ or $S_{5}$). As this leaves one of the real emissions
unresolved, we call this contribution single-unresolved. Since we perform the
integral over $[\mathrm{d}{}q_{5}]$ analytically, we do not need to further
regulate the $S_{5}$ and $C_{45}$ limits. Thus, we rearrange these terms such
that for the gluon emissions they read

$\displaystyle
2M_{H}\,{\mathrm{d}\Gamma}_{\mathrm{RR}}^{\mathrm{SU}}{({b\bar{b}{gg}})}={}$
$\displaystyle\bigl{\langle}\bigl{(}S_{5}(I-S_{45})(I-C_{45})+(I-S_{5})(I-S_{45})C_{45}$
$\displaystyle+S_{5}(I-S_{45})C_{45}\bigr{)}{}F_{LM}({b\bar{b}{gg}})\bigr{\rangle}$
$\displaystyle={}$
$\displaystyle\left\langle\bigl{(}S_{5}(I-S_{45})+(I-S_{45})C_{45}-S_{5}(I-S_{45})C_{45}\bigr{)}{}F_{LM}({b\bar{b}{gg}})\right\rangle$
$\displaystyle={}$
$\displaystyle\Bigl{\langle}g_{s}^{2}\left({\mathcal{P}^{(0)}_{gg,\mathrm{int}}-\mathcal{P}^{(0),\mathrm{soft}}_{gg,\mathrm{int}}+C_{A}(\mathcal{S}^{(0)}_{24,\mathrm{int}}+\mathcal{S}^{(0)}_{34,\mathrm{int}}-\mathcal{S}^{(0)}_{23,\mathrm{int}})}\right.$
$\displaystyle\left.{-C_{F}(\mathcal{S}^{(0)}_{22,\mathrm{int}}-2\mathcal{S}^{(0)}_{23,\mathrm{int}}+\mathcal{S}^{(0)}_{33,\mathrm{int}})}\right)(\xi_{1}/2)^{-2\epsilon}{(I-S_{4})}F_{LM}({b\bar{b}{g}})\Bigr{\rangle}\,.$
(4.6)

For the $q\bar{q}$ emission, we obtain

$\displaystyle
2M_{H}\,\langle{\mathrm{d}\Gamma}_{\mathrm{RR}}^{\mathrm{SU}}{({b\bar{b}{q\bar{q}}})}\rangle={}$
$\displaystyle
2n_{l}\,\langle{(I-S_{45})C_{45}}{}F_{LM}({b\bar{b}{q\bar{q}}})\rangle$
$\displaystyle={}$ $\displaystyle 2n_{l}\left\langle
g_{s}^{2}\mathcal{P}^{(0)}_{q\bar{q},\mathrm{int}}(\xi_{1}/2)^{-2\epsilon}{(I-S_{4})}F_{LM}({b\bar{b}{g}})\right\rangle\,,$
(4.7)

where we use the fact that in the collinear limit any infrared-safe observable
must be symmetric under $q\leftrightarrow\bar{q}$, which leads to the factor
$2$ in Section 4.1. To get from the first to the second lines of Sections 4.1
and 4.1, we use the factorisation formulae given in Appendix C and integrate
them over the unresolved phase space $[\mathrm{d}{}q_{5}]$, which yields the
integrated splitting functions $\mathcal{P}^{(0)}_{ij,\mathrm{int}}$ and the
integrated eikonal factors $\mathcal{S}^{(0)}_{ij,\mathrm{int}}$. Their
calculation is described given in Appendix D. The soft-regulated matrix
element ${(I-S_{4})}F_{LM}({b\bar{b}{g}})$ has the same form as the real
finite contribution at NLO, Eq. 3.4.

Finally, the last term in Section 4.1 is evaluated in the double-soft limit.
In this case, both partons are unresolved and we call this contribution
double-unresolved. We find

$\displaystyle
2M_{H}\,\langle{\mathrm{d}\Gamma}_{\mathrm{RR}}^{\mathrm{DU}}{({b\bar{b}{gg}})}\rangle={}$
$\displaystyle\langle{S_{45}}{}F_{LM}({b\bar{b}{gg}})\rangle$
$\displaystyle={}$ $\displaystyle\left\langle
g_{s}^{4}\mathrm{D}{}\mathrm{Soft}^{(0)}_{gg,\mathrm{int}}(q_{2},q_{3}){}F_{LM}({b\bar{b}{}})\right\rangle\,,$
(4.8) $\displaystyle
2M_{H}\,\langle{\mathrm{d}\Gamma}_{\mathrm{RR}}^{\mathrm{DU}}{({b\bar{b}{q\bar{q}}})}\rangle={}$
$\displaystyle 2n_{l}\,\langle{S_{45}}{}F_{LM}({b\bar{b}{q\bar{q}}})\rangle$
$\displaystyle={}$ $\displaystyle 2n_{l}\left\langle
g_{s}^{4}\mathrm{D}{}\mathrm{Soft}^{(0)}_{q\bar{q},\mathrm{int}}(q_{2},q_{3}){}F_{LM}({b\bar{b}{}})\right\rangle\,.$
(4.9)

The double-soft emissions decouple from the hard process and the whole term
can be written as a product of the double-soft function and the Born matrix
element. After performing the integration over the momenta of the soft
partons, we denote the integrated double-soft function as
$\mathrm{D}{}\mathrm{Soft}^{(0)}_{ij,\mathrm{int}}(q_{2},q_{3})$. Details of
this calculation are presented in Section D.4.

###  4.2 Real-virtual contribution

We need to consider the $H\to{b\bar{b}{g}}$ amplitude at one-loop level and
integrate its product with the tree-level amplitude for the
$H\to{b\bar{b}{g}}$ process over the phase space
$\mathrm{d}\Phi_{{b\bar{b}{g}}}$. In the following, we use the shorthand

$\displaystyle{}F_{LV}({b\bar{b}{g}})$
$\displaystyle=\mathrm{d}\Phi_{{b\bar{b}{g}}}(q_{1})\left({\frac{\alpha_{s}}{4\pi}}\right)2\operatorname{Re}\langle\mathcal{M}^{(0)}_{{b\bar{b}{g}}}|\mathcal{M}^{(1)}_{{b\bar{b}{g}}}\rangle\mathcal{F}_{\textrm{kin}}({b\bar{b}{g}})\,.$
(4.10)

As for the real contribution at NLO, we regulate the soft singularity using
the soft limit operator $S_{4}$ defined in Eq. 3.2 and find

$\displaystyle{}F_{LV}({b\bar{b}{g}})$
$\displaystyle={(I-S_{4})}F_{LV}({b\bar{b}{g}})+{S_{4}}F_{LV}({b\bar{b}{g}})\,.$
(4.11)

We split the one-loop amplitude into a singular part, containing all
$\epsilon^{-1}$ poles from loop integrals, and a finite remainder, see Section
2.3. We write

$\displaystyle
2\operatorname{Re}\langle\mathcal{M}^{(0)}_{{b\bar{b}{g}}}|\mathcal{M}^{(1)}_{{b\bar{b}{g}}}\rangle$
$\displaystyle=2\operatorname{Re}(\mathbf{Z}^{(1)}_{{b\bar{b}{g}}})\langle\mathcal{M}^{(0)}_{{b\bar{b}{g}}}|\mathcal{M}^{(0)}_{{b\bar{b}{g}}}\rangle+2\operatorname{Re}\langle\mathcal{M}^{(0)}_{{b\bar{b}{g}}}|\mathcal{F}^{(1)}_{{b\bar{b}{g}}}\rangle\,,$
(4.12)

where the coefficient of the $\mathbf{Z}$ operator is given in Section 2.3 and
the second term contains the finite remainder of the one-loop amplitude
defined in Eq. 2.10.

In total, we split the calculation into four contributions, which need to be
combined with the corresponding terms from the double-real and double-virtual
contributions to show pole cancellation. The first one is given by the soft-
regulated terms containing the finite remainder of the $H\to{b\bar{b}{g}}$
matrix element. It is free of $\epsilon^{-1}$ poles and thus we refer to it as
the real-virtual finite contribution. It reads

$\displaystyle
2M_{H}\,\langle{\mathrm{d}\Gamma}_{\mathrm{RV}}^{\mathrm{F}}{}\rangle={}$
$\displaystyle\langle{(I-S_{4})}F_{LV}^{\mathrm{fin}}({b\bar{b}{g}})\rangle$
$\displaystyle={}$
$\displaystyle\Bigl{\langle}{}F_{LV}^{\mathrm{fin}}({b\bar{b}{g}})$
$\displaystyle+g_{s}^{2}C_{F}(\mathcal{S}^{(0)}_{22,4}-2\mathcal{S}^{(0)}_{23,4}+\mathcal{S}^{(0)}_{33,4})$
$\displaystyle\hphantom{-}\times\left({S_{4}[\mathrm{d}{}q_{4}]{}F_{LV}^{\mathrm{fin}}({b\bar{b}{}})+\left({\frac{\alpha_{s}}{4\pi}}\right)\left[{\mathcal{R}^{(1)}_{23,4}+Z_{A}^{(1)}+Z_{\alpha_{s}}^{(1)}}\right]_{\epsilon^{0}}S_{4}[\mathrm{d}{}q_{4}]{}F_{LM}({b\bar{b}{}})}\right)\Bigr{\rangle}\,,$
(4.13)

where the necessary formulae for the factorisation of the one-loop matrix
element in the soft limit are given in Appendix C and the the renormalisation
constants $Z_{A}$ and $Z_{\alpha_{s}}$ are given in Appendix A. The notation
$[\ldots]_{\epsilon^{0}}$ indicates that the
$\operatorname{\mathcal{O}}\left(\epsilon^{0}\right)$ term of the expression
between the brackets should be taken. We use
${}F_{LV}^{\mathrm{fin}}({b\bar{b}{g}})$ in analogy to the definition in Eq.
3.7.

The corresponding term containing the soft-regulated $\mathbf{Z}$ operator for
the $H\to{b\bar{b}{g}}$ process is called real-virtual single-unresolved
contribution and is given by

$\displaystyle
2M_{H}\,\langle{\mathrm{d}\Gamma}_{\mathrm{RV}}^{\mathrm{SU}}{}\rangle={}$
$\displaystyle\left\langle\left({\frac{\alpha_{s}}{4\pi}}\right)(I-S_{4})\,2\operatorname{Re}(\mathbf{Z}^{(1)}_{{b\bar{b}{g}}}){}F_{LM}({b\bar{b}{g}})\right\rangle$
$\displaystyle={}$
$\displaystyle\Bigl{\langle}\left({\frac{\alpha_{s}}{4\pi}}\right)2\operatorname{Re}(\mathbf{Z}^{(1)}_{{b\bar{b}{g}}}){}F_{LM}({b\bar{b}{g}})$
$\displaystyle+g_{s}^{2}\left({\frac{\alpha_{s}}{4\pi}}\right)C_{F}(\mathcal{S}^{(0)}_{22,4}-2\mathcal{S}^{(0)}_{23,4}+\mathcal{S}^{(0)}_{33,4})$
$\displaystyle\hphantom{-}\times\left({\left[{\mathcal{R}^{(1)}_{23,4}+Z_{\alpha_{s}}^{(1)}+Z_{A}^{(1)}}\right]_{\mathrm{poles}}+2\operatorname{Re}(\mathbf{Z}^{(1)}_{{b\bar{b}{}}})}\right)S_{4}[\mathrm{d}{}q_{4}]{}F_{LM}({b\bar{b}{}})\Bigr{\rangle}\,,$
(4.14)

where the subscript “poles” refers to taking only the pole terms up to
$\operatorname{\mathcal{O}}\left(\epsilon^{-1}\right)$ into account.

The remaining two contributions are evaluated in the soft limit and are
integrated over the unresolved phase space $[\mathrm{d}{}q_{4}]$. The
integrated subtraction term containing the finite remainder carries the
superscript $\mathrm{FR}$ and reads

$\displaystyle
2M_{H}\,\langle{\mathrm{d}\Gamma}_{\mathrm{RV}}^{\mathrm{FR}}{}\rangle={}$
$\displaystyle\langle{S_{4}}F_{LV}^{\mathrm{fin}}({b\bar{b}{g}})\rangle$
$\displaystyle={}$ $\displaystyle-\left\langle
g_{s}^{2}C_{F}\left({\mathcal{S}^{(0)}_{22,\mathrm{int}}-2\mathcal{S}^{(0)}_{23,\mathrm{int}}+\mathcal{S}^{(0)}_{33,\mathrm{int}}}\right){}F_{LV}^{\mathrm{fin}}({b\bar{b}{}})\right\rangle\,.$
(4.15)

Again, the $\mathcal{S}^{(0)}_{ij,\mathrm{int}}$ correspond to integrated
eikonal factors discussed in Section D.2.

Finally, the contribution containing the soft limit of the $\mathbf{Z}$
operator integrated over the unresolved phase space is referred to as the
real-virtual double-unresolved contribution. We obtain

$\displaystyle
2M_{H}\,\langle{\mathrm{d}\Gamma}_{\mathrm{RV}}^{\mathrm{DU}}{}\rangle={}$
$\displaystyle\left\langle\left({\frac{\alpha_{s}}{4\pi}}\right)S_{4}\,2\operatorname{Re}(\mathbf{Z}^{(1)}_{{b\bar{b}{g}}}){}F_{LM}({b\bar{b}{g}})\right\rangle$
$\displaystyle={}$
$\displaystyle\Bigl{\langle}g_{s}^{2}\left({\frac{\alpha_{s}}{4\pi}}\right)C_{F}\left({\mathcal{R}^{(1)}_{\mathrm{int}}-\left[{Z_{\alpha_{s}}^{(1)}+Z_{A}^{(1)}+2\operatorname{Re}(\mathbf{Z}^{(1)}_{{b\bar{b}{}}})}\right]}\right.$
$\displaystyle\left.{\hphantom{g_{s}^{2}\left({\frac{\alpha_{s}}{4\pi}}\right)C_{F}}\times(\mathcal{S}^{(0)}_{22,\mathrm{int}}-2\mathcal{S}^{(0)}_{23,\mathrm{int}}+\mathcal{S}^{(0)}_{33,\mathrm{int}})}\right){}F_{LM}({b\bar{b}{}})\Bigr{\rangle}\,.$
(4.16)

The integrated one-loop soft function $\mathcal{R}^{(1)}_{\mathrm{int}}$ is
defined and calculated in Section D.3.

###  4.3 Double-virtual contribution

For the $H\to{b\bar{b}{}}$ process, the virtual corrections are described by a
single form factor,

$\displaystyle|\mathcal{M}_{{b\bar{b}{}}}\rangle$
$\displaystyle=F_{s}(M_{H}^{2},m_{b}^{2},\mu_{R}^{2})|\mathcal{M}^{(0)}_{{b\bar{b}{}}}\rangle\,.$
(4.17)

In Eq. 4.17, $|\mathcal{M}^{(0)}_{{b\bar{b}{}}}\rangle$ is the tree-level
$H\to{b\bar{b}{}}$ amplitude, and $F_{s}$ is the scalar heavy-quark form
factor, which can be computed perturbatively as

$\displaystyle F_{s}(M_{H}^{2},m_{b}^{2},\mu_{R}^{2})={}$ $\displaystyle
1+\left({\frac{\alpha_{s}(\mu_{R}^{2})}{4\pi}}\right)F_{s}^{(1)}(M_{H}^{2},m_{b}^{2},\mu_{R}^{2})$
$\displaystyle\hphantom{1}+\left({\frac{\alpha_{s}(\mu_{R}^{2})}{4\pi}}\right)^{2}F_{s}^{(2)}(M_{H}^{2},m_{b}^{2},\mu_{R}^{2})+\operatorname{\mathcal{O}}\left(\alpha_{s}^{3}\right)\,.$
(4.18)

The expansion coefficients $F_{s}^{(k)}$ depend on Higgs and $b$-quark masses,
the renormalisation scale $\mu_{R}$ and the regulator $\epsilon$. The two-loop
heavy-quark form factor was computed in Ref. [65] up to $\epsilon^{0}$ terms
and to higher orders in $\epsilon$ in Refs. [68] and [64]. In our calculation
we use the results of Ref. [64].

We again use the $\mathbf{Z}$ operator to split the double-virtual correction
into divergent terms, that contain all $\epsilon^{-1}$ poles, and a finite
remainder. In this case, given the functions
$F_{s}^{(1)}(M_{H}^{2},m_{b}^{2},\mu_{R}^{2})$ and
$F_{s}^{(2)}(M_{H}^{2},m_{b}^{2},\mu_{R}^{2})$ we reconstruct the finite
remainders $|\mathcal{F}^{(1)}_{{b\bar{b}{}}}\rangle$ and
$|\mathcal{F}^{(2)}_{{b\bar{b}{}}}\rangle$ using the definition in Eq. 2.10.
We also checked that all $\epsilon^{-1}$ poles of the form factor are in
agreement with the $\mathbf{Z}$ operator prediction of Eqs. 2.13 and 2.3. The
finite double-virtual contribution reads

$\displaystyle
2M_{H}\,\langle{\mathrm{d}\Gamma}_{\mathrm{VV}}^{\mathrm{F}}{({b\bar{b}{}})}\rangle={}$
$\displaystyle\Bigl{\langle}\left({\frac{\alpha_{s}}{4\pi}}\right)^{2}\mathrm{d}\Phi_{{b\bar{b}{}}}\left[{2\operatorname{Re}\langle\mathcal{M}^{(0)}_{{b\bar{b}{}}}|\mathcal{F}^{(2)}_{{b\bar{b}{}}}\rangle+\langle\mathcal{F}^{(1)}_{{b\bar{b}{}}}|\mathcal{F}^{(1)}_{{b\bar{b}{}}}\rangle}\right]\mathcal{F}_{\textrm{kin}}({b\bar{b}{}})\Bigr{\rangle}\,.$
(4.19)

Furthermore, we distinguish two contributions that contain all $\epsilon^{-1}$
poles. First, we have a term that is proportional to the finite remainder of
the one-loop $H\to{b\bar{b}{}}$ matrix element. It is given by

$\displaystyle
2M_{H}\,\langle{\mathrm{d}\Gamma}_{\mathrm{VV}}^{\mathrm{FR}}{}\rangle={}$
$\displaystyle\Bigl{\langle}\left({\frac{\alpha_{s}}{4\pi}}\right)2\operatorname{Re}(\mathbf{Z}^{(1)}_{{b\bar{b}{}}}){}F_{LV}^{\mathrm{fin}}({b\bar{b}{}})\Bigr{\rangle}\,.$
(4.20)

The other contribution containing $\epsilon^{-1}$ poles comes with the Born
$H\to{b\bar{b}{}}$ matrix element; it reads

$\displaystyle
2M_{H}\,\langle{\mathrm{d}\Gamma}_{\mathrm{VV}}^{\mathrm{DU}}{}\rangle={}$
$\displaystyle\Bigl{\langle}\left({\frac{\alpha_{s}}{4\pi}}\right)^{2}\left({2\operatorname{Re}(\mathbf{Z}^{(2)}_{{b\bar{b}{}}})+{\mathbf{Z}^{(1)}_{{b\bar{b}{}}}}^{\dagger}\mathbf{Z}^{(1)}_{{b\bar{b}{}}}}\right){}F_{LM}({b\bar{b}{}})\Bigr{\rangle}\,.$
(4.21)

The expansion coefficients of the $\mathbf{Z}$ operator are given in Eqs. 2.13
and 2.3.

###  4.4 Pole cancellation

After collecting formulae for all NNLO contributions in the preceeding
subsections, we demonstrate pole cancellation. We begin by considering
${\mathrm{d}\Gamma}^{\mathrm{FR}}{}={\mathrm{d}\Gamma}_{\mathrm{RV}}^{\mathrm{FR}}{}+{\mathrm{d}\Gamma}_{\mathrm{VV}}^{\mathrm{FR}}{}$.
From the real-virtual contribution, Section 4.2, we obtain the pole term

$\displaystyle
2M_{H}\,\langle{\mathrm{d}\Gamma}_{\mathrm{RV}}^{\mathrm{FR}}{({b\bar{b}{g}})}\rangle={}$
$\displaystyle\frac{1}{\epsilon}\left[{\left({\frac{\alpha_{s}}{4\pi}}\right)4C_{F}\left[{1+\frac{1+\beta^{2}}{2\beta}\log\left({\frac{1-\beta}{1+\beta}}\right)}\right]\langle{}F_{LV}^{\mathrm{fin}}({b\bar{b}{}})\rangle}\right]+\operatorname{\mathcal{O}}\left(\epsilon^{0}\right)\,,$
(4.22)

and the explicit expansion of the double-virtual contribution, Eq. 4.20,
yields

$\displaystyle
2M_{H}\,\langle{\mathrm{d}\Gamma}_{\mathrm{VV}}^{\mathrm{FR}}{({b\bar{b}{}})}\rangle={}$
$\displaystyle-\frac{1}{\epsilon}\left[{\left({\frac{\alpha_{s}}{4\pi}}\right)4C_{F}\left[{1+\frac{1+\beta^{2}}{2\beta}\log\left({\frac{1-\beta}{1+\beta}}\right)}\right]\langle{}F_{LV}^{\mathrm{fin}}({b\bar{b}{}})\rangle}\right]\,,$
(4.23)

which obviously cancels in the sum with Eq. 4.22.

For the single-unresolved contribution we have
${\mathrm{d}\Gamma}^{\mathrm{SU}}{}={\mathrm{d}\Gamma}_{\mathrm{RR}}^{\mathrm{SU}}{}+{\mathrm{d}\Gamma}_{\mathrm{RV}}^{\mathrm{SU}}{}$.
The double-real contribution, Sections 4.1 and 4.1, has an expansion which can
be cast into the form

$\displaystyle{2M_{H}\,\langle{\mathrm{d}\Gamma}_{\mathrm{RR}}^{\mathrm{SU}}{({b\bar{b}{gg}}+{b\bar{b}{q\bar{q}}})}\rangle}$
$\displaystyle={}$
$\displaystyle\Biggl{\langle}\left({\frac{\alpha_{s}}{4\pi}}\right)\Biggl{[}\frac{2C_{A}}{\epsilon^{2}}+\frac{1}{\epsilon}\biggl{[}4C_{F}+\beta_{0}(n_{l})+\frac{C_{A}-2C_{F}}{v_{23}}\log\left({\frac{1+v_{23}}{1-v_{23}}}\right)$
$\displaystyle+2C_{A}\log\left({\frac{m_{b}\mu_{R}}{2({q_{2}}\cdot{q_{4}})}}\right)+2C_{A}\log\left({\frac{m_{b}\mu_{R}}{2({q_{3}}\cdot{q_{4}})}}\right)\biggr{]}\Biggr{]}{(I-S_{4})}F_{LM}({b\bar{b}{g}})\Biggr{\rangle}+\operatorname{\mathcal{O}}\left(\epsilon^{0}\right)\,.$
(4.24)

A similar expansion holds for the real-virtual single-unresolved contribution,
Section 4.2,

$\displaystyle{2M_{H}\,\langle{\mathrm{d}\Gamma}_{\mathrm{RV}}^{\mathrm{SU}}{({b\bar{b}{g}})}\rangle}$
$\displaystyle={}$
$\displaystyle\Biggl{\langle}\left({\frac{\alpha_{s}}{4\pi}}\right)\Biggl{[}-\frac{2C_{A}}{\epsilon^{2}}-\frac{1}{\epsilon}\biggl{[}4C_{F}+\beta_{0}(n_{l})+\frac{C_{A}-2C_{F}}{v_{23}}\log\left({\frac{1+v_{23}}{1-v_{23}}}\right)$
$\displaystyle+2C_{A}\log\left({\frac{m_{b}\mu_{R}}{2({q_{2}}\cdot{q_{4}})}}\right)+2C_{A}\log\left({\frac{m_{b}\mu_{R}}{2({q_{3}}\cdot{q_{4}})}}\right)\biggr{]}\Biggr{]}{(I-S_{4})}F_{LM}({b\bar{b}{g}})\Biggr{\rangle}\,.$
(4.25)

As before, the poles of the two contributions cancel.

Finally, we turn to the double-unresolved contribution, given by
${\mathrm{d}\Gamma}^{\mathrm{DU}}{}={\mathrm{d}\Gamma}_{\mathrm{RR}}^{\mathrm{DU}}{}+{\mathrm{d}\Gamma}_{\mathrm{RV}}^{\mathrm{DU}}{}+{\mathrm{d}\Gamma}_{\mathrm{VV}}^{\mathrm{DU}}{}$.
The real-virtual and double-virtual contributions are known analytically, but
we only have numerical results for the double-real contribution at our
disposal. Therefore, we demonstrate pole cancellation numerically. We write
the pole terms as

$\displaystyle\langle{\mathrm{d}\Gamma}^{\mathrm{DU}}{}\rangle={}$
$\displaystyle\left({\frac{\alpha_{s}}{4\pi}}\right)^{2}\left[{\frac{\mathcal{C}^{\mathrm{DU},(-3)}}{\epsilon^{3}}+\frac{\mathcal{C}^{\mathrm{DU},(-2)}}{\epsilon^{2}}+\frac{\mathcal{C}^{\mathrm{DU},(-1)}}{\epsilon}}\right]\langle{}{}F_{LM}({b\bar{b}{}})\rangle+\operatorname{\mathcal{O}}\left(\epsilon^{0}\right)\,,$
(4.26)

and further subdivide the pole coefficients according to colour factors,

$\displaystyle\mathcal{C}^{\mathrm{DU},(k)}={}$ $\displaystyle
C_{F}C_{A}\,\mathcal{C}_{C_{F}C_{A}}^{\mathrm{DU},(k)}+C_{F}^{\,2}\,\mathcal{C}_{C_{F}^{2}}^{\mathrm{DU},(k)}+C_{F}T_{F}n_{l}\,\mathcal{C}_{C_{F}T_{F}n_{l}}^{\mathrm{DU},(k)}+C_{F}T_{F}\,\mathcal{C}_{C_{F}T_{F}}^{\mathrm{DU},(k)}\,.$
(4.27)

We list numerical values for these pole coefficients in Table 1. The pole
cancellation reported in the last row occurs to at least 7 to 8 significant
digits which proves pole cancellation in the double-unresolved term.

| $\mathcal{C}_{C_{F}C_{A}}^{\mathrm{DU},(-3)}$ | $\mathcal{C}_{C_{F}C_{A}}^{\mathrm{DU},(-2)}$ | $\mathcal{C}_{C_{F}^{2}}^{\mathrm{DU},(-2)}$ | $\mathcal{C}_{C_{F}T_{F}n_{l}}^{\mathrm{DU},(-2)}$ | $\mathcal{C}_{C_{F}C_{A}}^{\mathrm{DU},(-1)}$ | $\mathcal{C}_{C_{F}^{2}}^{\mathrm{DU},(-1)}$ | $\mathcal{C}_{C_{F}T_{F}n_{l}}^{\mathrm{DU},(-1)}$ | $\mathcal{C}_{C_{F}T_{F}}^{\mathrm{DU},(-1)}$
---|---|---|---|---|---|---|---|---
$\mathrm{RR}$ | $-22.11$ | $-279.75$ | $+244.32$ | $+14.74$ | $-1777.55$ | $+2672.58$ | $+185.77$ | $0$
$\mathrm{RV}$ | $+22.11$ | $+320.28$ | $-488.64$ | $-29.47$ | $+1732.44$ | $-2672.58$ | $-161.20$ | $+257.20$
$\mathrm{VV}$ | $0$ | $-40.53$ | $+244.32$ | $+14.74$ | $+45.11$ | $0$ | $-24.56$ | $-257.20$
Sum | $10^{-13}$ | $10^{-10}$ | $10^{-8}$ | $10^{-11}$ | $10^{-6}$ | $10^{-6}$ | $10^{-5}$ | $0$
Rel. canc. | $10^{-14}$ | $10^{-13}$ | $10^{-11}$ | $10^{-13}$ | $10^{-10}$ | $10^{-9}$ | $10^{-7}$ | $0$

Table 1: Numerical values of the pole coefficients of the double-unresovled
term as defined in Eq. 4.27. The numerical values correspond to
$m_{b}=4.78~{}\mathrm{GeV}$, $M_{H}=125.09~{}\mathrm{GeV}$ and the
renormalisation scale is $\mu_{R}=3M_{H}$. Each column corresponds to a
particular colour structure of a given $\epsilon$ pole. The three rows
correspond to the double-real, real-virtual, and double-virtual contributions.
In the last two rows, we report the absolute and relative level of
cancellation after adding up $\mathrm{RR}+\mathrm{RV}+\mathrm{VV}$
contributions. The last row is normalised to the largest value of each column.

###  4.5 Top-quark contribution to the $H\to{b\bar{b}{}}$ decay

An additional contribution that enters the $H\to{b\bar{b}{}}$ decay process at
$\mathcal{O}(\alpha_{s}^{2})$ is related to diagrams where the Higgs boson
couples to two gluons via a top-quark loop and the final state $b$-quarks are
generated by a gluon splitting. The corresponding diagrams are shown in Figs.
1(a) and 1(b) and their contributions to the decay rate are separately finite.

(a) Two-loop $H\to{b\bar{b}{}}$ diagram

(b) One-loop $H\to{b\bar{b}{g}}$ diagram

Figure 1: The top-Yukawa contributions to the $H\to{b\bar{b}{}}$ and
$H\to{b\bar{b}{g}}$ amplitudes. The solid thick lines represent a top quark
while the thin solid lines denote the external $b$-quarks.

This contribution to the total $H\to{b\bar{b}{}}$ width has been computed in
Refs. [11, 12] as an expansion in powers of $x_{t}=(M_{H}^{2}/m_{t}^{2})$. An
exact result, incorporated into a fully differential calculation, has been
recently published in Ref. [13]. At the level of the total decay width, the
differences between the approximate and the exact results are very small for
realistic Higgs and top-quark masses. Since we aim at testing our calculation
against results reported in Ref. [20], we will use the result of Ref. [11] to
include the top-quark contribution.

Note that this contribution should naively scale with one power of $y_{b}$ and
one power of $y_{t}$, hence we refer to this as the $y_{b}y_{t}$ contribution.
However, the amplitudes of Fig. 1 can interfere with the respective Born
amplitudes only if a helicity flip occurs on one of the $b$-quark lines. This
mechanism provides an additional power of $m_{b}$. Moreover, the top-quark
loop leads to a suppression factor of $1/m_{t}$, which makes the overall
scaling of the $y_{b}y_{t}$ contribution similar to all other terms, i.e.
proportional to a square of the $b$-quark mass.

We write the $y_{b}y_{t}$ contribution to the total decay width as [11]

$\displaystyle\Gamma^{y_{b}y_{t}}={}$
$\displaystyle\Gamma_{{\mathrm{LO}}}\,\left({\frac{\alpha_{s}}{\pi}}\right)^{2}\gamma_{2}^{y_{b}y_{t}}\,,$
(4.28)

where

$\displaystyle\gamma_{2}^{y_{b}y_{t}}={}$
$\displaystyle\beta^{-3}\left({f_{2}^{S,0}+x_{t}f_{2}^{S,1}+x_{t}^{\,2}f_{2}^{S,2}+x_{t}^{\,3}f_{2}^{S,3}+\mathcal{O}(x_{t}^{\,4})}\right)$
(4.29)

and the coefficients $f_{2}^{S,k}$ are given in Eq. (2) of Ref. [11]. The
factor of $\beta^{-3}$ appears since the normalisation factor in Eq. 4.28
involves the LO width for the massive $b$-quarks, in contrast to Ref. [11].

Since the contribution in Eq. 4.28 starts only at
$\mathcal{O}(\alpha_{s}^{2})$, changes due to decoupling, Eq. A.12, or due to
translating the on-shell Yukawa coupling to the $\overline{\textrm{MS}}$
scheme, Eq. 2.33, are of higher order in $\alpha_{s}$ and are discarded.

The real-virtual amplitude related to the diagram in Fig. 1(b) is computed
using the same techniques as the remaining one-loop amplitudes. The double-
virtual amplitude of Fig. 1(a) is obtained from the result of Ref. [11], see
Eq. 4.28, by subtracting the integrated real-virtual contribution of Fig.
1(b).

##  5 Results

In this section, we summarise our calculation and present results for the
total decay width of the $H\to{b\bar{b}{}}$ decay and a selection of jet
rates. We compare these predictions to results already available in the
literature. The main goal is to scrutinise our calculation as much as possible
to ensure its correctness.

The value of the strong coupling is set to $\alpha_{s}(M_{Z})=0.1181$ with
$M_{Z}=91.1876~{}\mathrm{GeV}$ [58]. The evolution is performed at two-loop
order with five active flavours using the package RunDec [69, 70]. We use a
value for the Higgs boson mass of $M_{H}=125.09~{}\mathrm{GeV}$ [1]. As a
starting point for the $b$-quark mass, we use
$\overline{m}_{b}(\mu=\overline{m}_{b})=4.18~{}\mathrm{GeV}$ [58]. From this
we calculate the pole mass $m_{b}=4.78~{}\mathrm{GeV}$, using the two-loop
matching formula at a matching scale of
$\mu=\overline{m}_{b}(\overline{m}_{b})$ implemented in the package RunDec.
For the value of the Fermi constant we use $G_{F}=1.166378\times
10^{-5}~{}\mathrm{GeV}^{-2}$ [58]. The central renormalisation scale is taken
to be equal to the mass of the Higgs boson, $\mu_{R}=M_{H}$, and for the scale
uncertainty estimation we vary it by a factor of $1/2$ and $2$.

###  5.1 Overview of the calculation

We combine all contributions discussed in Sections 3 and 4 to obtain the full
NLO and NNLO results. At leading order in QCD we have

$\displaystyle
2M_{H}\,\langle{\mathrm{d}\Gamma}_{{\mathrm{LO}}}^{{b\bar{b}{}}}{}\rangle={}$
$\displaystyle\langle{}{}F_{LM}({b\bar{b}{}})\rangle\,.$ (5.1)

At NLO we combine real and virtual corrections

$\displaystyle\langle{\mathrm{d}\Gamma}_{\delta{\mathrm{NLO}}}^{{b\bar{b}{}}}{}\rangle={}$
$\displaystyle\langle{\mathrm{d}\Gamma}_{\mathrm{R}}{}\rangle+\langle{\mathrm{d}\Gamma}_{\mathrm{V}}{}\rangle\,,$
(5.2)

which contain finite and unresolved parts as discussed in Section 3. The
situation is slightly more complicated at NNLO where there are more
ingredients entering the final result. The result reads

$\displaystyle\langle{\mathrm{d}\Gamma}_{\delta{\mathrm{NNLO}}}^{{b\bar{b}{}}}{}\rangle={}$
$\displaystyle\langle{\mathrm{d}\Gamma}_{\mathrm{VV}}{}\rangle+\langle{\mathrm{d}\Gamma}_{\mathrm{RV}}{}\rangle+\langle{\mathrm{d}\Gamma}_{\mathrm{RR}}{}\rangle+\langle{\mathrm{d}\Gamma}_{{b\bar{b}b\bar{b}}}{}\rangle+\langle{\mathrm{d}\Gamma}_{{b\bar{b}{}}}^{y_{b}y_{t}}{}\rangle\,.$
(5.3)

The double-virtual, real-virtual and double-real contributions consist of
finite, single- and double-unresolved parts as outlined in Section 4.

###  5.2 Total width of the $H\to{b\bar{b}{}}$ decay at NNLO

We start by considering the total width of the decay $H\to{b\bar{b}{}}$; this
implies using $\mathcal{F}_{\textrm{kin}}({b\bar{b}{X}})=1$. We first present
the results using the on-shell Yukawa coupling to compare the NNLO decay width
to a result in the large Higgs mass limit from Ref. [10]. Afterwards, we
present the results expressed in terms of the $\overline{\textrm{MS}}$ Yukawa
coupling.

#### Total width using the on-shell Yukawa coupling

We compare our results for the total decay width, Eq. 2.4, to the predictions
of Ref. [10], which has been obtained as an expansion in
$m_{b}^{2}/M_{H}^{2}$. This result was derived for the scenario of a decay of
a heavy scalar boson into a pair of top quarks. Nevertheless, it can
immediately be translated into the $H\to{b\bar{b}{}}$ decay rate if we neglect
the top-loop mediated contribution, discussed in Section 4.5. Since the
results of Ref. [10] are presented in terms of the on-shell Yukawa coupling,
we perform the comparison in this scheme.

To scrutinise our results, we split the NLO and NNLO coefficients into
independent colour structures

$\displaystyle\gamma_{1}^{{b\bar{b}{}}}={}$ $\displaystyle
C_{F}\gamma_{1}^{C_{F}}\,,$ (5.4) $\displaystyle\gamma_{2}^{{b\bar{b}{}}}={}$
$\displaystyle
C_{F}^{2}\gamma_{2}^{C_{F}^{2}}+C_{F}C_{A}\gamma_{2}^{C_{F}C_{A}}+C_{F}T_{F}n_{l}\gamma_{2}^{C_{F}T_{F}n_{l}}+C_{F}T_{F}\gamma_{2}^{C_{F}T_{F}}\,.$
(5.5)

Our findings are summarised in Table 2. The first row lists predictions of
Ref. [10] which are compared to our predictions in the second row. We see a
remarkable consistency between the two; the discrepancies are at most at the
level of our numerical errors. This good agreement is related to the smallness
of the expansion parameter in case of a Higgs decay to $b$-quarks,
$(m_{b}^{2}/M_{H}^{2})\approx 0.00146$.

| $\gamma_{1}^{C_{F}}$ | $\gamma_{2}^{C_{F}^{2}}$ | $\gamma_{2}^{C_{F}C_{A}}$ | $\gamma_{2}^{C_{F}T_{F}n_{l}}$ | $\gamma_{2}^{C_{F}T_{F}}$
---|---|---|---|---|---
Ref. [10] | $-7.446648$ | $+19.4192$ | $-53.5558$ | $+18.6286$ | $+14.7946$
Our res. | $-7.446648(7)$ | $+19.4199(10)$ | $-53.5557(20)$ | $+18.6283(2)$ | $+14.7945(1)$

Table 2: The results for the NLO and NNLO coefficients of the total decay
width split into independent colour structures. The renormalisation scale is
set to $\mu_{R}=M_{H}$. The uncertainties quoted for our results correspond to
errors from numerical integration. We note that the results do not include the
$y_{b}y_{t}$ contribution.

Note that the numerical values of the NNLO coefficients are much larger than
those of the NLO coefficients. Since this is only partially compensated by the
additional power of the strong coupling, $(\alpha_{s}/\pi)$, the NNLO
correction in this scheme still amounts to a sizeable change of the total
decay width. This behaviour is closely related to large quasi-collinear
logarithms, $L=\log(m_{b}/M_{H})$, discussed in Section 2.5. As already
indicated, this issue can be mitigated by reexpressing the results in terms of
the $\overline{\textrm{MS}}$ Yukawa coupling.

#### Total width using the $\overline{\textrm{MS}}$ Yukawa coupling

We now express the total decay width in terms of the $\overline{\textrm{MS}}$
coupling, as defined in Eq. 2.34. We use
$\overline{m}_{b}(\mu=\overline{m}_{b})=4.18~{}\mathrm{GeV}$ as an input
parameter [58] and evolve the $\overline{\textrm{MS}}$ mass with the RunDec
package [69, 70]. We obtain

$\displaystyle\overline{m}_{b}(\tfrac{1}{2}M_{H})={}$ $\displaystyle
2.9814~{}\mathrm{GeV}\,,$ $\displaystyle\overline{m}_{b}(M_{H})={}$
$\displaystyle 2.8095~{}\mathrm{GeV}\,,$
$\displaystyle\overline{m}_{b}(2M_{H})={}$ $\displaystyle
2.6641~{}\mathrm{GeV}\,.$ (5.6)

The evolution is performed with $n_{f}=5$ active flavours at two-loop order.
We use these values to evaluate the $\overline{\textrm{MS}}$ Yukawa coupling.
The NLO and NNLO coefficients together with a prediction for a total decay
width are presented in Table 3. We first discuss the results without the top-
quark contribution, described in Section 4.5. For comparison, we also include
the NLO and NNLO coefficients obtained analytically in the limit of massless
$b$-quarks given in Ref. [71].

$\mu_{R}$ | $\tfrac{1}{2}M_{H}$ | $M_{H}$ | $2M_{H}$
---|---|---|---
$\overline{\gamma}_{1}^{{b\bar{b}{}}}$ (our res.) | $+3.023597(10)$ | $+5.796203(15)$ | $+8.568783(11)$
$\overline{\gamma}_{1}^{{b\bar{b}{}}}$ (Ref. [20]) | $+3.024$ | $+5.798$ | $+8.569$
$\overline{\gamma}_{1}^{{b\bar{b}{}}}$ (Ref. [71], $m_{b}=0$) | $+2.8941$ | $+5.6667$ | $+8.4393$
$\overline{\gamma}_{2}^{{b\bar{b}{}}}$ (our res., w/o $y_{b}y_{t}$) | $-3.2466(31)$ | $+30.4376(33)$ | $+79.1755(38)$
$\overline{\gamma}_{2}^{{b\bar{b}{}}}$ (our res., with $y_{b}y_{t}$) | $+3.7123(31)$ | $+37.3965(33)$ | $+86.1345(38)$
$\overline{\gamma}_{2}^{{b\bar{b}{}}}$ (Ref. [20], with $y_{b}y_{t}$) | $+3.685$ | $+37.371$ | $+86.112$
$\overline{\gamma}_{2}^{{b\bar{b}{}}}$ (Ref. [71], $m_{b}=0$) | $-3.8368$ | $+29.1467$ | $+77.1844$
$\overline{\Gamma}_{{\mathrm{LO}}}^{{b\bar{b}{}}}~{}[\mathrm{MeV}]$ | $+2.17005$ | $+1.92702$ | $+1.73274$
$\overline{\Gamma}_{{\mathrm{NLO}}}^{{b\bar{b}{}}}~{}[\mathrm{MeV}]$ | $+2.43161$ | $+2.32781$ | $+2.21731$
$\overline{\Gamma}_{{\mathrm{NNLO}}}^{{b\bar{b}{}}}~{}[\mathrm{MeV}]$ (w/o $y_{b}y_{t}$) | $+2.42041(1)$ | $+2.40333(1)$ | $+2.36344(1)$
$\overline{\Gamma}_{{\mathrm{NNLO}}}^{{b\bar{b}{}}}~{}[\mathrm{MeV}]$ (with $y_{b}y_{t}$) | $+2.44441(1)$ | $+2.42059(1)$ | $+2.37628(1)$

Table 3:  The results for the LO, NLO and NNLO total decay width. The total
width is calculated using our results for the expansion coefficients,
$\overline{\gamma}_{1}^{{b\bar{b}{}}}$ and
$\overline{\gamma}_{2}^{{b\bar{b}{}}}$. For comparison we include
corresponding results from Ref. [20]. We also provide results in the limit of
massless $b$-quarks from Ref. [71], which do not contain the $y_{b}y_{t}$
contribution. The uncertainties quoted for our results correspond to errors
from numerical integration.

We see a reasonably good perturbative convergence of the predictions for the
total decay width, $\overline{\Gamma}^{{b\bar{b}{}}}$. At the central scale,
the NNLO corrections change the NLO result by a few percent and the NNLO
prediction stays within the NLO scale uncertainties. Since the large quasi-
collinear logarithms are removed by switching to the $\overline{\textrm{MS}}$
Yukawa coupling, the mass corrections are not large. The difference between
the NNLO coefficient of the massless and massive predictions is at the level
of about $4\%$ at the central scale, $\mu_{R}=M_{H}$.

Nevertheless, the fully massive treatment of $b$-quarks is desirable because
the corrections related to the top-quark Yukawa coupling cannot be
incorporated into a fully differential calculation with massless $b$-quarks
[17]. Furthermore, it is also important to study the impact of mass effects on
the kinematical distributions related to $H\to{b\bar{b}{}}$ decay.

#### Total decay width including the top-quark contribution

Finally, we incorporate the top-quark contribution into our predictions,
according to the discussion presented in Section 4.5. We set the top-quark
pole mass to $m_{t}=173.34~{}\mathrm{GeV}$ [20]. The NNLO coefficient
corresponding to top-quark induced contribution, Eq. 4.29, is then
$\gamma_{2}^{y_{b}y_{t}}=6.95895$ and is independent of the renormalisation
scale. Our predictions are included in Table 3.

We see that inclusion of the $y_{b}y_{t}$ contribution increases the total
width by about $1.7\%$ at the central renormalisation scale. A major part of
this correction, about $85\%$, comes from the real-virtual diagrams of the
$y_{b}y_{t}$ contribution, see Fig. 1(b).

Comparing our results for the NLO and NNLO coefficients of the total width to
those given in Ref. [20], we find excellent agreement at NLO. The NNLO
coefficients agree at the level of at least $0.7\%$ for $\mu_{R}=M_{H}/2$ and
at the sub-permill level for the other scales. Note that for
$\overline{\gamma}_{2}^{{b\bar{b}{}}}(\mu_{R}=M_{H}/2)$ there are large
cancellations when converting from the on-shell to the
$\overline{\textrm{MS}}$-scheme Yukawa coupling.

###  5.3 Jet rates for $H\to{b\bar{b}{}}$ at NNLO

We now turn to a discussion of jet rates which provides another stress test of
our calculation. We employ the Durham jet algorithm [72] with the default
recombination scheme of the parton momenta, $k_{(ij)}=k_{i}+k_{j}$. We use the
FastJet [73] implementation. We consider two cases of the clustering sequence
with $y_{\mathrm{cut}}=0.01$ and $y_{\mathrm{cut}}=0.05$, to facilitate a
numerical comparison against Ref. [20].

Similar to Eq. 2.34, we define an expansion of the differential quantities as

$\displaystyle\overline{\Gamma}^{{b\bar{b}{}}}(\mathrm{obs})={}$
$\displaystyle\overline{\Gamma}_{{\mathrm{LO}}}^{{b\bar{b}{}}}\,\left[{\overline{\gamma}_{0}^{{b\bar{b}{}}}(\mathrm{obs})+\left({\frac{\alpha_{s}}{\pi}}\right)\overline{\gamma}_{1}^{{b\bar{b}{}}}(\mathrm{obs})+\left({\frac{\alpha_{s}}{\pi}}\right)^{2}\overline{\gamma}_{2}^{{b\bar{b}{}}}(\mathrm{obs})+\mathcal{O}(\alpha_{s}^{3})}\right]\,,$
(5.7)

where “$\mathrm{obs}$” denotes a generic observable and the
$\overline{\textrm{MS}}$ quantities are related to their counterparts in the
pole scheme via

$\displaystyle\overline{\gamma}_{0}^{{b\bar{b}{}}}(\mathrm{obs})={}$
$\displaystyle\gamma_{0}^{{b\bar{b}{}}}(\mathrm{obs})\,,$ (5.8)
$\displaystyle\overline{\gamma}_{1}^{{b\bar{b}{}}}(\mathrm{obs})={}$
$\displaystyle\gamma_{1}^{{b\bar{b}{}}}(\mathrm{obs})+r_{1}\gamma_{0}^{{b\bar{b}{}}}(\mathrm{obs})\,,$
(5.9) $\displaystyle\overline{\gamma}_{2}^{{b\bar{b}{}}}(\mathrm{obs})={}$
$\displaystyle\gamma_{2}^{{b\bar{b}{}}}(\mathrm{obs})+r_{1}\gamma_{1}^{{b\bar{b}{}}}(\mathrm{obs})+r_{2}\gamma_{0}^{{b\bar{b}{}}}(\mathrm{obs})\,.$
(5.10)

Note that the sum of expansion coefficients of all relevant jet
multiplicities, $\overline{\gamma}_{i}^{{b\bar{b}{}}}(n\mathrm{jet})$, yields
the expansion coefficient of the total decay width
$\overline{\gamma}_{i}^{{b\bar{b}{}}}$.

We keep using the setup presented in Section 5.2 and report results for the
two-, three- and four-jet rates in Table 4. We present these results with and
without top-quark contributions whenever relevant, i.e. for the two- and
three-jet rate NNLO coefficients.

| $y_{\mathrm{cut}}=0.01$ | $y_{\mathrm{cut}}=0.05$
---|---|---
$\mu_{R}$ | $\tfrac{1}{2}M_{H}$ | $M_{H}$ | $2M_{H}$ | $\tfrac{1}{2}M_{H}$ | $M_{H}$ | $2M_{H}$
$\overline{\gamma}_{1}^{{b\bar{b}{}}}(\mathrm{2jet})$ | $-5.0559(1)$ | $-2.2832(1)$ | $+0.4894(1)$ | $+0.2903(1)$ | $+3.0629(1)$ | $+5.8355(1)$
$\overline{\gamma}_{1}^{{b\bar{b}{}}}(\mathrm{2jet})$ (Ref. [20]) | $-5.055$ | $-2.282$ | $+0.490$ | $+0.291$ | $+3.063$ | $+5.836$
$\overline{\gamma}_{2}^{{b\bar{b}{}}}(\mathrm{2jet})$ (w/o $y_{b}y_{t}$) | $-60.50(3)$ | $-70.68(1)$ | $-65.83(2)$ | $-25.42(1)$ | $-6.59(1)$ | $+27.31(1)$
$\overline{\gamma}_{2}^{{b\bar{b}{}}}(\mathrm{2jet})$ (with $y_{b}y_{t}$) | $-56.40(3)$ | $-66.58(1)$ | $-61.73(2)$ | $-19.52(1)$ | $-0.69(1)$ | $+33.21(1)$
$\overline{\gamma}_{2}^{{b\bar{b}{}}}(\mathrm{2jet})$ (Ref. [20]) | $-56.351$ | $-66.532$ | $-61.658$ | $-19.496$ | $-0.650$ | $+33.250$
$\overline{\gamma}_{1}^{{b\bar{b}{}}}(\mathrm{3jet})$ | $+8.0794(1)$ | $+8.0794(1)$ | $+8.0794(1)$ | $+2.7333(1)$ | $+2.7333(1)$ | $+2.7333(1)$
$\overline{\gamma}_{1}^{{b\bar{b}{}}}(\mathrm{3jet})$ (Ref. [20]) | $+8.079$ | $+8.079$ | $+8.079$ | $+2.733$ | $+2.733$ | $+2.733$
$\overline{\gamma}_{2}^{{b\bar{b}{}}}(\mathrm{3jet})$ (w/o $y_{b}y_{t}$) | $+34.09(3)$ | $+77.96(1)$ | $+121.84(2)$ | $+21.25(1)$ | $+36.09(1)$ | $+50.94(1)$
$\overline{\gamma}_{2}^{{b\bar{b}{}}}(\mathrm{3jet})$ (with $y_{b}y_{t}$) | $+36.95(3)$ | $+80.82(1)$ | $+124.70(2)$ | $+22.30(1)$ | $+37.15(1)$ | $+51.99(1)$
$\overline{\gamma}_{2}^{{b\bar{b}{}}}(\mathrm{3jet})$ (Ref. [20]) | $+36.873$ | $+80.741$ | $+124.609$ | $+22.256$ | $+37.096$ | $+51.937$
$\overline{\gamma}_{2}^{{b\bar{b}{}}}(\mathrm{4jet})$ | $+23.164(1)$ | $+23.163(1)$ | $+23.163(1)$ | $+0.9323(1)$ | $+0.9322(1)$ | $+0.9322(1)$
$\overline{\gamma}_{2}^{{b\bar{b}{}}}(\mathrm{4jet})$ (Ref. [20]) | $+23.163$ | $+23.163$ | $+23.163$ | $+0.926$ | $+0.926$ | $+0.926$

Table 4: The jet rates expansion coefficients
$\overline{\gamma}_{i}^{{b\bar{b}{}}}(\mathrm{obs})$ as defined in Eq. 5.7 and
computed using the Durham clustering algorithm with $y_{\mathrm{cut}}=0.01$
and $y_{\mathrm{cut}}=0.05$ for various choices of the renormalisation scale.
Whenever necessary, we report results without and with top-quark
contributions.

We see that the NLO coefficient of the three-jet rate and the NNLO coefficient
of the four-jet rate are scale-independent since they involve only tree-level
subprocesses. We also observe a migration of events from higher to lower jet
multiplicities when jet-cut parameter is increased.

For the NLO and NNLO coefficients of the jet rates we see agreement between
our result and that of Ref. [20] at the level of $0.1\%$ to $0.2\%$ for all
quantities considered for $y_{\mathrm{cut}}=0.01$. Similar agreement is
observed for $y_{\mathrm{cut}}=0.05$ with the exception of
$\overline{\gamma}_{2}^{{b\bar{b}{}}}(2\mathrm{jet},\mu_{R}=M_{H})$, which
however suffers from large cancellations in the course of the on-shell to
$\overline{\textrm{MS}}$-scheme conversion, and the four-jet rates.

$y_{\mathrm{cut}}$ | $0.01$ | $0.05$
---|---|---
$\overline{\Gamma}^{{b\bar{b}{}}}(\mathrm{2jet})~{}[\mathrm{MeV}]$ | $1.60395(4)$ | $2.13710(2)$
$\overline{\Gamma}^{{b\bar{b}{}}}(\mathrm{3jet})~{}[\mathrm{MeV}]$ | $0.75917(4)$ | $0.28118(2)$
$\overline{\Gamma}^{{b\bar{b}{}}}(\mathrm{4jet})~{}[\mathrm{MeV}]$ | $0.05747$ | $0.0023$
Table 5: The total jet rates at NNLO for $\mu_{R}=M_{H}$ using the Durham
clustering algorithm with $y_{\mathrm{cut}}=0.01$ and $y_{\mathrm{cut}}=0.05$,
including the $y_{b}y_{t}$ contribution.

Finally, in Table 5 we present the total jet rates for the central choice of
the renormalisation scale, $\mu_{R}=M_{H}$. We see a clear hierarchy of the
jet rates, with the two-jet rate being the largest and the four-jet rate the
smallest for both clustering parameters. This is more pronounced for the
larger jet cut, $y_{\mathrm{cut}}=0.05$, as expected.

To conclude, we note that the main objective of this section was to put our
NNLO calculation under a thorough examination. We performed a series of
comparisons of our calculation against various results available in the
literature. The positive outcome of these checks assures us of the validity of
our approach.

##  6 Conclusions

In this paper, we presented an independent calculation of the NNLO QCD
corrections to the Higgs boson decay into massive $b$-quarks. We worked in the
framework of the nested soft-collinear subtraction scheme introduced in Refs.
[25, 17, 26, 27].

A complete discussion of all necessary NNLO contributions was presented in
Section 4. In particular, we demonstrated cancellation of all $\epsilon^{-1}$
poles, related to soft and collinear singularities of QCD amplitudes. The
cancellation was obtained pointwise in phase space, without referring to a
specific form of the matrix elements. Furthermore, a full treatment of the
$b$-quark mass allowed for an inclusion of the additional contribution which
originates from a direct interaction of top quarks with the Higgs boson, at a
differential level.

Our fully differential calculation was implemented in a computer program. We
carefully tested it by performing a number of cross-checks with results
available in the literature both for the total decay width [10, 20] and jet
rates [20].

We note that the calculation presented in this paper is an important step
towards a broader phenomenological goal, namely, combining a description of
Higgs boson production with its decay into massive $b$-quarks. In the context
of associated Higgs production, where large radiative corrections for
important observables have been reported [21, 17, 18], a thorough study of
$b$-quark mass effects is an interesting topic for future research.

## Acknowledgements

We wish to thank Kirill Melnikov for stimulating discussions on the topics
treated here and very valuable comments on the manuscript. We would like to
thank Peter Marquard and Matthias Steinhauser for useful conversations about
renormalisation and decoupling. We are also grateful to the authors of Ref.
[20] for useful correspondence regarding the numerical cross-checks. This
research was supported by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) under grant 396021762 - TRR 257. The work of A.B. was in
part supported by BMBF grant 05H18VKCC1.

##  Appendix A Renormalisation

The subtraction scheme that we apply in this calculation is formulated in
terms of UV-renormalised amplitudes. Here, we describe the renormalisation
prescriptions used in our calculation.

We employ a hybrid scheme in which we renormalise the quark and gluon fields,
the quark masses and the Yukawa coupling in an on-shell scheme, whereas we use
the $\overline{\textrm{MS}}{}$ scheme with five active flavours for the strong
coupling constant. We write

$\displaystyle\alpha_{s,b}$
$\displaystyle=(\mu_{R}^{2})^{\epsilon}S_{\epsilon}Z_{\alpha_{s}}\alpha_{s}\,,$
$\displaystyle\psi_{b}^{0}$ $\displaystyle=\sqrt{Z_{\psi}}\psi_{b}\,,$
$\displaystyle G_{\mu}^{0,a}$ $\displaystyle=\sqrt{Z_{A}}G_{\mu}^{a}\,,$
$\displaystyle m_{b}^{0}$ $\displaystyle=Z_{m}m_{b}\,,$ $\displaystyle
y_{b}^{0}$ $\displaystyle=Z_{m}y_{b}\,,$ (A.1)

where
$S_{\epsilon}=(4\pi)^{-\epsilon}\exp(\epsilon\gamma_{\mathrm{E}})$.444This
agrees with the renormalisation scheme in [64], but differs from [65] where a
different choice for the renormalisation of the strong coupling is used. Both
schemes yield the same result for physical observables, but differ by terms
proportional to $\zeta_{2}$ in IR-divergent intermediate steps; this
difference has already been discussed in Ref. [20]. Due to the choice to
enforce the relation $y_{b}=m_{b}(2\sqrt{2}G_{F})^{1/2}$, through the
renormalisation condition of the Yukawa coupling, the mass and Yukawa coupling
share the same renormalisation constant. The remaining renormalisation
constants (for the massless quarks, ghosts and the gauge parameter) do not
explicitly appear in our calculation. The renormalisation of the strong
coupling constant requires [74, 75]

$\displaystyle Z_{\alpha_{s}}$
$\displaystyle=1+\frac{\alpha_{s}}{4\pi}Z_{\alpha_{s}}^{(1)}+\operatorname{\mathcal{O}}\left(\alpha_{s}^{2}\right)=1-\frac{\alpha_{s}}{4\pi}\frac{\beta_{0}(n_{l}+1)}{\epsilon}+\operatorname{\mathcal{O}}\left(\alpha_{s}^{2}\right)\,,$
(A.2)

with $\beta_{0}(n_{l}+1)=\frac{11}{3}C_{A}-\frac{4}{3}T_{F}(n_{l}+1)$. The
mass and field renormalisation constants for the massive quark, $Z_{\psi}$ and
$Z_{m}$, are given by [76]

$\displaystyle Z_{m}={}$ $\displaystyle
1+\frac{\alpha_{s,b}}{4\pi}Z_{m}^{(1)}+\left(\frac{\alpha_{s,b}}{4\pi}\right)^{2}Z_{m}^{(2)}+\operatorname{\mathcal{O}}\left(\alpha_{s,b}^{3}\right)$
$\displaystyle={}$ $\displaystyle
1+\frac{\alpha_{s,b}}{4\pi}(m^{2})^{-\epsilon}S_{\epsilon}^{-1}C_{F}\left(-\frac{3}{\epsilon}-4-\epsilon\left(8+\frac{3}{2}\zeta_{2}\right)\right)$
$\displaystyle+\left(\frac{\alpha_{s,b}}{4\pi}\right)^{2}(m^{2})^{-2\epsilon}S_{\epsilon}^{-2}\left[C_{F}^{2}\left(\frac{9}{2\epsilon^{2}}+\frac{45}{4\epsilon}+\frac{199}{8}+\left(-\frac{51}{2}+48\ln(2)\right)\zeta_{2}-12\zeta_{3}\right)\right.$
$\displaystyle+C_{A}C_{F}\left(-\frac{11}{2\epsilon^{2}}-\frac{91}{4\epsilon}-\frac{605}{8}+\left(\frac{5}{2}-24\ln(2)\right)\zeta_{2}+6\zeta_{3}\right)$
$\displaystyle\left.+C_{F}T_{F}\left(\frac{2}{\epsilon^{2}}+\frac{7}{\epsilon}+\frac{69}{2}-14\zeta_{2}\right)+C_{F}T_{F}n_{l}\left(\frac{2}{\epsilon^{2}}+\frac{7}{\epsilon}+\frac{45}{2}+10\zeta_{2}\right)\right]+\operatorname{\mathcal{O}}\left(\alpha_{s,b}^{3}\right)\,,$
(A.3) $\displaystyle Z_{\psi}={}$ $\displaystyle
1+\frac{\alpha_{s,b}}{4\pi}Z_{\psi}^{(1)}+\left(\frac{\alpha_{s,b}}{4\pi}\right)^{2}Z_{\psi}^{(2)}+\operatorname{\mathcal{O}}\left(\alpha_{s,b}^{3}\right)$
$\displaystyle={}$ $\displaystyle
1+\frac{\alpha_{s,b}}{4\pi}(m^{2})^{-\epsilon}S_{\epsilon}^{-1}C_{F}\left(-\frac{3}{\epsilon}-4-\epsilon\left(8+\frac{3}{2}\zeta_{2}\right)\right)$
$\displaystyle+\left(\frac{\alpha_{s,b}}{4\pi}\right)^{2}(m^{2})^{-2\epsilon}S_{\epsilon}^{-2}\left[C_{F}^{2}\left(\frac{9}{2\epsilon^{2}}+\frac{51}{4\epsilon}+\frac{433}{8}+\left(-\frac{147}{2}+96\ln(2)\right)\zeta_{2}-24\zeta_{3}\right)\right.$
$\displaystyle+C_{A}C_{F}\left(-\frac{11}{2\epsilon^{2}}-\frac{101}{4\epsilon}-\frac{803}{8}+\left(\frac{49}{2}-48\ln(2)\right)\zeta_{2}+12\zeta_{3}\right)$
$\displaystyle\left.+C_{F}T_{F}\left(\frac{4}{\epsilon^{2}}+\frac{19}{3\epsilon}+\frac{1139}{18}-28\zeta_{2}\right)+C_{F}T_{F}n_{l}\left(\frac{2}{\epsilon^{2}}+\frac{9}{\epsilon}+\frac{59}{2}+10\zeta_{2}\right)\right]+\operatorname{\mathcal{O}}\left(\alpha_{s,b}^{3}\right)\,.$
(A.4)

The gluon field renormalisation constant is non-vanishing in the on-shell
scheme due to heavy-quark loops and reads [77, 78]

$\displaystyle Z_{A}$
$\displaystyle=1+\frac{\alpha_{s,b}}{4\pi}Z_{A}^{(1)}+\operatorname{\mathcal{O}}\left(\alpha_{s,b}^{2}\right)=1+\frac{\alpha_{s,b}}{4\pi}\frac{\beta_{0,Q}}{\epsilon}\left(\frac{\mu_{R}^{2}}{m_{b}^{2}}\right)^{\epsilon}\Gamma(1+\epsilon)e^{\epsilon\gamma_{\mathrm{E}}}+\operatorname{\mathcal{O}}\left(\alpha_{s,b}^{2}\right)\,,$
(A.5)

with $\beta_{0,Q}=-\frac{4}{3}T_{F}$.

For single-virtual and double-virtual amplitudes of the $H\to{b\bar{b}{}}$
process we start with the unrenormalised amplitude calculated from all one-
particle-irreducible Feynman diagrams. Expanded in terms of the bare strong
coupling, we have

$\displaystyle|\hat{\mathcal{M}}_{{b\bar{b}{}}}\rangle$
$\displaystyle=|\hat{\mathcal{M}}^{(0)}_{{b\bar{b}{}}}\rangle+\frac{\alpha_{s,b}}{4\pi}|\hat{\mathcal{M}}^{(1)}_{{b\bar{b}{}}}\rangle+\left(\frac{\alpha_{s,b}}{4\pi}\right)^{2}|\hat{\mathcal{M}}^{(2)}_{{b\bar{b}{}}}\rangle+\operatorname{\mathcal{O}}\left(\alpha_{s,b}^{3}\right)$
(A.6)

and together with the LSZ factors we obtain the renormalised amplitude as

$\displaystyle|\mathcal{M}_{{b\bar{b}{}}}\rangle={}$ $\displaystyle
Z_{m}Z_{\psi}|\hat{\mathcal{M}}_{{b\bar{b}{}}}\rangle+\left(\frac{\alpha_{s,b}}{4\pi}\right)^{2}Z_{m}^{(1)}|\hat{\mathcal{C}}^{(1)}_{m,{b\bar{b}{}}}\rangle+\operatorname{\mathcal{O}}\left(\alpha_{s,b}^{3}\right)$
$\displaystyle={}$
$\displaystyle|\hat{\mathcal{M}}^{(0)}_{{b\bar{b}{}}}\rangle+\frac{\alpha_{s}}{4\pi}(\mu_{R}^{2})^{\epsilon}S_{\epsilon}\left[|\hat{\mathcal{M}}^{(1)}_{{b\bar{b}{}}}\rangle+|\hat{\mathcal{M}}^{(0)}_{{b\bar{b}{}}}\rangle(Z_{m}^{(1)}+Z_{\psi}^{(1)})\right]$
$\displaystyle+\left(\frac{\alpha_{s}}{4\pi}\right)^{2}\left[(\mu_{R}^{2})^{\epsilon}S_{\epsilon}Z_{\alpha_{s}}^{(1)}\left(|\hat{\mathcal{M}}^{(1)}_{{b\bar{b}{}}}\rangle+|\hat{\mathcal{M}}^{(0)}_{{b\bar{b}{}}}\rangle(Z_{m}^{(1)}+Z_{\psi}^{(1)})\right)\right.$
$\displaystyle+(\mu_{R}^{2})^{2\epsilon}S_{\epsilon}^{2}\left(|\hat{\mathcal{M}}^{(2)}_{{b\bar{b}{}}}\rangle+(Z_{m}^{(1)}+Z_{\psi}^{(1)})|\hat{\mathcal{M}}^{(1)}_{{b\bar{b}{}}}\rangle\right.$
$\displaystyle\left.\left.+(Z_{m}^{(2)}+Z_{m}^{(1)}Z_{\psi}^{(1)}+Z_{\psi}^{(2)})|\hat{\mathcal{M}}^{(0)}_{{b\bar{b}{}}}\rangle+Z_{m}^{(1)}|\hat{\mathcal{C}}^{(1)}_{m,{b\bar{b}{}}}\rangle\right)\right]+\operatorname{\mathcal{O}}\left(\alpha_{s}^{3}\right)$
$\displaystyle={}$
$\displaystyle|\mathcal{M}^{(0)}_{{b\bar{b}{}}}\rangle+\frac{\alpha_{s}}{4\pi}|\mathcal{M}^{(1)}_{{b\bar{b}{}}}\rangle+\left(\frac{\alpha_{s}}{4\pi}\right)^{2}|\mathcal{M}^{(2)}_{{b\bar{b}{}}}\rangle+\operatorname{\mathcal{O}}\left(\alpha_{s}^{3}\right)\,.$
(A.7)

Here, $|\hat{\mathcal{C}}^{(1)}_{m,{b\bar{b}{}}}\rangle$ denotes the amplitude
of the mass counterterm diagrams and
$|\mathcal{M}^{(k)}_{{b\bar{b}{}}}\rangle$ denotes the renormalised amplitude
of order $\alpha_{s}^{k}$.

In contrast to the Born process, the real-emission amplitude for
$H\to{b\bar{b}{g}}$ starts at
$\operatorname{\mathcal{O}}\left(\alpha_{s}\right)$, which means that coupling
renormalisation already affects the one-loop correction. In accordance with
Section 2.3 we leave powers of the strong coupling constant related to real
emissions implicit in the amplitude and write the perturbative expansion of
the bare $H\to{b\bar{b}{g}}$ amplitude as

$\displaystyle|\hat{\mathcal{M}}_{{b\bar{b}{g}}}\rangle$
$\displaystyle=|\hat{\mathcal{M}}^{(0)}_{{b\bar{b}{g}}}\rangle+\frac{\alpha_{s,b}}{4\pi}|\hat{\mathcal{M}}^{(1)}_{{b\bar{b}{g}}}\rangle+\operatorname{\mathcal{O}}\left(\alpha_{s,b}^{3}\right)\,,$
(A.8)

In analogy to Appendix A, the $H\to{b\bar{b}{g}}$ process is renormalised via

$\displaystyle|\mathcal{M}_{{b\bar{b}{g}}}\rangle={}$ $\displaystyle
Z_{m}Z_{\psi}\sqrt{Z_{A}}|\hat{\mathcal{M}}_{{b\bar{b}{g}}}\rangle+\frac{\alpha_{s,b}}{4\pi}Z_{m}^{(1)}|\hat{\mathcal{C}}^{(1)}_{m,{b\bar{b}{g}}}\rangle+\operatorname{\mathcal{O}}\left(\alpha_{s,b}^{3}\right)$
$\displaystyle={}$
$\displaystyle(\mu_{R}^{2})^{\epsilon}S_{\epsilon}|\hat{\mathcal{M}}^{(0)}_{{b\bar{b}{g}}}\rangle+\frac{\alpha_{s}}{4\pi}\left[(\mu_{R}^{2})^{\epsilon}S_{\epsilon}Z_{\alpha_{s}}^{(1)}|\hat{\mathcal{M}}^{(0)}_{{b\bar{b}{g}}}\rangle\right.$
$\displaystyle\left.+(\mu_{R}^{2})^{2\epsilon}S_{\epsilon}^{2}\left(|\hat{\mathcal{M}}^{(1)}_{{b\bar{b}{g}}}\rangle+\left(\frac{Z_{A}^{(1)}}{2}+Z_{m}^{(1)}+Z_{\psi}^{(1)}\right)|\hat{\mathcal{M}}^{(0)}_{{b\bar{b}{g}}}\rangle+Z_{m}^{(1)}|\hat{\mathcal{C}}^{(1)}_{m,{b\bar{b}{g}}}\rangle\right)\right]+\operatorname{\mathcal{O}}\left(\alpha_{s}^{3}\right)$
$\displaystyle={}$
$\displaystyle(\mu_{R}^{2})^{\epsilon}S_{\epsilon}\left(|\mathcal{M}^{(0)}_{{b\bar{b}{g}}}\rangle+\frac{\alpha_{s}}{4\pi}|\mathcal{M}^{(1)}_{{b\bar{b}{g}}}\rangle\right)+\operatorname{\mathcal{O}}\left(\alpha_{s}^{3}\right)\,.$
(A.9)

Note that for the expansion of the renormalised amplitude, we pull out a
global factor of $(\mu_{R}^{2})^{\epsilon}S_{\epsilon}$ and we move it to the
normalisation of the phase-space measure, see Eq. 2.16. Analogously, the
renormalised double-real amplitudes are written as

$\displaystyle|\mathcal{M}_{{b\bar{b}{gg}}}\rangle$
$\displaystyle=(\mu_{R}^{2})^{2\epsilon}S_{\epsilon}^{2}|\mathcal{M}^{(0)}_{{b\bar{b}{gg}}}\rangle+\operatorname{\mathcal{O}}\left(\alpha_{s}^{3}\right)\,,$
(A.10) $\displaystyle|\mathcal{M}_{{b\bar{b}{q\bar{q}}}}\rangle$
$\displaystyle=(\mu_{R}^{2})^{2\epsilon}S_{\epsilon}^{2}|\mathcal{M}^{(0)}_{{b\bar{b}{q\bar{q}}}}\rangle+\operatorname{\mathcal{O}}\left(\alpha_{s}^{3}\right)\,.$
(A.11)

Both for the calculation of the amplitudes and the application of the
subtraction scheme, we choose to work in a theory with $n_{f}=n_{l}+1$ active
flavours. This amounts to renormalising the strong coupling with
$\beta_{0}(n_{l}+1)$ and the gluon field with a non-trivial renormalisation
constant ($\sqrt{Z_{A}}\neq 1$). This procedure takes care of insertions of
heavy-quark loops on gluon into propagators and external gluon legs.

A possible alternative [30] is to formulate the subtraction scheme in a theory
where the strong coupling evolves with $n_{l}$ active flavours. If we
calculate the amplitudes in the renormalisation scheme with $n_{l}+1$ active
flavours described above, we can use the decoupling relation [79, 80, 81, 82,
83, 84, 85]

$\displaystyle\alpha_{s}^{(n_{l})}$
$\displaystyle=\zeta_{\alpha_{s}}\alpha_{s}^{(n_{l}+1)}$ (A.12)

to absorb the effects of the heavy-quark loops into the running of
$\alpha_{s}$. We stress that in this case, the higher-order terms in
$\epsilon$ of the decoupling constant $\zeta_{\alpha_{s}}$ need to be taken
into account. To $\operatorname{\mathcal{O}}\left(\alpha_{s}\right)$ we need
[86]

$\displaystyle\zeta_{\alpha_{s}}$
$\displaystyle=1+\frac{\alpha_{s}^{(n_{l})}}{4\pi}\frac{\beta_{0,Q}}{\epsilon}\left[\left(\frac{\mu_{R}^{2}}{m_{b}^{2}}\right)^{\epsilon}\Gamma(1+\epsilon)e^{\epsilon\gamma_{\mathrm{E}}}-1\right]+\operatorname{\mathcal{O}}\left((\alpha_{s}^{(n_{l})})^{2}\right)\,.$
(A.13)

The subtraction scheme then operates on amplitudes in a theory with $n_{l}$
active flavours. Finally, after all IR poles cancel we can use Eq. A.12 in the
opposite direction to go back to the theory with $n_{l}+1$ active flavours.
The results of this procedure are identical to the results obtained when
working with $n_{l}+1$ active flavours throughout the whole calculation.

##  Appendix B Useful formulae

For the convenience of the reader, in this appendix, we collect useful
formulae that are available in the literature and are used in our calculation.

###  B.1 Anomalous dimensions for IR factorisation

Below, we give explicit expressions for the anomalous dimensions that appear
in the $\mathbf{Z}$ operator. All of these coefficients were taken from Ref.
[30]; they can further be traced to Refs. [45, 50].

We write an expansion of the anomalous dimensions in terms of the strong
coupling constant,

$\displaystyle\gamma_{i}(\alpha_{s})={}$
$\displaystyle\sum_{n=0}^{\infty}\gamma_{i}^{(n)}\left({\frac{\alpha_{s}}{4\pi}}\right)^{n+1}\,,$
(B.1)

where $i$ stands for the type of the anomalous dimension. We report formulae
for the massless and massive cusp anomalous dimensions, heavy-quark and gluon
anomalous dimensions.

The expansion coefficients of the massless cusp anomalous dimension are given
by

$\displaystyle\gamma_{\mathrm{cusp}}^{(0)}$ $\displaystyle=4\,,$
$\displaystyle\gamma_{\mathrm{cusp}}^{(1)}$
$\displaystyle=\left(\frac{268}{9}-\frac{4}{3}\pi^{2}\right)C_{A}-\frac{80}{9}T_{F}n_{l}\,,$
(B.2)

with $n_{l}$ being a number of massless quark flavours. For the cusp anomalous
dimension of massive emitters the expansion coefficients depend on

$\displaystyle v_{ij}$
$\displaystyle=\sqrt{1-\frac{m_{i}^{2}m_{j}^{2}}{({q_{i}}\cdot{q_{j}})^{2}}}$
(B.3)

and read

$\displaystyle\gamma_{\mathrm{cusp,Q}}^{(0)}(v)={}$
$\displaystyle\gamma_{\mathrm{cusp}}^{(0)}\frac{1}{v}\left[\frac{1}{2}\ln\left(\frac{1+v}{1-v}\right)-i\pi\right]\,,$
(B.4) $\displaystyle\gamma_{\mathrm{cusp,Q}}^{(1)}(v)={}$
$\displaystyle\gamma_{\mathrm{cusp}}^{(1)}\frac{1}{v}\left[\frac{1}{2}\ln\left(\frac{1+v}{1-v}\right)-i\pi\right]+8C_{A}\Bigg{\\{}\zeta_{3}-\frac{5}{6}\pi^{2}+\frac{1}{4}\ln^{2}\left(\frac{1+v}{1-v}\right)$
$\displaystyle+\frac{1}{v^{2}}\left[\frac{1}{24}\ln^{3}\left(\frac{1+v}{1-v}\right)+\ln\left(\frac{1+v}{1-v}\right)\left(\frac{1}{2}\mathrm{Li}_{2}\left(\frac{1-v}{1+v}\right)-\frac{5\pi^{2}}{12}\right)\right.$
$\displaystyle\left.+\mathrm{Li}_{3}\left(\frac{1-v}{1+v}\right)-\zeta_{3}\right]+\frac{1}{v}\left[\frac{5}{6}\pi^{2}+\frac{5}{12}\pi^{2}\ln\left(\frac{1+v}{1-v}\right)-\frac{1}{4}\ln^{2}\left(\frac{1+v}{1-v}\right)\right.$
$\displaystyle\left.-\ln\left(\frac{2v}{1+v}\right)\ln\left(\frac{1+v}{1-v}\right)-\frac{1}{24}\ln^{3}\left(\frac{1+v}{1-v}\right)+\mathrm{Li}_{2}\left(\frac{1-v}{1+v}\right)\right]$
$\displaystyle+i\pi\left\\{\frac{1}{v^{2}}\left[\frac{\pi^{2}}{6}-\frac{1}{4}\ln^{2}\left(\frac{1+v}{1-v}\right)-\mathrm{Li}_{2}\left(\frac{1-v}{1+v}\right)\right]\right.$
$\displaystyle\left.+\frac{1}{v}\left[-\frac{\pi^{2}}{6}+2\ln\left(\frac{2v}{1+v}\right)+\ln\left(\frac{1+v}{1-v}\right)+\frac{1}{4}\ln^{2}\left(\frac{1+v}{1-v}\right)\right]\right.$
$\displaystyle\left.-\ln\left(\frac{1+v}{1-v}\right)\right\\}\Bigg{\\}}$ (B.5)

Furthermore, the heavy-quark anomalous dimensions are

$\displaystyle\gamma_{Q}^{(0)}$ $\displaystyle=-2C_{F}\,,$
$\displaystyle\gamma_{Q}^{(1)}$
$\displaystyle=C_{F}C_{A}\left(\frac{2}{3}\pi^{2}-\frac{98}{9}-4\zeta_{3}\right)+\frac{40}{9}C_{F}T_{F}n_{l}\,,$
(B.6)

while for gluons we have

$\displaystyle\gamma_{g}^{(0)}$
$\displaystyle=-\beta_{0}(n_{l})=-\frac{11}{3}C_{A}+\frac{4}{3}T_{F}n_{l}\,.$
(B.7)

###  B.2 Coefficients for on-shell to $\overline{\textrm{MS}}$-scheme
conversion

The coefficients for the conversion relation between the on-shell and
$\overline{\textrm{MS}}$ Yukawa coupling, defined in Eq. 2.33, are given by
[87, 20]

$\displaystyle r_{1}={}$ $\displaystyle-2d_{1}\,,$ $\displaystyle{}r_{2}={}$
$\displaystyle 3d_{1}^{\,2}-2d_{2}\,,$ (B.8)

and

$\displaystyle d_{1}(m_{b},\mu)={}$ $\displaystyle-
C_{F}\left({1+\frac{3}{4}L}\right)\,,$ (B.9) $\displaystyle
d_{2}(m_{b},\mu)={}$ $\displaystyle
C_{F}^{2}\left[{\frac{7}{128}-\frac{3}{4}{\zeta_{3}}+3\ln(2){\zeta_{2}}-\frac{15}{8}{\zeta_{2}}+\frac{21}{32}L+\frac{9}{32}L^{2}}\right]$
$\displaystyle+C_{A}C_{F}\left[{-\frac{1111}{384}+\frac{3}{8}{\zeta_{3}}+\frac{1}{2}{\zeta_{2}}-\frac{3}{2}\ln(2){\zeta_{2}}-\frac{185}{96}L-\frac{11}{32}L^{2}}\right]$
$\displaystyle+C_{F}T_{F}n_{l}\left[{\frac{71}{96}+\frac{1}{2}{\zeta_{2}}+\frac{13}{24}L+\frac{1}{8}L^{2}}\right]+C_{F}T_{F}\left[{\frac{143}{96}-{\zeta_{2}}+\frac{13}{24}L+\frac{1}{8}L^{2}}\right]\,,$
(B.10)

with the abbreviation $L=\ln\left({\mu^{2}/m_{b}^{2}}\right)$.

##  Appendix C Factorisation formulae

Here, we collect the factorisation formulae which are necessary to evaluate
the singular limits of the squared matrix elements. We specialise to the
$H\to{b\bar{b}{g}}$, $H\to{b\bar{b}{gg}}$ and $H\to{b\bar{b}{q\bar{q}}}$
processes. For a useful collection of factorisation formulae for general NNLO
QCD processes, we refer the reader to Ref. [30].

###  C.1 Single-collinear factorisation

In order to discuss the single-collinear limit of the $H\to{b\bar{b}{gg}}$ and
$H\to{b\bar{b}{q\bar{q}}}$ matrix element, we have to define a perpendicular
direction $k_{\perp}^{\mu}$ which determines how the collinear limit is
approached. We choose

$\displaystyle k_{\perp}^{\mu}={}$ $\displaystyle\lim_{\eta\to
0}\frac{\hat{q}_{5}-\hat{q}_{4}}{\|\hat{q}_{5}-\hat{q}_{4}\|}\,,$ (C.1)

where $\eta$ refers to the phase-space parametrisation in Section 2.4. Then
the factorisation formula in the single-collinear limit reads

$\displaystyle C_{45}|\mathcal{M}^{(0)}_{{b\bar{b}{ij}}}|^{2}={}$
$\displaystyle\frac{g_{s}^{2}}{({q_{4}}\cdot{q_{5}})}\langle\mathcal{M}^{(0)}_{{b\bar{b}{g}},\mu}|\mathcal{P}^{(0),\mu\nu}_{ij}(z,k_{\perp};\epsilon)|\mathcal{M}^{(0)}_{{b\bar{b}{g}},\nu}\rangle\,,$
(C.2)

where $ij$ stands for either the $gg$ or $q\bar{q}$ channel and
$|\mathcal{M}^{(0)}_{{b\bar{b}{g}},\nu}\rangle$ is the spin-correlated
$H\to{b\bar{b}{g}}$ amplitude with the gluon polarisation vector removed, i.e.

$\displaystyle|\mathcal{M}^{(0)}_{{b\bar{b}{g}}}\rangle={}$
$\displaystyle(\varepsilon^{\mu}(q_{4},\lambda))^{\ast}|\mathcal{M}^{(0)}_{{b\bar{b}{g}},\mu}\rangle\,.$
(C.3)

The splitting functions $\mathcal{P}^{(0),\mu\nu}_{ij}$ carry Lorentz indices
which are contracted with the corresponding index of the gluon. For the $gg$
and $q\bar{q}{}$ channels they read

$\displaystyle\mathcal{P}^{(0),\mu\nu}_{gg}(z,k_{\perp};\epsilon)={}$
$\displaystyle
2C_{A}\left[{-g^{\mu\nu}\left({\frac{z}{1-z}+\frac{1-z}{z}}\right)-2(1-\epsilon)z(1-z)\frac{k_{\perp}^{\mu}k_{\perp}^{\nu}}{k_{\perp}^{\,2}}}\right]\,,$
(C.4)
$\displaystyle\mathcal{P}^{(0),\mu\nu}_{q\bar{q}}(z,k_{\perp};\epsilon)={}$
$\displaystyle
T_{F}\left[{-g^{\mu\nu}+4z(1-z)\frac{k_{\perp}^{\mu}k_{\perp}^{\nu}}{k_{\perp}^{\,2}}}\right]\,,$
(C.5)

with

$\displaystyle z={}$
$\displaystyle\frac{E_{4}}{E_{4}+E_{5}}=1-\frac{\xi_{2}}{2}\,.$ (C.6)

As the factorisation formula above indicates, we have to evaluate spin-
correlated matrix elements. In the term where they are contracted with the
metric tensor $g^{\mu\nu}$, the result corresponds to an uncorrelated matrix
element via the polarisation sum. This leads to

$\displaystyle\langle\mathcal{M}^{(0)}_{{b\bar{b}{g}},\mu}|(-g^{\mu\nu})|\mathcal{M}^{(0)}_{{b\bar{b}{g}},\nu}\rangle={}$
$\displaystyle|\mathcal{M}^{(0)}_{{b\bar{b}{g}}}|^{2}\,.$ (C.7)

For four-dimensional matrix elements, the product
$k_{\perp}^{\mu}|\mathcal{M}^{(0)}_{{b\bar{b}{g}},\nu}\rangle$ can be linked
to helicity amplitudes using the fact that the perpendicular vector
$k_{\perp}^{\mu}$ may be decomposed as

$\displaystyle k_{\perp}^{\mu}={}$
$\displaystyle-\sum_{\lambda}({\varepsilon(q_{4},\lambda)}\cdot{k_{\perp}})(\varepsilon^{\mu}(q_{4},\lambda))^{\ast}\,.$
(C.8)

This is true since, by construction, $k_{\perp}^{\mu}$ lies in the plane
perpendicular to the momentum of a gluon, $q_{4}$, which is spanned by the
polarisation vectors $\varepsilon^{\mu}(q_{4},\lambda)$.

###  C.2 Single-soft factorisation (tree-level)

In the single-soft limits, the factorisation formulae for the
$H\to{b\bar{b}{g}}$ and $H\to{b\bar{b}{gg}}$ squared matrix elements are given
by

$\displaystyle S_{4}|\mathcal{M}^{(0)}_{{b\bar{b}{g}}}|^{2}={}$
$\displaystyle-
g_{s}^{\,2}\,C_{F}\left({\mathcal{S}^{(0)}_{22,4}-2\mathcal{S}^{(0)}_{23,4}+\mathcal{S}^{(0)}_{33,4}}\right)|\mathcal{M}^{(0)}_{{b\bar{b}{}}}|^{2}\,,$
(C.9) $\displaystyle S_{5}|\mathcal{M}^{(0)}_{{b\bar{b}{gg}}}|^{2}={}$
$\displaystyle-
g_{s}^{\,2}\,\left[{C_{F}\left({\mathcal{S}^{(0)}_{22,5}-2\mathcal{S}^{(0)}_{23,5}+\mathcal{S}^{(0)}_{33,5}}\right)-C_{A}\left({\mathcal{S}^{(0)}_{24,5}+\mathcal{S}^{(0)}_{34,5}-\mathcal{S}^{(0)}_{23,5}}\right)}\right]|\mathcal{M}^{(0)}_{{b\bar{b}{g}}}|^{2}\,,$
(C.10)

where $\mathcal{S}^{(0)}_{ij,k}$ is the usual tree-level eikonal factor

$\displaystyle\mathcal{S}^{(0)}_{ij,k}={}$
$\displaystyle\frac{({q_{i}}\cdot{q_{j}})}{({q_{i}}\cdot{q_{k}})({q_{j}}\cdot{q_{k}})}\,.$
(C.11)

###  C.3 Single-soft factorisation (one-loop)

The soft limit of the one-loop amplitudes with massive quarks has been studied
in Refs. [88, 89]. In our case, it can be written as

$\displaystyle
S_{4}\,2\operatorname{Re}\langle\mathcal{M}^{(0)}_{{b\bar{b}{g}}}|\mathcal{M}^{(1)}_{{b\bar{b}{g}}}\rangle={}$
$\displaystyle-
g_{s}^{\,2}C_{F}\left({\mathcal{S}^{(0)}_{22,4}-2\mathcal{S}^{(0)}_{23,4}+\mathcal{S}^{(0)}_{33,4}}\right)$
$\displaystyle\times\left[{2\operatorname{Re}\langle\mathcal{M}^{(0)}_{{b\bar{b}{}}}|\mathcal{M}^{(1)}_{{b\bar{b}{}}}\rangle+\left({\mathcal{R}^{(1)}_{23,4}+Z_{\alpha_{s}}^{(1)}+Z_{A}^{(1)}}\right)|\mathcal{M}^{(0)}_{{b\bar{b}{}}}|^{2}}\right]\,,$
(C.12)

where in the second line we single out the contribution that comes from the
renormalisation procedure, i.e. it involves terms resulting from the strong-
coupling renormalisation, $Z_{\alpha_{s}}^{(1)}$, as well as a term from the
gluon wave-function renormalisation, $Z_{A}^{(1)}$, see Appendix A.
Furthermore, the functions $\mathcal{R}^{(1)}_{ij,4}$ denote the one-loop
eikonal factor that can be expanded in an $\epsilon$ series as

$\displaystyle\mathcal{R}^{(1)}_{ij,4}={}$ $\displaystyle
4C_{A}\left({\tfrac{1}{2}\mu_{R}^{2}\mathcal{S}^{(0)}_{ij,4}}\right)^{\epsilon}\sum_{k=-2}^{1}\epsilon^{k}R_{k}(q_{i},q_{j};q_{4})\,,$
(C.13)

where functions $R_{k}(q_{i},q_{j};q_{4})$ have been calculated in Ref. [89],
and further simplified in Ref. [90]. In our case we use the formulae from Eq.
(4) of Ref. [90]; we emphasise that the expressions therein correspond to the
unrenormalised one-loop soft-gluon current. This is particularly convenient,
since we perform the renormalisation in a hybrid scheme, as outlined in
Appendix A.

Finally, we note that, for our calculation, we split the one-loop matrix
elements in the factorisation formula, Section C.3, into finite terms and
terms containing poles in $\epsilon^{-1}$ using the $\mathbf{Z}$ operator.
This is reflected in Sections 4.2 and 4.2.

###  C.4 Double-soft factorisation

The relevant factorisation formulae for amplitudes that involve massive
particles can be obtained from Ref. [91, 29], see also Ref. [30]. In general,
the double-soft limit requires single- and double-eikonal factors and colour-
correlated matrix elements. For the $H\to{b\bar{b}{gg}}$ and
$H\to{b\bar{b}{q\bar{q}}}$ matrix elements, the colour-correlated matrix
elements can be expressed explicitly through $\mathrm{SU}(3)$ colour factors
so that we arrive at

$\displaystyle S_{45}|\mathcal{M}^{(0)}_{{b\bar{b}{ij}}}|^{2}={}$
$\displaystyle
g_{s}^{\,4}\,\mathrm{D}{}\mathrm{Soft}^{(0)}_{ij}(q_{2},q_{3};q_{4},q_{5})|\mathcal{M}^{(0)}_{{b\bar{b}{}}}|^{2}\,,$
(C.14)

where $\mathrm{D}{}\mathrm{Soft}^{(0)}_{ij}(q_{2},q_{3};q_{4},q_{5})$ denotes
the double-soft functions for the partons $ij\in\\{gg,q\bar{q}{}\\}$. They are
given by

$\displaystyle\mathrm{D}{}\mathrm{Soft}^{(0)}_{gg}(q_{2},q_{3};q_{4},q_{5})={}$
$\displaystyle
C_{F}^{2}\left({\mathcal{S}^{(0)}_{22,4}-2\mathcal{S}^{(0)}_{23,4}+\mathcal{S}^{(0)}_{33,4}}\right)\left({\mathcal{S}^{(0)}_{22,5}-2\mathcal{S}^{(0)}_{23,5}+\mathcal{S}^{(0)}_{33,5}}\right)$
$\displaystyle-C_{A}C_{F}\
\left({\mathcal{S}^{gg}_{22,45}-\mathcal{S}^{gg}_{23,45}-\mathcal{S}^{gg}_{32,45}+\mathcal{S}^{gg}_{33,45}}\right)\,,$
(C.15)
$\displaystyle\mathrm{D}{}\mathrm{Soft}^{(0)}_{q\bar{q}}(q_{2},q_{3};q_{4},q_{5})={}$
$\displaystyle
C_{F}T_{F}\left({\mathcal{S}^{q\bar{q}{}}_{22,45}-\mathcal{S}^{q\bar{q}{}}_{23,45}-\mathcal{S}^{q\bar{q}{}}_{32,45}+\mathcal{S}^{q\bar{q}{}}_{33,45}}\right)\,,$
(C.16)

where $\mathcal{S}^{(0)}_{ij,k}$ is the usual single-eikonal factor, defined
in Eq. C.11, and $\mathcal{S}^{gg}_{ij,45}$ and
$\mathcal{S}^{q\bar{q}}_{ij,45}$ are the double-eikonal factors. For the $gg$
emission case, we have

$\displaystyle\mathcal{S}^{gg}_{ij,45}={}$
$\displaystyle\mathcal{S}^{(m=0)}_{ij,45}+m_{i}^{\,2}\mathcal{S}^{(m\neq
0)}_{ij,45}+m_{j}^{\,2}\mathcal{S}^{(m\neq 0)}_{ji,45}\,,$ (C.17)

where [91]

$\displaystyle\mathcal{S}^{(m=0)}_{ij,45}={}$
$\displaystyle\frac{(1-\epsilon)}{({q_{4}}\cdot{q_{5}})^{2}}\frac{({q_{i}}\cdot{q_{4}})({q_{j}}\cdot{q_{5}})+({q_{i}}\cdot{q_{5}})({q_{j}}\cdot{q_{4}})}{({q_{i}}\cdot{q_{45}})({q_{j}}\cdot{q_{45}})}$
$\displaystyle-\frac{1}{2}\frac{({q_{i}}\cdot{q_{j}})}{({q_{i}}\cdot{q_{4}})({q_{j}}\cdot{q_{4}})}\frac{({q_{i}}\cdot{q_{j}})}{({q_{i}}\cdot{q_{5}})({q_{j}}\cdot{q_{5}})}\left[{2-\frac{({q_{i}}\cdot{q_{4}})({q_{j}}\cdot{q_{5}})+({q_{i}}\cdot{q_{5}})({q_{j}}\cdot{q_{4}})}{({q_{i}}\cdot{q_{45}})({q_{j}}\cdot{q_{45}})}}\right]$
$\displaystyle+\frac{1}{({q_{4}}\cdot{q_{5}})}\left[{\frac{({q_{i}}\cdot{q_{j}})}{({q_{i}}\cdot{q_{4}})({q_{j}}\cdot{q_{5}})}+\frac{({q_{i}}\cdot{q_{j}})}{({q_{i}}\cdot{q_{5}})({q_{j}}\cdot{q_{4}})}}\right.$
$\displaystyle\left.{-\frac{({q_{i}}\cdot{q_{j}})}{({q_{i}}\cdot{q_{45}})({q_{j}}\cdot{q_{45}})}\left({2+\frac{\left[{({q_{i}}\cdot{q_{4}})({q_{j}}\cdot{q_{5}})+({q_{i}}\cdot{q_{5}})({q_{j}}\cdot{q_{4}})}\right]^{2}}{2({q_{i}}\cdot{q_{4}})({q_{j}}\cdot{q_{5}})({q_{i}}\cdot{q_{5}})({q_{j}}\cdot{q_{4}})}}\right)}\right]$
(C.18)

and [29]

$\displaystyle\mathcal{S}^{(m\neq 0)}_{ij,45}={}$
$\displaystyle-\frac{1}{4({q_{i}}\cdot{q_{4}})({q_{i}}\cdot{q_{5}})({q_{4}}\cdot{q_{5}})}+\frac{1}{2}\frac{({q_{i}}\cdot{q_{j}})}{({q_{i}}\cdot{q_{4}})({q_{i}}\cdot{q_{5}})({q_{j}}\cdot{q_{4}})({q_{j}}\cdot{q_{5}})}\frac{({q_{j}}\cdot{q_{45}})}{({q_{i}}\cdot{q_{45}})}$
$\displaystyle-\frac{1}{2}\frac{1}{({q_{4}}\cdot{q_{5}})({q_{i}}\cdot{q_{45}})({q_{j}}\cdot{q_{45}})}\left({\frac{({q_{j}}\cdot{q_{4}})^{2}}{({q_{i}}\cdot{q_{4}})({q_{j}}\cdot{q_{5}})}+\frac{({q_{j}}\cdot{q_{5}})^{2}}{({q_{i}}\cdot{q_{5}})({q_{j}}\cdot{q_{4}})}}\right)\,.$
(C.19)

In Sections C.4 and C.4 the shorthand $q_{45}=(q_{4}+q_{5})$ is used. For the
case of a soft $q\bar{q}{}$ pair emission, the double-eikonal factor is given
by [91]

$\displaystyle\mathcal{S}^{q\bar{q}{}}_{ij,45}={}$
$\displaystyle\frac{1}{({q_{4}}\cdot{q_{5}})^{2}}\frac{({q_{i}}\cdot{q_{4}})({q_{j}}\cdot{q_{5}})+({q_{i}}\cdot{q_{5}})({q_{j}}\cdot{q_{4}})-({q_{i}}\cdot{q_{j}})({q_{4}}\cdot{q_{5}})}{({q_{i}}\cdot{q_{45}})({q_{j}}\cdot{q_{45}})}\,.$
(C.20)

In the strongly-ordered double-soft limit, $S_{5}S_{45}$, where we take both
$\xi_{1}\to 0$ and $\xi_{2}\to 0$, the double-soft function simplifies to

$\displaystyle
S_{5}\mathrm{D}{}\mathrm{Soft}^{(0)}_{gg}(q_{2},q_{3};q_{4},q_{5})={}$
$\displaystyle
C_{F}^{2}\left({\mathcal{S}^{(0)}_{22,4}-2\mathcal{S}^{(0)}_{23,4}+\mathcal{S}^{(0)}_{33,4}}\right)\left({\mathcal{S}^{(0)}_{22,5}-2\mathcal{S}^{(0)}_{23,5}+\mathcal{S}^{(0)}_{33,5}}\right)$
$\displaystyle+C_{F}C_{A}\left({\mathcal{S}^{(0)}_{22,4}-2\mathcal{S}^{(0)}_{23,4}+\mathcal{S}^{(0)}_{33,4}}\right)\left({\mathcal{S}^{(0)}_{23,5}-\mathcal{S}^{(0)}_{24,5}-\mathcal{S}^{(0)}_{34,5}}\right)\,.$
(C.21)

The corresponding limit for the case of $q\bar{q}$ emission is regular.

Moreover, we need the double-soft single-collinear limit, $S_{45}C_{45}$,
which can be obtained by taking the collinear limit of the double-soft
function. However, a simpler expression arises if we use an iterated
factorisation formula, taking first the collinear ($q_{4}\parallel q_{5}$) and
then the soft limit of the parent parton of the splitting ($q_{45}^{0}\to 0$).
We obtain

$\displaystyle S_{45}C_{45}|\mathcal{M}^{(0)}_{{b\bar{b}{ij}}}|^{2}$
$\displaystyle=\frac{g_{s}^{4}}{({q_{4}}\cdot{q_{5}})}\mathcal{P}^{(0),\mu\nu}_{ij}C_{F}\left({\mathcal{S}^{(0)}_{22,(45),\mu\nu}-2\mathcal{S}^{(0)}_{23,(45),\mu\nu}+\mathcal{S}^{(0)}_{33,(45),\mu\nu}}\right)|\mathcal{M}^{(0)}_{{b\bar{b}{}}}|^{2}\,,$
(C.22)

where we use the shorthand notation

$\displaystyle\mathcal{S}^{(0)}_{ij,(45),\mu\nu}$
$\displaystyle=\frac{q_{i,\mu}\,q_{j,\nu}}{({q_{i}}\cdot{q_{45}})({q_{45}}\cdot{q_{j}})}\,.$
(C.23)

##  Appendix D Integrated subtraction terms

In this section we report formulae for integrated subtraction terms that we
use throughout the calculation. For the convenience of the reader, we also
include results available in the literature.

###  D.1 Single-collinear subtraction terms

The relevant factorisation formula for the single-collinear limit is given in
Eq. C.2. For the integrated subtraction terms, we integrate the splitting
function over the unresolved phase space in $d$ dimensions.

We recall that the splitting functions in Eqs. C.4 and C.5 contain a term
proportional to $k_{\perp}^{\mu}k_{\perp}^{\nu}/k_{\perp}^{\,2}$ which is
contracted with the spin-correlated matrix element
$\langle\mathcal{M}^{(0)}_{{b\bar{b}{g}},\mu}|\mathcal{M}^{(0)}_{{b\bar{b}{g}},\nu}\rangle$.
In contrast to the single-collinear subtraction terms, where the spin
correlations are required to make the subtraction local, the integrated
subtraction term can be averaged over the azimuthal directions of momentum
$q_{5}$. The reduced matrix element depends only on the sum of momenta,
$q_{45}=(q_{4}+q_{5})$, which is independent of the azimuthal direction of
$q_{5}$. Therefore, the integral over $\mathrm{d}\Omega_{5}^{(2-2\epsilon)}$
decouples and yields

$\displaystyle\left(\int\mathrm{d}\Omega_{5}^{(2-2\epsilon)}\right)^{-1}\int\mathrm{d}\Omega_{5}^{(2-2\epsilon)}\frac{k_{\perp}^{\mu}k_{\perp}^{\nu}}{k_{\perp}^{\,2}}={}$
$\displaystyle\frac{1}{2(1-\epsilon)}\left[{g^{\mu\nu}-\frac{n^{\mu}\bar{n}^{\nu}+\bar{n}^{\mu}n^{\nu}}{({n}\cdot{\bar{n}})}}\right]\,,$
(D.1)

with $\bar{n}^{\mu}=n_{\nu}$ and $n^{\mu}=q_{4}^{\mu}/q_{4}^{0}$ in the
collinear limit. To derive this result, we use the fact that the right-hand
side of Eq. D.1 is invariant under rotations in the
$(2-2\epsilon)$-dimensional sphere, that it has to be orthogonal to
$q_{4}^{\mu}$ and that by definition the time components ($\mu=0$ and $\nu=0$)
need to vanish [30]. Note that, due to the Ward identity
$q_{4}^{\mu}|\mathcal{M}^{(0)}_{{b\bar{b}{g}},\mu}\rangle={}0$, the
$n^{\mu}\bar{n}^{\nu}$ and $\bar{n}^{\mu}n^{\nu}$ terms in Eq. D.1 drop out
when contracted with the spin-correlated squared matrix elements.555Note that
this argument is not necessarily valid if more that one gluon is spin-
correlated in the reduced amplitude.

This means that after azimuthal averaging and using Eq. C.7 we replace

$\displaystyle\frac{k_{\perp}^{\mu}k_{\perp}^{\nu}}{k_{\perp}^{\,2}}\langle\mathcal{M}^{(0)}_{{b\bar{b}{g}},\mu}|\mathcal{M}^{(0)}_{{b\bar{b}{g}},\nu}\rangle\longrightarrow{}$
$\displaystyle-\frac{1}{2(1-\epsilon)}|\mathcal{M}^{(0)}_{{b\bar{b}{g}}}|^{2}\,.$
(D.2)

As a result the factorisation formula of Eq. C.2 simplifies to

$\displaystyle C_{45}|\mathcal{M}^{(0)}_{{b\bar{b}{ij}}}|^{2}={}$
$\displaystyle\frac{g_{s}^{2}}{({q_{4}}\cdot{q_{5}})}\langle\mathcal{P}^{(0)}_{ij}(z;\epsilon)\rangle|\mathcal{M}^{(0)}_{{b\bar{b}{g}}}|^{2}\,,$
(D.3)

where $\langle\mathcal{P}^{(0)}_{ij}(z;\epsilon)\rangle$ are the azimuthally
averaged splitting functions and read

$\displaystyle\langle\mathcal{P}^{(0)}_{gg}(z;\epsilon)\rangle={}$
$\displaystyle 2C_{A}\left[{\frac{z}{1-z}+\frac{1-z}{z}+z(1-z)}\right]\,,$
(D.4) $\displaystyle\langle\mathcal{P}^{(0)}_{q\bar{q}}(z;\epsilon)\rangle={}$
$\displaystyle T_{F}\left[{1-\frac{2z(1-z)}{1-\epsilon}}\right]\,.$ (D.5)

To obtain the integrated collinear subtraction term, we perform the integral
over the unresolved phase space, $[\mathrm{d}{}q_{5}]$, using the
parametrisation of Section 2.4 in the collinear limit. We arrive at

$\displaystyle
C_{45}\int\mathrm{d}\Phi_{{b\bar{b}{gg}}}^{E_{4}>E_{5}}(q_{1})|\mathcal{M}^{(0)}_{{b\bar{b}{ij}}}|^{2}={}$
$\displaystyle
g_{s}^{2}\mathcal{P}^{(0)}_{ij,\mathrm{int}}\int\mathrm{d}\Phi_{{b\bar{b}{g}}}(q_{1})\,(\xi_{1}/2)^{-2\epsilon}|\mathcal{M}^{(0)}_{{b\bar{b}{g}}}|^{2}\,,$
(D.6)

with the integrated splitting function $\mathcal{P}^{(0)}_{ij,\mathrm{int}}$
defined as

$\displaystyle\mathcal{P}^{(0)}_{ij,\mathrm{int}}={}$
$\displaystyle(\mu_{R}^{2})^{\epsilon}S_{\epsilon}E_{\mathrm{max}}^{2-2\epsilon}\left(\frac{\xi_{1}}{2}\right)^{2}\int\frac{\mathrm{d}\Omega_{5}^{(2-2\epsilon)}}{2(2\pi)^{3-2\epsilon}}\int_{0}^{1}\frac{\mathrm{d}\eta}{\eta^{\epsilon}}\int_{0}^{1}\mathrm{d}\xi_{2}(\xi_{2}(2-\xi_{2}))^{1-2\epsilon}\frac{\langle\mathcal{P}^{(0)}_{ij}(1-\tfrac{\xi_{2}}{2};\epsilon)\rangle}{({q_{4}}\cdot{q_{5}})}$
$\displaystyle={}$
$\displaystyle\frac{1}{2}(\mu_{R}^{2})^{\epsilon}S_{\epsilon}E_{\mathrm{max}}^{-2\epsilon}\int\frac{\mathrm{d}\Omega_{5}^{(2-2\epsilon)}}{2(2\pi)^{3-2\epsilon}}\int_{0}^{1}\frac{\mathrm{d}\eta}{\eta^{1+\epsilon}}\int_{0}^{1}\mathrm{d}\xi_{2}(\xi_{2}(2-\xi_{2}))^{-2\epsilon}\langle\mathcal{P}^{(0)}_{ij}(1-\tfrac{\xi_{2}}{2};\epsilon)\rangle\,.$
(D.7)

To perform the integral in Section D.1 we use the fact that the integrand is
symmetric under $\xi_{2}\leftrightarrow(2-\xi_{2})$ exchange and hence we can
extend the integration region to $0<\xi_{2}<2$ at at cost of introducing a
factor of $1/2$. Then we obtain

$\displaystyle\mathcal{P}^{(0)}_{gg,\mathrm{int}}={}$
$\displaystyle\left({\frac{\mu_{R}^{2}}{16E_{\mathrm{max}}^{2}}}\right)^{\epsilon}\,\frac{C_{A}}{(4\pi)^{2}}\,\frac{6(2-3\epsilon)}{\epsilon^{2}}\,\frac{e^{\epsilon\gamma_{\mathrm{E}}}\Gamma^{2}(2-2\epsilon)}{\Gamma(4-4\epsilon)\Gamma(1-\epsilon)}\,,$
(D.8) $\displaystyle\mathcal{P}^{(0)}_{q\bar{q},\mathrm{int}}={}$
$\displaystyle\left({\frac{\mu_{R}^{2}}{16E_{\mathrm{max}}^{2}}}\right)^{\epsilon}\,\frac{T_{F}}{(4\pi)^{2}}\,\frac{\epsilon(10-8\epsilon)-4}{\epsilon}\,\frac{e^{\epsilon\gamma_{\mathrm{E}}}\Gamma(1-2\epsilon)\Gamma(2-2\epsilon)}{\Gamma(4-4\epsilon)\Gamma(2-\epsilon)}\,.$
(D.9)

In the case of $\mathcal{P}^{(0)}_{gg,\mathrm{int}}$, it is necessary to also
calculate the integrated splitting function in the $\xi_{2}\to 0$ limit. We
find

$\displaystyle\mathcal{P}^{(0),\mathrm{soft}}_{gg,\mathrm{int}}={}$
$\displaystyle\left({\frac{\mu_{R}^{2}}{4E_{\mathrm{max}}^{2}}}\right)^{\epsilon}\,\frac{C_{A}}{(4\pi)^{2}}\,\frac{2}{\epsilon^{2}}\,\frac{e^{\epsilon\gamma_{\mathrm{E}}}}{\Gamma(1-\epsilon)}\,.$
(D.10)

###  D.2 Single-soft subtraction terms (tree-level)

When considering a single-soft emission of a gluon with momentum $q_{k}$ that
involves partons with momenta $q_{i}$ and $q_{j}$, we encounter the tree-level
eikonal factor $\mathcal{S}^{(0)}_{ij,k}$, defined in Eq. C.11. In the
following we will use the notation

$\displaystyle\hat{q}^{\mu}$ $\displaystyle=\frac{q^{\mu}}{q^{0}}$ (D.11)

to denote a momentum rescaled by its energy component.

The soft function is integrated over the unresolved phase space using the
parametrisation in Section 2.4, where the Born phase space decouples in the
soft limit. Thus, the gluon energy is unconstrained unless we insert some
bound by hand. For simplicity, we keep the integration domain of $\xi_{1}$
unchanged, which corresponds to an upper limit of the energy of
$E_{\mathrm{max}}$. We arrive at

$\displaystyle\mathcal{S}^{(0)}_{ij,\mathrm{int}}={}$
$\displaystyle\int[\mathrm{d}{}q_{4}]\mathcal{S}^{(0)}_{ij,4}$
$\displaystyle={}$
$\displaystyle\left(\frac{\mu_{R}^{2}}{E_{\mathrm{max}}^{2}}\right)^{\epsilon}S_{\epsilon}\int_{0}^{1}\frac{\mathrm{d}\xi_{1}}{\xi_{1}^{1+2\epsilon}}\int\frac{\mathrm{d}\Omega_{4}^{(3-2\epsilon)}}{2(2\pi)^{3-2\epsilon}}\frac{({q_{i}}\cdot{q_{j}})}{({q_{i}}\cdot{\hat{q}_{4}})({\hat{q}_{4}}\cdot{q_{j}})}$
$\displaystyle={}$
$\displaystyle\left(\frac{\mu_{R}^{2}}{E_{\mathrm{max}}^{2}}\right)^{\epsilon}S_{\epsilon}\left[-\frac{1}{2\epsilon}\right]\int\frac{\mathrm{d}\Omega_{4}^{(3-2\epsilon)}}{2(2\pi)^{3-2\epsilon}}\frac{({q_{i}}\cdot{q_{j}})}{({q_{i}}\cdot{\hat{q}_{4}})({\hat{q}_{4}}\cdot{q_{j}})}\,,$
(D.12)