Higgs boson decay into b-quarks at NNLO accuracy

We compute the fully differential decay rate of the Standard Model Higgs boson into b-quarks at next-to-next-to-leading order (NNLO) accuracy in alpha_S. We employ a general subtraction scheme developed for the calculation of higher order perturbative corrections to QCD jet cross sections, which is based on the universal infrared factorization properties of QCD squared matrix elements. We show that the subtractions render the various contributions to the NNLO correction finite. In particular, we demonstrate analytically that the sum of integrated subtraction terms correctly reproduces the infrared poles of the two-loop double virtual contribution to this process. We present illustrative differential distributions obtained by implementing the method in a parton level Monte Carlo program. The basic ingredients of our subtraction scheme, used here for the first time to compute a physical observable, are universal and can be employed for the computation of more involved processes.


Introduction
In run I, the ATLAS and CMS collaborations of the Large Hadron Collider (LHC) discovered a new particle [1,2] with quantum numbers corresponding to those of the Higgs boson in the Standard Model (SM) within the experimental accuracy of the measurements [3][4][5][6]. Thus by now it is widely accepted that the new particle is the Higgs boson of the SM. Nevertheless, further more precise measurements are being prepared for the upcoming run II. In particular, a lot of emphasis is put on the precise determination of the couplings of the Higgs boson to the heavy fermions to check whether the fermion masses are consistent with fermion mass generation in the SM.
Since the b-quark is quite light (its mass is only about 2 % of the vacuum expectation value of the Higgs field), the rate of associated production of a b-quark pair with a Higgs boson is rather low. This fact, together with the overwhelming number of background events coming from direct QCD b-quark pair production makes the determination of the b-quark Yukawa coupling through Hbb production impossible. A better option that gives direct access to the Hbb Yukawa coupling is to measure the H → bb decay in the associated production of a Higgs boson with a W or a Z boson in a boosted or semi-boosted regime [7]. In this scenario it is possible to use the kinematic and topological properties of the final states to isolate the H → bb decay. In this respect, first measurements have been performed by the CMS [8] and ATLAS [9] collaborations. Such search strategies may be aided by accurate modeling of QCD radiation in the H → bb decay, which motivates the computation of the fully differential decay rate at next-to-next-toleading order (NNLO) accuracy in QCD perturbation theory. Computing fully differential cross sections and decay rates at NNLO turns out to be rather involved, however the last decade has witnessed substantial development  leading to a number of differential results for specific processes .
The first computation of the fully differential decay rate of the SM Higgs boson into b-quarks at NNLO accuracy was published in ref. [47]. That computation was performed with the method of sector decomposition based on non-linear mappings [13]. Here we offer a different approach based on the numerical implementation of the general subtraction scheme developed in a series of papers for the computation of QCD jet cross sections at NNLO accuracy [31][32][33][34][35][36][37][38][39][40][41]. This method, which is used for the first time in this paper to compute a physical observable at NNLO, employs the universal infrared factorization of QCD squared matrix elements to define local subtraction terms for regulating the singularities emerging in unresolved real radiation.
Specifically, we can write the NNLO correction to the cross section of a generic m-jet process as a sum of three contributions, the tree level double real radiation, the one-loop plus a single radiation, and the two-loop double virtual terms of the basic process under consideration, and rearrange it as follows, where, The subscripts on the integral signs are simply reminders that the integration is over the phase space of n = m, m + 1 or m + 2 final state particles. Above J n denotes the value of some infrared-safe observable J evaluated on an n parton final state. The right-hand sides of eqs. (1.3) and (1.4) are integrable in four dimensions by construction [31][32][33][34], while the integrability of eq. (1.5) in four dimensions is ensured by the Kinoshita-Lee-Nauenberg (KLN) theorem on infrared-safe quantities, provided that our subtraction scheme is well defined.
The counterterms which contribute to dσ NNLO m+2 and to dσ NNLO m+1 were introduced in refs. [33] and [34]. The integration of the real-virtual counterterms (the last two terms of eq. (1.5)) was performed in refs. [35,36,38]. The integral of the iterated single unresolved counterterm (the third term of eq. (1.5)) was computed in ref. [39]. The integration of the collinear-type contributions to the double unresolved counterterm (the second term of eq. (1.5)) was performed in ref. [40]. The softtype contributions to the same counterterm were presented in ref. [41]. Most of these results were given as expansions in whose coefficients were computed numerically. Here we present the relevant integrals with pole coefficients evaluated analytically, while the finite parts are given numerically. The final test on the consistency of our subtraction scheme is then to verify that eq. (1.5) is free of singularities, as prescribed by the KLN theorem. In this paper, we perform that check analytically for the first time by computing the fully differential decay rate 1 of the Higgs boson into b-quarks at NNLO.
The present work is the first physical application of this method, therefore in order to facilitate reading we present the full computation as implemented in a parton level Monte Carlo program in detail. As usual in such codes, the jet function J is computed from generated momenta in d = 4 dimensions, therefore, the implementation of any infrared-safe physical quantity is straightforward as demonstrated here.
The paper is organized as follows: in section 2, the notation and conventions we use are introduced; in sections 3 and 4, we show the decay width at leading order and next-to-leading order (NLO) accuracy in α s ; in section 5, we display the counterterms and the insertion operators which are necessary to define the double real (1.3) and the real-virtual (1.4) contributions to the decay width, and we show that the double virtual contribution (1.5) is free of singularities; in section 6, we show a selection of illustrative results; we draw our conclusions in section 7. The two appendices provide details on the matrix elements we use, as well as on the insertion operator used in the NLO computation.

