NNLO QCD corrections to associated $WH$ production and $H \to b \bar b$ decay

We present a computation of the next-to-next-to-leading order (NNLO) QCD corrections to the production of a Higgs boson in association with a W boson at the LHC and the subsequent decay of the Higgs boson into a b-bbar pair, treating the b-quarks as massless. We consider various kinematic distributions and find significant corrections to observables that resolve the Higgs decay products. We also find that a cut on the transverse momentum of the W boson, important for experimental analyses, may have a significant impact on kinematic distributions and radiative corrections. We show that some of these effects can be adequately described by simulating QCD radiation in Higgs boson decays to b-quarks using parton showers. We also describe contributions to Higgs decay to a b-bbar pair that first appear at NNLO and that were not considered in previous fully-differential computations. The calculation of NNLO QCD corrections to production and decay sub-processes is carried out within the nested soft-collinear subtraction scheme presented by some of us earlier this year. We demonstrate that this subtraction scheme performs very well, allowing a computation of the coefficient of the second order QCD corrections at the level of a few per mill.


I. INTRODUCTION
Production of the Higgs boson in association with the W boson pp → W H plays an important role in Higgs physics explorations at the LHC [1][2][3][4]. For example, it provides direct access to the HW W coupling, which is completely fixed by the gauge symmetry of the Standard Model (SM) but may receive new contributions in its extensions. The W H associated production is known to provide important constraints on such anomalous couplings, see e.g. Ref. [5]. Furthermore, as was pointed out in Ref. [6], by selecting Higgs bosons with relatively high transverse momenta, it is possible to identify and study the decay of a Higgs boson into a bb pair with high efficiency. The associated W H production then becomes sensitive to the value of the bottom quark Yukawa coupling which currently is only constrained to within a factor of two relative to its SM value [4,7].
The importance of associated W H production inspired a large number of computations of higher-order QCD and electroweak (EW) corrections to this process. The next-to-leading order (NLO) QCD and EW corrections to pp → W H were computed in Refs. [8] and [9,10], respectively. NLO QCD and EW fixed order computations were subsequently matched to parton showers in Refs. [11,12]. The inclusive next-to-next-to-leading order (NNLO) QCD corrections to pp → W H can be deduced [13] from the NNLO QCD corrections to the Drell-Yan process pp → W * computed in Refs. [14,15]. Additional NNLO QCD effects that distinguish associated production from the Drell-Yan process originate from diagrams where the Higgs boson is emitted by loops of top quarks; these effects were computed in Ref. [16] in the large top mass approximation. The numerical program VH@NNLO, which allows highprecision computations of the inclusive cross section of associated Higgs boson production, was developed in Ref. [17].
Fully-differential NNLO QCD results for associated W H production were obtained in Refs. [18,19] using slicing techniques. The NNLO calculation of Ref. [18] was matched to a parton shower in Ref. [20]. NLO QCD corrections to H → bb decay were combined with NNLO QCD corrections to the pp → W H production process in Ref. [21], in the limit of a vanishing b-quark mass, and in Ref. [19], retaining the full m b dependence. Recently, the computation of Ref. [21] was extended [22] to include the NNLO QCD corrections to H → bb decay computed earlier in Ref. [23] (see also Ref. [24]), in the limit m b → 0. Very large effects, apparently caused by an improved treatment of radiative corrections in the decay H → bb, were found for some kinematic distributions.
The purpose of this paper is to repeat the computation of Ref. [22]. There are several reasons for doing so. First, it is important to check the appearance of large effects when QCD corrections to decays are included. Also, we note that some peculiar contributions to Higgs decay to a bb pair that appear at NNLO QCD for the first time were not considered in the computations of Refs. [23,24] and we discuss them here.
Second, the type of distributions for which large QCD corrections were found in Ref. [22] are typically pathological at leading order. For example, kinematic requirements can result in certain regions of phase space only being populated at NLO. In these kinematic regions, the NNLO computations provide next-to-leading order corrections so that moderately large effects are not too surprising. In addition, severe cuts on the final state particles imply the appearance of kinematic boundaries that may cause genuine large effects that signal poor convergence of perturbation theory. In general, many of these effects are driven by parton emissions and may be properly described by parton showers. It is then interesting to check to what extent the large radiative corrections found in Ref. [22] can be described by a parton shower applied to H → bb decay.
Finally, we perform the computation using the local subtraction scheme described recently in Ref. [25]. This scheme is an extension of the original sector-improved residue subtraction scheme developed in Refs. [26,27]. As we already mentioned, all previous computations of W H production at NNLO QCD were performed using variants of the slicing method and it is interesting to perform the computation using fully local subtractions.
The remainder of the paper is organized as follows. In Section II we briefly review the computational scheme of Ref. [25] with an eye on its application to the production process pp → W H. In Section III we illustrate the performance of the subtraction scheme by showing numerical results for NNLO QCD corrections to the pp → W H process, treating the H → bb decay at leading order. In Section IV we discuss the generalization of the scheme of Ref. [25] to the decay process H → bb, and point out differences between the production and decay cases. We also present numerical results for the NNLO QCD corrections to the H → bb decay process, to illustrate the performance of the subtraction scheme in this case as well. Finally, we discuss the phenomenology of the process pp → W (lν)H(bb), consistently including NNLO QCD corrections to both production and decay. We present numerical results for cross sections and selected distributions in Section V and compare them with the approximate treatment of QCD corrections to H → bb decay using a parton shower in Section VI. We conclude in Section VII.

