The fully differential decay rate of a Higgs boson to bottom-quarks at NNLO in QCD

The decay of a light Higgs boson to bottom quarks is dominant and can be exploited for the discovery of the Higgs particle and the measurement of its properties at the LHC and future collider experiments. We perform a first computation of the fully differential decay at next-next-to-leading order in perturbative QCD. We employ a novel method of non-linear mappings for the treatment of singularities in the radiative processes which contribute to the decay width. This constitutes the first physical application of the method.


Introduction
Direct searches and electroweak precision tests suggest that a Standard Model Higgs boson is light, and that it should decay predominantly into a bottom quark (bb) pair. Inclusive searches at the LHC of a Higgs resonance in the bottom-pair invariant mass distribution are unfortunately hampered by a very large QCD background. Direct searches are focusing on other decay channels, such as H → γγ, with a much smaller branching ratio. In view of this, perturbative corrections to the H → bb decay are most interesting for the inclusive rate, due to its contribution to the total decay width and the branching ratios of other decays. A remarkable effort has been made for a precise calculation of Γ H→bb and now QCD corrections are known to O(α 4 s ) [1][2][3][4][5][6][7][8][9]. Recently, confidence has grown that the H → bb decay can be observed at the LHC directly in events where the Higgs boson is produced in association with a massive vector boson (W, Z) [10] or a tt pair [11]. Backgrounds from tt,V +jj and multi-jet production are challenging. However, the excellent b−jet tagging of the ATLAS and CMS detectors as well as sophisticated selections of jet events described in [10,11], render these channels hopeful, as explicitly demonstrated by the ATLAS collaboration in a full detector simulation analysis for the ZH case [12]. For this channel, further suppression of the background can be obtained with the sub-jet algorithms of [13].
These search strategies rely on a selection of phase space corners which are rich in potential Higgs events. An accurate modeling of QCD radiation is necessary in order to assess the efficiency of these selections. This motivates the computation of the fully differential H → bb decay rate at next-to-next-to-leading-order (NNLO) in perturbative QCD.
In Ref. [74] we introduced a new method to extract the singularities of multi-dimensional loop and phase space integrals which emerge at NNLO. The method employs non-linear transformations for the factorization of overlapping singularities which in turn allow for an efficient numerical evaluation. In this article, we demonstrate the method in a realistic physical example. We compute for the first time the fully differential decay rate for H → bb at NNLO.
The paper is organized as follows. In Section 2 we fix our notation and set up the calculation. In Section 3 we present the matrix elements for the processes which contribute to the decay width through NNLO. In Section 4 we describe the parameterization of the phase space for the final state partons. In Section 5 and Section 6 we explain how our method leads to a simple factorizable form for the singularities in the integrals which emerge at NNLO. In Section 7 we present the inclusive width and examples of differential observables with our numerical code.

Notation
We compute the partial width Γ H→bbX [J ] through NNLO in perturbative QCD for any infrared safe observable J , such as an appropriate jet-algorithm, in the decay of a Higgs boson to a pair of bottom-quarks, We consider that the Higgs boson couples to a bottom-quark pair with an unrenormalized Yukawa coupling y B b . Due to the largely different mass of the Higgs boson, m H , and the small bottom-quark, m b , we treat the latter as independent of the Yukawa coupling and set it to zero, m b = 0, in matrix elements and phase space integrals. In addition, we neglect top-quarks both in virtual corrections and in the renormalization procedure and consider n f = 5 light quark flavors. For the renormalization of ultraviolet (UV) divergences we work in the MS-scheme where the bare couplings, y B b , and g B s = 4πα B s , are related to their renormalized counterparts through the following relations: The renormalized strong coupling and Higgs Yukawa coupling depend implicitly on the renormalization scale µ, The SU(N) Casimirs are given by 5) where N = 3 is the number of quark colors. Through NNLO, the partial decay width receives contributions from the following partonic processes: where in parenthesis we indicate the four-momenta of the corresponding particles. The width is the sum In the above, we integrate the appropriate matrix elements computed in conventional dimensional regularization (CDR) over the phase space in D = 4 − 2ǫ dimensions of the final-state partons and the function J defining an infrared safe observable. The phase space measure reads Lorentz invariants are defined as Regarding the analytic continuations of complex functions presented in this paper the positive "epsilon prescription" must be used, since all occurring mass scales are external,i.e.
where ε is a small positive parameter.

Matrix elements
In this section we present the square of the scattering amplitudes needed for the computation of the H → bb decay width through O(α s (µ) 2 ). These can be easily derived with modern methods for computing Feynman diagrams and have been ingredients of many prior calculations. For example, the same matrix elements enter the calculation of the NNLO total cross-section bb → H [75]. We have made an independent computation and present them here for completeness. For the generation of the matrix elements we used QGRAF [76]. For further symbolic manipulations, such as color and Dirac algebra, we used FORM [77] and MAPLE [78]. We reduced one and two-loop amplitudes to master integrals using the method of integration by parts [79,80] and the Laporta algorithm [81] with AIR [82].