Notation
We consider the partial decay width Γ H→bb [J] of the Higgs boson into a b-quark pair, for any infrared-safe observable J. Through NNLO in QCD, this decay width receives contributions from the following partonic subprocesses: two-loop where we show also the four-momenta of the particles in parentheses. We report the matrix elements corresponding to all subprocesses up to the required loop level in appendix A. We use the colour and spin space notation of ref. [77] where the matrix element for a given subprocess, |M n , is a vector in color and spin space, normalized such that the squared matrix element summed over colours and spins is given by where n is the number of particles in the final state. The matrix element has the following formal loop expansion with the dots denoting higher-loop contributions. We will always consider matrix elements computed in conventional dimensional regularization (CDR) with MS subtraction. We will also use the following ⊗ product notation to indicate the insertion of colour charge operators between M ( 1) | and |M ( 2) : We use the customary normalization of T R = 1/2 for the colour-charge operators, thus the quadratic The b-quark mass is much smaller than the scale of the problem that is the Higgs boson mass, therefore, we treat the b-quarks as massless, both in the matrix elements and phase space integrals, retaining the b-quark mass only in the Yukawa coupling. We neglect the t-quark throughout and consider n f = 5 light quark flavours.
In QCD the renormalized amplitudes are obtained from the unrenormalized ones by replacing the bare couplings y B b and α B s with their renormalized counterparts evaluated at the renormalization scale µ and S MS = (4π) exp(− γ E ) corresponds to MS subtraction. Although the factor (4π) exp(− γ E ) is often abbreviated as S in the literature, we reserve the latter to denote On the right-hand side, y b ≡ y b (µ) and α s ≡ α s (µ) are the dimensionless renormalized couplings in the MS scheme evaluated at the renormalization scale µ. The n particle massless phase space measure reads Throughout the paper, we will use y ik to denote twice the dot-product of two momenta, scaled by the total momentum squared Q 2 . For example, We also introduce the combination for later convenience.

Leading order
Let us denote the Born differential decay rate by, Then the leading order decay width is, Here J is an infrared-safe observable whose value evaluated on a kinematic configuration with two partons is J 2 . For the inclusive decay width (J ≡ 1) at leading order we have where the expression on the right-hand side is the four-dimensional result.
4 Next-to-leading order

Real emission contribution
The real emission contribution to the differential decay width reads dΓ R 3 is divergent when the radiated gluon becomes unresolved (soft, or collinear with one of the b-quarks). In order to regularize it, we subtract an approximate decay rate, where the counterterm for processes with m + 1 partons in the final state is given by [32,33], In eq. (4.3) the functions C (0,0) ir and S (0,0) r appearing in the right-hand side correspond to counterterms which regularize the p i ||p r collinear limit and the p r → 0 soft limit. In order to avoid double counting in the overlapping soft-collinear region, we must add back a soft-collinear counterterm, C ir S (0,0) r . The precise definitions of these subtractions are given in refs. [32,33]. In our convention the indices of C Since the sums over i and r in eq. (4.3) are likewise not ordered, the factor of 1 2 assures that we count each collinear limit precisely once. Finally, the superscript ( 1 , 2 ) means that the corresponding counterterm involves the product (in colour or spin space) of an 1 -loop unresolved kernel (an Altarelli-Parisi splitting function or a soft eikonal current) with an 2 -loop squared matrix element. Thus, (0, 0) means that we consider a tree level collinear or soft function acting on a tree level reduced matrix element. Such superscripts will appear also for other counterterms throughout the paper. For definitiveness, we spell out eq. (4.3) explicitly for H → bbg (m = 2) below, where the b,b and gluon carry the labels 1, 2 and 3.
With the counterterms given in refs. [32,33] it is straightforward to check that the difference is integrable in all kinematic limits. Then, the regularized real contribution to the decay rate, is finite in four dimensions for any infrared-safe observable. An explicit calculation for the contribution to the total decay width from the real emission part plus subtractions yields (4.7)

