Analytic results for decays of color singlets to gg and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q \bar{q}$$\end{document}qq¯ final states at NNLO QCD with the nested soft-collinear subtraction scheme

We present compact analytic formulas that describe the decay of colorless particles to both \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q \bar{q}$$\end{document}qq¯ and gg final states through next-to-next-to-leading order in perturbative QCD in the context of the nested soft-collinear subtraction scheme. In addition to their relevance for the description of decays like \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V \rightarrow q \bar{q}'$$\end{document}V→qq¯′, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V= Z,W$$\end{document}V=Z,W, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H \rightarrow b \bar{b}$$\end{document}H→bb¯ and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H \rightarrow gg$$\end{document}H→gg, these results provide an important building block for calculating NNLO QCD corrections to arbitrary processes at colliders within the nested soft-collinear subtraction scheme.


Introduction
The development of an efficient and physically transparent subtraction scheme for next-to-next-to-leading order (NNLO) computations in QCD is an important problem in theoretical particle physics that attracted a lot of attention recently [26][27][28][29][30]. However, among the many subtraction schemes that have been proposed, there is not a single one that is generic, fully local and fully analytic (in a sense that all the integrated subtraction terms are available in an analytic form). Given the impressive practical successes of many subtraction schemes in describing physical processes, it is unclear whether or not locality and analyticity are truly essential. However, we believe that it is useful to develop a scheme that is general, physically transparent and efficient, especially in view of the need to extend the functionality of existing subtraction schemes beyond 2 → 2 processes for forthcoming LHC applications.
In Ref. [30], we introduced the nested soft-collinear subtraction scheme. It is based on the idea of sector decomposition [10,11] but it relies heavily on the phenomenon of color coherence in constructing soft and collinear approximations to matrix elements. This subtraction scheme is local by construction; however, initially, some subtraction terms were not known analytically. Recently, this problem was solved for both the double-soft [31] and triple-collinear [32] subtraction terms so that analytic results for all double-unresolved subtraction terms are now available. Building on that, in Ref. [33] we presented analytic results for the production of a color-singlet final state in hadron collisions obtained within the nested soft-collinear subtraction scheme. In addition to their phenomenological relevance, we view these results as building blocks that should, eventually, allow us to describe arbitrary hard processes at hadron colliders through NNLO QCD. Typically, these building blocks are obtained by partitioning the phase space for a particular process in such a way that only emissions off two hard particles at a time lead to infra-red and collinear singularities when integration over the phase space is attempted. These hard emittors can be both in the initial or in the final state or one of them can be in the initial and the other one in the final state. When looking at the problem of constructing a subtraction scheme from this perspective, the results presented in Ref. [33] should facilitate the description of the two initial-state emittors.
The goal of this paper is to take one further step towards the application of the nested soft-collinear subtraction scheme to the description of generic LHC processes by considering a situation when the hard emittors are in the final state. An important physical example of this situation is decays of colorless particles into a qq or gg final state. The NNLO QCD results for the qq final state have already been used by us in Ref. [34] to describe the decay of the Higgs boson into a massless bb pair; however, we did not provide analytic formulas for this final state in that reference. The goal of this paper is to provide such formulas and to supplement them with the analytic results for decays of a color singlet into a gg final state.
Although, conceptually, the computation of NNLO QCD corrections to the production and decay of a color singlet are very similar, there are a few differences between the two that are worth pointing out.
• In the case of the double-real corrections to the gg final state we need to carefully separate unresolved gluons from the resolved ones. This issue does not appear in case of production where incoming particles are always the hardest ones and their momenta are fixed. • The computation of the integrated collinear counterterms requires modifications since, in the initial-state case, the integrated collinear subtraction terms are functions of fractions of the initial energy that a hard parton carries into the hard process, while in case of the finalstate emissions one has to integrate over fractions of energies that are shared by partons in the collinear splitting. • Construction of the double-collinear phase space, i.e. the phase space appropriate for the description of a kinematic situation where singularities occur when each unresolved parton is emitted by a different emittor, is straightforward in the production and non-trivial in the decay cases. • Obviously, no renormalization of parton distribution functions is needed to describe decay processes; for this reason, cancellation of infra-red and collinear singularities works differently in the production and decay cases.
The rest of the paper is organized as follows. In Sect. 2 we set the stage for the calculations described in the following sections and introduce our notation. We then discuss in detail the calculation of QCD corrections to H → gg decay to explain our approach. In particular, in Sect. 3, we present the computation of the NLO QCD corrections to the decay rate H → gg. In Sect. 4 we discuss how to set up the calculation of NNLO QCD corrections to H → gg decay and then consider the H → 4g channel in detail. We present our final results for the NNLO QCD corrections to the decay of a color singlet to two gluons in Sect. 4 and to a qq final state in Sect. 5. We discuss the validation of our results in Sect. 6, and conclude in Sect. 7. Many useful formulas and intermediate results are collected in several appendices.