COLLINEAR SUBTRACTION FRAMEWORK
The goal of this Section is to review the subtraction scheme for NNLO QCD calculations [25].
We consider the collision of two partons and ask for the fiducial volume cross section defined by an infra-red and collinear-safe observable O. The fiducial cross section is schematically written as where dLips is the Lorentz-invariant phase space and M is the amplitude for the process X. In Eq. (1), final states X of increasing multiplicity have to be included to arrive at a high-order result for σ f . In our case, the leading-order computation includes partonic processes of the type qq → W H followed by the decays H → bb and W → lν. Both the production and the H → bb decay processes are affected by QCD corrections. In this Section we focus on the QCD corrections to W H associated production and consider Higgs decay in the leading-order approximation.
We note that the NNLO QCD corrections to inclusive W H production are known since long ago [13,16,17]. The challenge for an exclusive computation is to extract the soft and collinear divergences from, say, a matrix element squared with two additional final state partons relative to the leading-order matrix element, while avoiding integration over momenta of partons that can get resolved.
At next-to-leading order, an understanding of how to do this in full generality using both slicing and subtraction methods was achieved more than twenty years ago [28][29][30][31]. Unfortunately, the generalization of these methods to NNLO proved to be difficult and required significant effort. This effort started to pay off in the past two to three years, and a large number of fully-differential NNLO QCD results for important LHC processes has been obtained using different computational methods [26,27,[32][33][34][35][36][37][38].
One of these methods, the so-called sector improved residue-subtraction scheme, was developed in Refs. [26,27] (see also [33,39] for related work). Recently, it was shown [25] how to modify the original formulation of the method by exploiting the fact that in QCD soft and collinear singularities are not entangled. This allows one to closely follow the so-called FKS subtraction scheme [30,31], developed for NLO QCD computations, and perform the required soft and collinear subtractions in a nested way [25]. As a consequence, the computational framework becomes very transparent and, as we show below, numerically efficient.
We will illustrate the main idea of Ref. [25] by considering the double-real emission contribution, taking the process q(p 1 )q (p 2 ) → W H + g(p 4 )g(p 5 ) as an example. Final states with lower multiplicity can be treated along the same lines although the details can be slightly different. Schematically, we write the corresponding cross section as where and is the phase-space element for a gluon, supplemented with a θ-function that ensures that the gluon energy is bounded from above. Note that we introduced the energy ordering of gluons in Eq. (2) to remove the 1/2! identical particles factor.
Our goal is to extract singularities from Eq. (2). These singularities can occur in several ways. For example, the so-called double-soft singularity arises if the energies of the two gluons vanish simultaneously. A single-soft singularity appears if E 5 vanishes at fixed E 4 .
Note that due to the energy ordering in Eq. (2) the opposite limit (E 4 → 0 at fixed E 5 ) cannot occur. In addition to these soft singularities, there are also collinear singularities that occur when the gluons are emitted along the direction of the incoming quark, incoming anti-quark or if they are emitted collinear to each other.
We need to extract all these singularities in an unambiguous way. We begin with soft singularities. We write where S S is an operator that extracts the double-soft 1 singularity from F LM . When the operator S S acts on F LM , it removes the four-momenta of the gluons from both the energymomentum conserving δ-function inside dLips and the observable O, and extracts the leading singular behavior from the matrix element squared. The result is well known where Eik (1,2,4,5) is the square of the eikonal factor derived in Ref. [36]. It is also given in Ref. [25] using the same notation as we use in this paper.
We deal with the two terms on the right-hand side of Eq. (5) in different ways. In the first term the hard matrix element decouples thanks to Eq. (6) and only the eikonal factor needs to be integrated over the two-gluon phase-space. This integral was performed numerically in Ref. [25]. The second term in Eq. (5) has its double-soft divergences regularized. However, both the E 5 → 0 divergence at fixed E 4 as well as the collinear divergences are still present there. To take care of them, we repeat the procedure and subtract the E 5 → 0 singularities at fixed E 4 . We call the corresponding operator S 5 and write The operator S 5 acting on F LM (1, 2, 4, 5) removes the gluon g 5 from the phase space and the observable and extracts the leading singularity We use the notation ρ ij = 1 − cos θ ij in Eq. (8), where θ ij is the relative angle between partons i and j. Among the two terms on the right hand side in Eq. (7), the first has only collinear divergences and the second has a simplified (i.e. independent of g 5 ) matrix element. Therefore, the integration over the energy and emission angles of the gluon g 5 can be performed in this term. The remaining matrix element for the process qq → W H + g 4 can then be treated similarly to a normal NLO computation.
The procedure continues with collinear subtractions that are again applied to the terms on the right hand side in Eq. (7) on top of the soft subtractions shown there. However, an additional step, similar to the energy ordering in Eq. (2), is required. Indeed, we need to 1 Here we define the double-soft limit as further split the phase space into sectors such that in each of them only a particular type of collinear singularity can occur.
There are two major ingredients to this phase space splitting. First, we partition the phase space into two double-collinear partitions and two triple-collinear partitions. In the two double-collinear partitions, the gluons can only have singularities if p 4 || p 1 , p 5 || p 2 , or if p 4 || p 2 , p 5 || p 1 , respectively. In the two triple-collinear partitions, singularities appear if p 1 || p 4 || p 5 or if p 2 || p 4 || p 5 , respectively. We note that in the two latter cases the singularities can also appear if p 4 || p 5 .
The contributions of the double-collinear partitions can be computed right away since all singular limits are uniquely established. The situation is more complex for the triple-collinear partitioning where this is not the case. Indeed, in triple-collinear configurations we need to consider the two cases of the gluons being either close or well-separated in rapidity. To this end, we further partition the phase space into four sectors. Taking as an example the p 1 || p 4 || p 5 partitioning, we introduce four sectors according to the following formula Note that this splitting is largely arbitrary. The important point is that in each of the four sectors only a well-defined type of singular collinear limit can occur; by choosing an appropriate parametrization, these singularities can be resolved and isolated. The nested subtraction of these collinear limits can then be performed, similar to what we discussed in connection with the soft limits. A convenient phase-space parametrization for each of the four sectors can be found in Ref. [26].
A detailed discussion of this approach can be found in Ref. [25] which an interested reader should consult. Below, we list a few aspects of the current computation that go beyond that reference.
• We extend the computation of Ref. [25] by including the qg → W H + q g partonic channel. The difference with the quark-anti-quark annihilation channel is that the quark-gluon channel appears for the first time at NLO and, therefore, to obtain the result relevant for the NNLO computation, we only need to include one-loop corrections to this channel and consider one additional gluon in the final state. There are no conceptual differences with the computations described in Ref. [25] and, similarly to that reference, compact formulas are obtained for the NNLO contribution of the quarkgluon channel. In our implementation, we used a slightly different parametrization of the phase space compared to Ref. [25] making use of the fact that there are no singlesoft singularities related to quark emission.
• We compute all the channels with an additional quark-anti-quark pair in the final state qq → W H + q 1q2 . If the quark-anti-quark pair comes from gluon splitting, the corresponding process has a double-soft singularity that is different from the one described above; the integral of the respective eikonal factor has to be computed anew.
• We compute the contribution of the gg → W H + qq channel. This channel has simple collinear divergences and their extraction is straightforward.
• We include the NNLO contributions to the associated production pp → W H where the Higgs boson is emitted from a loop of virtual top quarks. These include two-loop corrections to the qq → W H partonic process, as well as one-loop corrections to the qq → W H + g process. These finite contributions were first computed in Ref. [16], where they were referred to as V I and R I , respectively. We take the amplitude for V I in an expansion in 1/m top and the amplitude for R I with the exact mass dependence from Refs. [16,19].  Fermi constant is G F = 1.16639 × 10 −5 GeV −2 , and we take sin 2 θ W = 0.2226459 as the sine squared of the weak mixing angle. We also approximate the CKM matrix by an identity matrix. For the decay of the Higgs boson, we take the b-quark Yukawa coupling to be We report our results for W + and W − production in Table I. NLO QCD corrections increase the leading order cross section by about 15%; the NNLO QCD corrections increases the NLO cross sections by an additional 2%. We note that the scale dependence of the NNLO cross sections is below a percent. Therefore, it is both completely negligible and unlikely to be a reliable estimate of the actual theory uncertainty. This issue has been discussed at length in Ref. [41], and we do not comment on it here. Ratios of W + and W − cross sections stay close to 1.57 − 1.58, independent of both the order of perturbation theory and the choice of (5)   the factorization and the renormalization scales. We have cross-checked all these numbers against VH@NNLO [17] and found perfect agreement.
As the next step, we study the NNLO corrections in more detail, focusing on the case of  (7), 0.2193 (7), 0.2464(7)} fb, for the renormalization and factorizations scales Note that the numerical integration error on the NNLO coefficients is just a few per mill.
Also in this case, full agreement with Ref. [17] was found. The fact that our computational method is capable of delivering results at this level of numerical precision for the NNLO QCD coefficients has already been noticed in the calculation reported in Ref. [25]. However, since the calculation of Ref. [25] was performed for a simplified case, it is gratifying to see that this feature persists in a more complex situation where all the different partonic channels are included in the calculation, and significant numerical cancellations between their contributions occur.
NNLO QCD corrections to kinematic distributions can also be computed with a high degree of numerical stability. In Fig. 1 we display the Higgs boson rapidity and transverse momentum distributions in consecutive orders of QCD perturbation theory, for the case of W − H production. In the lower panels of Fig. 1 we also display ratios of NLO to LO and NNLO to 2 Results for W + H production show a similar qualitative behavior.

