Differential decay rates of CP-even and CP-odd Higgs bosons to top and bottom quarks at NNLO QCD

We consider the decay of a neutral Higgs boson of arbitrary CP nature to a massive quark antiquark pair at next-to-next-to-leading order in perturbative QCD. Our analysis is made at the differential level using the antenna subtraction framework. We apply our general set-up to the decays of a CP-even and CP-odd heavy Higgs boson to a top-quark top-antiquark pair and to the decay of the 125 GeV Higgs boson to a massive bottom-quark bottom-antiquark pair. In the latter case we calculate, in particular, the two-jet, three-jet, and four-jet decay rates and, for two-jet events, the energy distribution of the leading jet.


Introduction
The detailed investigation of the production and decay modes of the 125 GeV Higgs boson and the search for new Higgs resonances remain among the major research topics at the Large Hadron Collider (LHC). The interactions of the 125 GeV Higgs boson, as investigated by the ATLAS and CMS experiments, are in accord with the predictions of the Standard Model (SM), and considerable improvement of the experimental precision is expected with the future high luminosity LHC program. (For a recent overview see, for instance, [1].) Precision Higgs physics and searches for new spin-zero resonances are also key issues for putting forward plans for future linear or circular e + e − colliders [2][3][4][5][6][7].
A plethora of theoretical investigations have been made on non-standard Higgs bosons and their interactions within SM extensions and, specifically, on the production and the decay modes of the SM and of non-SM Higgs bosons including higher-order QCD and electroweak radiative corrections. In this paper we are concerned with the decays of heavy neutral non-SM Higgs bosons and of the 125 GeV Higgs boson -these states will be generically denoted by h in the following -into a massive quark antiquark pair (QQ) at next-to-next-to-leading (NNLO) order in the QCD coupling. Let us briefly recapitulate previous work on neutral Higgs-boson decays to quarks. The order α s QCD corrections to the decay width of h → QQ and several distributions including the full quark-mass dependence were computed both for scalar and pseudoscalar Higgs bosons by [8][9][10][11][12][13][14]. Higher order QCD corrections to the hadronic decay widths of neutral Higgs bosons, in particular the decay widths into quark antiquark pairs, were determined from the imaginary part JHEP07(2018)159 of respective two-point current correlation functions. The corrections to order α 2 s were computed in [15] for massless quarks, while quadratic and higher-order quark-mass corrections in the m Q /m h expansion were calculated both for scalar and pseudoscalar Higgs bosons by [16][17][18][19] and by [20][21][22]. The order α 2 s scalar and pseudoscalar current correlators were determined by [23] for arbitrary values of m Q /m h using Padé approximations. The hadronic decay width of a scalar Higgs boson was computed for massless quarks to order α 3 s and α 4 s in [24,25] and [26][27][28][29], respectively. Partial results for the decay of a pseudoscalar Higgs boson to three massless partons at NNLO QCD were presented in [30]. Electroweak corrections were determined in [31,32]. For Higgs-boson decays into quark antiquark pairs, differential distributions were computed by [33,34] for the decay of a scalar Higgs boson into massless b quarks at order α 2 s . We analyze in this paper the decay of a neutral Higgs boson h into a massive quark antiquark pair, h → QQX , Q = t, b , (1.1) at order α 2 s , i.e. at NNLO QCD, and to lowest order in the Yukawa and electroweak interactions. The quark-mass dependence of the matrix elements is taken into account without any approximation. 1 We consider the general case of a neutral Higgs boson with both scalar and pseudoscalar couplings to quarks. Its Yukawa interaction is where Q and h denote bare fields, y 0,Q is the bare SM Yukawa coupling, Because we consider in this work only QCD corrections, the reduced Yukawa couplings a Q , b Q do not get renormalized. Our set-up applies to the general case (1.2) where h is not a CP eigenstate. For the decay of a heavy Higgs boson into tt we analyze the cases where h is CP-even (i.e., a scalar where b t = 0) and CP-odd (a pseudoscalar where a t = 0), respectively. Because we consider unpolarized Q,Q our NNLO QCD results can be combined and applied to the case of h being a CP mixture (see below). In our analysis of h(125GeV) → bbX we assume that h is CP-even.
Our paper is organized as follows. We outline in the next section the calculation of the differential decay rate and distributions of IR-safe observables of (1.1) at order α 2 s using the