General considerations
We begin by describing common features of QCD corrections to color singlet decays and by introducing notations that we will use throughout the paper. We consider decays of a color-singlet particle Q to quarks and gluons. Our goal is to provide formulas that describe NNLO QCD corrections to these decays at a fully-differential level. Specifically, we study the decay process Q → f i f j + X , where { f i , f j } can be either {g, g} or {q,q}. We first discuss the decays into the gg final state since, compared to Q → qq, the singularity structure of the decay Q → gg is more complex. Therefore, once the calculation of the NNLO QCD corrections to Q → gg is understood, NNLO QCD corrections to Q → qq are easily established.
We write the perturbative expansion of the differential decay rate as The different contributions in Eq. (2.1) are obtained by integrating various matrix elements squared over the phase space of final state particles. To describe this integration in a compact way, we introduce the notation analogous to our earlier papers [30,33] and define where N is a symmetry factor for identical final-state particles, d = 4 − 2 is the space-time dimensionality, is the phase-space element for a parton f i , M tree (1 f 1 , . . . , n f n ) is the matrix element for the process (2.4) and O is a function that depends on partons' energies and angles. Furthermore, E max is an auxiliary parameter with the dimension of energy that should be large enough to accommodate all events that are allowed by the energy-momentum conservation constraints. Its relevance will become clear in what follows. In the rest of this paper, we will use E max = m H /2. We note that the explicit constraint on the energy in Eq. (2.3) breaks Lorentz invariance at intermediate stages of the calculation; for this reason all energies in this paper are defined in the rest frame of the decaying particle Q.
To describe contributions of loop-corrected processes, we introduce similar quantities 1 and Finally, we define where X = LM, LV, LVV and the sum runs over all allowed final states. Using these notations, the three contributions to the differential width Eq. (2.1) are written as The symbol · · · δ indicates that the integration over the momenta of partons that are explicitly shown as arguments of a function F X is not performed, so that the right hand side 1 We note that in this paper we always work with UV-renormalized amplitudes.
of Eq. (2.8) provides a fully-differential description of the decay rate. Starting from next-to-leading order, the individual terms appearing on the right hand sides of Eq. (2.8) are infra-red divergent and cannot be integrated in four dimensions when taken separately. The goal of a subtraction scheme is to rearrange them in the following way are finite in four dimensions and contain contributions from final states with at most i partons. In Refs. [30,33] we explained how this can be done for hadroproduction of color-singlet states. We now use a very similar procedure to discuss color singlet decays.
Since the required computations are often quite similar, we do not describe the calculational details if the results for the decay follow easily from the ones for the production. To this end, we note that a detailed introduction to our subtraction scheme can be found in Refs. [30,33] and we extensively refer to these papers in what follows. In this paper, we highlight differences between the computations required for the production and decay cases and present formulas for color singlet decay to either gg or qq final states. We begin with the discussion of the NLO QCD corrections to H → gg.