Decay to two partons
For the process we require up to 2-loop corrections: where (let ℜ(z) shall denote the real part of z) The O(α s (µ) 2 ) contribution can be split into two pieces corresponding to the one-loop squared amplitude and the two-loop amplitude interfered with the tree amplitude where

Decay to three partons
For the process we require up to 1-loop corrections: The tree contribution can be expressed as is the quark gluon splitting kernel. The one-loop correction can be expressed as There are only two different master integrals which appear here. The bubble as well as the following box integral which can be expressed to all orders in ǫ in terms of hypergeometric functions

Decay to four partons
For the decay to four partons we need The H → bbqq amplitude can be expressed as while for two bb pairs in the final state we get The form factors A

Integration over phase space
The matrix elements of Section 3 are integrated over the phase space of the final-state partons. For convenience, we perform our calculation in the rest frame of the Higgs boson.
Our results can be trivially extended to any frame of reference. The matrix-element for the 1 → 2 process H → bb is independent of any phase space integration variables, and it is multiplied by the phase space volume: The matrix elements for the 1 → 3 processes are integrated over the corresponding phase space with volume: (4. 2) The invariants take the simple form: In the above, we have introduced the shorthand notation We parameterize the phase space for the 1 → 4 processes as follows: where The invariants in this parameterization are and This parameterization is very similar to the ones used in Refs [37,57].
In the following sections we shall describe the method of non-linear mappings that we have employed in order to perform the non-trivial phase space integrations over the matrix elements at NNLO. Specifically, we shall discuss the "real-virtual" contributions which require both an integration over the three-parton phase space and an one-loop amplitude, and the "double-real" contributions which require the integration of tree matrix elements over the four-parton phase space.

Real-Virtual NNLO contribution
In the one loop amplitude for H → bbg we will face integrals like Since the numerator in this case is not well defined for s 23 , s 13 → 0, we can not simply apply a naive "plus-prescription" subtraction, to extract the singularities. To tackle this difficulty we use Euler's integral representation of the hypergeometric function, and apply a non-linear mapping to disentangle an overlapping singularity of the integration variable in the hypergeometric representation and the phase space variables which generate s 13 , s 23 . Defining we apply the transformation, which yields This representation of the hypergeometric function has finite values at the points where we need to make a subtraction, and we can obtain the Laurent expansion in ǫ of the real-virtual contribution easily with a "plus-prescription" subtraction. For non-special values of the kinematic variables the ǫ−expansion of hypergeometric function reads: The mapping of Eq. 5.4 simply re-derives a well known identity and there would be no need of it had we started with the following representation of the box master integral: On the other hand, more complicated master integrals for one-loop amplitudes do not always have representations in terms of hypergeometric functions with known expansions in ǫ and variable limits. However, it is always possible to derive a Feynman representation or a hypergeometric representation (via a Mellin-Barnes representation) for them and attempt a direct numerical evaluation of the combined phase space and virtual multiple integral using non-linear transformations to factorize the singularities. We have compared the two approaches: 1. Evaluate numerically the coefficients of the Laurent expansion for the two-dimensional phase space integral after we express the loop amplitude and its hypergeometric functions in terms of polylogarithms, 2. Evaluate numerically the coefficients of the Laurent expansion for the three-dimensional combined phase space and loop integral.
We find that the three-dimensional integration is as fast as the two-dimensional integration where we require the evaluation of polylogarithmic functions.

