The parton-level structure of Higgs decays to hadrons at N$^3$LO

We present the quantum chromodynamics (QCD) corrections for Higgs boson decays to hadronic final states at next-to-next-to-next-to-leading order (N$^3$LO) in the strong coupling constant $\alpha_s$. In particular, we consider the Higgs boson decay to massless bottom quarks and the Higgs boson decay to a pair of gluons in the limit of a heavy top quark. The tree-level five-parton, the one-loop four-parton, the two-loop three-parton, and the three-loop two-parton matrix elements are integrated separately over the inclusive phase space and classified by partons appearing in the final state and by colour structure. As a check, we reproduce known results for the hadronic $R$-ratios at N$^3$LO. We study patterns of infrared singularity cancellation within the colour layers of the integrated expressions and observe an agreement in the highest trascendental weight terms in the decay of different colour singlets to quarks. We anticipate that our result will be an essential ingredient for the formulation of N$^3$LO subtraction schemes.


Introduction
The Higgs boson plays a special role in the Standard Model. Its discovery in 2012 [1,2] through the clean decay modes H → γγ, H → 4l and H → 2l2ν marked a major milestone in particle physics. Amongst the remaining decay modes, the decay to bottom quarks has the largest branching ratio and allows probing the Higgs boson coupling to third generation fermions. On the other hand, the loop-induced decays of a Higgs boson to gauge bosons (gg, γγ or γZ) are key to an indirect determination of the couplings to Z, W and t. A summary of the current experimental evidence for the Higgs production and decay modes can be found in reviews by ATLAS [3] and CMS [4]. The importance of Higgs decays to quarks and gluons extends beyond phenomenology. Thanks to the simplicity of the 1 → 2 Born kinematics, the computation of radiative QCD corrections to these processes has been achieved to very high orders in perturbation theory. The H → gg decay rate in an effective theory with the top quark integrated out [5][6][7] is known up to N 4 LO [8][9][10][11][12][13]. The H → bb decay rate in a massless approximation where only the bottom Yukawa coupling is kept different from zero is also known up to N 4 LO [8,[14][15][16][17]. The correction to the total decay rate of the Higgs boson to hadrons, including also contributions induced by the effective Higgs-bottom quark interaction where only the dominant m 2 b mass terms are retained, has been calculated at order N 3 LO in [18] and at order N 4 LO in [19]. We refer the reader to a review [20] for further discussion on terms suppressed by the top-quark and bottom-quark mass and electroweak effects.
All the above results follow a common approach based on the optical theorem, which amounts to the computation of the imaginary part of the Higgs self-energy. This technique is completely agnostic to the number and species of particles in the final state. Hence it does not reveal the infrared structure of the result and obscures the interplay between real radiation and virtual corrections in the different final-state partonic channels.
In this paper, we extend this analysis by separately integrating all the possible physical cuts of the four-loop QCD correction to the Higgs two-point function. In other words, we analytically integrate the tree-level five-parton, the one-loop four-parton, the two-loop three-parton, and the three-loop two-parton matrix elements over the respective phase space: where M l n denotes the l-loop matrix element for the decay of the Higgs into n final state QCD particles. We refer to the four terms in (1.1) as the triple-real (RRR), double-realvirtual (VRR), double-virtual-real (VVR), and triple-virtual (VVV) layers of the calculation. For completeness we recompute also the next-to-leading order (NLO) and next-tonext-to-leading order (NNLO) corrections. Our method was first described in [21] in the context of the decay of a virtual photon into hadrons and leveraged the reverse unitarity relation [22][23][24][25][26] to gain access to modern multi-loop techniques.
The paper is organised as follows. In Section 2, we introduce the notation and we briefly describe our method. In Section 3, we analyse our results which are explicitly reported in Appendices A and B. We conclude in Section 4 with an outlook of future applications.

Method
We consider the H → gg and H → bb decays. In the first case, we work in the heavy-top effective theory, with the QCD Lagrangian supplemented with an effective Lagrangian given by with G µν a the renormalised gluon field-strength, H the Higgs field and λ 0 the bare effective coupling, obtained by matching the effective theory to the full Standard Model [10,64,65]. In the second case, we implement the standard vertex between the Higgs field and a fermion line with ψ the bottom quark field and y b 0 the bare Yukawa coupling. In our calculation, the bottom quark is treated as massless but the Yukawa coupling is non-vanishing.
We present results for the integration of renormalised squared amplitudes in the MS scheme. We replace the bare coupling α 0 with the renormalised coupling α s according to [66] where α 0 is the bare coupling, µ 2 0 is the mass parameter introduced in dimensional regularisation to maintain a dimensionless coupling in the bare QCD Lagrangian density. We fix the renormalisation scale µ 2 to be the invariant mass of the decaying particle q 2 .
In the calculation of H → bb, the renormalisation of the bare Yukawa coupling y b 0 is needed, which is done through the replacement y b 0 = Z y y b , with Z y as in [67] In the calculation of H → gg, we additionally need to renormalise the effective coupling λ 0 = Z λ λ [68] according to [66] We follow the strategy outlined in [21], where the relevant decay diagrams are generated with QGRAF [69] as self-energies of the Higgs boson with cut internal propagators. They are matched onto the integral families reported in [21] using Reduze2 [70] and the Feynman rules are inserted and evaluated in FORM [71]. The integrals appearing in the matrix elements have up to eleven propagators in the denominator and a maximum of five scalar products in the numerator, compared to four scalar products in the photon decay. The integrals are reduced with the help of Reduze2 [70] to a set of 22, 27, 35 and 31 master integrals for the four terms of (1.1), respectively. The master integrals required for the NNLO calculation can be found in [72] and have been extended up to weight 6 in [21,66,73]. The integrals required for the N 3 LO calculation were computed in [73,74].

Results
We illustrate the general structure of the different partonic contributions by adopting a notation similar to [21] for ease of reference. Given a set of n final-state particles denoted by I, we can generically write the associated amplitude for the H → ij process as We denote the integration over the respective phase space of the matrix element ⟨M ij |M ij ⟩ I summed over spins, colours and quark flavours as and for where ij = gg, qq. The label k denotes the perturbative order: contributions with the same k sum to the N k LO result for the total cross section. The long explicit expressions for T ij,(3,[ℓ 1 ×ℓ 2 ]) I are provided in Appendix A, while in Appendix B we report the lower-order results expanded up to transcendental weight six. We denote the coefficient of each colour factor C as T ij,(k,[ℓ 1 ×ℓ 2 ]) I C , and we omit the superscript [ℓ 1 × ℓ 2 ] in case of no ambiguity. All results are also provided in the ancillary files in computer-readable format, with the notation T All expressions are renormalised and in time-like kinematics. Higher-order results are normalised to T gg,(0) for the decay to gluons and for the decay to bottom quarks, with N the number of colours and C F = (N 2 − 1)/(2N ). T gg,(0) gg and T qq,(0) qq are the respective Born-level cross sections. Note that the factor 2C F in (3.6) is included in the normalisation of H → bb as it appears in all colour layers starting from NLO. Finally, P 2 is the volume of the two-particle phase space, The structure or (non-)appearance of certain colour factors is key in understanding the cancellation patterns between real and virtual corrections. We therefore summarize the colour factors appearing in the various final-state configurations in Appendix C. In the final state with two quark lines, we explicitly separate the configurations with same-or different-flavour quark pairs. Note that for all colour factors C in T ij,(k) qqq ′q′ (g) , we have that .

(3.8)
For H → gg, in the two-particle final state, the terms proportional to N k F for k = 1, 2, 3 in T gg,(k) gg are introduced as part of renormalisation and as such are not present in the finite part or beyond. Note that the only allowed two-particle final state in the H → gg corrections is a pair of gluons, even at higher loops. This is explained by the fact that for a scalar particle to decay into a quark-antiquark pair, a chirality flip along the fermionic line is needed due to spin conservation. Therefore, all the diagrams contributing to the Higgs decay to a massless quark-antiquark pair via QCD interactions vanish.
We perform several checks on our results. First, we can directly compare our expressions for the two-particle final states to the calculations of the quark and gluon form factors up to three loops [66,67]. Second, we observe the complete cancellation of infrared (IR) poles in the sum of all the partonic final states at each perturbative order, both in H → gg and H → bb. Third, the finite parts agree with the known results for the total cross sections e.g. in [8]. Indeed, the sum of (A.1)-(A. 46 (3.10)

Infrared structure of the deepest poles
It is interesting to compare the colour and singularity structure of H → bb to that of γ * → qq derived in our previous work [21], both representing decays of a colour-singlet state into a quark-antiquark pair. Once we normalise with respect to the respective Bornlevel cross sections, we observe the presence of identical colour factors at all perturbative orders, except for the singlet contributions which are only present in γ * → qq. For the Higgs case, the singlet contribution vanishes because the Yukawa interaction in (2.2) introduces a chirality flip which would break chirality conservation along a closed massless fermionic loop. We notice that for all layers and orders, the two deepest poles appearing in any given colour factor of the integrated renormalised matrix elements coincide between the H → bb and γ * → qq processes. To clarify, we refer to the two deepest non-vanishing poles in each partonic configuration and colour layer, not simply to the ϵ −2k and ϵ −2k+1 poles at order k. For the two-particle final states, the infrared singularity structure is predicted through universal IR factorisation formulae [75,76]. In particular, the deepest poles can be interpreted in a completely process-independent way in terms of the I (1) qq insertion operator [77]. For purely real corrections up to NNLO, we can build on a recent algorithmic construction of idealised NNLO antenna functions [78]. In this work, the deepest poles of M 0 4 in the γ * → qq decay are identified with structures coming from the integrations of double-unresolved limits. Defining the operator P n [·] which extracts from an expression the n deepest non-vanishing poles in ϵ, the relations read Dsof t γγ + 2T col q h γγ , where Ssof t, Scol, Dsof t and T col are given in Appendix A and B of [78], with the superscript h in Scol and T col indicating the hard radiator. The quantities in the last line of (3.11) only exhibit a single pole and C 0 4 , which encapsulates the triple collinear limit of a q ∥q ∥ q configuration, is given in (5.57) of [78]. Our results confirm that (3.11) holds for the H → bb process as well. Beyond the two deepest poles, the finite parts of lower-order integrated matrix elements contribute to the singularities, resulting in differences between H → bb and γ * → qq. An analogous study for single-real emission at one loop will allow for a physically motivated description of the singularities of the real-virtual contribution [79]. The fact that the exact correspondence of the deepest non-vanishing poles (and not just the ϵ −2k and ϵ −2k+1 poles at order k) between the two processes extends also to N 3 LO in all the layers, colour factors and partonic final states is notable. Therefore, we expect that a universal and process-independent description exists for unresolved radiation between a hard quark-antiquark pair also at this perturbative order.
We observe that for every power of ϵ, the terms with the highest trascendental weight always coincide in the γ * → qq and the H → bb process. This is reflected in the N 3 LO decay rates: the weight 4 contribution in the N 2 F colour layer is absent in both processes; the weight 5 contribution in the N F /N colour layer is identically 5 ζ 5 ; the weight 5 contribution in the N F N colour layer is identically −10/3 ζ 5 ; the weight 6 contributions are absent in both processes. In fact, the correspondence between the highest weight contribution in every colour factor extends to the N 4 LO decay rates of γ * → qq and H → bb [8,80]. Moreover, we observe that the weight 2k contribution vanishes in the decay rate of γ * → qq, H → bb and H → gg all the way up to N 4 LO. Since the decay rate at N k LO is given by the absorptive part of a (k + 1)-loop propagator, this is directly related to the maximal trascendental weight appering in the two-point function. It would be interesting to analyse our observations in the framework of [81].
For H → gg, no analogous processes have been computed, so it is not possible to perform a comparison as done for H → bb. From [78], it is possible to argue that for the real radiation corrections up to NNLO in the purely gluonic final state, the following relations hold: On the other hand, the dependence on lower orders is manifest already at the seconddeepest poles in the qqgg and qqq ′q′ partonic channels [78]. Nonetheless, one can expect a weaker statement to hold also for H → gg, namely the universality of the ϵ −2k and ϵ −2k+1 poles at order k across all the partonic final states and colour factors.

Cancellation of infrared singularities
The pattern of cancellation of infrared singularities among the four layers contributing at N 3 LO is highly non-trivial due to the interplay of different soft currents and splitting functions in the real-emission contributions. Despite recently becoming available in their unintegrated form [49][50][51][52][53][54][55][56][57][58][59][60][61][62][63], such structures are yet to be fully understood at the integrated level. One can expect the cancellation pattern of the deepest poles to be easier to explain in terms of universal factorization properties of QCD. Moreover, the deeper the pole, the fewer infrared-divergent partonic configurations contribute to it, simplifying the analysis.
To facilitate the inspection of the deepest singularities, in Tables 1 and 2 we display the coefficients of the ϵ −6 and ϵ −5 poles for different colour layers and partonic final states. Trivially, the sum of the coefficients in each column vanishes. In the following, we make some basic observations, postponing a thorough analysis, which would require detailed knowledge of the integrated structures contributing at each layer.
For H → bb in Table 1, the easiest colour factor to inspect is N −2 , because it only receives contributions from the emission of abelian gluons from the hard pair of quarks.
For both the virtual and real corrections, the coefficients of the ϵ −6 and ϵ −5 poles are proportional to (I (1) qq ) 3 . Moreover, we notice a pair-wise cancellation between the poles of the triple-virtual and triple-real, and the double-real-virtual and double-virtual-real contributions, which can be explained by the following argument. Due to a complete factorization in the deepest poles, each virtual (or real) abelian emission from a hard quark-antiquark pair contributes to the poles independently and with a factor I (1) qq ). Hence at N k LO, the structure in every layer is proportional to (I (1) qq ) k and the relative value of the coefficients between the layers is entirely of combinatorial origin. Namely if the deepest poles of the virtual layer are captured by N (I (1) qq ) k , then the poles of the r-fold real, (k − r)-fold virtual correction are given by since there are k r ways to cut open r of the k loops and substitute a real emission for a virtual correction. In (3.13), N is an overall normalisation factor common to all layers at order k.
At NLO, (3.13) Table 1, justifying the pair-wise cancellation. In general, the cancellation of the deepest abelian poles at any order k is guaranteed by 14) The second term in the previous equation reflects the exponentiation of multiple abelian emissions. In other words, the cancellation of infrared singularities proceeds independently for each emission. With a similar combinatorial logic, one can predict the highest poles of the individual ℓ 1 × ℓ 2 terms within each layer with ℓ = ℓ 1 + ℓ 2 loops. One can show that in the purely abelian case the coefficient is split among the loop configurations according to The coefficients of the remaining colour factors for purely gluonic emissions, N 2 and N 0 , are more complicated due to the appearance of non-abelian effects. For the triplevirtual contribution, the only colour structure present in the ϵ −6 coefficients is C 2 F , as  Blank cells indicate vanishing coefficients for these poles. Note that an overall factor of 2C F is factored out as indicated in (3.6).
indicated in the first line of Table 1, in accordance with [67] once we account for the extra colour factor of 2C F due to our normalisation in (3.6). For the other layers, the colour factors C F C A and C 2 A also contribute to the deepest singularities, indicating that the description of the poles associated with real emission contains ingredients different from I (1) qq already at ϵ −6 . An equivalent observation at NNLO was the basis of the analysis of the infrared structure performed in [82]. Additionally, the values of the coefficients in the N 2 contributions to the three-and four-particle final state suggest the presence of an integrated structure cancelling between these two layers. For the N F N −1 colour factor, we see a complete cancellation between the double-real-virtual and triple-real layers, explained by a single soft gluon emission on top of a four-quark final state.
For H → gg in Table 2, the rows representing final states containing four quarks are left blank since they only exhibit poles starting from ϵ −4 . This happens because the most infrared-divergent behaviour at the triple-real level with four quarks in the final state is given by the emission of two collinear quark-antiquark pairs with one of the pairs also soft, or alternatively by the emission of two collinear quark-antiquark pairs in association with a single soft gluon emission. Both of these configurations yield at most ϵ −4 poles.  the four-and five-particle final states. On the other hand, the N F N 2 coefficients exhibit a highly non-trivial interplay between real quark-antiquark pair emissions and fermionic loops.

Conclusions and outlook
In this paper, we carried out the analytical integration over the inclusive phase space for all the possible partonic channels for the decay of a Higgs boson into a pair of gluons and a bottom quark-antiquark pair up to N 3 LO. This work is a natural continuation of [21], where an analogous calculation was performed for the decay of a virtual photon into hadrons. Remarkably, the deepest IR poles of photon and Higgs decay to quarks coincide completely for any partonic final state. We observe a similar agreement also in the highest transcendental weight numbers in the coefficients of all powers of ϵ.
These results constitute an important contribution to the study of unresolved QCD radiation at the integrated level. In the context of the antenna subtraction method, they are necessary for the extraction of N 3 LO gluon-gluon antenna functions in final-final kinematics.
The set of integrated antenna functions at N 3 LO ought to be completed with analogous results for a quark-gluon pair of hard radiators. At NNLO, the expressions were extracted from QCD corrections to a neutralino decay [83]. We foresee extending this analysis to N 3 LO in a forthcoming publication.   Table 3: Colour factors appearing in the Higgs decay to bottom quarks, organised by finalstate particles I, perturbative order k and loop configuration ℓ 1 × ℓ 2 , in case of ambiguity. Note that an overall factor of 2C F is factored out as indicated in (3.6).