Higgs decay to gluons: a NLO computation
We consider the NLO QCD contribution to the differential decay rate of the Higgs boson to two gluons, H → gg. 2 We use the notations introduced in the previous section to write where n f is the number of massless quarks. We consider the three terms in Eq. (3.1) separately, starting with the realemission contribution F LM (1 g , 2 g , 3 g ). The first step is to identify all possible singularities that may appear in the computation of that contribution and to partition the phase space in such a way that for each partition only a small subset of singularities is present. An important consequence of any partitioning is the fact that certain partons are identified as "hard". This means that, for a given partition, we should know exactly which partons cannot produce infra-red singularities. Although there are many ways to construct partitions, we find it convenient to use scalar products of the gluons' four-momenta s i j = 2 p i · p j and the energy-momentum conservation We then write where we have introduced the notations i j ≡ s i j /m 2 H . We can use the symmetry of the matrix element and the phase space to rewrite this equation as (3.4) Thanks to the prefactors 12 , gluons g 1 and g 2 on the righthand side of Eq. (3.4) must be hard or "resolved" and the only potentially unresolved parton is the gluon g 3 . This means that the right-hand side of Eq. (3.4) is singular when g 3 is soft and when g 3 is collinear to either g 1 or g 2 ; it is, however, not singular when either g 1 or g 2 is soft or when g 1 and g 2 are collinear to each other. We can follow the approach described in the context of color singlet production [30,33] to extract singularities from the right hand side of Eq. (3.4). We begin by considering the soft contribution that arises when energy of the gluon g 3 , E 3 , becomes small. We find |M tree | 2 (1 g , 2 g ), (3.5) where C A = N c = 3 is the SU (3) color factor and g s,b is the bare strong coupling. The factorization formula Eq. (3.5) allows us to extract contributions of soft singularities from the decay rate. To do so, we introduce the soft operator S 3 that extracts the most singular contributions in the soft limit from the matrix element squared and the relevant phase space: , (3.6) Note that the function θ(E max − E 3 ) prevents the integral over E 3 from becoming unbounded from above. We rewrite Eq. (3.6) as where we defined 3 and (3.10) We note that in the H → gg decay discussed here η 12 = 1; however, we do not use this fact right away and write Eq. (3.8) in a more general way. The calculation that we just described allows us to remove the soft singularity. We obtain (3.11) We note that, since the reduced matrix element does not require further regularization, all singularities in the first term on the r.h.s. of Eq. (3.11) are explicit. The second term there is free of soft singularities, but it still contains collinear ones; these occur when η 31 = (1 − cos θ 31 )/2 or η 32 = (1 − cos θ 32 )/2 vanish. To isolate these singularities, we partition the phase space in such a way that only one of them can occur at a time. To this end, we introduce the partition of unity  32 , ω 32 = η 31 η 31 + η 32 . (3.13) Introducing this angular partitioning, we write (3.14) and consider the two terms in the sum separately. We start with the i = 1 term. Similarly to the soft case, we introduce a C 31 operator which extracts the corresponding collinear singularity, and apply it to ω 31 (I − S 3 )s 12 F LM (1 g , 2 g , 3 g ) . We define C 31 in such a way that it extracts the leading η 31 → 0 singularity from F LM (. . .) δ without acting on the phasespace elements [d f 1,.,3 ], see Ref. [33] for more details. We find In Eq. (3.15), we defined and denoted an on-shell gluon with momentum p 13 as 13 g . The function P gg in Eq. (3.15) stands for the g * → gg splitting function and we used the ⊗-sign to indicate its spin-correlated product with the matrix element squared, see Refs. [30,33] for details. In these references, we explicitly showed that at NLO spin correlations disappear after azimuthal averaging. As the result, Eq. (3.15) becomes where P gg is the spin-averaged g * → gg splitting function The term on the second line of Eq. (3.17) is very similar to F LM (1 g , 2 g ) δ . To make this similarity explicit, we change integration variables from E 1 and E 3 to E 13 and z = E 1 /E 13 . We obtain We also rename f 13 back to f 1 and obtain where E 1 = m H /2 and we used the fact that the integration over z starts at z = z min = min{0, 1− E max /E 1 }. Since E max must be chosen in such a way that the whole phase space is covered, E max should be larger than E 1 , E max > E 1 , for all E 1 . This implies z min = 0. Repeating these steps for the soft-collinear term S 3 C 31 , we find We use these results to write where C 31 (I − S 3 ) follows from Eqs. (3.20, 3.21). We find where γ 22 z,g→gg is a particular case of a general anomalous dimension defined as follows (3. 24) We note that in the first term on the right hand side in Eq. (3.22) all singularities are manifest and the reduced matrix element does not require regularization, whereas the second term is free of both soft and collinear singularities so that it can be immediately integrated in four dimensions.
We deal with the ω 32 term in the partition of unity Eq. (3.12) in a similar way. We obtain We combine Eqs. (3.11, 3.23, 3.25) and obtain the following result for the three-gluon contribution to Higgs boson decay We note that, thanks to Bose symmetry, the two terms in the sum in the last line in Eq. (3.27) are the same. Hence, we write (3. 28) This discussion implies that Bose symmetry can be efficiently used to partition the phase space in such a way that identical kinematic configurations of the three-gluon final states are accounted for only once in the calculation; this removes the original 1/3! symmetry factor.
Before combining this result with virtual corrections, we consider the other real-emission term in Eq. (3.1), n f F LM (1 q , 2 g , 3q ) , that describes the decay H → (g * → qq)g. Because (in this section) the qq pair does not directly couple to the Higgs boson, the singularity in this case is produced by the collinear splitting g * → qq. For this reason, we do not need any partitioning. We repeat steps that led to Eq. (3.22) and obtain 5 We can now combine the H → ggg and H → qgq decay channels and write the total real-emission contribution to d NLO , up to higher orders in , as where we have defined The two quantities γ g and γ g are given in Eq. (A.7).
It remains to combine Eq. (3.31) with virtual corrections. We follow Ref. [35] to separate the divergent and finite parts of the one-loop amplitude and define where F fin LV (1 g , 2 g ) δ is a finite remainder of the one-loop H → gg amplitude, see Appendix A in Ref. [33] for details. We combine Eq. (3.31) and Eq. (3.33), use η 12 = 1 and obtain a very simple result for the NLO QCD corrections to H → gg decay. It reads where the two contributions are defined in Eq. (2.9). We conclude this section by reminding the reader that the NLO construction we just described is identical to the FKS subtraction scheme [36,37]. In the next sections, we will show how to generalize the FKS scheme to NNLO.