JHEP07(2018)159
antenna subtraction method [35][36][37][38]. First we use renormalized matrix elements where the QCD coupling is defined in the MS scheme and the quark masses and the Yukawa coupling are defined in the on-shell scheme. In section 3 we express the on-shell Yukawa coupling in terms of the MS Yukawa coupling and derive a generic formula for the resulting differential decay widths of (1.1) to order α 2 s . In section 4 we apply this set-up to the decays of a scalar and a pseudoscalar Higgs boson h → ttX. We calculate the respective decay widths for a sequence of Higgs-boson masses and we compare with results that are obtained by expansion in (m t /m h ) 2 to fourth order [21]. In order to demonstrate the usefulness of our differential approach we compute the top-quark energy distribution in the Higgs-boson rest frame and the tt invariant mass distribution at NNLO QCD both for a scalar and a pseudoscalar Higgs boson of mass m h = 500 GeV. In section 5 we analyze h(125) → bbX at NNLO QCD for massive b quarks, using the formulas of sections 2 and 3, where we assume the 125 GeV Higgs boson to be purely CP-even. We compute the inclusive decay width and the decay rates into two, three, and four jets using the Durham cluster algorithm. Moreover, we compute for two-jet events the distribution of the energy of the leading jet. In addition we compare with results from the literature obtained at order α 2 s for massless b quarks. We conclude in section 6.
2 The differential decay rate at NNLO QCD In this section we briefly outline the salient features of computing the decay rate of h → QQX to order α 2 s at the differential level using the antenna subtraction method. Here Q denotes a massive quark, for instance, the b or t quark. The antenna method with massive quarks in the final state that we use has been described in detail in [39] for the process e + e − → QQX. Because the strategy and formulas presented in [39] can be applied in analogous fashion to the reaction at hand we will delineate in the following only the most important features of computing the various infrared-finite contributions to the differential Higgs-boson decay rate to a pair of massive quarks at order α 2 s and refer for details to [39]. We work in QCD with n f massless quarks q and one massive quark 2 Q. All matrix elements in this and in the following sections refer to renormalized matrix elements. We define the renormalized mass of Q in the on-shell scheme and denote it by m Q , while the QCD coupling α s is defined in the MS scheme. Dimensional regularization is used to handle IR singularities that appear in intermediate steps of our calculation. As to the renormalization of the Yukawa coupling y 0,Q , we proceed as follows. First, we renormalize it in the on-shell scheme and then convert it to the MS scheme. Details will be given in section 3. We use the following notation for the differential rate of the decay (1.1): where dΓ LO is the leading-order contribution and dΓ 1 , dΓ 2 are the contributions of order α s and α 2 s . The first order correction receives contributions from the interference of the Born and one-loop amplitude of h → QQ and the squared Born amplitude of h → QQg.

JHEP07(2018)159
In any subtraction scheme for handling the IR divergences the NLO QCD correction can be written as follows: where = (4 − D)/2 is the dimensional regularization parameter and the subscripts Φ n denote n-particle phase-space integrals. The second term in the first and second square bracket of (2.2) is the unintegrated and integrated subtraction term that renders the difference, respectively the sum of the terms in the square brackets finite in D = 4 dimensions.
Here and in the following the symbol n indicates the analytic integration over the phase space of n unresolved partons in D = 4 dimensions. The NLO subtraction terms labeled with the superscript S in (2.2) were computed within the antenna framework in [40]. The second-order term dΓ 2 in the expansion (2.1) of the differential decay rate receives the following contributions: i) the double real radiation contribution dΓ RR NNLO associated with the squared Born amplitudes of h → QQgg and h → QQqq (where q denotes a massless quark), and above the 4Q threshold, the squared Born matrix element of h → QQQQ, ii) the real-virtual cross section dΓ RV NNLO associated with the matrix element of h → QQg to order α 2 s (1-loop times Born), and iii) the double virtual correction dΓ V V NNLO associated with the matrix element of h → QQ to order α 2 s (i.e., 2-loop times Born and 1-loop squared). Apart from the QQQQ contribution, which is IR finite, the terms i), ii), iii) are IR divergent. When a subtraction method is used the second order correction dΓ 2 , where the different pieces are separately finite in D = 4 dimensions, is constructed schematically as follows: 3) The symbols dΓ S NNLO and dΓ T NNLO denote the double-real subtraction terms (for h → QQqq and h → QQgg) and the real-virtual subtraction term (for h → QQg), respectively. We shortly discuss in turn the various terms in (2.3).

