Associated production of a Higgs boson at NNLO

In this paper we present a Next-to-Next-to Leading Order (NNLO) calculation of the production of a Higgs boson in association with a massive vector boson. We include the decays of the unstable Higgs and vector bosons, resulting in a fully flexible parton-level Monte Carlo implementation. We also include all $\mathcal{O}(\alpha_s^2)$ contributions that occur in production for these processes: those mediated by the exchange of a single off-shell vector boson in the $s$-channel, and those which arise from the coupling of the Higgs boson to a closed loop of fermions. We study final states of interest for Run II phenomenology, namely $H\rightarrow b\bar{b}$, $\gamma\gamma$ and $WW^*$. The treatment of the $H\rightarrow b\bar{b}$ decay includes QCD corrections at NLO. We use the recently developed $N$-jettiness regularization procedure, and study its viability in the presence of a large final-state phase space by studying $pp\rightarrow V(H\rightarrow WW^*) \rightarrow$ leptons.

1 Introduction Run II of the LHC promises to shed new light on the mysteries behind the breaking of the electroweak symmetry. The standout result from Run I of the LHC was the discovery of a Higgs boson [1,2]. One of the principal physics goals of Run II is to pin down the precise nature of the Higgs boson and in particular how it interacts with the other particles of the Standard Model (SM). In order to do this a range of Higgs production and decay processes must be studied in greater detail than ever before. A significant improvement that is expected in Run II analyses is their ability to study the Higgs differentially for a wider range of processes. One such fascinating process is the production of a Higgs boson in association with a W or Z electroweak vector boson, i.e. pp → V H where V denotes the vector boson. At LHC energies these processes are the third (V = W ) and fourth (V = Z) largest production channels. V H production is somewhat special, in that it proceeds at Leading Order (LO) through an s-channel Feynman diagram. This results in the opportunity to probe the V V H vertex at high momentum transfer while keeping the final state vector and Higgs bosons on-shell, for instance by looking in the region of large m V H . This is an interesting region to study, since contributions arising from physics beyond the Standard Model (BSM) may induce a momentum dependent term in the V V H vertex [3][4][5]. New physics at TeV scales would modify the SM cross-section at the level of a few percent. Accordingly, it is essential that the SM cross-section is known at this level or better.
A second distinguishing feature of the pp → V H process is the ability to study the decay of the Higgs boson to a pair of bottom quarks, H → bb. Such a decay is extremely difficult to measure in inclusive Higgs boson production, given the small rate of gg → H → bb compared to the QCD production of the same final state. It is essential that the decay H → bb is measured experimentally since it provides a direct measurement of the coupling of the Higgs boson to fermions. Moreover, since it dominates the total width, the uncertainty on this branching ratio feeds into other searches, for instance, measurements of the Higgs invisible branching ratio. The presence of the vector boson in the final state allows experimental analyses to have manageable backgrounds, in particular when the Higgs is highly boosted and the two b-quarks reside inside a fat jet [6]. Again it is essential that accurate theoretical predictions are available, with the ability to apply intricate final state phase space selection requirements.
Given its importance, the pp → V H processes have been extensively studied in the theoretical literature. At LO the topology is essentially the same as that of Drell-Yan (DY) production, and this was utilized to obtain the first Next-to-Next-to Leading Order (NNLO) predictions for on-shell bosons in ref. [7]. However at O(α 2 s ) a second type of diagram appears, in which instead of coupling to the vector boson, the Higgs is radiated from a closed loop of heavy fermions. These "y t " pieces 1 were computed for on-shell vector bosons in ref. [8]. A fully differential calculation, including the decays of the bosons, was presented for the DY parts (i.e. neglecting the y t terms) of W H in refs. [9,10] and of ZH in ref. [11]. A subset of the y t diagrams, corresponding to those which are initiated by a pair of gluons gg → HZ, was also included in the calculation of ref. [11]. A primary motivation of this paper is to extend the calculations of refs. [9][10][11] to fully account for the contributions discussed in ref. [8] in a flexible Monte Carlo code. We will also extend the range of Higgs boson decays beyond the two-body ones presented previously. Electroweak corrections were calculated in ref. [12,13] while resummation effects have been studied in refs. [14][15][16]. There has also been significant progress in matching fixed order calculations to parton shower Monte Carlos, allowing for full event simulation. An implementation of the V H process in the POWHEG formalism was presented in [17] and extended to merge with the V H+jet process in ref. [18]. A SHERPA implementation that merges the V H and V H+jet processes was also presented recently in ref. [19].
The historical bottleneck for NNLO computations was in the construction of regularization schemes to handle the InfraRed (IR) singularities. These singularities are ubiquitous in a NNLO calculation, since they occur in the two-loop (double virtual), one-loop × real (real-virtual) and the real-real part of the calculation. The situation is made more complicated by the different dimensionality of the phase space in each part. The double-virtual has the same dimension as the Born, and IR singularities manifest themselves as poles in an expansion (where d = 4 − 2 , with parameterizing excursions from four dimensions). The real-virtual has one additional parton in the final state, and possesses IR singularities which manifest themselves as poles, and when the emitted parton becomes unresolved. Finally the real-real piece corresponds to the emission of two additional partons and its IR singularities correspond to when one, or both partons become unresolved. Constructing a scheme to regulate these divergences has been a ongoing task for many years [20][21][22]. Recently a new regularization scheme, based upon N -jettiness [23] has been proposed [24,25].
Here the idea is similar to that used in q T subtraction [21], and a calculation of the top quark decay at NNLO, based on Soft Collinear Effective Field Theory (SCET) methods [26]. These methods introduce a variable which separates the singly unresolved regions from the doubly unresolved ones. If an all-orders formulation (i.e. a factorization theorem) is known for the doubly unresolved region, then an expansion can be performed to a fixed order in the coupling. The singly unresolved region corresponds to the NLO calculation of the process with an additional parton, which can be evaluated using traditional means. For q T subtraction, applicable to production of colour neutral final states, the separation is obtained via a q T cut. If q T > q cut T then the electroweak (EW) system recoils against a parton, and only single unresolved limits can occur (i.e. the NLO calculation of the EW final state together with one additional parton). For q T < q cut T the all-orders factorization of Collins, Soper and Sterman [27] can be used. In refs. [24,25] N -jettiness [23] was proposed as the separation-cut i.e. τ N > τ cut N defines a NLO calculation. When τ N < τ cut N SCET [28][29][30][31][32][33] provides a factorization theorem [23,33] that can be used to compute the cross section. An advantage of this method is that it can be applied to coloured final states with jets. The recent advances in (a variety of) NNLO regularization schemes has led to a veritable explosion in the number of phenomenological predictions at NNLO for 2 → 2 scattering [24,[34][35][36][37][38][39][40][41][42] The aim of this paper is twofold. Our chief goal is to provide the first NNLO calculation including both the DY and y t contributions, with full flexibility in the boson decays for both W H and ZH processes. Second, we will apply the recently-developed SCET formalism to a detailed phenomenological study, including the process V H → V W W →leptons. Such decays have not previously been included in NNLO codes, but are studied experimentally. Given the large and intricate final state phase space (22 dimensions for the double-real part) this is a particularly good example to test the feasibility of the SCET regularization to provide NNLO predictions for complicated phenomenological applications. Our results are implemented in MCFM [43][44][45], and are available in MCFM 8.0 [46]. This paper proceeds as follows. In section 2 we present an overview of the component pieces needed to complete the calculation of V H at NNLO. Phenomenological results for the LHC Run II are then presented in section 3. We draw our conclusions in section 4. We present a detailed discussion of the helicity amplitudes needed in the computation of the NNLO correction in Appendices A and B.

Calculation
In this section we describe the details of our NNLO calculation for V H production and its implementation into a fully flexible Monte Carlo code. The aim of this section is to provide an overview of the calculation, and its subsequent implementation in MCFM. Technical details regarding the calculation of the amplitudes are presented in Appendices A-B. At NNLO the production cross-section dσ (2) pp→ 1 2 H can be written as the sum of two terms, dσ (2) Here the first term represents the contributions which have the same structure as single vector boson production, the second term represents a new type of contribution that occurs first at O(α 2 s ). These pieces arise from terms in which the Higgs boson couples directly to a heavy quark (predominantly a top-quark). In the following sections we first describe these two contributions in more detail, and then discuss our handling of the decays of the Higgs boson.

Drell-Yan type contributions
At LO and NLO the production cross-section dσ (i) pp→ 1 2 H (where i = 0, 1) has the same structure as the calculation of single vector boson production. At LO only qq initial states contribute, while the NLO corrections consist of virtual (one-loop) corrections to this process, and real-radiation in which the underlying matrix elements contain a qq pair and a gluon. Representative Feynman diagrams for these pieces are illustrated in Fig. 1 where, for simplicity, we have suppressed the decays of the vector and Higgs bosons. At NLO IR singularities are isolated using dimensional regularization and handled using Catani-Seymour dipole subtraction [47].
At NNLO the production cross-section receives contributions from the V H + 0, 1 and 2 parton phase spaces. Representative Feynman diagrams for each of these terms are presented in Figure 2. Utilizing the similarities with the NNLO calculation of the Drell-Yan process [48], cross-sections for inclusive on-shell V H production were presented at this order in ref. [7]. At O(α 2 s ) dσ   IR origin. In order to regularize these we use the recently developed N -jettiness slicing procedure [24,25,37,41]. This procedure uses the N -jettiness variable (τ N ) to divide the NNLO calculation into two pieces based on the value of τ N . Below the τ -cutoff parameter the technology of SCET [23,33,[49][50][51] is used to provide a factorization theorem. Above the τ N -cutoff the calculation reduces to a NLO computation of the (V H + j) process, and can be evaluated using traditional techniques. In MCFM the IR regularization of the NLO V Hj processes is obtained via the Catani-Seymour dipole formalism [47]. Since the SCET formalism below the τ N -cutoff is approximate and subject to power corrections, the value of τ N should be taken as small as possible. A check of the implementation is thus obtained by checking the cancellation of the logarithmic pieces above and below the cut. For our process, which does not contain any final state jets in the Born phase space, the τ N -cutoff procedure is similar to the q T subtraction technique [21] used in previous calculations [9]. A detailed study of the N -jettiness regularization for colour singlet final states and their implementation in MCFM is presented in ref. [46], to which we refer the interested reader for more details. MCFM also contains implementations of one-jet production in association with a Higgs [37], W [24] or Z [41] boson, and diphoton production [52]. We stress that in MCFM we only use N -jettiness slicing to calculate the coefficient of the O(α 2 s ) term in the perturbative expansion.
In order to implement the DY pieces in MCFM we need the two-loop virtual amplitude [48] interpreted in terms of the hard function of SCET [53,54], and the NLO implementation of V H + j. The results for the two-loop virtual amplitude are readily available in the literature [48]. We have calculated the NLO corrections to the V H + j process and implemented them in MCFM. Details of the relevant calculational ingredients are presented in Appendix A.

Top Yukawa contributions
A new type of process opens up at O(α 2 s ) and corresponds to diagrams in which the Higgs boson does not couple directly to the vector boson, but instead couples to a massive quark. Since the top quark has by far the largest Yukawa coupling, these contributions are dominated by the top-quark loops. These y t diagrams further sub-divide into two categories. Diagrams of the first kind, representatives of which are presented in Fig. 3, contain a closed loop of heavy quarks which does not radiate the vector boson. The second kind, illustrated   in Fig. 4, contains diagrams that include a closed loop of fermions which radiates both the Higgs and the vector boson. Charge conservation mandates that the latter examples are forbidden if the radiated boson is a W . Therefore the first topologies ( Fig. 3) occur for both W H and ZH production and the latter topologies occur only in the ZH case. Both sets of topologies can have two-loop qq topologies, which interfere with the LO amplitude, and one-loop qqg topologies, which interfere with the qqgV H tree amplitude. These pieces have been computed for on-shell final state particles in ref. [8] and we follow the nomenclature introduced in that paper. We refer to the two-loop diagrams by the label V and the one-loop diagrams by R. The sub-topologies of these sets are further distinguished by I (for diagrams that occur for both W H and ZH production) and II (ZH only). In ref. [8] these pieces were computed and found to contribute around 1-3% of the total NNLO crosssection. Whilst this may appear to be a small contribution that can safely be neglected, the total NNLO correction from the DY-type diagrams discussed previously is itself of the same order. Therefore in order to obtain a reliable prediction at O(α 2 s ) it is crucial to include both contributions. Hence a primary aim of this paper is to implement the corrections in this way in a fully flexible Monte Carlo code.
Finally we observe that for ZH production a gluon-initiated loop arises that interferes with itself at O(α 2 s ). Example diagrams are depicted in Fig. 5. Due to the enhancement from the gluon parton distribution function (pdf) these contributions represent a large part of the NNLO correction, particularly in the boosted regime that is defined by high vector boson, or Higgs boson, transverse momentum [55]. Throughout this paper we include the gg → ZH diagrams in the y t contribution. This is a slight abuse of nomenclature since, as can be seen in Fig. 5, there is a triangle diagram corresponding to gg → Z * → ZH, where the virtual Z boson radiates a Higgs boson. The gg → HZ contribution was included in the on-shell prediction of [7] and in the more differential calculation of ref. [11]. NLO corrections have been considered in the heavy top limit in ref. [56] and further improved through soft gluon resummation in ref. [57].
The one-loop processes R I , R II and gg → HZ can all be calculated in the full theory in which the top mass is retained. The calculation of the two-loop V I and V II contributions in the full theory is much more complicated and at present, the master integrals are not fully known. Therefore in our calculation and that of ref. [8] an asymptotic expansion in m t is performed. Since both V and R pieces are separately finite, there is some freedom in how the top quark is treated in each part of the calculation. Our strategy is to include the full top mass effects where possible and to perform an asymptotic expansion only when needed. Finally we note that R II was found to have a very small effect on the total prediction in ref. [8] so we do not include it in our calculation of ZH. Technical details regarding our implementation of these pieces in MCFM are presented in Appendix B.

Decays of unstable bosons
The aim of this paper is to present a fully flexible Monte Carlo code for the associated production of a Higgs and vector boson. A crucial element of this flexibility is to ensure that the relevant decays of the Higgs and vector bosons are included. Decays of vector bosons to leptons represent the cleanest experimental signature of these processes, so in this paper we focus on the decays W → ν and Z → + − . For the Higgs boson the bb decay is the most useful, primarily due to its high yield rather than its experimental cleanliness [58,59]. However, H → W W * decays also provide a viable experimental signature with the current data set [60,61]. The high luminosity runs of the LHC may also be able to study rarer channels, such as V H → γγ. These channels are much cleaner, since the dominant irreducible backgrounds from V γγ are much smaller, but the small H → γγ branching ratio makes detailed studies of this process impractical at present.
The decays discussed above are easily incorporated in Monte Carlo codes. However, when considering higher-order QCD corrections, the decay H → bb requires further discussion since radiation can occur in both production and decay stages. At NLO (O(α s )) the conservation of colour ensures a complete factorization between production and decay processes. The situation is more complicated at NNLO (O(α 2 s )) since here contributions exist which connect the initial-and final-state partons. It has been shown [62,63] that these these (non-factorizable) pieces contribute to the total rate at order Γ H /m bb where Γ H is the Higgs boson width. The situation for differential distributions is less precise, but it is plausible that for distributions where b andb are not distinguished the non-factorizable contributions should be similarly small. In our calculations we therefore neglect such effects and follow a factorized approach in which the Higgs decay is included to NLO accuracy. We follow the procedure outlined in refs. [10,11] and define, In the above equation dσ term in the perturbative expansion for the production of a Higgs boson and a pair of leptons. dΓ i H→bb represents the differential partial width at O(α i s ) for the H → bb decay, whilst Γ i H→bb represents the integrated partial width for these decays. In order to study the effect of the pure NNLO corrections it is also useful to define, defines the prediction that treats both radiation in production and decay stages at the NLO level.
Radiative corrections to the decay H → bb were first computed over thirty years ago [64]. It was shown that there are large differences between the partial width in a "massless" theory, in which the b-quark mass is kept in the Yukawa coupling but dropped in the matrix element and phase space, and the full theory in which a non-zero bottom quark mass is retained throughout. These large differences are the result of logarithms of the form log (m 2 b /m 2 H ) that can be absorbed into a re-definition of m b in the Yukawa coupling. As a result, if the running bottom quark mass is used then the massless and massive predictions are very similar. In our MCFM implementation we keep the mass of the b-quark in full, and do not run the b-quark mass in the LO partial width. In order to ensure that the parts of the cross-section that are only exposed to a LO partial width are not susceptible to the running mass corrections, we divide out the partial width and normalize to the branching ratio, BR(H → bb). In this way we can also take advantage of advanced theoretical predictions for this branching ratio, which is now known to O(α 4 s ) [65]. In our MCFM implementation we use the value obtained from the HDECAY code [66]. We note that, although we do not currently include effects beyond NLO in the decay, differential calculations for these quantities have been presented in the massless theory [67,68].

LHC Phenomenology
In this section we study the phenomenology of the V H processes at NNLO for the LHC Run II. For the larger rate H → bb decay we present results which can be compared to data collected with the current operating energy of √ s = 13 TeV. For the rarer H → W W * and H → γγ processes we instead focus on predictions which may be compared with a larger data set obtained in a future √ s = 14 TeV run.
Our predictions are obtained using the default MCFM EW scheme, which corresponds to the following parameter choices: m W = 80.398 GeV, m Z = 91.1876 GeV, Γ W = 2.1054 GeV, Γ Z = 2.4952 GeV, G F = 1.16639 × 10 −5 GeV −2 and m t = 172 GeV. These are sufficient to determine the remaining EW parameters. We use m b = 4.75 GeV and set the CKM matrix elements V ud = 0.975 and V cs = 0.222. Jets are clustered using the antik T jet algorithm with distance parameter R = 0.4. Unless otherwise stated, we use the CT14 pdf sets [69] matched to the appropriate order in perturbation theory. Our default renormalization and factorization scale choice is µ R = µ F = µ 0 with µ 0 = m V + m H . At NNLO the dependence on the unphysical renormalization and factorization scales is rather mild, especially for qq initiated processes such as those under consideration here. As a result any prescription for estimating the scale uncertainty, for instance by varying both scales in the same direction or varying them in opposite directions, yields similar results. However since the H → bb decay is computed at NLO, a larger scale dependence for this decay is observed, with the largest deviations at NNLO (for W H) arising from the case where the scales are varied in opposite directions. We will therefore present results obtained with µ R = kµ 0 and µ F = µ 0 /k, with k = 1/2 and k = 2.
In general our results will show that, once the scale uncertainties discussed above are taken into account, the results for NNLO cross sections still do not overlap those for NLO. This is consistent with other NNLO studies of processes that only receive contributions through qq initial states at LO, since the gluon parton distribution is dominant at the relevant partonic energy fractions of the LHC. At LO there is only a very mild scale dependence which is completely induced by the factorization scale in the parton distribution functions. At NLO the cross section becomes sensitive to the renormalization scale, but typically there is an accidental cancellation between the renormalization and factorization scales resulting in a weak scale dependence even at NLO [44]. Therefore interpreting the scale variation as indicative of the total theoretical error is unwise at NLO. It is difficult, without knowledge of the N 3 LO cross section, to predict whether the scale variation at NNLO will incorporate higher order predictions. However there is reason to believe this may be the case. Firstly the process has access to all initial state configurations, so there will be no new partonic channels at N 3 LO. Secondly the recently-reported calculation of the Higgs cross section at N 3 LO [70,71], is within the scale variation of the NNLO cross section for the first time in the perturbative expansion. Therefore we are reasonably confident that our scale variation can be interpreted as an indicator of theoretical uncertainty. As we look at more exclusive quantities, such as cross sections differential in the number of associated jets, this argument begins to break down and scale variation should not be taken as a rigorous estimate of the theoretical uncertainty.
A detailed study of color-singlet production (including V H processes) using N -jettiness slicing in MCFM 8 is presented in ref. [46]. We refer the interested reader to the detailed discussion of the methodology in that paper and instead briefly summarize the checks here.
Starting from the 0-jettiness of a parton k with momentum p k , where E a , E b are the energies of the beams [23], we define the 0-jettiness as the sum over all the M final state parton jettiness values, M takes the values 0, 1 and 2, depending on the particular phase space component of the NNLO calculation. We then define the above cut region as τ 0 > τ cut and the below cut region as τ 0 < τ cut . The terms above and below the cut combine to leave a residual dependence on τ cut that takes the form, in the limit that τ cut /Q → 0. Here Q defines a hard scale in the LO process (for us m V + m H ) [46] and c 2 and c 3 are coefficients that can be fitted numerically if desired. As a test of our implementation, we have validated our calculation in the absence of any cuts on the final state particles by comparison with the public code vh@nnlo [7,72]. In the limit τ cut → 0 the two are in perfect agreement. We also note that we have checked the calculation of the y t contributions for on-shell bosons with vh@nnlo, also finding perfect agreement. As a detailed discussion of the vh@nnlo checks are provided in ref. [46] we instead focus on similar fits for the phenomenologically relevant processes, in which bosonic decays are included, in the following section. The routines for decaying the Higgs boson in MCFM are well established and have been checked against calculations of branching ratios available in the literature.

Results: H → bb
In this section we present our results for LHC phenomenology for the H → bb decay. We will study a variety of phase space selection criteria, with cuts inspired by the ATLAS [58] and CMS [59] experiments. We define the following set of basic cuts, Before proceeding we first examine the dependence on the τ cut parameter in the context of the cuts specified in Eqs. (3.4)-(3.7). The τ cut -dependence of the NNLO coefficient in the expansion of the cross-section under these cuts (∆σ N N LO ), computed for the 14 TeV LHC, is shown in Figure 6. The remaining dependence on τ cut is a result of power corrections, Figure 6. The τ -dependence of the NNLO coefficient for the V H processes, under the H → bb cuts of Section 3.1. whose form is given in Eq. (3.3) in the previous section. The dashed lines indicate the fitting errors on the asymptotic result for τ cut → 0. It is clear that for W H production ∆σ N N LO is independent of τ cut at the level of a few percent for τ cut 0.01. For ZH production the same value of τ cut yields an accuracy of about 0.5% in the coefficient, where, as we will see shortly, the improvement is due to the fact that this process receives a larger y t contribution that does not depend on τ cut . The accuracy of the prediction for the NNLO cross-section can be assessed by combining this information with the order-by-order results that are shown in Table 1. It is thus clear that choosing τ cut = 0.01 is sufficient for permille accuracy in the full NNLO prediction for all processes. We shall make this choice henceforth.
In Fig. 7, we present the cross-section as a function of the LHC operating energy given the basic selection cuts described above. The left-hand plot illustrates the total rate at NLO (dashed) and NNLO (solid) for W H and ZH production. We plot the cross-section for W + (→ + ν)H and W − (→ − ν)H separately. At the LHC the production of W + H is dominant, since the ud initial state configuration has a larger flux than du for pp collisions.  to the total NNLO coefficient. For W H production the top induced pieces, after cuts, make up around 30-50% of the total O(α 2 s ) correction. The Z boson has a much smaller branching ratio to a single family of leptons compared to the W , so the cross-sections presented in Fig. 7 for + − bb are smaller than the corresponding W induced ones (for instance, compared to the inclusive results presented in ref. [46]). The NNLO corrections are also much more important for ZH production than for W H. This is due to the large contribution from the gg → ZH pieces. The importance of the gluon flux at the LHC can help to offset the α S suppression, resulting in a NNLO correction whose impact is more comparable to a NLO effect. This is clearly visible in the lowest plot on the right-hand side of Fig. 7  . Of the σ (2),yt ZH contribution the dominant effect is induced by the gg diagrams, although it is not possible to separate them from the V i and R i pieces at this order in perturbation theory.
The cuts described above result in a fairly inclusive selection. In order to reduce the backgrounds from top, diboson and V + jets processes, cuts on the transverse momentum of the vector boson are usually employed in experimental analyses. We therefore investigate the total cross section as a function of the minimum transverse momentum of the vector boson in Figs. 8 (W H) and 9 (ZH). We focus on the LHC operating at √ s = 13 TeV. To produce these results we apply the basic cuts described above. The results of the previous figures are also manifest in these plots: the initial impact of higher order corrections for W − H is slightly larger (at NLO), but the impact of the NNLO corrections is similar for both charges in W ± H. It is also clear that ZH has much larger NNLO corrections than W H. Particularly rich signal bins in the experimental analysis correspond to p V T > 120 GeV and p V T > 160 GeV. For these choices the signal cross-section is around 30-40% and 15−20% of the p V T -inclusive result, respectively. The impact of NLO is a mild enhancement in the tail of the p V T distribution for all process. For W H production the NNLO corrections are reasonably flat in p V T , while the NNLO corrections to the ZH process become more pronounced in the high p V T region. This is due almost exclusively to the y t correction,  which hardens the spectrum as can be clearly seen in the middle panel of Fig. 9.
We now turn our attention to jet-based observables. In Figure 10 we present the cross- section as a function of the total number of jets (i.e. b-jets plus light jets). The plot on the left-hand side has only the basic lepton cuts applied, while on the right p V T > 120 GeV is required in addition to the basic lepton cuts. Since the Higgs boson is a resonance decaying to massive quarks, a well-defined cross-section can be computed without any requirement on the number of b-jets present. An N n LO prediction can then have between 0 and (n + 2) jets in the final state, with the (n + 2)-jet bin corresponding to the LO prediction for V H + n jets. In Fig. 10 we present NLO and NNLO predictions, with NLO (solid) and LO (dashed) H → bb decays. As expected the largest scale variation occurs in the four-jet bin (NNLO) and three-jet bin (NLO) predictions, since these are LO predictions in this observable. Including the decay has a significant impact on the jet counting, particularly in the two-and three-jet bins where it changes the predicted rate by O(20%). The higher p V T selection has relatively more three-jet events than the more inclusive selection, which arises from the kinematic favorability of balancing a high p T vector boson with a jet and the Higgs boson. The Higgs decay at LO has no dependence on α S and therefore no scale dependence. When we introduce the decay at NLO we include α S for the first time, and acquire a larger dependence on the choice of scale. Reducing the scale dependence further requires consideration of NNLO effects in the decay stage [67,68].
The impact of including the NLO decay is shown differentially in Figs. 12 (W H) and Figs. 13 (ZH). We present differential distributions for the hardest b-jet (left) and p bb T (right), applying the basic lepton cuts, demanding two b-jets and enforcing p V T > 120 GeV. The general impact of including the higher order corrections in the decay is immediately apparent as a general softening of both spectra. This is easily understood from the decay kinematics. The invariant mass of the system is constrained to be very close to m 2 H . When there are three particles present in the decay to share the energy, the result is a softer spectrum. For p bb T the impact of higher order corrections in production and decay are  particularly important. At LO, p bb T = p V T so that the cut on the vector boson momentum is also a cut on p bb T . At NLO this is no longer necessarily the case, since the real corrections allow for an unclustered parton to balance the total momentum. Therefore the region p bb T < 120 GeV is first accessible at NLO. Since in this region of phase space the total transverse momentum of the bb is by definition relatively small, the resulting transverse momentum of the b-quark pair is also relatively soft. As a result the region of phase space where p b T (hard) < 120 GeV also has large higher order corrections. This is highlighted in the middle panel of the figures which presents the impact of higher order corrections in production (for NLO decays). Going from LO to NLO there are large corrections to the p T spectrum of the hardest b quark, however the NNLO prediction is relatively stable illustrating that the perturbative expansion is well-behaved beyond LO.
For p bb T there is a strong feature at the edge of the phase space for the NLO decays that is not present for the LO decay option. This is due to the phase space boundary at the LO threshold in the decay phase space [73]. The virtual decay corrections reside in the p bb T > 120 GeV region of phase space, whilst the real corrections H → bbg can fill the region both above and below this value. However in the real phase space when p bb T = 120 GeV there is a restriction on the phase space for soft gluon emission and a large logarithm arises. Boundary problems such as these occur frequently in perturbation theory [73] and have  been observed for this specific process in previous calculations [10,11]. For ZH production the spectrum is smoothed-out somewhat by the turn-on of the gg → ZH contribution, which resides in the Born phase space.

Results
The decay H → W W * represents a significant fraction, about 20%, of the total decay rate. Although requiring leptonic decays of the W bosons reduces the rate further, the signal cross sections are large enough to have warranted experimental investigation in Run I of the LHC [74]. Going forward into Run II, triboson signatures represent a fresh environment in which to search for new physics. From the technical point of view the 2 → 8 phase space is large, corresponding to a 22-dimensional Monte Carlo integration in the double-real part of the calculation. This therefore represents a demanding application of the N -jettiness slicing technique that tests its suitability for complex phenomenological studies at NNLO. We study the production of V H(→ W W * ) at the 14 TeV LHC with the following simple phase space selection criteria,

MET
: Results for the cross-sections obtained under these cuts are presented in Table 2. We present results for a single family of leptons for each V decay. Including all possibilities of electron and muon configurations would thus increase the total rates presented in the table by a factor of eight. We do not consider any interference between the leptons arising from the decay of the associated vector boson and those arising from the Higgs boson decay. The primary aim of this section is to address the feasibility of producing NNLO results using the N -jettiness slicing method for a high-dimensional final state configuration. We therefore postpone the treatment of interference effects for a future study. An important variable when considering the decay H → W W * is the transverse mass of the electroweak final state. For the W H process it is defined by, The equivalent definition for ZH production (m ZH T ) is obtained by making the replacement 3 → 4 . The transverse mass is important since it can be used as a proxy for m V H , which is not experimentally observable for H → W W * decays. The study of this variable is interesting due to its sensitivity to high-energy structures that may be present in the HV V vertex in BSM scenarios. For instance, treating the SM as an EFT introduces six-(and higher-) dimensional operators that induce momentum dependent couplings between the Higgs and the vector bosons 2 . In general these operators will induce small deviations in the tails of the transverse mass distribution, so it is crucial to have control over percent-level effects from the SM in this region.
Our predictions for this observable are presented in Fig. 14 for W + H (left) and ZH (right). The difference in shape between m W H T (W + H left) and m ZH T (ZH right) is apparent. Since the transverse mass for ZH is closer in definition to m V H , it has a much harder spectrum, whereas m W H T is much softer. The higher-order corrections for ZH are much larger, primarily due to the significant gg → ZH contribution. These pieces are sensitive to the 2m t thresholds present in the loop integrals, and as a result turn on at m ZH T ∼ 300 GeV. 2 A study of these operators at NLO (+PS) in the MCFM framework was presented recently in ref. [75].  In this region the NNLO corrections are much larger and the shape of m ZH T is significantly altered. It is therefore essential to include NNLO predictions for this observable in order to avoid attributing any observed change in shape to the presence of BSM physics. Fig. 15 we present the cross-section as a function of the number of additional jets, where the basic jet definition is used from the previous section, p j T > 25 GeV and |η j | < 2.5. The n j distribution for these decays are different from those studied previously in the H → bb section, since now the jets are only produced through initial state radiation, with no contamination from jets arising from the decay. For the W H process, around 40% of the events have one or more jets in the final state. For ZH production the percentage  Table 3. Cross-sections for V γγ production at the 14 TeV LHC, for a single family of leptons. The corresponding phase space selection criteria are described in the text.

Finally in
drops to around 35% due to the presence of the gg → ZH contribution that only populates the 0-jet bin.

Results: H → γγ
The decay H → γγ provides a relatively clean experimental signature, at the cost of a very small branching ratio. However during Run II enough data should be collected to allow experimental studies of this channel. In Table 3 we collect cross-sections for V H(→ γγ) processes at the LHC operating at 14 TeV, after application of the following basic selection criteria: Photons : Jets : p j T > 25 GeV, |η j | < 2.5 (3.14) Leptons : We set the central renormalization and factorization scale as µ 0 = m 1 2 γγ and, as in the previous section, vary the central scale by a factor of two in opposite directions. Although the cross-sections for this process are rather small, around 0.05 fb, the advantage they possess over H → bb decays is a much smaller irreducible background. We illustrate this by including in Table 3 the background cross-sections in the Higgs resonance region (σ Back LO,mγγ ) that is defined by, 120 GeV < m γγ < 130 GeV (3.17) These cross-sections are for illustration only and are therefore evaluated at LO. The results of Table 3 clearly demonstrate the potential for an excellent signal-to-background ratio in this channel, although we note that experimental analyses would also have to cope with a large reducible background from V +jets that may contaminate this significantly. Unsurprisingly the impact of the NNLO corrections for this decay channel are essentially the same as those discussed in detail in the previous sections. We demonstrate the impact on a canonical differential observables in Fig. 16, which depicts the p γγ T spectrum for both W + H and ZH production at the 14 TeV LHC.

Conclusions
In this paper we have presented a NNLO calculation of V H production and its implementation into a fully flexible Monte Carlo code. These processes can provide a useful handle on the coupling of the Higgs boson to bottom quarks, in addition to serving as sensitive probes of anomalous interactions between the Higgs, W and Z bosons. At NNLO the calculation of these processes includes both Drell-Yan-like contributions and corrections where the Higgs boson is radiated from a heavy quark loop rather than from the vector boson. These two contributions are comparable in size for W H production. For the ZH process the latter corrections involve gg-initiated diagrams that dominate the NNLO correction. Including both contributions is therefore imperative and our calculation enables the combined effects to be studied in a differential manner for the first time. Analytic results for these amplitudes can be found in the appendix and the distributed MCFM code.
The V H processes are among the most interesting Higgs production modes for phenomenological studies in Run II. We have therefore studied a number of decay modes of the Higgs boson in some detail. In the case of the decay to a pair of bottom quarks, H → bb, we have also included the effect of radiation in the decay at NLO accuracy. This has a considerable impact on, for instance, the transverse momentum of the bb pair. Our results suggest that including QCD corrections in the decay to NNLO is important, although it is beyond the scope of this paper. We have investigated the phenomenology of other Higgs decay channels that might be explored more fully in Run II of the LHC, namely the decays H → W W * → leptons and H → γγ.
From the theoretical point of view, presenting a NNLO calculation for the 2 → 6 process, pp → V H → V W W * → leptons+ / E T is technically challenging due to the large final state phase space. In the double-real part of the calculation the phase space is 22dimensional, and provides a challenging environment in which to test the jettiness-based approach to NNLO calculations. The H → γγ decay results in very small cross-sections, but has the advantage that the irreducible background in the neighborhood of the Higgs boson mass is also small, resulting in comparable signal and background rates.
The results of this paper have been implemented into MCFM, and are publicly available.

A Amplitudes for DY contributions
In this appendix we present analytic expressions for the parts of the calculation that are closely-related to the Drell-Yan process. We present results for contributions above and below τ cut separately. Since the Higgs boson is a scalar particle the decay amplitude factorizes from the production amplitude, modulo O(α 2 s ) corrections that we neglect in our calculation. Therefore the results in this section are presented for an on-shell Higgs boson and a V boson that decays to leptons are included. We have used the Mathematica package S@M [76] frequently in our calculations.

A.1 Below τ cut
Below τ cut the calculation requires the soft, beam, and hard functions of the SCET formalism. The necessary two-loop soft and beam functions were computed in refs. [49,50] and [51] respectively. The process-dependent hard-function can be extracted from the two-loop virtual form-factor as in refs [53,54]. We have repeated the calculation and for completeness we reproduce the 1-and 2−loop results below, where L = log −s 12 /µ 2 and M (0) represents the LO amplitude. Note that the above amplitudes are defined in four-dimensions and as such the LO amplitude (for W H production) can be defined in terms of helicity amplitudes as follows, In the above equation P represents the propagator function, The tree-level helicity amplitude is then defined as follows 3 , It is instructive to re-write this amplitude as, Modulo the overall factor of s 12 , which is a result of the internal W propagator function, the holomorphic piece of the above expression (the first term) corresponds exactly to the MHV tree-level amplitude for qq ν. The second, non-holomorphic, term is a correction that vanishes in the soft Higgs boson limit. Since the amplitudes for W + 3 and W + 4 partons have been calculated analytically [78,79] the above decomposition provides a useful check of all of our calculated amplitudes. The helicity breakdown for the ZH process requires a summation over the left and right-handed couplings The fermionic (quark or lepton) coupling to the Z boson is given by v The sign in the v q L term is determined by whether the quark is up (+) or down (−) type. The helicity amplitudes are then obtained from the equivalent amplitudes for the W H process, by applying line reversal symmetries as necessary.

A.2 Above τ cut
Above τ cut the calculation corresponds to a NLO one for the V H+ jet process. We have calculated helicity amplitudes for this process which, to the best of our knowledge, have not been presented in the literature before. The LO amplitude for W Hj can be written as follows, The tree-level MHV helicity amplitude is defined as, 34 15 25 (A.10) where P 25 = p 2 + p 5 . The helicity amplitude for h 5 = −1 can be readily obtained from the above by the conjugation operation. At NLO we require the one-loop amplitude for W Hj and the tree-level amplitudes for W Hjj. The one-loop amplitude for W Hj can be written as follows, As discussed in the previous subsection, the amplitudes for V Hj can be decomposed in terms of those for V j and a piece which vanishes in the soft limit (where it is understood that momentum conservation is altered accordingly in the V j amplitude). We use this decomposition on our one-loop amplitudes, defining, where A (0) is understood to be given by eq. (A.10). V i and F i then correspond exactly to the amplitudes obtained in the calculation of the W + 3 parton one-loop amplitude, and S H i corresponds to the missing piece which vanishes in the soft Higgs limit. For brevity we do not reproduce the results for V i and F i here, they can be found in the literature [79] or in the distributed MCFM code. The leading colour contribution S H 1 has the following form,  [21] (A.14) Here F 1m 4F represents the finite part of the one-mass box integral and is given by, F 1m 4F (s, t; P 2 ) = −2 Li 2 1 − P 2 s + Li 2 1 − P 2 t + 1 2 log 2 s t + π 2 6 (A. 15) and the auxiliary function L 0 (s, t) is, L 0 (s, t) = log s/t 1 − s/t (A.16) The amplitudes for ZH production can be obtained in similar fashion to the tree-level discussion in the previous sub-section, i.e. by appropriate dressing of the helicity amplitudes by the left-and right-handed couplings and modification of the electroweak pre-factor accordingly.

B Amplitudes for O(α 2 s ) y t contributions
In this section we present formulae for some of the amplitudes that contribute to the term dσ (2),yt . There are in principle five such terms, labelled in ref. [8] as V I , V II , R I , R II and gg-initiated pieces. In this section we follow closely the formalism and nomenclature used in that reference, which presents these contributions for on-shell bosons. Here we do not include any results for R II since its effects are tiny and not accounted-for in our calculation. We also do not present our analytic results for the gg → ZH contribution, which are too lengthy to include here but may be inspected in the distributed MCFM code.

B.1 V I pieces
We begin by considering the V I pieces, which occur for both W and Z associated production (cf. Fig. 3, right). Using the method of asymptotic expansions 4 it was shown in ref. [8] that the leading terms in the expansion correspond to replacing the top quark loop by the effective ggH vertex. Therefore we can obtain V I by calculating the results in the effective field theory. We write the one-loop amplitude as follows,  [43] (B. 2) The finite part of the two-mass easy box is defined as F 2me 4F (s, t; P 2 , Q 2 ) = −2 Li 2 1 −

B.2 V II pieces
Next we consider the V II diagrams, in which the vector boson also couples to the closed fermion loop. (cf. Fig. 4, left). These contributions therefore only exist for ZH production.
The results of ref. [8] show that the leading terms in the m t asymptotic expansion correspond to the qqZH effective vertex multiplying a two-loop massless tadpole diagram. The leading term in the expansion thus has the form of a tree-level amplitude The (v t L − v t R ) factor arises since only the axial part of the two-loop massive tadpole contributes to the amplitude.

B.3 R I pieces
Finally we present the results for the R I contribution (cf. Fig. 3, left). These contributions are universal and occur for either W H or ZH production. They correspond to one-loop calculations and can be performed either including the top quark mass in full, or in the effective field theory. The general structure of the result is as follows, R I (1 − q , 2 + q , 3 − , 4 + , 5 h 5 g , p H )I F T /EF T (B.6) The helicity amplitude is universal, but the prefactor J depends on whether one is in the full or the effective field theory, Finally, the function I is either equal to unity in the effective field theory or is a combination of scalar loop integrals in the full theory,