Higgs decay to gluons: a NNLO computation
In this section we generalize the discussion of the NLO QCD corrections to the decay of a color singlet to the NNLO case. We will follow Refs. [30,33] and perform subtractions of soft and collinear divergences in an iterated manner, starting from the soft ones. Many technical details are similar to the production case described at length in the above references and we do not discuss them here. Instead, we focus on the peculiarities of the decay.

Double-real contribution
There are four different partonic final states that we have to consider. They are a) 4 gluons, b) 2 gluons, 2 quarks, c) two quark pairs of different flavors and d) two quark pairs of the same flavor. We write In full analogy to the NLO case, we partition the phase space in such a way that only a subset of partons are allowed to become unresolved. In case of the NNLO contributions, two partons can become unresolved simultaneously; we will systematically rename partons so that, eventually, the unresolved partons are always referred to as f 3 and f 4 . We first consider the four-gluon channel, and introduce a partition of unity following what has already been done at NLO 1 =s 12 +s 13 +s 14 +s 23 +s 24 +s 34 .
We insert this partition inside the integrand for F LM (1 g , 2 g , 3 g , 4 g ) δ , use the symmetry of the phase space and the matrix element and arrive at 6 The prefactors 12 ensures that no singularity arises in the products 12 |M(1 g , 2 g , 3 g , 4 g )| 2 when gluons 1 and 2 become either soft or collinear to each other. To proceed further, we introduce an energy ordering for potentiallyunresolved gluons g 3 and g 4 , use g 3 ↔ g 4 symmetry and write ) . (4.5) 6 In this subsection, the "tree" superscript on M is always assumed.
We now consider the 2q2g final state. In principle, it contains fewer singularities than the four-gluon final state. Therefore, one may use a simpler partition of unity to single out the potentially unresolved partons. However, to streamline the bookkeeping, we find it convenient to use identical partitioning for all final states. Our starting point is then where the partition of unity Eq. (4.3) has already been employed. We note that the amplitude is symmetric with respect to permutations of the two gluons, so that Furthermore, since in this amplitude the quark-antiquark pair arises from gluon splitting, the amplitude squared summed over quark and anti-quark polarizations satisfies We can use these symmetries of the amplitude squared as well as the symmetry of the phase space to re-write Eq. (4.6) in the following way (4.8) To proceed further, we introduce the energy ordering for the two potentially unresolved partons f 3,4 and use symmetries of the amplitude to remove the factor 1/2 in the above equation. In cases when f 3 and f 4 are partons of a different type, this requires us to combine the different contributions in a particular way. As an example, consider the second and the third term in Eq. (4.8). Relabelling parton momenta where appropriate, we write (4.9) Using these transformations, we obtain which we can write as (4.11) We note that the six terms in Eq. (4.11) have very different singularity structures. For example, all the terms in Eq. (4.11) that contain gluon g 4 give rise to single soft singularities that arise when E 4 → 0. In the remaining three terms, the energy E 4 is associated with an anti-quark and, therefore, these terms are not singular in the single-soft limit. Similarly, the collinear limit C 41 corresponds to an (anti)quark and a gluon becoming collinear in the first, second, fifth and sixth terms in Eq. (4.11). However, the same limit describes a kinematic configuration with two collinear gluons in the third term in Eq. (4.11). Clearly, the two limiting cases result in different splitting functions and different reduced matrix elements. Finally, we turn to the four-quark channels, where we need to make a further distinction between cases when quarks have same or different flavors. If they are different, i.e. q = q , we write If the flavors are identical, we can use the same amplitude M as for the different-flavor case, accounting for a permutation of two identical particles. We write (4.13) We denote the interference term as (1 q , 2 q , 4q , 3q ) , (4.14) and write the complete four-quark contribution to the decay rate, including both different and identical flavors, as The interference term in Eq. (4.15) is not singular and can be evaluated in four dimensions; for this reason we keep it as it is. Moreover, the first term in that equation only produces singularities when either one or two qq pairs become collinear. Despite this simplicity, we find it convenient to treat the four-quark contributions Eq. (4.15) in the same way as the two other channels that we discussed previously. To this end, we insert the partition of unity Eq. (4.3) into the integrands in Eq. (4.15), re-label partonic momenta, use the symmetry of the amplitude squared and obtain Int(1, 2, 3, 4). (4.17) The prospective unresolved partons are f 3,4 . Similar to other channels, we introduce the energy ordering E 3 > E 4 and again use the symmetry of the amplitude squared to simplify the result. We obtain where we have defined Upon combining all the channels, we obtain the final result for the double-real contribution to the decay width. It reads (4.20) To illustrate how soft and collinear singularities are extracted from the double-real emission contribution Eq. (4.20), we focus on the four-gluon final state F LM (1 g , 2 g , 3 g , 4 g ). This contribution possesses the richest singularity structure yet, at the same time, it is one of the simplest as far as the bookkeeping is concerned. After explaining how the singularities are extracted in this case, we present the results for all channels in Sect. 4.4.

