P3H-19-047, TTP19-041 Higgs decay into massive b-quarks at NNLO QCD in the nested soft-collinear subtraction scheme

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. ar X iv :1 91 1. 11 52 4v 1 [ he pph ] 2 6 N ov 2 01 9


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 → bb branching fraction. Although such measurements are often very challenging, due to overwhelming QCD backgrounds, the H → bb 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 → bb 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 O(α 4 s ) in the limit of massless b-quarks [9]. Mass effects at O(α 2 s ) have been estimated using a large momentum expansion [10]. Furthermore, the impact of a separate class of corrections from diagrams arising at O(α 2 s ) 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 → bb decay have been presented [15][16][17][18] with first N 3 LO 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 → bb 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 → HV → bbV . NNLO QCD corrections to this process have already been studied in the limit of massless b-quarks [17,18,21], 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 [17,[25][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 MS Yukawa couplings. In Section 3 and Section 4 we review the NLO QCD calculation of the H → bb 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.

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 [17,[25][26][27]].

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 H(q 1 ) −→ b(q 2 ) +b(q 3 ) , (2.1) with q 2 1 = M 2 H and q 2 2 = q 2 3 = m 2 b . The leading order (LO) partial decay width of the Higgs boson into a pair of b-quarks is where dΦ bb (q 1 ) is the two-particle phase-space volume element for the production of two b-quarks with total momentum q 1 , and |M (0) bb | 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 where N c = 3 is the number of colours, β = 1 − 4m 2 b /M 2 H and y b stands for the b-quark Yukawa coupling, 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 In Eq. (2.4), α s is the 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 Lorentzinvariant phase space of particles involved in a particular (sub)process. Similar to the notation in Ref. [25], we define F LM (bbX) = dΦ bbX (q 1 )|M (0) bbX | 2 F kin (bbX) , (2.5) where bbX denotes the constituents of the final state of a considered subprocess, dΦ bbX (q 1 ) and |M (0) bbX | 2 stand for the Lorentz-invariant phase-space measure and the squared amplitude of the H → bbX process, respectively. The momentum q 1 refers to the initial state Higgs boson, while F kin is an infrared-safe observable that depends on the kinematical configuration of the particles involved in the process. We will use the notation A = dΦA (2.6) to denote the integration of some quantity A over the phase space, dΦ.

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 = (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 crosssection, i.e., the combination of squared matrix elements and phase-space measure, scales as E −1+a i dE i in soft limits and as θ −1+b ij dθ ij in collinear limits of massless partons, where E i is the energy of the soft parton, θ ij is the angle between the collinear partons and a, b ∈ R. Thus, 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 −1 and to regulate the integrals in such a way that the expansion in 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 [17,[25][26][27] which is closely related to the sectorimproved 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 where O is an operator which extracts the asymptotic behaviour of F and the phase space in a singular limit. The term OF , which is integrated over the unresolved phase space in d dimensions, then carries explicit poles in −1 , while the regulated term F − OF is expanded in 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. Introduce subtraction terms for the soft limits. At NNLO this involves up to two single-soft limits and the double-soft limit.
2. Introduce a partition of unity for the phase space that isolates the collinear limits, 1 = {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. In each collinear partition, choose a suitable phase-space parametrisation in terms of angles and energies of the partons that can become unresolved.
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 → bb 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 −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 phasespace 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 → bb process in greater detail.

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 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 −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 where |M is the UV-renormalised amplitude, Z is an operator in colour space which contains poles in −1 and |F is a finite remainder which does not contain any poles. Expanding all pieces in the strong coupling, we find (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 Eq. (2.9) as a prediction in order to check the −1 poles of the UV-renormalised amplitudes. On the other hand, we can also use Eq. (2.9) to define the finite remainders, i.e.
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 −1 poles in the definition of the Z operator, cf. Refs. [45,50]. In general, Z is an operator acting on vectors in colour space, which expresses nontrivial correlations between different colour configurations [50,51]. 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 → bb and H → bbg). The expansion coefficients Z (k) bb for the H → bb process are given by where γ (i) cusp and γ (0) cusp,Q (v 23 ) are the massless and massive cusp anomalous dimensions, is the zeroth-order coefficient of the QCD β-function with n l massless flavours, β 0 (n l ) = 11 3 C A − 4 3 T F n l and β 0, Q denote the expansion coefficients of the anomalous dimensions of the massive quark and µ R is the renormalisation scale. We collect the necessary formulae in Appendix B.1. For the H → bbg process, where we only need the one-loop amplitude, we find The gluon anomalous dimension γ (0) g is also given in Appendix B.1.

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 → bb) and we include final states with one additional gluon (H → bbg) as well as two additional massless partons (H → bbgg or H → bbqq) 1 . 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 space-time dimensions reads where we denote the parton momentum by q i and its energy by E i . Note that we also include a global factor (µ 2 R ) S that originates from the strong coupling renormalisation, see discussion below Eq. (A.9). We do not introduce an upper bound on the energy of the emitted gluon since it naturally appears due to the energy-momentum conserving δ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 (2.17) and discuss its phase-space parametrisation. The phase-space measure reads where dΦ bb (Q) stands for the Born phase space of the two b-quarks with total momentum Q, andq 4 determines the direction of the gluon momentum, We further parametrise the gluon energy as This finally leads us to (2.21) The dΩ 4 element denotes the angular integral over the direction ofq 4 . The remaining angular integrals can be performed using For a single gluon emission, the only unresolved limit is the single-soft one, i.e. ξ 1 → 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 dΦ bb (q 1 − q 4 ) by dΦ bb (q 1 ) in Eq. (2.21).
Double-emission phase space We now focus on the parametrisation of the H → bbgg 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 qq emission phase space is nearly identical. We comment on the differences where necessary.
We consider the process We denote the sum of the gluon momenta by q 45 = q 4 + q 5 . The phase space measure then reads where the vectorsq 4 andq 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 which leads to the split 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 ↔ q 5 . 2 Hence, we parametrise the gluon energies as [31] where E 45,max is to be chosen such that the integration ranges ξ 1 ∈ [0, 1] and ξ 2 ∈ [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

28)
2 Note that thanks to the symmetry between the gluons we have dΦ bbgg (q 1 ) = 2 dΦ E4>E5 bbgg (q 1 ). However when considering a qq emission such a simplification may only be used if the symmetry q ↔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.
whereq 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 (ξ 1 = E 45 /E 45,max ) and their ratio (ξ 2 = 2E 5 /E 45 ).
We also explicitly parametrise the angle θ 45 between the two gluons as follows The first step in the phase-space construction is to choose a direction forq 4 . Here, we explicitly parametrise angle θ 4 between the emission and theẑ-axis. Next, we fix the direction ofq 5 relative toq 4 using the angle θ 45 between them and the angle φ which is the azimuthal angle ofq 5 around the direction ofq 4 . Given the directionsq 4 andq 5 as well as ξ 2 , we calculate the vectorq 45 viā 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 The unparametrised angles in dΩ 4 and dΩ 5 can be integrated in the end using Eq. (2.22). The parametrisation shown in Eq. (2.31) achieves the desired decoupling of the emission phase space in unresolved limits. In the collinear (η → 0) and single-soft (ξ 2 → 0) limits the energy bound E 45,max simplifies to E max = 1 2 β 2 M H and becomes independent of η and ξ 2 . In the single-soft limit the momentum q 5 decouples from the energy-momentum conserving δ-function and the integrations over η and ξ 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 ξ 2 . Again, the integrations over η and ξ 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 ξ 1 , ξ 2 and η can be carried out for a fixed Born configuration.

Pole vs. 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 1. It has also been shown there that these large logarithms can be avoided by reexpressing the result in terms of the MS Yukawa coupling evaluated at the renormalisation scale µ = M H .
The Yukawa coupling in the MS scheme y b (µ) is related to the MS mass m b (µ) via where µ is the renormalisation scale. Thus, the relation between the on-shell and MS Yukawa couplings can be deduced from the corresponding relation between the masses, which we need up to O (α 2 s ) [52,53]. It reads where the coefficients r i (m b , µ) are presented in Appendix 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 MS Yukawa coupling are denoted with a bar, i.e.
where the expansion coefficients in the two schemes are related by As discussed before, the large logarithmic corrections to the total decay width can be mitigated by reexpressing the result in terms of the running 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.

H → bb decay at NLO
In this section we briefly describe the calculation of the NLO QCD corrections to the H → bb 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.

Real contribution
At NLO we consider one real emission in addition to the Born process, which means that we need to integrate the function over the bbg 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 where ξ 1 refers to the parametrisation of Eq. (2.21). We define the operator to act on all quantities to the right of the S 4 symbol, extracting the leading asymptotic behaviour in ξ 1 of the quantity A if the actual limit does not exist. Denoting the identity operation by I, we can immediately write 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 ξ 1 in S 4 F LM (bbg) is the leading behaviour ξ −1−2

1
. Thus, the ξ 1 integration becomes trivial and the soft singularity manifests itself as an explicit −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 Here, we use the factorisation formula for the soft limit as discussed in Appendix C.2, where also the eikonal factors S (0) ij,k are defined. Moreover, we have the unresolved contribution, which contains the integrated subtraction term and reads with the integrated eikonal factors S (0) ij,int given in Appendix 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.

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 −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 where the term with the Z operator contains an explicit −1 pole, while the second term is finite. As a shorthand we introduce Accordingly, we define two contributions to the virtual correction: the virtual finite contribution and the virtual unresolved contribution The expansion coefficients of the Z operator are given in Eq. (2.13).

Pole cancellation
At this point we can combine all contributions that enter the NLO calculation. We remind the reader that the dΓ F R (bbg) and dΓ F V (bb) terms are free of −1 poles, while the dΓ U R (bbg) and dΓ U V (bb) terms feature −1 poles. Expanding the explicit results for the real unresolved contribution from Eq. (3.10) Analogously, the virtual unresolved contribution, defined in Eq. (3.9), yields Note that the constant term in is absent since our definition of the Z operator of Eq. (2.8) includes only −1 poles. Obviously, the poles of the two contributions cancel. A crucial part of the argument is to notice that the dΓ U R (bbg) part, after integrating out the real emission, is a function of a Born-like phase-space configuration, as is the dΓ U V (bb) contribution. This allows us to demonstrate the pole cancellation for any point of the dΦ bb phase space and without specifying the explicit form of the LO matrix element.

H → bb decay at NNLO
We now consider the NNLO QCD corrections to the H → bb 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 → bbgg and H → bbqq) or an additional pair of b-quarks (H → bbbb); • the real-virtual contribution (RV) -where we consider one-loop virtual corrections to the process H → bbg; • the double-virtual contribution (VV) -where we consider two-loop virtual corrections to the Born process H → bb.
Except for the H → bbbb 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 dimensions. Apart from these divergent contributions we distinguish the subprocess H → bbbb which enters the calculation at O (α 2 s ). 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 dΦ bbbb . 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 → bb decay that appears at second order of perturbation theory is related to Feynman diagrams where the H → bb or H → bbg 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 spinorhelicity formalism. 3 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 −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.

Double-real contribution
In this section we only focus on the H → bbgg and H → bbqq subprocesses. The H → bbbb process is completely finite on its own and does not require any regularisation procedure. For the record we write 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. ξ 2 → 0; • the double-soft limit(S 45 ) -where the energies of both additional partons vanish at a similar rate, i.e. ξ 1 → 0; • the collinear limit (C 45 ) -where the momenta of the two additional partons become collinear to each other, i.e. η → 0.
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 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, To derive Eq. (4.3), 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 Eq. (4.3) we start with the collinear limit, where appropriate, having in mind the simplicity of the corresponding factorisation formulae. We subdivide Eq. (4.3) 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 Eq. (4.3) represents the fully regulated doublereal contribution. Therefore, it can be evaluated in four dimensions using standard numerical techniques. We write and similarly for the qq emission 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 qq emission we do not subtract the single-soft limit (S 5 ) since this limit is not singular in case of a g → qq 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 Eq. (4.3) are regulated in the doublesoft 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 [dq 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 For the qq emission, we obtain where we use the fact that in the collinear limit any infrared-safe observable must be symmetric under q ↔q, which leads to the factor 2 in Eq. (4.7). To get from the first to the second lines of Eqs. (4.6) and (4.7), we use the factorisation formulae given in Appendix C and integrate them over the unresolved phase space [dq 5 ], which yields the integrated splitting functions P Finally, the last term in Eq. (4.3) is evaluated in the double-soft limit. In this case, both partons are unresolved and we call this contribution double-unresolved. We find 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 DSoft

Real-virtual contribution
We need to consider the H → bbg amplitude at one-loop level and integrate its product with the tree-level amplitude for the H → bbg process over the phase space dΦ bbg . In the following, we use the shorthand 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 (4.11) We split the one-loop amplitude into a singular part, containing all −1 poles from loop integrals, and a finite remainder, see Section 2.3. We write where the coefficient of the Z operator is given in Eq. (2.15) 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 → bbg matrix element. It is free of −1 poles and thus we refer to it as the real-virtual finite contribution. It reads 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 αs are given in Appendix A. The notation [. . .] 0 indicates that the O ( 0 ) term of the expression between the brackets should be taken. We use F fin LV (bbg) in analogy to the definition in Eq. (3.7).
The corresponding term containing the soft-regulated Z operator for the H → bbg process is called real-virtual single-unresolved contribution and is given by where the subscript "poles" refers to taking only the pole terms up to O ( −1 ) into account. The remaining two contributions are evaluated in the soft limit and are integrated over the unresolved phase space [dq 4 ]. The integrated subtraction term containing the finite remainder carries the superscript FR and reads ij,int correspond to integrated eikonal factors discussed in Appendix D.2. Finally, the contribution containing the soft limit of the Z operator integrated over the unresolved phase space is referred to as the real-virtual double-unresolved contribution.
We obtain The integrated one-loop soft function R int is defined and calculated in Appendix D.3.

Double-virtual contribution
For the H → bb process, the virtual corrections are described by a single form factor, In Eq. (4.17), |M (0) bb is the tree-level H → bb amplitude, and F s is the scalar heavy-quark form factor, which can be computed perturbatively as The expansion coefficients F (k) s depend on Higgs and b-quark masses, the renormalisation scale µ R and the regulator . The two-loop heavy-quark form factor was computed in Ref. [65] up to 0 terms and to higher orders in in Refs. [68] and [64]. In our calculation we use the results of Ref. [64].
We again use the Z operator to split the double-virtual correction into divergent terms, that contain all −1 poles, and a finite remainder. In this case, given the functions F Furthermore, we distinguish two contributions that contain all −1 poles. First, we have a term that is proportional to the finite remainder of the one-loop H → bb matrix element. It is given by The other contribution containing −1 poles comes with the Born H → bb matrix element; it reads The expansion coefficients of the Z operator are given in Eqs. (2.13) and (2.14).

Pole cancellation
After collecting formulae for all NNLO contributions in the preceeding subsections, we demonstrate pole cancellation. We begin by considering dΓ FR = dΓ FR RV + dΓ FR VV . From the real-virtual contribution, Eq. (4.15), we obtain the pole term 22) and the explicit expansion of the double-virtual contribution, Eq. (4.20), yields

(4.24)
A similar expansion holds for the real-virtual single-unresolved contribution, Eq. (4.14), As before, the poles of the two contributions cancel.  Finally, we turn to the double-unresolved contribution, given by dΓ DU = dΓ DU RR +dΓ DU RV + dΓ DU VV . 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 and further subdivide the pole coefficients according to colour factors, 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.

Top-quark contribution to the H → bb decay
An additional contribution that enters the H → bb decay process at O(α 2 s ) 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. 1a and 1b and their contributions to the decay rate are separately finite. This contribution to the total H → bb width has been computed in Refs. [11,12] as an expansion in powers of x t = (M 2 H /m 2 t ). 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] and the coefficients f S,k 2 are given in Eq. (2) of Ref. [11]. The factor of β −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 O(α 2 s ), changes due to decoupling, Eq. (A.12), or due to translating the on-shell Yukawa coupling to the MS scheme, Eq. (2.33), are of higher order in α s and are discarded.
The real-virtual amplitude related to the diagram in Fig. 1b is computed using the same techniques as the remaining one-loop amplitudes. The double-virtual amplitude of Fig. 1a is obtained from the result of Ref. [11], see Eq.

Results
In this section, we summarise our calculation and present results for the total decay width of the H → bb 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 α s (M Z ) = 0.1181 with M Z = 91.1876 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 GeV [1]. As a starting point for the b-quark mass, we use m b (µ = m b ) = 4.18 GeV [58]. From this we calculate the pole mass m b = 4.78 GeV, using the two-loop matching formula at a matching scale of µ = m b (m b ) implemented in the package RunDec. For the value of the Fermi constant we use G F = 1.166378 × 10 −5 GeV −2 [58]. The central renormalisation scale is taken to be equal to the mass of the Higgs boson, µ R = M H , and for the scale uncertainty estimation we vary it by a factor of 1/2 and 2.

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 At NLO we combine real and virtual corrections 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 dΓ bb δNNLO = dΓ VV + dΓ RV + dΓ RR + dΓ bbbb + dΓ y b yt bb .

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

Total width of the H → bb decay at NNLO
We start by considering the total width of the decay H → bb; this implies using F kin (bbX) = 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 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 2 b /M 2 H . 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 → bb decay rate if we neglect the top-loop mediated contribution, discussed in Ref. [  The results for the NLO and NNLO coefficients of the total decay width split into independent colour structures. The renormalisation scale is set to µ 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. 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 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 2 b /M 2 H ) ≈ 0.00146. 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, (α s /π), 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 MS Yukawa coupling.
Total width using the MS Yukawa coupling We now express the total decay width in terms of the MS coupling, as defined in Eq. (2.34). We use m b (µ = m b ) = 4.18 GeV as an input parameter [58] and evolve the MS mass with the RunDec package [69,70]. We obtain The evolution is performed with n f = 5 active flavours at two-loop order. We use these values to evaluate the 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].   [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, Γ bb . 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 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, µ 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 → bb 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 GeV [20]. The NNLO coefficient corresponding to top-quark induced contribution, Eq. (4.29), is then γ y b yt 2 = 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. 1b.
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 µ R = M H /2 and at the sub-permill level for the other scales. Note that for γ bb 2 (µ R = M H /2) there are large cancellations when converting from the on-shell to the MS-scheme Yukawa coupling.

Jet rates for H → bb 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 cut = 0.01 and y 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 where "obs" denotes a generic observable and the MS quantities are related to their counterparts in the pole scheme via γ bb 2 (obs) = γ bb 2 (obs) + r 1 γ bb 1 (obs) + r 2 γ bb 0 (obs) . (5.10) Note that the sum of expansion coefficients of all relevant jet multiplicities, γ bb i (njet), yields the expansion coefficient of the total decay width γ bb i . We keep using the setup presented in Section 5.2 and report results for the two-, threeand 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. 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 cut = 0.01. Similar agreement is observed for y cut = 0.05 with the exception of γ bb 2 (2jet, µ R = M H ), which however suffers from large cancellations in the course of the on-shell to MS-scheme conversion, and the four-jet rates.
Finally, in Table 5 we present the total jet rates for the central choice of the renormalisation scale, µ 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 cut = 0.05, as expected.   Table 5: The total jet rates at NNLO for µ R = M H using the Durham clustering algorithm with y cut = 0.01 and y cut = 0.05, including the y b y t contribution.
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.

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. [17,[25][26][27].
A complete discussion of all necessary NNLO contributions was presented in Section 4.
In particular, we demonstrated cancellation of all −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 [17,18,21], a thorough study of b-quark mass effects is an interesting topic for future research.

A Renormalisation
The subtraction scheme that we apply in this calculation is formulated in terms of UVrenormalised 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 MS scheme with five active flavours for the strong coupling constant. We write 4 Due to the choice to enforce the relation y b = m b (2 √ 2G 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] with β 0 (n l + 1) = 11 3 C A − 4 3 T F (n l + 1). The mass and field renormalisation constants for the massive quark, Z ψ and Z m , are given by [76] The gluon field renormalisation constant is non-vanishing in the on-shell scheme due to heavy-quark loops and reads [77,78] with β 0,Q = − 4 3 T F . For single-virtual and double-virtual amplitudes of the H → bb 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 and together with the LSZ factors we obtain the renormalised amplitude as Here, |Ĉ m,bb denotes the amplitude of the mass counterterm diagrams and |M (k) bb denotes the renormalised amplitude of order α k s . In contrast to the Born process, the real-emission amplitude for H → bbg starts at O (α s ), which means that coupling renormalisation already affects the one-loop correction. In accordance with Eq. (2.9) 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 → bbg amplitude as In analogy to Eq. (A.7), the H → bbg process is renormalised via Note that for the expansion of the renormalised amplitude, we pull out a global factor of (µ 2 R ) S 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 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 β 0 (n l +1) and the gluon field with a non-trivial renormalisation constant ( √ Z A = 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] α (n l ) s = ζ αs α (n l +1) s (A. 12) to absorb the effects of the heavy-quark loops into the running of α s . We stress that in this case, the higher-order terms in of the decoupling constant ζ αs need to be taken into account. To O (α s ) we need [86] 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.

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 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, 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 with n l being a number of massless quark flavours. For the cusp anomalous dimension of massive emitters the expansion coefficients depend on and read Furthermore, the heavy-quark anomalous dimensions are while for gluons we have

B.2 Coefficients for on-shell to MS-scheme conversion
The coefficients for the conversion relation between the on-shell and MS Yukawa coupling, defined in Eq. (2.33), are given by [20,87] and with the abbreviation L = ln (µ 2 /m 2 b ).

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 → bbg, H → bbgg and H → bbqq 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 → bbgg and H → bbqq matrix element, we have to define a perpendicular direction k µ ⊥ which determines how the collinear limit is approached. We choose where η refers to the phase-space parametrisation in Eq. (2.31). Then the factorisation formula in the single-collinear limit reads where ij stands for either the gg or qq channel and |M bbg,ν is the spin-correlated H → bbg amplitude with the gluon polarisation vector removed, i.e.
The splitting functions P (0),µν ij carry Lorentz indices which are contracted with the corresponding index of the gluon. For the gg and qq channels they read 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 µν , the result corresponds to an uncorrelated matrix element via the polarisation sum. This leads to For four-dimensional matrix elements, the product k µ ⊥ |M bbg,ν can be linked to helicity amplitudes using the fact that the perpendicular vector k µ ⊥ may be decomposed as This is true since, by construction, k µ ⊥ lies in the plane perpendicular to the momentum of a gluon, q 4 , which is spanned by the polarisation vectors ε µ (q 4 , λ).

C.2 Single-soft factorisation (tree-level)
In the single-soft limits, the factorisation formulae for the H → bbg and H → bbgg squared matrix elements are given by ij,k is the usual tree-level eikonal factor . (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 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 αs , as well as a term from the gluon wave-function renormalisation, Z A , see Appendix A. Furthermore, the functions R (1) ij,4 denote the one-loop eikonal factor that can be expanded in an series as 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, Eq. (C.12), into finite terms and terms containing poles in −1 using the Z operator. This is reflected in Eqs. (4.13) and (4.14).

C.4 Double-soft factorisation
The relevant factorisation formulae for amplitudes that involve massive particles can be obtained from Ref. [29,91], see also Ref. [30]. In general, the double-soft limit requires single-and double-eikonal factors and colour-correlated matrix elements. For the H → bbgg and H → bbqq matrix elements, the colour-correlated matrix elements can be expressed explicitly through SU(3) colour factors so that we arrive at where DSoft (0) ij (q 2 , q 3 ; q 4 , q 5 ) denotes the double-soft functions for the partons ij ∈ {gg, qq}. They are given by and [29] S (m =0) ij,45 .
(C. 19) In Eqs. (C. 18) and (C.19) the shorthand q 45 = (q 4 + q 5 ) is used. For the case of a soft qq pair emission, the double-eikonal factor is given by [91] S qq ij,45 = In the strongly-ordered double-soft limit, S 5 S 45 , where we take both ξ 1 → 0 and ξ 2 → 0, the double-soft function simplifies to The corresponding limit for the case of qq 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 q 5 ) and then the soft limit of the parent parton of the splitting (q 0 45 → 0). We obtain where we use the shorthand notation ij, (45),µν = q i,µ q j,ν (q i · q 45 )(q 45 · q j ) .
(C. 23) 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 µ ⊥ k ν ⊥ /k 2 ⊥ which is contracted with the spin-correlated matrix element M bbg,ν . 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 dΩ withn µ = n ν and n µ = q µ 4 /q 0 4 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 )dimensional sphere, that it has to be orthogonal to q µ 4 and that by definition the time components (µ = 0 and ν = 0) need to vanish [30]. Note that, due to the Ward identity q µ 4 |M (0) bbg,µ = 0, the n µnν andn µ n ν terms in Eq. (D.1) drop out when contracted with the spin-correlated squared matrix elements. 5 This means that after azimuthal averaging and using Eq. (C.7) we replace As a result the factorisation formula of Eq. (C.2) simplifies to where P (0) ij (z; ) are the azimuthally averaged splitting functions and read To obtain the integrated collinear subtraction term, we perform the integral over the unresolved phase space, [dq 5 ], using the parametrisation of Eq. (2.31) in the collinear limit. We arrive at with the integrated splitting function P ij,int defined as To perform the integral in Eq. (D.7) we use the fact that the integrand is symmetric under ξ 2 ↔ (2 − ξ 2 ) exchange and hence we can extend the integration region to 0 < ξ 2 < 2 at at cost of introducing a factor of 1/2. Then we obtain In the case of P (0) gg,int , it is necessary to also calculate the integrated splitting function in the ξ 2 → 0 limit. We find . (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 S ij,k , defined in Eq. (C.11). In the following we will use the notation to denote a momentum rescaled by its energy component. The soft function is integrated over the unresolved phase space using the parametrisation in Eq. (2.21), 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 ξ 1 unchanged, which corresponds to an upper limit of the energy of E max . We arrive at where the pole in the last line arises from performing the energy integral. We are left with angular integrals only, for which we write Here, we introduce the auxiliary function . (D.14) We write it as a Laurent series in , i.e.
where the coefficients I (k) (q i , q j ) depend on the momenta q i and q j , in particular on whether they are massless or massive. All necessary coefficients for the single-soft integrated subtraction terms, except for the O ( 2 ) terms, have been obtained in Ref. [92], see Appendix A therein. For completeness, we collect below the formulae for those I (k) (q i , q j ) which are relevant to our calculation. We write the full integrated single-soft eikonal factor as We have to distinguish the case where the emitters are two massive particles, which we denote by a subscript MM, and the case where emitter i is massive and emitter j is massless, which we denote by M0. For our calculation, we need the coefficients I (k) MM (q i , q j ) and I (k) M0 (q i , q j ) for k ∈ {−1, 0, 1} when the directions of the momenta q i and q j are arbitrary. Moreover, we need I (k) MM (q i , q j ) for k ∈ {−1, 0, 1, 2} in two special cases: for the case where the two momenta are equal, q i = q j and for the case where q i and q j are in a back-to-back configuration.
Two massive emitters The reported formulae correspond to the result outlined in Eqs. (A.41) to (A.51) of Ref. [92]. We consider two time-like momenta q 2 i = m 2 i and q 2 j = m 2 j . With the definition of v ij from Eq. (B.3) and the notation we introduce the shorthand notations Furthermore, we need the following arguments which will be used in the function With these abbreviations, the coefficients of the massive-massive angular integral read One massive and one massless emitter The reported formulae correspond to the result outlined in Eqs. (A.22) to (A.24) of Ref. [92]. We consider one time-like momentum q 2 i = m i 2 and one light-like momentum q 2 j = 0. We define the symbol Then the coefficients of the massive-massless angular integral are Two massive back-to-back emitters We consider the special case of two massive emitters with the same mass, q 2 i = m 2 and q 2 j = m 2 , arranged in a back-to-back configuration, q i = − q j . We then have E i = E j = M H /2 and, therefore, κ = β. The expansion coefficients are given by MM,b2b (q i , q j ) = where the subscript "b2b" indicates a back-to-back configuration of the emitters. Formulae for I MM,b2b was calculated independently; this O( 2 ) term is needed in the soft limit of the real-virtual contribution.
Self-correlated massive emitter We consider the special case of a self-correlated massive emitter, i = j, for a time-like momentum q i with q 2 i = m 2 . Then we have The O ( 2 ) term was calculated independently.

D.3 Single-soft subtraction terms (one-loop)
We consider the soft limit of the one-loop H → bbg amplitude which is given in Eq. (C.12). This factorisation involves both tree-level and one-loop soft functions. The integration of the tree-level eikonal factors is discussed in Appendix D.2, while in this section we focus on the integration of those terms in the one-loop soft function of Eq. (C.12) which contain R (1) ij,4 . We integrate the one-loop eikonal factors over the soft-gluon phase space, i.e. appears, we restrict ourselves to the case of two emitters with the same non-vanishing mass in a back-to-back configuration. The soft-gluon phase space is parametrised using Eq. (2.21). Note that even though the soft-gluon momentum factorises from the energy-momentum conserving δ-function, we restrict ourselves to the same integration region as stated in Eq. (2.21); this is in accordance with the choice made for the integrated tree-level eikonal factor. The only non-trivial integral is the integration over the angle θ between the momentum of the soft gluon q 4 and the b-quark momentum q i . We reexpress this angle in terms of the variable λ = cos θ and use the symmetry under λ → −λ to restrict the domain of integration to λ ∈ [0, 1]. We arrive at The energy dependence factorises from the integrand, which means that the integral over dξ can be solved trivially. The angular integral dΩ is performed using Eq. (2.22). The expressions for the expansion coefficients R k (q i , q j ; q 4 ) in Ref. [90] are given in terms of the variables α i = m 2 i 2 (q j · q 4 ) (q i · q 4 )(q i · q j ) = (1 − β 2 )(1 + βλ) 2(1 + β 2 )(1 − βλ) , (D.38) α j = m 2 j 2 (q i · q 4 ) (q j · q 4 )(q i · q j ) = (1 − β 2 )(1 − βλ) 2(1 + β 2 )(1 + βλ) , (D. 39) and contain at most classical polylogarithms with arguments composed of these variables. We rewrite those special functions in terms of iterated integrals of argument λ over the alphabet As the rational coefficients in front of the iterated integrals also only contain these letters, the integration over λ can again be performed in terms of iterated integrals over the same alphabet. These iterated integrals are then evaluated at 1 and we rewrite them in terms of harmonic polylogarithms [66] of argument β with up to weight four. All manipulations of the iterated integrals are performed with the Mathematica package HarmonicSums [66,[93][94][95][96][97][98][99][100][101][102]. Note that the term R int is again a Laurent series in . Therefore, we write The analytical expressions can be found in computer-readable form in an ancillary file supplied with this article.

D.4 Double-soft subtraction terms
For the integrated subtraction terms of the double-real contribution we have to integrate the double-soft function, discussed in Appendix C.4, over the unresolved phase space in d dimensions. We define the integrated double-soft function as  We use the phase-space parametrisation outlined in Eq. (2.31), taking into account that the gluon momenta decouple from the overall energy-momentum conserving δ-function. Note that the decoupling of soft particles means that, in principle, their energies are unbounded unless we introduce some constraint by hand. For simplicity, we choose to keep the bound that explicitly appears in the formulation of the phase-space measure, i.e. we keep E 45,max unchanged while the dξ 1 dξ 2 integration region still covers the unit square.
It is particularly useful to keep the Born configuration fixed so that the b-quark points along theẑ-axis, i.e. q 2 = 1 2 M H t µ + βẑ µ and q 3 = 1 2 M H t µ − βẑ µ , (D. 44) wheret µ andẑ µ are unit vectors along time andẑ axes, respectively. After expressing the integrand using the phase-space parametrisation, the ξ 1 dependence factorises and can be integrated analytically. The remaining integrals are performed numerically. Note that the highest pole that occurs in the gg channel is −3 . This arises from taking all singular limits in the integrand, i.e. the strongly-ordered, collinear and the double-soft limits. In the qq case we only have an −2 pole since the strongly-ordered soft limit, S 5 , is regular. Furthermore, we split the C (k) ij coefficients into separate colour structures, We present numerical values for the integrated double-soft function coefficients, evaluated for m b = 4.78 GeV and M H = 125.09 GeV, in Table 6.