Double real radiation corrections
The first term on the right-hand side of (2.3) receives the following contributions: The last term on the right-hand side of this equation is the unsubtracted, IR-finite differential decay rate associated with the 4Q final state, while the first two terms denote the subtracted, IR-finite differential decay rates associated with the QQqq and QQgg final states. Within the antenna subtraction framework they are defined as follows: In the case of ij = qq a sum over all massless quark-antiquark pairs is implicit on the right-hand side of (2.5 The antenna subtraction terms required in (2.5) were determined in [41] and [42] for the QQqq and QQgg final states, respectively. They involve redefined on-shell momenta constructed from the final-state parton momenta. For the construction in the case where two massive quarks are involved, see for instance [39].
There is a subtlety associated with angular correlations in the unsubtracted squared matrix element of QQij, caused by the splitting into a pair of massless partons i, j of the virtual gluon radiated off the Q orQ. The subtraction term dΓ S,b,2,QQij NNLO contains these angular correlations, too, but the subtraction terms dΓ S,a,QQij NNLO and dΓ S,b,1,QQij NNLO do not [39]. This precludes a local cancellation of the IR singularities. A very efficient IR cancellation can be achieved by averaging out these angular correlations. This can be done as follows [43][44][45][46]. Denoting the final-state four-momenta by k 1 , k 2 , k 3 , k 4 where k 3 , k 4 are the momenta of the massless partons i, j, one evaluates the two terms in (2.5) that contain the angular correlations for each set of phase-space momenta k 1 , k 2 , k 3 , k 4 and k 1 , k 2 , k 3r , k 4r and takes the average. The four-momenta k 3r , k 4r are obtained by rotating the spatial parts of k 3 , k 4 by an angle π/2 around the collinear axis k = k 3 + k 4 . We will apply this procedure in our calculation of the inclusive cross sections h → QQX (Q = t, b) and distributions associated with the massive quark. In section 5.2 we compute, in addition, the two-jet, three-jet, and four-jet rates for the decay h → bbX. In this case this averaging procedure is not applicable because the rotated momenta k 3r , k 4r are located, in general, outside of the respective n-jet phase space. Nevertheless, straightforward subtraction as prescribed in (2.5) leads to a numerically stable IR cancellation in the calculation of the n-jet cross sections.

Real-virtual corrections
Within the antenna subtraction framework, the second term on the right-hand side of (2.3) receives the following contributions: .
(2.6) The term dΓ RV NNLO involves the interference of the tree-level and renormalized one-loop amplitude associated with h → QQg. This term contains explicit IR poles, i.e., single and doubles poles in . These poles are canceled by the explicit IR poles of the second term on the right-hand side of (2.6) that is determined by The integrals in (2.7) denote the analytic integration of the a-type subtraction terms appearing in (2.5) over the phase space of one unresolved massless parton in D = 4 dimensions. The result of the integration involves integrated tree-level antenna functions determined

JHEP07(2018)159
in [40,47]. Apart from its explicit IR poles dΓ RV NNLO is IR-singular also in the limit where the external gluon becomes soft. The term dΓ T,b,QQg NNLO is constructed such that it cancels those terms of dΓ RV NNLO which become singular in the limit k gluon → 0. This subtraction term involves massive one-loop antenna functions that were computed in [48].
In certain regions of phase space, the subtraction terms (2.7) and dΓ T,b,QQg NNLO develop IR singularities that do not coincide with the IR singularities of dΓ RV,QQg NNLO . For removing these spurious singularities one must introduce an additional subtraction term given by where the integrands are subtraction terms that appear in (2.5). The analytic integration in D = 4 dimensions over the phase space of one unresolved parton was performed in [48]. By construction the expression (2.6) is finite in D = 4 dimensions when integrated over the three-parton phase space, or arbitrary sections thereof when computing differential distributions. We recall that the terms dΓ T,a,QQg NNLO and dΓ T,c,QQg NNLO are counterbalanced by the double-real subtraction terms dΓ S,a,QQij NNLO and dΓ S,b,1,QQij NNLO (ij = gg, qq), respectively, that were introduced in (2.5). Hence, only the integrated form of dΓ T,b,QQg NNLO has to be added back to the contribution from the double virtual corrections.

Double virtual corrections
Recalling the subtraction terms introduced in (2.5) and (2.6) we see that those that are not yet counterbalanced are dΓ S,b,2,QQij NNLO (ij = qq, gg) and dΓ T,b,QQg NNLO . They have to be integrated over the unresolved two-parton, respectively unresolved one-parton phase space and then serve as counterterms for the subtraction of the IR poles of the double virtual correction. Thus the last term on the right-hand side of (2.3) is given by The term dΓ V V,QQ NNLO is the order α 2 s correction computed from the interference of the renormalized two-loop and Born matrix element and the squared one-loop matrix element of h → QQ. In the sum (2.9) all IR poles cancel in analytic fashion. One can then take the limit → 0 and perform the integration over the two-parton phase space in four dimensions.

Decay rate in terms of the MS Yukawa coupling
After having outlined our subtraction formalism we discuss a few details of our computation of the differential decay rate of (1.1) for general Yukawa couplings (1.2) and unpolarized Q,Q. For constructing the subtraction terms we need the (un)integrated antenna functions determined at NLO QCD in [40,47] and at NNLO QCD in [41,42,48] and, in addition, the squared tree level matrix elements of h → QQ and h → QQg, and the O(α s ) interference JHEP07(2018)159 matrix element (tree level and renormalized one-loop amplitude) of h → QQ. These matrix elements and the antenna subtraction terms [40,47] determine also the differential decay rate dΓ 1 at NLO. Because we sum over spins, dΓ 1 does not contain a CP-odd term proportional to a Q b Q .
As to the unsubtracted matrix elements that contribute to dΓ 2 . Computation of the squared tree-level matrix elements of h → QQgg, h → QQqq, and h → QQQQ is straightforward. Summation over the spins of the partons in the final state eliminates the CP-odd terms proportional to a Q b Q . That is, these squared matrix elements are given by the incoherent sum of the respective squared matrix elements for the decay of a CP-even and CP-odd Higgs boson h. This holds true also for the contributions to dΓ 2 associated with three-parton and two-parton final states.
The renormalized one-loop amplitude for h → QQg required in (2.6) is computed with standard one-loop perturbation theory methods. As to the O(α 2 s ) contribution to h → QQ. The renormalized (but IR-divergent) non-singlet two-loop scalar and pseudoscalar hQQ form factors (and the one-loop form factors including the terms of order ) that enter the respective matrix element were computed in [49]. The two-loop singlet scalar form factor is contained in the results of [49] while the singlet pseudoscalar form factor was determined in [50]. The singlet form factors are UV-and IR-finite. The scalar and pseudoscalar twoloop form factors were recently computed also in [51]. While the unrenormalized non-singlet form factors of [51] agree with those of [49] the renormalized ones differ by terms proportional to ζ(2) = π 2 /6. This is because the renormalization constant for the QCD coupling in the MS scheme used in [51] is not identical to the one used 3 in [49]; they differ at order 2 by a term proportional to ζ(2). This does not affect the UV renormalization of Green functions with at most a single UV pole per loop, but it affects the two-loop quark wave function and mass renormalization constants in the on-shell scheme that have double IR -poles and the renormalized non-singlet form factors. However, we emphasize that these differences are unphysical infrared artefacts. Our NNLO subtraction antenna subtraction terms [48] were computed with the same MS coupling renormalization convention as used in [49,52,53]. Therefore these IR artefacts cancel in the computation of observables. A further remark pertains to the order α 2 s flavor-singlet contributions to the scalar and pseudoscalar h → QQ and h → QQg amplitude, depicted in figure 1a and 1b, respectively. These amplitudes are ultraviolet-and infrared-finite. Let us first discuss h → tt. The case where h has a scalar coupling to the internal top-quark triangle (that is, equal masses associated with the inner and outer fermion line) has been taken into account in [49], while the amplitude associated with a pseudoscalar coupling of h to the internal top-quark loop was determined in [50]. Here we neglect the masses of the five quarks b, . . . , u in the propagators of the respective fermion triangle. Therefore, they do not contribute to the flavor-singlet scalar and pseudoscalar amplitudes h → tt and h → ttg because the respective massless quark triangle is proportional to the trace of the product of an odd number of γ matrices. The contribution to h → ttg with a top quark circulating in the triangle is computed with standard methods. In our calculation of h(125GeV) → bbX we assume h to be CP-even, take the b quark to be massive and u, d, c, s to be massless. As just mentioned, concerning the h → bb amplitude, the flavor-singlet scalar equal-mass case with an internal and external massive b quark is contained in the results given in [49]. The massive b-quark triangle contribution to h → bbg is again computed with standard methods. The scalar flavor-singlet contribution to h → bbX with a top quark circulating in the triangle -which does not decouple -and an external massive b quark, that is, the sum of the contributions shown in figure 1a and 1b, where Q = t and Q = b, was computed in [20] to leading order in the b-quark mass as an expansion in inverse powers of the top-quark mass. (The first four terms of this expansion are given in [20].) For our computation of the two-jet and three-jet rates for h → bbX in section we will need the contributions figure 1a and 1b with a top quark (and b quark) circulating in the respective quark triangle separately. We compute the contribution to the differential decay rate of h → bbg with a top-quark in the triangle, and by subtraction of this term from the result of [20] we obtain the top-quark triangle contribution to h → bb.

JHEP07(2018)159
As already mentioned above, we first compute the differential decay width in terms of the on-shell mass m Q of Q and in terms of the on-shell Yukawa coupling y Q = m Q /v. It reads schematically to order α 2 s : where the MS coupling α s (µ) is defined in QCD with n l massless and one massive quark, µ is the renormalization scale, and m h denotes the Higgs-boson mass. To lowest order in α s we havê where β Q = 1 − 4m 2 Q /m 2 h and N c denotes the number of colors.

JHEP07(2018)159
It has been known for a long time [8] that it is not a good idea to express the decay rate in terms of the on-shell Yukawa coupling when m Q /m h 1, because then the decay rate contains already at NLO a large logarithm in m Q /m h that arises from the onshell renormalization. This logarithm can be absorbed by using the MS Yukawa coupling y Q (µ) = m Q (µ)/v, where m Q (µ) is the running MS mass, and choosing µ = m h . In the following we convert in (3.1) the on-shell Yukawa coupling y Q to y Q (µ) using the relation between these couplings to order α 2 s . One may express the coefficients dΓ QQ i defined in (3.2) also in terms of the MS mass m Q (µ). However, we keep their dependence on the on-shell mass m Q because one may argue that this yields a somewhat more realistic description of the decay kinematics when computing distributions. We notice that this hybrid renormalization prescription also avoids large logarithms in the case m Q /m h 1. In our analysis of h → ttX in section 4 we consider a range of Higgs-boson masses for which the ratio m t /m h is not very small. Nevertheless it turns out that using the MS Yukawa coupling y t is appropriate also in this case, see section 4.
For expressing y Q in terms of y Q we need the relation between the on-shell mass m Q and the MS mass m Q (µ) to order α 2 s . It reads in QCD with n l massless and one massive quark [54,55]: where and Here Inserting (3.4) into the squared Yukawa coupling in the on-shell scheme, y 2 Q = m 2 Q /v 2 , and expanding to order α 2 s we get

JHEP07(2018)159
where y 2 Q (µ) is the squared Yukawa coupling in the MS scheme, Let us now apply (3.9) to the differential decay width of h → QQX. Inserting (3.9) into (3.1) and expanding to order α 2 s we obtain schematically: (3.11) Because of the truncation of the perturbation series, equation (3.11) is not identical to (3.1); therefore we use the bar notation. The decay width can be cast into the form The MS Yukawa coupling makes the LO decay width y 2 Q (µ)Γ QQ 0 dependent on µ. In our application of (3.11) and (3.12) to h → ttX we work in 6-flavor QCD with α s = α (6) s and n l = 5 while in the application to h → bbX we evaluate these formulas with α s = α (3.14) i.e., the logarithm L Q in u 1 , u 2 is replaced bȳ .

(3.16)
For computing the MS top and bottom Yukawa couplings at arbitrary scales µ we use the solution of the renormalization group equation for m Q (µ) at two-loops. It reads 17) The coefficients are, putting N c = 3 and n f = n l + 1:

JHEP07(2018)159
and α s in (3.17) is the coupling in n f -flavor QCD. Finally, a remark concerning the relation between the on-shell and MS Yukawa coupling of Q. When the iterated relation (3.14) instead of (3.4) is used to express y 2 Q in terms of y 2 Q to order α 2 s , the coefficients r 1 and r 2 in (3.9) change; in particular, they depend then on the MS mass of Q. In cases where the ratio m Q /m h is not very small, one may consider it to be a matter of convention which one of these two variants one uses. In cases where m Q /m h tends to zero the formulas (3.11) and (3.12) are to be used if the matrix elements are computed in terms of the on-shell mass m Q , as we do. This is because for small m Q the dΓ QQ i (i = 1, 2) develop large logarithms in the on-shell mass which are compensated by the large logarithms in m Q that appear in r 1 , r 2 . We will show in section 5.1 that by choosing a small b-quark mass, eq. (3.12) allows to recover the known h → bbX decay width at order α 2 s for massless b quarks.

Decays of non-standard Higgs bosons to ttX
In this section we analyze the decay of a non-standard heavy neutral Higgs boson of arbitrary CP nature into tt for a set of Higgs-boson masses. As emphasized before we consider unpolarized final states. Then, as mentioned in the previous section, the differential decay rate to order α 2 s and to lowest order in the electroweak couplings decomposes into the incoherent sum where dΓ tt S (dΓ tt P ) denotes the differential decay rate of a CP-even (CP-odd) Higgs boson with reduced Yukawa couplings a t = 1, b t = 0 (a t = 0, b t = 1). In the following we compute Γ tt S , Γ tt P and the respective top-quark energy distribution in the Higgs-boson rest frame and the tt invariant mass distribution using the NNLO antenna subtraction framework outlined in section 2 and the formulas (3.11), (3.12).
We use m t = 173.34 GeV for the on-shell mass of the top quark. With the relation (3.4) we get m t (µ = m t ) = 163.46 GeV for the MS mass at µ = m t . This value is used as an input for computing with (3.17) the running top-quark mass m t (µ) with which the MS Yukawa coupling y t (µ) is determined. The renormalization scale is chosen to be µ = m h . The effect of scale variations is assessed by varying µ between µ = m h /2 and µ = 2m h . For the QCD coupling we use α   Higgs boson with mass m h = 500 GeV and µ = m h the α 2 s correction increases the NLO width by 9.2% (5.8%) while for m h = 680 GeV the increase amounts to 5.3% (3.4%).
If the ratio m t /m h is not significantly smaller than one, one may argue that the decay widths can also be computed using the on-shell top Yukawa coupling instead of the MS Yukawa coupling. In table 1 we list, for two Higgs-boson masses, our results for the tt decay widths at NNLO QCD obtained with (3.1) and (3.12). The decay widths Γ tt X (X = S, P ) are slightly smaller than the corresponding widths Γ tt X (X = S, P ). Using the on-shell top Yukawa coupling we note that for m h 600 GeV the order α 2 s contribution to the respective decay width becomes of the same order of magnitude than the order α s contribution, while figure 2 shows that the perturbation expansion in α s works well when the MS top Yukawa coupling is used.
Moreover, we can compare our NNLO QCD results for the tt decay widths, which are exact in m t , with those of 4 [21] that were obtained by expanding the scalar and pseudoscalar current correlation functions to fourth order in (m t /m h ) 2 . Let us denote the order α 2 2,S /Γ tt 2,S = 1.007 and Γ tt,HS 2,P /Γ tt 2,P = 1.026. As expected, for smaller ratios m t /m h the agreement becomes even better. For m h = 680 GeV and µ = m h we find Γ tt,HS 2,S /Γ tt 2,S = 0.9997 and Γ tt,HS 2,P /Γ tt 2,P = 1.0003. Next we turn to differential distributions. We consider the distribution of the normalized top-quark energy x t = 2E t /m h in the rest frame of the Higgs boson and the tt invariant mass distribution. We recall that the distributions expressed in terms of the MS Yukawa coupling are computed with (3.11). For definiteness we choose m h = 500 GeV.
The left plot of figure 3 shows dΓ  the SM Higgs boson, that is, we assume it to be CP-even with a b = 1. Moreover, we use the Higgs-boson mass m h = 125.09 GeV [56].
For the b-quark mass we use m b (µ = m b ) = 4.18 GeV [56] in n f = 5 flavor QCD as input in (3.17) for computing the running MS b-quark mass and the Yukawa coupling y b (µ) at arbitrary scales µ. As in section 4 we use α (5) s (m Z ) = 0.118. From the two-loop relation (3.14) we obtain m b = 4.78 GeV. This value of the on-shell b-quark mass is used in computing the coefficients r 1 , r 2 , the dΓ bb i , and the γ bb i defined in eq. (3.10), (3.11), and (3.13), respectively. (Higher order corrections to (3.14) shift this value of m b to a somewhat higher value [57], but for the sake of consistency of the order of the perturbative expansion, we stick to m b = 4.78 GeV.) With the above value of m b (µ = m b ) we get, using (3.17), the values m b (µ) = 2.97 GeV, 2.80 GeV, and 2.65 GeV at µ = m h /2, m h , and 2m h , respectively. With these values we compute the MS Yukawa couplings y b (µ). As in section 4 we use α (5) s (m Z ) = 0.118. The displayed digits of our numerical results given below are not affected by our numerical integration errors.