Double-soft contribution for H → gggg
Similar to the production case, we begin with the double-soft limit that occurs when E 3 , E 4 → 0. We follow Refs. [30,33] and introduce an operator SS that extracts the leading double-soft singularity from the product of the matrix element squared and the phase space, and write (4.21) The double-soft limit is computed in exactly the same way as in the production case [30,33]. We find 12 (4.24) The term on the second line in Eq. (4.24) does not contain double-soft singularities anymore but it still contains both single-soft and collinear ones. We discuss how to extract them in what follows.

Single-soft contribution
We need to extract the single-soft singularity from see Eq. (4.24). The soft limit of the amplitude squared reads and Equation (4.27) is free from soft singularities, but it still contains collinear ones; these arise when the gluon g 3 becomes collinear to gluon g 1 or gluon g 2 . We proceed as in the NLO computation. Namely, we introduce a partition of unity, use the symmetry of the process under the exchange of gluons 1 and 2 and write (4.30) All singularities are regulated in the second term on the r.h.s. of Eq. (4.30). We now consider the first term on the r.h.s. of Eq. (4.30). Taking the η 31 → 0 limit in J ggg S 4 , we obtain where we used Since we have to apply the C 31 operator to 3 J ggg S 4 (I − S 3 )s 12 F LM (1 g , 2 g , 3 g ) δ and since the limit of F LM (1 g , 2 g , 3 g ) is identical to what we already discussed in the NLO case, the computation proceeds similarly to the NLO case. Note that since the (E 2 3 ) − prefactor in J S 4 gives an extra factor (1 − z) −2 the relevant anomalous dimension in this case is γ 24 z,g→gg , c.f. Eq. (3.24). The result of the calculation reads (4.33) We combine the different contributions and obtain (4.34) In Eq. (4.34) the third and fourth lines are free from unregulated singularities whereas the second line contains unregulated collinear singularities that need to be extracted. We explain how to do that in the next section.

Collinear singularities: general structure
Having regulated all the soft singularities, we are left with only one contribution on the right hand side of Eq. (4.34), where w 3i,4 j are functions of the partons' emission angles. These functions are constructed in such a way that a product of w 3i,4 j with the matrix element squared has non-integrable collinear singularities if gluon g 3 is collinear to gluon g i or gluon g 4 is collinear to gluon g j . The singularities that arise when gluons g 3 and g 4 become collinear can only occur in the partitions w 31,41 and w 32,42 . Following Refs. [30,33], we refer to w 31,41 and w 32,42 as the triple-collinear partitions and w 31,42 and w 41,32 as the double-collinear partitions. A possible choice for these functions is given in Appendix A.
The double-collinear partitions can be dealt with in a relatively straightforward way since the collinear singularities are clearly isolated. The only issue that we need to address is the construction of a proper phase space for this contribution; we discuss how this can be done in Appendix B. For the triple-collinear partitions, we need to order the emission angles of gluons g 3 and g 4 and we refer to these orderings as "sectors" that we label as a, b, c, d, see Refs. [30,33] for details. Explicitly, we write (4.37) Once partitions and sectors are introduced, we can extract the collinear limits from the decay rates following the procedure already discussed for the production case [30,33]. Note, however, that similar to the NLO computations discussed in Sect. 3, we have to integrate the various splitting functions appearing in the calculation over energies to obtain (generalized) anomalous dimensions. We now summarize the relevant steps for the extraction of the collinear singularities, closely following the procedure and notation of Ref. [30,33]. We introduce the short-hand notation G(1, 2, 3, 4) ≡ 12(I − SS)(I − S 4 )s 12 (4.40) • the soft-regulated triple-and double-collinear contribution, defined as (4.41) • and, finally, the soft-regulated collinear-regulated term (4.42) We remind the reader that the notations in Eqs. (4.40, 4.41, 4.42) are such that collinear operators act on everything that appears to the right of them. In particular, the notation · · · C[d f i ] · · · implies that a particular collinear limit should be taken in the phase-space element of the parton f i . More details can be found in Refs. [30,33] where we show an explicit parametrization of the emission angles for gluons g 3 and g 4 and define the action of collinear operators in terms of this parametrization. We discuss the terms G s r ,c s (1, 2, 3, 4) and G s r ,c t (1, 2, 3, 4) in the next two subsections. The term G s r ,c r (1, 2, 3, 4) is finite and can be immediately computed in four dimensions. This point is again discussed in Refs. [30,33] in the context of color singlet production. Since there is no conceptual difference between how this contribution is computed in the production and decay cases, we won't repeat the discussion here.