IV. HIGGS DECAY TO A PAIR OF BOTTOM QUARKS
In this Section, we discuss a fully exclusive computation of NNLO QCD corrections to Higgs boson decay to a bb pair. Such computations were performed in Refs. [23,24]. Unfortunately, both of these references did not consider an interesting subtlety related to this decay that we will explain first.
We consider the Standard Model Lagrangian, integrate out the top quark and neglect the interaction of the Higgs boson with quarks of the first two generations. Interactions of the Higgs boson with hadronic constituents are then described by an effective Lagrangian The two terms in Eq. (11) refer to interactions of Higgs bosons with gluons and b-quarks, respectively. The first term originates from the Htt interaction and, therefore, is proportional to the Higgs-top Yukawa coupling; the second term is proportional to the Higgs-bottom Yukawa coupling.
The two constants C 1,2 in Eq. (11) are the Wilson coefficients of the corresponding operators.
Their perturbative expansions in the strong coupling constant -to the order relevant to us -read (see e.g. [42]) The computation of NNLO QCD corrections to Higgs boson decay to two b-quarks reported in Refs. [23,24] was performed under a tacit assumption C 1 = 0 and C 2 = 1. As we explain below, C 1 = 0 leads to additional contributions to Higgs decay to bb starting at NNLO.
In the limit of a small b-quark mass, these contributions scale like /v 2 , so they are parametrically indistinguishable from terms proportional to y 2 b coming from C 2 alone. As a consequence, they should be included in an NNLO computation. However, before discussing this point, we repeat the computation of the decay H → bb reported in Refs. [23,24] by setting C 1 = 0, C 2 = 1.
A. Higgs decay to a bb pair: contribution proportional to bottom Yukawa coupling squared In this subsection, we compute NNLO QCD corrections to the decay H → bb in the approximation C 1 = 0. We treat b-quarks as massless but with a non-vanishing Yukawa coupling to the Higgs boson. The generalization of the computational method described in Ref. [25] to this case is straightforward. Since the collinear renormalization of parton distribution functions is obviously not needed in this case, the computation is simpler and the final formulas for the NNLO QCD corrections are more compact. There are, however, a few subtleties, which we point out in this Section.
First, as we already mentioned, we work in the approximation of massless b-quarks. This means that the only place where the b-quark mass appears is in the Yukawa coupling. We renormalize the Yukawa coupling in the MS-scheme at the scale µ = m H . It is well-known from the computation of the inclusive rate that this choice of the renormalization scale reduces the magnitude of QCD radiative corrections that are very large otherwise [43].
Second, integrals of the double-soft eikonal factors are identical to the production case and can be re-used in the H → bb computation. Other numerical components of the computation, i.e. integrals of the triple-collinear splitting functions, are different from the production case but they actually become simpler. 3 Third, it turns out that the calculation of the double-collinear contributions is non-trivial for the decay kinematics. This is in stark contrast to the computation of the NNLO QCD corrections to the production case where the double-collinear contribution is among the simplest. The reason for this difference is as follows. The double-collinear contributions refer to sectors where collinear singularities appear if, say, the gluon g 4 is emitted collinear to the b-quark and the gluon g 5 is emitted collinear to theb-anti-quark. To extract collinear divergences in this case, it is convenient to choose cosines of the relative angles between p b and p 4 and between pb and p 5 as independent kinematic variables. For the decay case, we work in the rest frame of the Higgs boson. Hence, in contrast to the production case, the directions of p b and pb are not fixed. It then appears to be non-trivial to use the two angles as independent variables and to have the phase space properly simplify in soft and collinear limits, while also satisfying the constraint p H = p 4 + p 5 + p b + pb. Nevertheless, this can be done and we will present the corresponding formulas in a separate publication. Here, we only note that this complexity is a particular feature of the process at hand. Since the Born process H → bb involves too few particles, the momentum conservation constraint makes it difficult to find a parametrization in terms of the two angles discussed above. For more complicated decay processes, for example for Z decays to three jets, this issue is not present.
The last point concerns the contribution of the bbbb final state to the decay rate of the Higgs boson. This final state is different from everything that we considered so far because we cannot say a priori which of the two bb pairs comes from the Higgs vertex and which from the g * → bb splitting. Without this information, we cannot separate the phase space into a hard part and a radiation part, which is central for the method of Ref. [25]. To get around this problem, we use the symmetry of the process H → bbbb with respect to the permutations of the two b-quarks and the twob-anti-quarks and split the matrix element into a part that is equivalent to the singlet component H → bb + qq, q = b, and an identical quark interference contribution. In each of the interference contributions, there is either a quark line or an anti-quark line that always originates from the Higgs decay vertex. We assign this line to belong to the hard phase space. The remaining lines can originate either from the Higgs decay vertex or from the g * → bb splitting. Which line belongs to the hard phase space and which one to the radiative phase space is a matter of choice at this point.
The interference terms only contain a purely triple-collinear singularity. It corresponds to the interference term in the non-singlet triple-collinear splitting function [44] and can be easily extracted and integrated numerically.
We continue by presenting some numerical results of the calculation. Again, our goal in this Section is not to discuss phenomenology of the Higgs boson decay to a bb pair but to show that our method is capable of producing high-precision results.
The numerical computation yields the following result for the decay rate of the Higgs boson to a bb pair where . The value of the Yukawa coupling constant has already been discussed in the previous Section. The renormalization scale for the strong coupling constant is set to the mass of the Higgs boson.
It is instructive to compare Eq. (13) with the results of an analytic computation [45]. The analytically-known two-loop coefficient evaluates to 116.59... which is in better than per mill agreement with the result of the numerical computation shown in Eq. (13).
It is also interesting to compute jet rates in H → bb decay since such, more exclusive, calculations provide a stronger test of the numerical stability of the method. Similar to Ref. [24], we use the JADE clustering algorithm with y cut = 10 −2 to define jets. 4 We obtain Figure 3: Illustrative interference diagrams that contribute to the H → bb decay rate for C 1 = 0.
See text for details.
The sum of the exclusive jet rates in Eq. (14) gives the total decay rate; computing this sum, we obtain We mentioned above that a non-vanishing Wilson coefficient C 1 gives rise to additional contributions to H → bb decays starting at NNLO in QCD, which were not considered in previous fully-differential calculations [23,24]. We describe these contributions in more However, since the calculation in the previous subsection was performed with massless bquarks, we would like to compute the diagrams shown in Fig. 3 in the same approximation.
Unfortunately, doing so leads to problems. Indeed, if we factor out one power of m b caused by the helicity flip, the reduced matrix element has peculiar soft and collinear limits in the m b = 0 approximation, that are typically not present in QCD amplitudes at leading power.
For example, it develops a logarithmic singularity when a single b-quark becomes soft.
The validity of the massless approximation assumes that the logarithmic dependence on the b-quark mass cancels out in infra-red safe quantities. It is easy to see, however, that this cancellation does not take place for the interference contributions, and it is not possible to give a proper inclusive definition of this process in the massless approximation. Indeed, the logarithmic mass dependence cancels between the diagrams in Fig. 3 and similar diagrams with a b-mediated Higgs decay into gluons. One could try to circumvent this problem by regulating the collinear singularity related to g * → bb with a flavored jet algorithm, e.g. the one in Ref. [46]. This would trade the logarithmic dependence on the b-quark mass for a logarithmic sensitivity to a jet radius R. However, even this does not solve the problem completely as the single-soft quark singularity is not regulated by the jet algorithm of Ref. [46].