Virtual contribution
The virtual contribution to the differential decay width reads and is of course divergent in four dimensions. Its -expansion reads (see eq. (A.2)) where we have introduced the abbreviation L = ln µ 2 m 2 H . In eq. (4.9), dΓ B denotes the ddimensional Born decay rate as given in eq. (3.1). By the KLN theorem, the integral of the approximate decay rate precisely cancels the divergences of the virtual piece, so adding back what we have subtracted from the real correction, the virtual contribution becomes finite as well. We have performed the integration of the various subtraction terms analytically in ref. [32] and here we only quote the result, which can be written as, where the ⊗ product is defined in eq. (2.2) and the insertion operator is in general given by [32] 2 The variables y iQ and Y ik,Q were defined in eqs. (2.8) and (2.9) and Q µ is the total incoming momentum. The functions C (Y ik,Q ; ) have been computed as Laurent expansions in in ref. [32] and are recalled here up to finite terms in appendix B. We mention that there is no one-to-one correspondence between the unintegrated subtraction terms in eq. (4.3) and the kinematic functions that appear in eq. (4.11). The latter are obtained from the former after summing over all unobserved quantum numbers (colour and flavour) in addition to integrating over the unresolved momentum, and organizing the result in colour and flavour space. Loosely speaking, the integrated form of C  the integrated form of C ir S (0) r to either of the integrated counterterms and this final organization was performed differently in ref. [32] and in this paper. In ref. [32], the integrated form of C ir S (0) r was grouped into S (0),(i,k) 1 , while here we find it more convenient to group it into C (0) 1,i . For H → bb, with only two partons in the final state the colour connections factorize completely, Furthermore, momentum conservation implies that Thus, the insertion operator I 1 becomes, where, as indicated, we must evaluate all functions with arguments equal to one. The Laurent expansion of eq. (4.14) in is, where, for future reference, we have also provided the O( ) part in terms of rational numbers and known transcendental constants. The uncertainty of the O( 2 ) numerical result, as well as those of all other numerical results we show affect the last quoted digit, unless specifically stated otherwise. It is easy to check that the expression is free of -poles. Hence is finite in four dimensions for any infrared-safe observable. For the contribution to the total width from the virtual part plus integrated subtractions we find Combining eqs. (4.7) and (4.18), we obtain the full NLO correction to the total decay rate, As C F = 4 3 in the conventions used, we recover the well-known NLO result [78][79][80].
5 Next-to-next-to-leading order

Double real emission contribution
The double real emission contribution to the differential decay width is and its integral over the phase space is divergent in four dimensions due to kinematic singularities emerging in unresolved regions. In order to regularize the singularities of eq. (5.1) due to two unresolved partons, we subtract an approximate decay rate, where the double unresolved counterterm for processes with m + 2 partons in the final state is [33] ir;s and S (0,0) rs denote counterterms which regularize the p i ||p r ||p s triple collinear, the p i ||p r , p j ||p s double collinear, the p i ||p r , p s → 0 one collinear, one soft (collinear+soft) and the p r → 0, p s → 0 double soft limits. The rest of the counterterms which appear in eq. (5.3) account for the double or triple overlap of limits, their role is to make sure that no multiple subtractions are performed in overlapping double unresolved regions. Thus, for instance, C irs CS (0,0) ir;s accounts for the triple collinear limit of the collinear+soft counterterm, and the rest of the counterterms have a similar interpretation as suggested by the notation. The precise definitions of all functions appearing in eq. (5.3) were given in ref. [33]. As in our convention the collinear indices of counterterms and the sums over them in eq. (5.3) are not ordered, the factors of where A 1 has been defined in eq. (4.3). To avoid double subtraction in overlapping single and double unresolved regions of phase space, we must also consider dΓ RR,A 12 4 The general formula for the iterated single unresolved counterterm is where the three terms above are given by [33], The interpretation of the various terms in eqs. (5.8)-(5.10) are suggested by the notation: for instance, C kt C (0,0) ktr in eq. (5.8) accounts for the p k ||p t single collinear limit of the C (0,0) ktr triple collinear counterterm, while, for example, S t C (0,0) irt in eq. (5.9) represents the counterterm appropriate to the p t → 0 soft limit of C However, very importantly, it can also be shown [33] that A 12 |M (0) m+2 | 2 simultaneously cancels the double unresolved singularities of the single unresolved subtraction term A 1 |M (0) m+2 | 2 and so properly accounts for the overlap of single and double unresolved subtractions. All of the counterterms appearing in eqs. (5.8)-(5.10) were precisely defined in ref. [33]. As before, the collinear indices and sums over them in eqs. (5.7)-(5.10) are not ordered, hence the appearance of the factors of 1 2 at various instances. With these definitions, the difference can be shown to be integrable in all kinematic limits [33]. Thus, the regularized double real contribution to the decay rate This numerical value has been obtained by implementing eq. (5.12) in a fully differential parton level Monte Carlo program using four dimensional double real emission matrix elements and phase space. However, we have also reproduced the result by integrating the matrix elements and subtraction terms directly in d dimensions and then summing the separate contributions. We stress that this is a highly non-trivial cross check, as both calculations are very different conceptually and technically.