Soft-regulated single-collinear contribution
To obtain an expression for the soft-regulated single-collinear contribution G s r ,c s (1,2,3,4) in Eq. (4.40), we follow the same steps as in the production case [30,33]. After a tedious but otherwise straightforward computation we obtain 7 G s r ,c s (1,2,3,4) The anomalous dimensions γ i j z,g→gg , that appear in Eq. (4.43), are defined in Eq. (3.24) whereasγ g ( ),γ g ( , k ⊥ ) and δ g ( ) can be found in Refs. [30,33]. For completeness, we report them in Appendix A, see Eqs. (A.10, A.11, A.12). Finally, γ 24,r z,g→gg is defined as (4.44) Note that, as a consequence of spin correlations, the result in Eq. (4.43) contains a finite term r μ r ν F μν LM . This term should be understood as the corresponding matrix element squared where the polarization vector for the gluon g 3 is taken to be a particular four-vector r μ . The precise form of the vector r depends on the specific way in which the limit where gluons g 3 and g 4 become collinear is approached. Since we use the same parametrization of the triple-collinear phase space as in Refs. [30,33], the explicit form of the vector r μ can be taken from these references. As an example, consider the ω 31,41 partition, where p 3 is written as (4.45) Here, θ 31 is the relative angle between the momenta of g 1 and g 3 . Upon parametrizing the collinear limit of g 3 and g 4 as described in Refs. [30,33], we find the following expression for the vector r μ r μ = (0, − cos θ 31 cos ϕ 3 , − cos θ 31 sin ϕ 3 , sin θ 31 ). (4.46) Similar to Refs. [30,33], damping factors with tildes in Eq.

Soft-regulated triple-and double-collinear contribution
We now discuss the triple-and double-collinear contribution G s r ,c t (1,2,3,4) shown in Eq. (4.41). As indicated in the previous section, this term includes all the double-unresolved collinear contributions which arise when both gluons g 3,4 are collinear to either gluon g 1 or gluon g 2 , as well as singlecollinear contributions where gluons g 3 and g 4 are collinear to each other. This contribution requires a non-trivial integration of the triple-collinear splitting function over energies and angles of particles that participate in the splitting. The relevant computation was performed in Ref. [32]. Using the results presented there, we can write the final result for the soft-regulated tripleand double-collinear contribution as where the first term on the right hand side comes from doublecollinear configurations and the second one from the triplecollinear ones.

Real-virtual contribution
We now turn to the discussion of real-virtual contributions. Their calculation is similar to the NLO case discussed in Sect. 3. As in the previous section, we illustrate the most important steps of the real-virtual calculation for the three-gluon final state. Similar to NLO, we introduce a phase-space partitioning and write (4.52) We note that the 1 ↔ 2 symmetry was used to simplify Eq. (4.52). The first term on the right hand side of Eq. (4.52) is fully regulated. The terms on the second line are soft and collinear subtractions, which we now discuss. The starting point for the calculation of the soft subtraction contribution is the factorization property of the one-loop amplitude [38], that leads to with The appearance of the β 0 term in Eq. (4.53) is related to the fact that we work with UV-renormalized amplitudes. Starting from Eq. (4.53), we follow the discussion presented in Sect. 3 and obtain (4.55) Next we consider the collinear subtraction. At one-loop, the collinear factorization of one-loop amplitudes leads to [39] where s 13 = 2 p 1 · p 3 + i0. We remind the reader that the notation "13 g " indicates a gluon that has the same direction as the gluon g 1 but whose energy E 13 is given by E 13 = E 1 + E 3 . As in Sect. 3, the symbol ⊗ in Eq. (4.56) indicates a contraction of the one-loop spin-correlated splitting function P (1) gg with the relevant scattering amplitudes. The one-loop splitting function P (1) gg was computed in Ref. [39]; we report it in Appendix A for convenience. We note that, at variance with the production case, the splitting function P (1) gg is manifestly real for the decay kinematics. Following the same steps as in the NLO calculation described in Sect. 3, we obtain In Eq. (4.57), γ 1−loop z,g→gg is the one-loop anomalous dimension, analogous to γ z,g→gg , obtained by integrating P (1) gg in Eq. (4.56) over the energy fraction E 1 /E 13 . Its explicit expression is reported in Appendix A.
Finally, it is also convenient to explicitly extract the 1/poles from the F LV terms in Eqs. (4.52, 4.55, 4.57). Their structure is well-known [35] and we have already discussed it in Refs. [30,33] using our notations. For completeness, we report the relevant formulas below (4.58) In Eq. (4.58) the functions F fin LV are finite in four dimensions and s i j = 2 p i · p j > 0.