It is clear that a proper description of the interference contributions requires a computation
with fully massive b-quarks. In its absence, we estimate the order of magnitude of these effects by simply imposing restrictions on the phase space of the b-quarks and gluons that reproduce the leading logarithmic terms. We find that these contributions may change the NNLO corrections to the inclusive Higgs decay rate shown in Eq. (13) by up to O(30%).
Using the same setup in the fiducial region that will be discussed in the next Section, we find that this interference contribution is somewhat reduced. Since their impact on the decay rate appears to be limited, we will omit these terms from the phenomenological analysis in the next Section but we stress that it is important to understand them better. As we explained, this will require a fully-differential computation of the Higgs decay to massive bottom pairs at NNLO. We leave this for future investigations.

V. THE PHYSICAL PROCESS
We are now in position to discuss the physical process pp → W (lν)H(bb), including QCD corrections to both production and decay. Given the results of the preceding Sections, it is straightforward to do so. The only subtlety is how to treat the Higgs boson decay width that appears in the cross section in the narrow width approximation. We write We note that in the approximation of massless b-quarks, the Higgs boson decay rate to a bb pair and therefore the Higgs branching ratio to a bb-pair subtly depends on the definition of a b-quark. However, the effect on the total decay rate is relatively small, as discussed in the previous Section, and we set C 1 = 0 for the phenomenological studies in this paper. We use Br(H → bb) = 0.5824 [41] as a fixed quantity, not subject to an α s expansion.
To define an expansion of Eq. (16) in α s , we follow Ref. [22], write the production cross section and the decay width to bb as an expansion in α s and introduce Note that dγ (i) = 1, provided that the integration goes over the unrestricted phase space.
Using this notation, we define the physical cross sections computed through different orders in QCD perturbation theory In addition, for comparison with the previous computations of Refs. [19,21], it is convenient to introduce an approximate NNLO cross section that includes NNLO corrections to the production process but only NLO corrections to the decay. It reads We are now in a position to discuss the results of the computation. To define the W (lν)H(bb) final state, we reconstruct b-jets using the infra-red safe flavor-k t jet algorithm [46] 5 with ∆R = 0.5 and require that an event should contain at least one b-jet and oneb-jet with An identified light (non-b) jet is required to have a transverse momentum p ⊥ > 25 GeV as well but no pseudo-rapidity cut is applied in this case. In addition, we impose the following cuts on the pseudorapidity and transverse momentum of the charged lepton For the cross sections in Eqs. (23,24)  The results for the fiducial cross sections Eqs. (23,24) show that NLO QCD effects are larger if a transverse momentum cut is imposed on the W boson. This is expected since the W boson can evade this cut by recoiling against additional radiation which appears at NLO.
The approximate NNLO results, which include NNLO corrections to the production process only, show a similar effect: this cross section is about 3% higher than the NLO cross section without the p ⊥,W cut, but about 6% higher when this cut is imposed. Including Higgs decay at NNLO decreases the approximate cross section by about 9% without the p ⊥,W cut, and 7% in the presence of this cut. Therefore there are cancellations between corrections to the production and decay sub-processes, that make the size of the full NNLO QCD corrections quite sensitive to the value of the p ⊥,W cut.
We now turn to differential distributions. We begin by identifying the bb system comprised of a b-jet and ab-jet whose invariant mass best approximates the mass of the Higgs boson, and consider the invariant mass m bb distribution of this bb system. Since we work in the narrow width approximation, at leading order this distribution is described by a deltafunction δ(m 2 bb − m 2 H ). At next-to-leading order, this situation changes: a gluon emitted in the Higgs boson decay can decrease the invariant mass of the bb system while a gluon emitted in the production process can increase it. Hence, the m bb distribution has tails both above and below m bb = m H that start to appear if the next-to-leading order correction to either production or decay is included in the computation. In Fig. 4 we compare predictions for this observable obtained using full and approximate NNLO computations, defined in Eqs. (19,20), respectively. We study this observable both without (left) and with (right) the cut on the W boson transverse momentum p ⊥,W > 150 GeV. It is seen from Fig. 4 that the application of this cut affects the shape of m bb distribution in a minor way. For example, in both cases, full NNLO results deplete the distribution at m bb > m H and enhance the distribution at m bb < m H relative to approximate NNLO predictions. Since the full NNLO provides a better description of the radiation in the decay, compared to the approximate NNLO, and since radiation in the decay predominantly reduces m bb , this re-shaping is not unexpected.
However, the magnitude of this O(α 2 s ) effect -O(60%) correction at m bb ∼ 80 GeV and O(−15%) at m bb > m H -is somewhat surprising. We note that similarly large corrections have also been observed in Ref. [22].
To understand what causes these large effects, we split the difference between approximate NNLO and full NNLO into two terms -NNLO radiation in the decay (NNLO dec ) and NLO radiation in the production followed by the NLO radiation in the decay (NLO prod ×NLO dec ). We define such that dσ NNLO,approx W H(bb) + δ dec. + δ NLO×NLO = dσ NNLO W H(bb) . We display the two distributions in Fig. 5. As we said already, the radiation in the decay does not populate the m bb region to the right of m H , so that the O(−15%) correction at such values of the bb invariant mass comes exclusively from the NLO prod × NLO dec contribution. On the other hand, for m bb < m H the NNLO corrections to the decay play a dominant role, increasing the distribution by about 40%, as compared to the O(20%) increase from NLO prod × NLO dec .
Next, we consider the transverse momentum of the bb system whose invariant mass provides the best approximation to the Higgs boson mass. The NNLO and approximate NNLO distributions for this observable are compared in Fig. 6; the cut p ⊥,W > 150 GeV is applied to events displayed in the right pane. It follows from Fig. 6 that the cut on the W boson transverse momentum re-shapes the distribution, pushing its maximum to larger values.
Again, this is easily understood by observing that the p ⊥,W cut implies the requirement p ⊥,bb > 150 GeV at LO. In addition, if the cut on the W transverse momentum is applied, both the full and the approximate NNLO calculations develop a Sudakov shoulder below p ⊥,bb = p cut ⊥,W = 150 GeV. We note that this feature is somewhat less prominent in the full NNLO distribution.
To understand the relative impact of different contributions, we again split the full NNLO into two different parts, δ dec. and δ NLO×NLO , and display them separately in Fig. 7. For values of p ⊥,bb larger than p cut ⊥,W , the approximate NNLO is larger than the full NNLO by about O(5% − 10%), independent of whether or not the cut on the W boson transverse momentum is applied, due to the corrections from both NLO prod × NLO dec and the NNLO decay. When the p ⊥,W cut is imposed, the slight increase at low values of p ⊥,bb is the result of a cancellation between the somewhat larger contributions from the NNLO decay and the NLO prod × NLO dec . We also note that the NLO prod × NLO dec contribution smears the Sudakov shoulder.
It is also interesting to study the angular separation ∆R bb = ∆η 2 bb + ∆φ 2 bb of the band b-jets that are used to reconstruct the Higgs boson; the corresponding distributions without (left pane) and with (right pane) the p ⊥,W cut are shown in Fig. 8. The impact of the W boson transverse momentum cut on the angular separation of the jets is dramatic, as the comparison of left and right panes shows. The shift to lower values of ∆R bb is again expected, as imposing the p ⊥,W cut selects boosted Higgs kinematics whose decay products are closer together. Both with and without the p ⊥,W cut, the NLO corrections modify the shape of ∆R bb distributions significantly, while the NNLO corrections have a much smaller impact.
Another distribution that is subject to large modifications if the cut on the vector boson   transverse momentum is applied is the transverse momentum distribution of the hardest b-jet; it is shown in Fig. 9. In this case, large radiative corrections appear below the value of the transverse momentum where the distribution reaches its maximum. If the p ⊥W cut is not applied, large corrections at NLO are followed by moderate corrections at NNLO. On the contrary, if the p ⊥W cut is in place, both the NLO and NNLO corrections are very large and perturbation theory does not appear to converge (see the right pane in Fig. 9). Clearly, the situation is completely different at high values of p b ⊥ where NNLO effects are relatively small and the NNLO/NLO K-factor is flat and close to one.
As the last example, we show in Fig. 10    lepton that originates from the W decay. In this case, the cut on the W boson transverse momentum has a significant impact on the shape of the distribution, but the NLO and NNLO corrections to the two cases are very similar. In particular, the NNLO corrections in both cases are relatively small and do not change the shape of the respective NLO distributions.