Real-virtual contribution
The real-virtual contribution to the differential decay rate reads bbg , (5.14) which contains explicit -poles coming from the one-loop matrix element and furthermore it is divergent in phase space regions where the gluon becomes unresolved. The explicit poles are cancelled by the integral of the single unresolved subtraction term in the double real emission contribution to the full NNLO decay rate, where the real emission differential decay rate is dΓ R 3 is given by eq. (4.1), while the insertion operator I (0) 1 (p 1 , p 2 , p 3 ; ) is given by eq. (4.11). As there are only three partons in the final state, the colour connections that appear in the generic case in eq. (4.11) factorize completely, Thus, Using the expressions in appendix B, it is straightforward to check that which matches the kinematic singularity structure of dΓ RV 3 . The general definition of the real-virtual counterterm is [34], The basic structure of this subtraction in terms of unresolved limits is the same as the tree level single unresolved counterterm in eq. (4.3). However, in accordance with the form of infrared factorization of one-loop QCD matrix elements [81][82][83][84], in eq. (5.21) we have terms with tree level collinear or soft functions multiplying (in colour or spin space) one-loop matrix elements (those with the (0, 1) superscript), as well as terms with one-loop collinear or soft functions multiplying tree level matrix elements (denoted with the (1, 0) superscript). The precise definitions of the functions appearing in eq. (5.21) are given in ref. [34]. Then we consider the counterterm, . In general, the counterterm is given by [34], The organization of this subtraction in terms of unresolved limits is again identical to the tree level single unresolved counterterm in eq. (4.3). However, for each limit, we have two types of terms, labeled by the different superscripts. The reason is as follows. This counterterm is built from the infrared factorization formulae for the product of a QCD squared matrix element times the I 1 insertion operator of eq. (4.11). It turns out that these factorization formulae are sums of two pieces. Both of these involve the product of a tree level collinear or soft function times a tree level matrix element, but one piece is further multiplied by the I (0) 1 insertion operator appropriate to the reduced matrix element, while the other is multiplied with a well-defined remainder function R [34]. Hence the superscripts on the various terms in eq. (5.23).
It can be shown that the combination is both free of -poles and integrable in all kinematically singular limits [34]. Thus, the regularized real-virtual contribution to the decay rate