Inclusive decay width
First, we determine the inclusive decay width of h(125) → bbX at NNLO QCD using the antenna subtraction framework of section 2. We have to take into account also h → bbbb whose contribution is IR finite. We represent our result for the inclusive decay width at NNLO QCD in the form (3.12): where Γ bb LO = y 2 b (µ)Γ bb 0 , g 1 = γ bb 1 + r 1 , g 2 = γ bb 2 + r 1 γ bb  figure 1 constitute an important part of the order α 2 s corrections. This contribution, which is µ-independent, to g 2 listed in table 2 is g 2 (t) = 6.898.
As a check of our computational set-up in the case of h → bbX we determine the QCD correction coefficients g 1 and g 2 in the limit of small b-quark mass and compare with the results for massless b quarks, g 1 (m b = 0) = 5.6666 and g 2 (m b = 0) = 29.1467 [15,24]. We cannot perform the limit m b → 0 analytically. Therefore we choose a value for m b that is very small compared to m h , to wit, we choose m b = 0.5 GeV and compute the resulting g 1 and g 2 for µ = m h . For the order α s QCD coefficient we obtain g 1 (m b = 0.5GeV) = 5.6685. The α 2 s QCD coefficient given in [15,24] does not include the topquark triangle contribution to h → bb and h → bbg, respectively. In fact these contributions vanish for m b = 0 because of helicity mismatch. Omitting these contributions we obtain g 2 (m b = 0.5GeV) = 29.187 which is very close to the number for m b = 0. Using the same Yukawa coupling y b (m h ) as in the massive case, the resulting NNLO decay width for