VI. COMPARISON OF FIXED ORDER AND PARTON SHOWER PREDICTIONS
The goal of this Section is to compare fixed order QCD predictions for pp → W (lν)H(bb), described in the previous Section, with the results obtained when parton showers are used to account for QCD radiation in H → bb decays, as typically done in many experimental analyses. We use the publicly available HWJ generator [12] implemented in the POWHEG BOX framework [48][49][50] to compute the process pp → W (lν)H +j at NLO QCD accuracy. In order to be as close as possible to the NNLO calculation, and since the HWJ generator allows it, we run it with the improved MiNLO method [51,52]. This allows observables that are inclusive in the production of the color-neutral system, i.e. quantities in which the jet is unresolved, to be computed with NLO QCD accuracy. Thus, the difference between the NNLO fixed order calculation and the NLO parton shower simulation for the process pp → W (lν)H is formally due only to the missing two loop amplitudes in the HWJ generator. The decay of the Higgs boson to a bb pair and an arbitrary number of gluons is instead simulated with a parton shower using PYTHIA-8 [53] with the default tune. Since we want to compare the parton shower results with a fixed order calculation, we do not include any non-perturbative effects in the simulation, i.e. the hadronization and the multi-parton interactions are switched off. In the parton shower simulation we reconstruct jets using the anti-k t algorithm [54], and select b-jets according to Monte Carlo truth, in order to be as close as possible to experimental analyses. Following Ref. [22] and the analysis in the previous Section, we use Moreover, some of these regions, e.g. p ⊥,bb ∼ p cut ⊥,W or hardest p ⊥,b → 0, are close to kinematic boundaries where parton showers are known to accurately describe radiation effects. Other regions and observables, for example the case m bb < m H require a relatively hard gluon emission and it is unclear a priori if parton showers do a good job in describing them.
As in the previous Section, we study the b andb jets whose invariant mass m bb is closest to the Higgs mass. We show a comparison of the NNLO and parton shower predictions for the m bb distribution in Fig. 11, for the transverse momentum distribution of the bb system in Fig. 12, and for the hardest b (orb) jet p ⊥ distribution in Fig. 13. In all of these cases, the   Figure 12: Same as Fig. 11 but for the transverse momentum of the bb system that is used to reconstruct the Higgs boson. See text for further details. distributions are normalized to their inclusive result so that their shapes can be compared.
However, we note that, while the fixed order and parton shower results use the same jet radius, the former makes use of the flavor-k t jet algorithm while the latter uses the standard anti-k t algorithm, and therefore the comparison between the two is not straightforward. We will return to this point at the end of this Section.
For the m bb distribution, we observe that the parton shower does quite a good job in describing the NNLO corrections, although it predicts more events at both low and high values of m bb . Interestingly, the parton shower smears the peak at m bb = m H more significantly in the case where the p ⊥,W cut is not applied. When this cut is imposed, the parton shower predicts fewer events at the peak but the smearing effect is not as dramatic. See text for further details.
Turning to the p ⊥,bb distribution, we observe that the parton shower is able to describe the NNLO distributions quite well. When the p W ⊥ cut is not imposed, the parton shower prediction is in excellent agreement with the fixed order one, except in the very high transverse momentum region. However, there is a difference at low p ⊥,bb if the p W ⊥ cut is applied, with the parton shower predicting more events in this region than the fixed order calculation. As expected, the parton shower also removes the Sudakov shoulder in this distribution that was observed in both the approximate and the full NNLO distributions.
Next, in Fig. 13 we show the p ⊥ distribution of the b-(orb-) jet with largest transverse momentum. Without the cut on p W ⊥ , the NNLO and shower results are similar, although the latter predicts slightly more events at large p ⊥ . On the other hand, if the cut p W ⊥ > 150 GeV is imposed, the fixed order and shower calculations deviate significantly at small p ⊥ . Large shower effects in this region are expected, since as we have shown in Section V, the fixed order predictions are not reliable here.
Given the different jet algorithms used in the fixed order and parton shower calculations, it is interesting to investigate to what extent the details of the jet definition affect these results.
In Figs. 14 and 15, we show the invariant mass m bb and transverse momentum distribution p ⊥,bb , obtained from the parton shower simulation for different choices of the jet algorithm and radius. We compare the flavor-k t jet algorithm [46] with both the k t [55] and anti-k t [54] algorithms. For the invariant mass distribution, Fig. 14 shows that both with and without the p W ⊥ > 150 GeV cut the result is quite insensitive to the recombination algorithm, and it only depends on the choice of the jet radius: smaller values of R lead to more events below the Higgs peak. For the p ⊥,bb on the other hand, Fig. 15 shows that without the p W ⊥ cut all jet algorithms and radii lead to the same result, apart from the high p ⊥,bb tail where the flavor-k t jet algorithm [46] predicts fewer events compared to the k t and anti-k t cases. With the additional p W ⊥ > 150 GeV cut, a qualitative dependence on the jet radius similar to the one seen in the m bb distribution is observed: smaller values of R lead to a softer spectrum.