Double-virtual corrections
The double-virtual contribution is identical to those in the production case described in Refs. [30,33]. For convenience, we report the relevant formulas here. Following Ref. [35], we extract all the -poles from the loop amplitudes and write wherẽ with β 0 defined in Eq. (4.54) and T R n f , Finally, we note that F fin LV (1 g , 2 g ) is defined in Eq. (3.33), and F fin LVV (1 g , 2 g ) , F fin LV 2 (1 g , 2 g ) are finite remainders, see Appendix A in Ref. [33] for details.

Final result
The sum of the different contributions discussed in the previous sections gives a result that is finite in the → 0 limit. Repeating similar calculations for all the other partonic channels, we obtain the full NNLO QCD corrections to the decay H → gg. We write the result as the sum of contributions with different final state multiplicities, cf. Eq. (2.9) The contribution of the four-parton final state reads where F LM (1,2,3,4) is defined in Eq. (4.20). Similarly, the three-parton contribution reads wherê and J i jk = J (1) i jk + J (2) i jk . (4.66) The functions J (1,2) are defined as and J (2) i jk = γ i + γ j +γ k −ω 31,41 4||1 ln where C q = C F and C g = C A . The various constants and functions used in Eqs. (4.67, 4.68) can be found in where γ q is given in Eq. (A.8) and F fin LV (1 b , 2b) is a finite virtual remainder analogous to F LV (1 g , 2 g ) in Eq. (3.33), see Appendix A in Ref. [33] for its explicit definition.
At NNLO, we also do not require any additional partitioning except perhaps for the 4b final state that arises from the prompt decay of the Higgs boson. 9 We show in Appendix C that the contribution of this subprocess to the decay rate can be written as a sum of two terms: a term that coincides with the contribution of the decay H → bbqq, q = b, where only b andb can be prompt and the qq pair originates from gluon splitting, and an interference term. The first term can be treated without any partitioning since the hard partons are always the two b-quarks. The interference term has only a triple-collinear singularity that maps onto the corresponding splitting function. Its proper treatment is described in Appendix C.
The NNLO contribution to H → bb decay is then computed following the steps discussed in the previous section. We write the NNLO contribution as a sum of "fixedmultiplicity" terms d NNLO Q→i , i = 2, 3, 4. The four-parton contribution reads where now with F int LM defined in Appendix C. The three-parton contribution reads 9 To avoid confusion, we emphasize that in the previous section a 4q final state originating from the decay H → (g * → qq) (g * → qq) was discussed whereas in this section we consider prompt decays to fermions. For this reason, the 4b final state originates from e.g. (5.4) where the function J qqg is defined in Eq. (4.66) and γ k ⊥ ,g can be found in Appendix A. Finally, the two-parton contribution reads  (5.5) where i j are defined in Eqs. (4.48, 4.49).

Validation of results
In this section, we use the analytic formulas for the fullydifferential decay rates presented above to calculate the NNLO QCD corrections to decays H → gg and H → bb. 10 We compare these results with analytic formulas extracted from Refs. [42][43][44] to validate our calculations. 11 We begin with the decay process H → gg, which was discussed in Sect. 4. We consider a Higgs boson of mass m H = 125 GeV which couples to gluons through the effective Lagrangian L Hgg = −λ Hgg H G (a) μν G μν,(a) , (6.1) where in the MS scheme with α s = α s (μ) being the renormalized coupling in a theory with 5 massless flavors and v is the Higgs vacuum expectation value, see e.g. Ref. [65]. For the numerical results presented below, we use m t = 173.2 GeV. For numerical checks, we split the width for the H → gg decay into different color factors, which allows us to check different partonic channels separately. We write where the LO decay width that has been factored out is given by LO . The comparison between our results for the NLO and NNLO coefficients R (1,2) and those presented in Ref. [42] is given in Table 1. We present numerical results for a scale μ = 2m H , in order to avoid accidental cancellations between the renormalization scale μ and the Higgs mass m H that happen for μ = m H . We observe agreement well below the per mille level for all coefficients.
We turn now to the decay H → bb. Again, we consider a 125 GeV Higgs boson and five flavors of massless quarks, which allows us to use the results presented in Sect. 5. The Higgs couples to bottom quarks only, through a Yukawa interaction The comparison between the coefficients S (1,2) obtained from our numerical code and from the analytic formulas of Ref. [44] are displayed in Table 2. Again, we use the scale μ R = 2m H for this comparison. The agreement is consistently below the per mille level across all color structures. Finally, we compare exclusive jet rates for the H → bb decay with those reported in Ref. [44]. To do so, we use the JADE clustering algorithm with y cut = 0.01 and the distance measure defined as y i j = ( p i + p j ) 2  Clearly, the level of numerical precision achieved for the NNLO coefficients in our calculation is excessive since for phenomenological applications it is enough to know widths with sub-percent accuracy. We note that to achieve this level of numerical precision within our framework, one would typically require up to one CPU hour of computation time.