Two-jet, three-jet, and four-jet rates
Concerning the decay of h(125) to b quarks, other interesting observables for future experimental analyses include the decay rates into two, three, and four jets. For definiteness we use here the Durham jet algorithm [58] (and the recombination scheme k (ij) = k i + k j ) with jet resolution parameters y cut = 0.01 and y cut = 0.05. The n-jet rates can be represented, JHEP07(2018)159 y cut = 0.01 y cut = 0.05  Table 3. The coefficients g i (n jet) defined in eqs. (5.6)-(5.8) and computed with the Durham algorithm using y cut = 0.01 and y cut = 0.05 for three renormalization scales µ. They determine the n-jet rates (5.3)-(5.5).
in analogy to (3.11) and (3.12), to order α 2 s as follows: where g 1 (2 jet) = γ bb 1 (2 jet) + r 1 , g 2 (2 jet) = γ bb 2 (2 jet) + r 1 γ bb 1 (2 jet) + r 2 , (5.6) g 1 (3 jet) = γ bb 1 (3 jet) , g 2 (3 jet) = γ bb 2 (3 jet) + r 1 γ bb 1 (3 jet) , (5.7) g 2 (4 jet) = γ bb 2 (4 jet) . (5.8) Our results for the coefficients g i (n jet) are given in table 3 for µ = m h /2, m h , and 2m h . With decreasing resolution parameter y cut the two-jet rate decreases while the three-and four-jet rates increase. This is because a smaller y cut resolves more and more gluon radiation respectively massless quark radiation. The three-jet rate has a maximum below y cut = 0.01 and then decreases for smaller y cut , while the four-jet rate still increases for these resolution parameters till its maximum is reached. Notice that for each value of y cut and µ the sum of the g i (n jet) listed in table 3 yields the inclusive coefficient g i (i = 1, 2) given in table 2. The coefficients g i (n jet) were computed in [33] for µ = m h and massless b quarks using the JADE algorithm and y cut = 0.01.
Finally we analyze, for two-jet events, the distribution of the energy of the leading jet, i.e., the jet with the largest energy in the Higgs-boson rest frame. We use the dimensionless variable x max = max(E j 1 /m h , E j 2 /m h ) . Obviously, x max ≥ 0.5. We choose the Durham algorithm with y cut = 0.05 and the parameter set-up as specified at the beginning of this section. The distribution of this observable is shown in figure 5 at LO, NLO, and NNLO QCD. The figure shows that the perturbation series is well-behaved for this observable. This distribution was presented before for massless b quarks in [33] (where the JADE algorithm with y cut = 0.1 was used) and in [34] (where the Durham algorithm with y cut = 0.05 was used). Comparing figure 5 with figure 2 of [34] shows that taking the non-zero b-quark mass into account has only a minor impact on the shape of this distribution.

Conclusions
We have presented, within the antenna subtraction framework, the set-up for calculating the fully differential decay rate of a scalar and pseudoscalar Higgs boson to a massive quark antiquark pair at NNLO in the perturbation series in α s . Our approach is fully differential in the phase-space variables and can be used to compute any infrared-safe observable for these decays. So far, methods for calculating Higgs-boson decays at order α 2 s in QCD at the differential level were available only for massless quarks. We have applied our set-up to the decay of a scalar and pseudoscalar Higgs boson to t,t and to the decay of the h(125GeV) Higgs boson to massive b,b quarks. Apart from computing the respective decay widths at NNLO accuracy we have determined also several distributions, in order to show the use of our approach. We have presented our results in terms of the MS Yukawa coupling y Q , while the on-shell mass m Q is used in the matrix elements. For the NNLO decay width of h(125GeV) → bbX we have shown numerically that this formulation recovers, for m b /m h → 0, the known result for massless b quarks.
Our results can be used as a building block in a theoretical description of Higgs-boson production and decay to massive quarks at order α 2 s at the differential level in hadron or electron-positron collisions.