VII. CONCLUSIONS
In this paper we presented a computation of the NNLO QCD corrections to the associated production of the Higgs boson pp → W H at the LHC. We considered the H → bb decay of the Higgs boson and included radiative corrections to this decay through NNLO in perturbative QCD.
We pointed out an interesting contribution to Higgs decay to bb pairs that was ignored in previous fully-differential NNLO QCD computations to this process. This contribution is infrared-sensitive even after standard jets algorithms are applied and understanding it necessitates the computation of fully-differential NNLO corrections to the H → bb decay with massive bottom quarks. Although we argued that its numerical importance should be small, a more refined analysis is required to properly quantify these effects.
We found a number of kinematic distributions in the pp → W (lν)H(bb) process that receive large perturbative corrections if certain cuts on the final state, and especially a cut on the transverse momentum of the W boson, are applied. These findings are in accord with an earlier discussion given in Ref. [22].
We compared fixed order predictions for the pp → W (lν)H(bb) process with calculations where a parton shower is used to describe QCD radiation in H → bb decay. Parton showers confirm the existence of large effects observed in fixed order computations. Since at the moment fixed order NNLO QCD computations for H → bb are performed for massless bquarks, one has to use specially-tailored jet algorithms to describe flavored jets in fixed order computations [46]. Although we showed in Section VI that the results are largely insensitive to the jet recombination algorithm, it would be interesting to repeat the fixed order studies reported here for the setup used in experimental analyses. This requires the computation of the fully-differential decay H → bb through NNLO QCD keeping the full dependence on the b-quark mass. It would also be interesting to compare fixed order predictions to more advanced parton shower implementations, as described e.g. in Refs. [20,[56][57][58]. We leave these investigations for future work.