Double virtual contribution
The double virtual contribution to the differential decay rate reads which contains explicit -poles coming from the two-loop matrix element and the square of the one-loop matrix element: The finite part of dΓ VV 2 is also known exactly [85] which we recall in appendix A (see eqs. (A.3) and (A.4)). In order to regulate these poles, we add the integrals of the counterterms which have been subtracted in sections 5.1 and 5.2. The KLN theorem then ensures that, provided the physical observable we are to compute is infrared-safe and our subtraction scheme is internally consistent, the ensuing result will be free of infrared divergences. It is our task in this section to verify that this is indeed the case.
Let us begin with the integral of the double unresolved subtraction term, eq. (5.2), which can be written as, where the insertion operator has five contributions according to the possible colour structures, The kinematic functions in eq. (5.30) have been defined and computed as expansions in in refs. [40,41]. Again, there is no one-to-one correspondence between the unintegrated double unresolved subtraction terms in eq. (5.3) and the kinematic functions that appear in eq. (5.30). The latter are obtained from the former after integration over unresolved momenta and summation over unobserved colours and flavours. This remark applies to the rest of the insertion operators to be discussed below.
For H → bb, the colour connections that appear in eq. (5.30) are simply given by eq. (4.12), and the kinematic variables simplify as in eq. (4.13). Furthermore, when evaluating eq. (5.30) the coincidence of certain summation indices is allowed. In particular, i in the second line need not be distinct from j and l, while in the last line we only require that i and k as well as j and l are different, with no further restrictions, as shown in the formula. As a result, some indices of kinematic functions coincide once we explicitly write out eq. (5.30). Specifically, since in our case there are only two hard partons in the final state, only CS are absent from the sum, as those require at least three hard partons if all indices are different. In such cases we also simplify the list of arguments of the functions so that we do not display arguments that are the same or identically zero. For instance, in CS whose -expansion is The coefficients of the poles are all given in terms of rational numbers and known transcendental constants.
Next, we consider the integral of the iterated single unresolved subtraction term, eq. (5.6), which can be written as, The kinematic functions in eq. (5.34) have been defined and computed as expansions in in ref. [39]. The discussion below eq. (5.30) applies to eq. (5.34) as well, hence, using eqs. (4.12) and (4.13), we obtain the I whose -expansion is  has the following structure in colour and flavour space, The bare kinematic functions in eq. (5.39) have been defined and computed as expansions in in ref. [35]. Using eqs. (4.12) and (4.13), the unrenormalized I   whose -expansion is

Inclusive and differential results
In this section, we show that using the fully differential two-, three-and four-parton contributions of eqs. in agreement with the known analytic prediction [78][79][80]. In figure 1, we compute the inclusive decay rate at µ = m H /2 and µ = 2m H and compare it to the known analytic result for the scale dependence, finding excellent agreement.
To illustrate the impact of NNLO QCD corrections on differential distributions, we apply the Durham jet algorithm [86] with resolution parameter y cut = 0.05 to cluster final state partons and order the resulting jets in energy. In the top panel of figure 2 we show the energy distribution of the leading jet in the rest frame of the decaying Higgs boson for two-jet events. In ref. [47] the same distribution was computed for jets clustered according to the JADE algorithm with y cut = 0.1. We have repeated that calculation and found excellent agreement with the published results. However, for two-parton kinematics the energy of the leading jet is just E max = m H /2, so at leading order the leading jet energy distribution is a delta function. Furthermore, double unresolved subtractions for four parton matrix elements, as well as single unresolved subtractions for three parton matrix elements also contribute to this distribution only at E max = m H /2. Then, to show the subtraction method at work on an observable that has a non-trivial distribution already at leading order, we consider the absolute value of the pseudorapidity of the leading jet, |η 1 |, with respect to an arbitrary axis. The effect of higher order corrections on this distribution is shown on the bottom panel of figure 2. In this last illustrative example we note that going from the leading order to NNLO, the uncertainty bands shrink, and that the NNLO band falls within the NLO band, thereby showing the good convergence of the perturbative series.
We have shown that our subtractions render both the double real and real-virtual contributions to the NNLO correction integrable in four dimensions. We have also presented the integrated forms of our subtraction terms with pole coefficients evaluated analytically, while the finite parts were given numerically. We confirmed that the sum of the double virtual contribution and the integrated subtractions is free of infrared singularities as required by the KLN theorem. We have implemented our computation in a parton level Monte Carlo program and presented illustrative examples of differential distributions at NNLO.
The successful application of our subtraction scheme reported here opens the way to the computation of other, more involved processes and is also encouraging to further developments of the scheme to deal with initial state radiation. These directions of development are under way and will be the subject of further publications.

A Matrix elements
We present the matrix elements in the form used in our parton level Monte Carlo program. In particular, in our scheme we need the the four-parton tree level and the three-parton one-loop matrix elements only up to finite terms in . Higher order terms must of course be included when integrating the matrix elements and subtraction terms separately in d dimensions. When needed for our cross checks, we take these higher order terms directly from ref. [47].

A.1 Two partons
For H → bb at tree level we have We computed the one-loop correction and obtained We used the formula at two loops as given in ref. [47]: We checked that the poles of this expression satisfy the general formula given in ref. [87], while the finite part agrees with that in ref. [85]. The square of the one-loop matrix element is

A.3 Four partons
In our computation we need the H → four partons squared matrix elements at tree level in d = 4 dimensions. We checked our formulae, presented below, with GoSam [88,89].
For H → bbqq we have For H → bbbb we find