Conclusion
We presented analytic formulas that describe fully-differential decays of color-singlet particles to qq and gg final states through NNLO QCD. The results are obtained within the nested soft-collinear subtracted scheme that we proposed earlier in Ref. [30]. The results are remarkably compact and simple to implement in a numerical code. We have validated these results by computing the NNLO QCD corrections to the H → gg and H → bb decay rates and comparing them to independent numerical and analytic computations finding per mille level agreement for observables that are known analytically. In addition to their phenomenological relevance for decays of the Higgs boson and electroweak vector bosons, such as H → gg, H → bb and V → qq, these results provide an important building block for the extension of the nested soft-collinear subtraction scheme which will make it applicable for computations of NNLO QCD corrections to arbitrary processes at hadron colliders.
with κ ⊥ = k ⊥ / −k 2 ⊥ . The transversal metric tensor g μν ⊥ and the transversal vector k ⊥ are defined relative to the fourmomentum of the collinear gluon, in the standard way [38]. The d−dimensional spin averages of the splitting functions give We use these results to construct the spin-averaged splitting function P (1) gg ; we then integrate it over z to obtain the anomalous dimension γ 1−loop z,g→gg following a similar procedure to the one described for the tree-level splitting function, see the discussion leading to Eq. (3.24). We find γ Finally, the function K i j reads (B.12) The energy E 2 is then obtained from energy conservation. Finally, we write the phase space for parton f 4 in terms of the angle θ 13,4 [d f 4 ] = dE 4 E 1−2 Fig. 1 Diagrams for H → bbbb decay. The gluon emission off the bubble in Hbb vertex describes diagrams where the gluon is emitted from one of the outgoing quark lines There are four subamplitudes that contribute to this process; they are shown in Fig. 1. The difference between these amplitudes is in the fermion lines that originate from the Hbb vertex and the ones that originate from the gluon splitting, g * → bb. It is clear that b-quarks from the Hbb vertex are hard, in a sense that they cannot produce infra-red singularities, whereas b-quarks from gluon splitting can be soft. Since whether a given b-quark is hard or soft changes from diagram to diagram, the extraction of singularities becomes intricate.
To overcome this problem we make use of the symmetries of the H → bbbb decay. To this end, we write the matrix element as the sum of four subamplitudes shown in The H → bbbb decay width reads In the squares of amplitudes, the choice of (potentially) hard and soft fermions is unambiguous. We label the hard momenta as 1 and 2 and the soft momenta as 3 and 4. Using the symmetry of the phase space, we obtain where we have included a factor of 4 for the four diagrams and another factor of 2 for the energy ordering E 3 > E 4 . It is straightforward to extract the various singularities from this contribution; in fact the result is identical to the qq contribution to NNLO QCD corrections to H → bb decay. The interference terms in Eq. (C.4) are more involved since it is not possible to choose hard and soft momenta unambiguously. Before discussing this, we note that since helicities of massless quarks are conserved and since H → bb and g * → bb produce quarks with different (same) helicities, respectively, the interferences of diagrams (a) and (d) m ad as well as diagrams (b) and (c) m bc vanish. We then classify the possible collinear divergences in the remaining interference contributions. We find the following divergences in various interference terms:  Using these results, we can write 2b, 4 b , 1b) . (C.8) By construction, c.f. Fig. 2, the interference contributions are only singular in the limit when momenta of partons f 1 , f 3 and f 4 become collinear, whereas the non-interference term has multiple singularities, including the double-soft one that occurs when f 3 and f 4 become soft. We denote the interference term as and use this notation in the main text when we discuss the computation of NNLO QCD contribution to Higgs decay to two quarks in Sect. 5. The non-interference term in Eq. (C.8) is accounted for as part of the n f -dependent contributions in the NNLO QCD computation.