Integration of double-real contributions
phase space integrals for four-parton processes develop singularities at multiple soft and collinear kinematic configurations. These singularities are overlapping. In this section, we demonstrate how they can be factorized with non-linear transformations. We begin the discussion with an integral with four singular denominators, where the function J represents a non singular function of the partonic momenta, composed of numerators in the squared matrix elements and the infrared-safe observable J . We disentangle the overlapping singularities usinḡ We note that x → α(x, A, B) ⇒x → α(x, B, A). After this transformation, the expansion in ǫ is straightforward with simple subtractions. For the special case J = 1, we obtain Other integrals may be mapped to I 1 of Eq. 6.1 with a relabeling of the partonic momenta. For example, the integral which has an apparently more complicated singularity structure due to s 14 in the denominator, is transformed to the integral of Eq. 6.1 by exchanging p 2 and p 4 . A second class of integrals has the following denominator structure: Notice, that the integrand is singular not only at the boundaries of the integration variables λ i , but it also develops a "line singularity" inside the integration volume. We eliminate this by casting the integral in the form It is easy to see that Eq. 6.7 is equivalent to I 2 [J], by exchanging the momenta p 1 ↔ p 2 in the term with J(p 2 , p 1 , p 3 , p 4 ). The denominator s 13 s 24 +s 14 s 23 = m 4 H λ 1λ1λ2 2λ 3 λ 4λ4 + λ 2λ3 (λ 2 4 +λ 2 4 ) + 2 cos(πλ 5 )(1 − 2λ 4 ) λ 2 λ 3λ3 λ 4λ4 (6.8) is only singular at the boundaries of the integration region. Applying the following sequence of mappings leads to a factorized form. Numerically we then obtain . (6.11) The term in parenthesis in the last denominator is given by s 13 s 23 +s 14 s 24 = m 4 H λ 1λ1λ2 2λ 2λ3 λ 4λ4 + λ 3 (λ 2 4 +λ 2 4 ) + 2 cos(πλ 5 )(2λ 4 − 1) λ 2 λ 3λ3 λ 4λ4 (6.12) We let dI 3 = λ 2 dI 3 +λ 2 dI 3 and apply λ 2 → α(λ 2 , λ 3 , 1) λ 3,4 → α(λ 3,4 ,λ 2 , 1) (6.13) toλ 2 dI 3 and λ 4 → α(λ 4 , λ 3 , 1) (6.14) to λ 2 dI 3 . After these transformations, both integrals have only factorized singularities and can be evaluated numerically. We obtain: Other integrals with four denominators are mapped to I 1,2,3 [J] with simple relabeling of the partonic momenta. Integrals with fewer denominators are simpler and the non-linear mappings for them are obvious upon inspection, following the examples of Ref. [74]. One case which requires special attention arises when a denominator is raised to the second power in the squares of diagrams where a gluon splits into a qq or gg pair as has already been discussed in the literature (see, for example, [68,83] is the one which emerges in the physical matrix elements. After a non-linear mapping to factorize an overlapping singularity, the integrand is also free of quadratic singularities in the integration variables λ i .

Numerical results
In this section we present our numerical results for the decay rate Γ H→bbX [J ] for selected infrared safe observables J . We select the value of the renormalization scale µ = m H . First, we compute the inclusive decay width. Our numerical result is, This is in agreement with the known analytic expression [1] Γ We have checked that our results for the jet-rates add up to the inclusive rate for general y cut values. In Fig. 1, we present the leading jet energy, E max , in the rest frame of the Higgs boson, for events with two jets (y cut = 0.1). This is a new result which can only be obtained with a fully differential NNLO calculation and cannot be inferred from the inclusive decay width and NLO calculations. At leading order, E max = m H 2 . At higher orders, there is a range of jet energies which are kinematically allowed. We choose a value of α s (m Z ) = 0.118 at the Z boson mass and evolve consistently through LO, NLO and NNLO up to the Higgs boson mass, which we assume to be m H = 120GeV.
The above numerical results demonstrate the applicability of our method to physical processes. A number of phenomenological studies which are relevant to the searches of the Higgs boson can be made. It is easy to interface our numerical code with a Monte-Carlo code for the production of a Higgs boson. We we will present complete phenomenological studies pursuing this direction in the future.

Conclusions
In this paper, we present a first physical application of a new method for the factorization of overlapping singularities in phase space and loop integrations. We compute the fully differential decay width of a Higgs boson to a bottom-quark pair. We produce the required tree, one-loop and two-loop amplitudes with standard Feynman diagrammatic methods.
Our article focuses on the phase space integrations which emerge at NNLO. We apply non-linear mappings to factorize all overlapping singularities in all real-virtual and doublereal integrations. Consequently, we perform the expansion of all integrals in the dimension regulator ǫ with simple subtractions. The formalism allows for the computation of the decay rate for arbitrary physical observables.
We verify that we can reproduce the known results for the NNLO inclusive decay width and compute the differential two,three and four jet rates with the JADE algorithm. We also present the leading jet energy distributions, a new result that cannot be inferred from previous calculations. In the future, we will interface our NNLO Monte-Carlo to Monte-Carlo programs for the production of a Higgs boson.
We believe that our method for NNLO computations is powerful and suitable for complicated physical processes. We are looking forward to further apply our method on other interesting processes at hadron colliders.

Acknowledgements
We thank Andrea Banfi and Zoltan Kunszt for useful discussions. This research is supported by the ERC Starting Grant for the project "IterQCD" and the Swiss National Foundation under contract SNF 200020-126632.

A. Double-real matrix elements
The function B(p 1 , p 2 , p 3 , p 4 ) occurring in the H → bbbb squared amplitude can be expressed as