Next-to-next-to-leading order event generation for $VH$ production with $H\to b\bar{b}$ decay

We consider the Higgsstrahlung process in hadronic collisions and present the computation of next-to-next-to-leading order predictions matched to parton showers for both production and $H\to b\bar{b}$ decay employing the MiNNLO$_{\rm PS}$ method. We present predictions for $ZH$ and $W^{\pm}H$ production including spin correlations and off-shell effects by calculating the full processes $pp \to \ell^+\ell^-H \to \ell^+\ell^-b\bar{b}$, $pp \to \nu_\ell\bar\nu_\ell H \to \nu_\ell\bar\nu_\ell b\bar{b}$ and $pp \to \ell^\pm \nu_\ell H \to \ell^\pm \nu_\ell b\bar{b}$ in the narrow-width approximation for the Higgs boson. For the $W^{\pm}H$ process, NNLO+PS accuracy in production and decay is achieved for the first time. Our calculations are validated against earlier simulations in the NNLOPS approach that includes NNLO corrections via multi-differential reweighting. The new MiNNLO$_{\rm PS}$ generators for these processes, which evaluate NNLO corrections on-the-fly in the event generation, will supersede those earlier calculations. Our predictions are in good agreement with recent measurements of the Higgsstrahlung cross sections.


Introduction
The Higgs boson is a unique particle in the Standard Model of Particle Physics (SM). Although it interacts with all massive particles, which obtain their masses precisely through their interaction with the Higgs field, only with the high-energy proton-proton collisions of the Large Hadron Collider (LHC) the Higgs boson could be discovered due it its relatively large mass of 125 GeV. The discovery in 2012 was based on 5 fb −1 of data collected at 7 and 8 TeV [1,2]. A discovery in this early stage of the LHC was possible since the Higgs boson has a sufficiently large branching fraction to photons and Z bosons, which lead to experimentally very clean signatures with an excellent mass resolution. Since the discovery studying the properties of the observed scalar resonance has been one of the cornerstones of the rich physics programme at the LHC. In particular, the Higgs sector provides a window to new-physics phenomena that could for instance affect the Higgs couplings. With all Higgs properties appearing to be SM-like thus far and without any other clear signs of physics beyond the SM (BSM), precision tests of the Higgs sector have become a valuable path towards the observation of new physics.
Higgs-boson production at the LHC proceeds predominantly through gluon fusion and is induced by a top-quark loop, while the Higgs decay to bottom quarks (H → bb) has the largest branching ratio. Other Higgs-production modes are suppressed by roughly one order of magnitude or more with respect to gluon fusion, but they are indispensable in the study of the Higgs properties to determine as many of the Higgs couplings as possible. For instance, the associated production of a Higgs boson and a vector boson (V H production or "Higgsstrahlung") or vector-boson fusion (VBF) provide access to Higgs couplings to vector bosons, increasing the sensitivity to those couplings in addition to the Higgs decay modes to vector bosons. Even more importantly, the H → bb decay, which provides direct access to the bottom Yukawa coupling, could only be observed [3,4] by exploiting the V H production mode with a boosted Higgs boson. In fact, it had been thought for a long time that Higgsstrahlung would not be accessible at the LHC until it was proposed to consider V H events where the Higgs boson has large transverse momentum (i.e. it is boosted) and to apply substructure techniques to identify the bottom quarks and reconstruct the Higgs boson [5]. This renders the accurate modeling of both V H production and the H → bb decay at the level of fully exclusive events indispensable. While this was already crucial for the experiments to observe this decay mode, it will become even more critical for an accurate determination of the bottom Yukawa coupling in the future.
A remarkable progress in next-to-next-to-leading order (NNLO) QCD calculations for SM processes has been made in the past years, with all 2 → 2 colour-singlet scattering processes computed at this accuracy , and even first NNLO calculations for genuine 2 → 3 processes emerging in the recent past [40][41][42][43]. The Higgsstrahlung process is actually much simpler from a computational point of view: As far as QCD corrections are concerned the calculation corresponds to an extension of the Drell-Yan (DY) process in which the Higgs boson is emitted from the vector boson. Indeed, the inclusive NNLO cross section was computed almost twenty years ago [44][45][46], and about ten years later fully differential NNLO calculations appeared [6][7][8][9][47][48][49]. Part of the NNLO corrections to ZH production is the loop-induced gluon-fusion (gg) process, which is mediated by a top-quark loop, while other contributions mediated by a top-quark loop have been shown to be relatively small [45] and are thus often neglected in differential calculations. In fact, the effectively only leading-order (LO) accurate loop-induced gg process dominates the theoretical uncertainties in ZH analyses, which calls for the inclusion of NLO corrections to gg → ZH production, especially since first approximations suggested that these are quite large [50][51][52][53]. Recently, substantial progress was made in the calculation of the relevant two-loop amplitude to (on-shell) gg → ZH production [54][55][56], which eventually led to the first calculation at NLO QCD in a small mass expansion [53]. Strategies to isolate the gg → ZH component in experimental analyses have been suggested in ref. [57]. The analytic resummation of logarithms at threshold and due to a jet veto in V H production have been considered in refs. [51,[58][59][60]. Recently, also W ± H production in association with a jet was calculated at NNLO accuracy in QCD [61]. As far as NLO electroweak (EW) corrections are concerned, they have been available for a while [62][63][64].
Since V H measurements exploit a particularly differential strategy, requiring a boosted Higgs and substructure techniques, as discussed before, it is mandatory to combine state-of-the-art higher-order calculations with full-fledged event simulations through a parton shower (PS) Monte Carlo generator. NLO QCD corrections matched to parton showers are available as a merged calculation for V H + 0, 1-jet production within Powheg [65] using the MiNLO method [66,67]. Similarly, both NLO QCD and EW corrections were simultaneously matched to a parton shower in a merged V H+0,1-jet calculation in ref. [68]. Predictions at LO+PS for gg → ZH can be easily obtained with several tools, including with Powheg-Box [65]. Merged predictions at LO, up to gg → ZH+ 1-jet, have been considered in refs. [69,70], and a study on how different approximations compare has also been performed in ref. [71]. The effects from anomalous couplings in V H production have been considered within Powheg in ref. [72]. In recent years, a great effort has been devoted to the matching of NNLO QCD calculations and parton showers (NNLO+PS), and essentially four different methods have been introduced: The NNLOPS method [73] based on reweighting MiNLO events [67], the UNNLOPS method [74,75], the Geneva method [76], and, most recently, the MiNNLO PS approach [77,78]. At first these methods have been applied to the simplest colour singlet processes, namely Higgs-boson production [73,75,79] and Drell Yan [74,80,81]. Especially in the NNLOPS, MiNNLO PS and Geneva approaches a variety of colour-singlet production processes are now available at NNLO+PS, including ZH [82,83], W ± H [83,84], γγ [85], Zγ [86,87], W ± γ [88], ZZ [89,90] and W + W − [91,92] production. Very recently, also the NNLO+PS matching for the very first production process with coloured particles in the initial and final state has been achieved in ref. [93], which extents the MiNNLO PS method to top-quark pair production.
Similarly important in this context is an accurate and fully exclusive description of the H → bb decay. The total decay rate for massless bottom quarks is known to N 4 LO QCD by now [94][95][96][97][98][99][100][101][102]102]. Electroweak corrections were computed in refs. [103][104][105][106]. We refer to refs. [107,108] for a review of calculations for the inclusive H → bb branching ratio. Fully differential calculations are available at fixed-order through NLO QCD for massive bottom quarks since a long time [109][110][111][112]. NNLO QCD predictions have been first computed for massless bottom quarks [8,[113][114][115], and the corresponding results retaining massive bottom quarks were presented in refs. [116][117][118][119]. Recently, even the N 3 LO corrections were calculated in the massless approach [120,121]. The effect of top-quark loops have been studied in refs. [122,123]. A complete calculation at NLO QCD in the SM Effective Field Theory (SMEFT) was presented in ref. [124], while anomalous couplings between the vector bosons in the production of the V H final state, including the H → bb decay, were considered at NNLO QCD in ref. [125]. As far as Monte Carlo predictions are concerned, the simulation of H → bb decay events has been achieved at NNLO+PS in refs. [126][127][128].
In this paper, we present a fully exclusive simulation of ZH and W ± H production with H → bb decay in the MiNNLO PS framework. Both production and decay events are simulated at NNLO+PS accuracy and consistently combined to obtain accurate predictions for the + − bb, ν ν bb and ± ν bb final states, including off-shell effects and spin correlations of the vector bosons and treating the Higgs boson in the narrow-width approximation. NNLO+PS accuracy for the pp → ν ν H → ν ν bb and pp → ± ν H → ± ν bb processes is achieved for the first time. Our implementation has been done both within Powheg-compatible within numerical uncertainties. For the + − bb final state our MiNNLO PS predictions are compared against the NNLOPS results of ref. [126], while for ± ν H production we validate our calculation against the NNLOPS results presented in ref. [84] for on-shell Higgs bosons. The new MiNNLO PS implementation shall supersede those NNLOPS implementations in the future within the Powheg-Box framework [129]. Finally, we present a complete study of NNLO+PS events for both ZH and W ± H production with H → bb decay, finding agreement with a recent LHC measurement from the ATLAS collaboration [131].
This manuscript is organized as follows: in section 2 we introduce our calculation by discussing the processes (section 2.1), sketching the MiNNLO PS method (section 2.2), describing our practical implementation (section 2.3) and recalling the combination of production and decay events (section 2.4). An extensive validation of our MiNNLO PS generator is presented in section 3 against NNLOPS prediction for W ± H production without decay (section 3.1) and for ZH production with H → bb decay (section 3.2). In section 4 we present new results for both ZH and W ± H production including the H → bb decay: After introducing the setup (section 4.1), we discuss inclusive and fiducial cross sections (section 4.2), distributions in the fiducial phase space (section 4.3), the effect of different clustering algorithms in the definition of bottom-quark jets (section 4.4), and compare against data in (section 4.5). We conclude in section 5. While in the main text we focus on results for W + H and ZH production, in appendix A we present results for W − H production.

Description of the processes
We consider the processes  Figure 1: Sample Feynman diagrams for V H production in the H → bb decay channel with (a,c) V = Z and Z → + − and (b) V = W ± and W ± → ± ν , where by ν we indicate generally a neutrino or anti-neutrino. Panels (a) and (b) are tree-level diagrams in the quark-annihilation (qq) channel at LO, while panel (c) shows a loop-induced diagram in the gluon-fusion (gg) channel entering at NNLO QCD.
leptons throughout this paper. For simplicity, we will refer to the + − bb and ± ν bb final states as ZH and W ± H production in the following. Sample LO diagrams of the two processes are shown in figure 1 (a) and (b). They are qq-initiated with an s-channel vector boson emitting a Higgs boson and subsequent decays of the two bosons. Figure 1 (c) shows a diagram for the loop-induced process in the gg channel that due to charge conservation exists only for ZH, but not for W ± H production. This process enters the ZH cross section at O(α 2 s ) and is therefore a NNLO contribution relative to the qq-initiated process. Its dominant contribution proceeds via a top-quark loop, since all other quark loops are suppressed by their respective Yukawa couplings to the Higgs boson. In this work, we focus on the accurate simulation of the qq-initiated processes. The loop-induced gg process can be computed completely independently at LO+PS and added to the qq-initiated cross section. We only include the loop-induced gg contribution when presenting integrated cross sections in section 4.2 and when comparing to data in section 4.5. In fact, we would like to stress that very recently the gg → ZH two-loop amplitude for an on-shell Z boson has been calculated [54][55][56], which will enable NLO+PS predictions for the loop-induced gg process in the near future (with the approximation of the Z boson being on-shell in the two-loop contribution). Apart from the best predictions for the qq-initiated ZH production presented here, this calculation will be crucial to reduce the theoretical uncertainties in the experimental analyses of boosted ZH production at the LHC.
We note that there are other contributions to both W ± H and ZH production that are mediated by a closed top-quark loop radiating the Higgs boson, which were calculated in ref. [45]. We adopt the approximation that the virtual diagrams denoted as V I and V II in ref. [45] are neglected in our calculation, while the real contributions R I and R II , where the latter only contributes to ZH production, are included. In order to facilitate the comparison to previous results we have excluded those contributions to W ± H production, but included them for ZH production, as discussed below. We recall that the contribution of the R I/II and V I/II terms to the total cross section is positive. Both R I/II and V I/II contribute to about 0.6%, giving a total positive correction of about 1.2% with respect to the total qq-initiated NNLO result obtained neglecting these terms.

MiNNLO PS method
Our calculation of V H production at NNLO+PS relies on the MiNNLO PS method, which was formulated and applied to 2 → 1 processes in refs. [77,78] and later extended to generic colour-singlet processes in ref. [86] and to heavy-quark pair production in ref. [93]. The interested reader is referred to those publications for details, while here we only sketch the main ideas behind the MiNNLO PS procedure in a simplified notation.
The MiNNLO PS method that formulates a matching of NNLO corrections with a parton shower for a system F of colour-singlet particles proceeds essentially in three steps: In Step I, a Powheg calculation [129,[133][134][135] generates F in association with one light parton at NLO inclusively over the second radiation. In Step II, an appropriate Sudakov form factor and higher-order corrections are applied such that the simulation remains finite in the unresolved limit of the light partons and becomes NNLO accurate for inclusive F production.
Step III involves both the exclusive generation of the second radiated parton (accounted for inclusively in Step I) through the Powheg Sudakov and the generation of subsequent emissions through the parton shower. Since all emissions are correctly ordered in p T (when using p T -ordered showers) and the Sudakov in Step II matches the leading logarithms generated by the parton shower, the MiNNLO PS method preserves the (leading logarithmic) accuracy of the parton shower.
We write the fully differential MiNNLO PS cross section symbolically as which corresponds to a Powheg calculation for F plus one light parton (FJ), but including NNLO accuracy for F production through a modification of theB function. Indeed, ∆ pwg is the standard Powheg Sudakov form factor with a default cutoff of Λ pwg = 0.89 GeV, Φ rad and p T,rad are phase space and transverse momentum of the second radiation, respectively, and B FJ and R FJ denote the squared tree-level matrix elements for FJ and FJJ production, respectively, while Φ FJ is the FJ phase space. In the second line dσ (1,2) FJ are the firstand second-order term in the α s expansion of the differential FJ cross section, while e −S denotes the Sudakov form factor in the transverse momentum of F (p T ), with S (1) being the first term in the expansion of the exponent. Renormalization and factorization scales are evaluated as µ R ∼ µ F ∼ p T in MiNNLO PS and NNLO accuracy is achieved by the last term in theB function, which is of order α 3 s (p T ) and adds the relevant (singular) contributions [77], while regular contributions at this order are subleading. The function D is derived from the following formula that determines the relevant singular terms in p T dσ res which can be obtained by suitably manipulating a p T resummation formula, as explained in section 4 of ref. [77]. L is the luminosity factor up to NNLO, including the convolution of the collinear coefficient functions with the parton distribution functions (PDFs) and the squared hard-virtual matrix elements for F production. Instead of truncating eq. (2.3) at the third order, i.e. evaluating s ), as originally done in ref. [77], we keep all terms of O(α 4 s ) and higher, as suggested in ref. [78], in order to preserve the total derivative in eq. (2.4). Finally, the factor F corr in eq. (2.3) ensures that the Born-like NNLO corrections are consistently included in the FJ Powheg events by spreading them in FJ phase space.

Practical implementation
We have implemented MiNNLO PS generators for V H production both in the Powheg-Box-V2 [129] and the Powheg-Box-Res [130] framework. Both calculations yield fully compatible results within numerical uncertainties. Our Powheg-Box-V2 implementation uses the V H+jet generators developed in ref. [65], which were already used in the former NNLOPS calculations of refs. [82,84]. It then includes NNLO QCD corrections to V H production by exploiting the original implementation of the MiNNLO PS method for Drell-Yan production of refs. [77,78], which is suitable since the respective two-loop corrections for V H production correspond to those of the Drell-Yan (DY) process with the Higgs boson being emitted from the vector boson. Our MiNNLO PS V H generators in Powheg-Box-Res, on the other hand, are based on the V H+jet implementation of ref. [68]. The V H+jet generators have then been upgraded to include NNLO QCD accuracy for V H production through the general MiNNLO PS implementation for colour-singlet production developed for Powheg-Box-Res in ref. [86].
As far as the physical amplitudes are concerned, both calculations exploit the same results, documented in refs. [65,68]. The tree-level real and double-real matrix elements (i.e. for V H+1,2-jet production) are obtained using the automated interface [136] between Powheg-Box and MadGraph v4 [137]. The V H+jet one-loop amplitude is taken from the analytic calculation of ref. [138]. The one-loop and two-loop amplitudes have been obtained by applying the form-factor correction to the quark vertex function (equivalent to Drell-Yan production) to the Born amplitude. We note that the V H+jet generator of ref. [68] in Powheg-Box-Res in principle also allows us to use OpenLoops [139][140][141] for the evaluation of the relevant amplitudes, but we have not made use of that feature. Note that, as discussed before, our V H MiNNLO PS generators include the R I and R II real contributions, first calculated in ref. [45], while we neglect the V I and V II virtual amplitudes, which contribute below one percent.
The calculation of D in eq. (2.4), which is the central ingredient to include the NNLO QCD corrections through MiNNLO PS , involves the evaluation of several convolutions with the parton distribution functions (PDFs), which are performed through hoppet [142]. Moreover, the collinear coefficient functions require the computation of polylogarithms, for which we employ the hplog package [143].
For the calculation of the decay of the Higgs boson to bottom quarks at NNLO+PS we follow closely the procedure of ref. [126] using the NNLOPS method. Thus, pure H → bb events are generated, such that the inclusive H → bb decay width is NNLO QCD accurate and observables involving an extra jet are NLO QCD accurate. While for the generation of decay events we could also apply the MiNNLO PS method, we stick to the NNLOPS approach based on including NNLO corrections through an a-posteriori reweighting. This is because in the present case the reweighting is completely straightforward as one needs to reweight to one single number, i.e. to the NNLO decay width of the Higgs boson to bottom quarks. All details about the reweighting are given in section 2.2 of ref. [126]. In particular, we follow the procedure and settings recommended there, i.e. we use the three-jet resolution variable in the Cambridge algorithm (with y 3,pow = 2 and y 3ref = e −4 ) to smoothly turn off the reweighting of events accompanied by hard radiation, which are already NLO accurate. All relevant amplitudes have been computed in ref. [114].
Finally, we report the most relevant non-standard settings that we have used in this paper. To turn off spurios higher-order logarithmic terms at large transverse momenta we use a modified logarithm for p T > m V H /2, where m V H is the invariant mass of the Higgs and vector-boson system, which smoothly vanishes at p T = m V H , as introduced in ref. [93]. The hard scale in the logarithm is consistently changed to m V H /2 as explained in ref. [144]. At small transverse momenta, we use the standard MiNNLO PS scale setting given in eq. (14) of ref. [78] with Q 0 = 0 GeV. We regularize the Landau singularity by freezing the strong coupling and the PDFs for scales below 0.8 GeV. Furthermore, we activate the option largeptscales 1 to set the scales entering the cross section at large transverse momenta [78]. We turn on the Powheg-Box option doublefsr 1, which was introduced and discussed in detail in ref. [145]. In our interface to the parton-shower we have kept all settings to the standard ones, in particular for the recoil scheme.

Combining V H production and H → bb decay at NNLO+PS
In order to combine V H production at NNLO+PS with the decay of the Higgs boson to bottom quarks at NNLO+PS, we follow the procedure outlined in ref. [126], which we summarize here. First, + − H and ± ν H production events as well as H → bb decay events are generated, using the calculations outlined in the previous section.
Once available production and decay events are merged following a two-step procedure. First, one replaces, in each production event, the Higgs boson with the decay products of a Higgs boson in a decay event. This procedure is explained in all detail in section 3.1 of ref. [126]. Note that decay events are generated on-shell, while we include the Higgs width in the production events. The procedure performs a rescaling of the decay event in such a way that momentum is conserved exactly but small off-shell effects suppressed by Γ H /m H are neglected. Second, the weights of the combined events are set using the production and decay weights. In particular, when multiple weights are present (e.g. to probe scale variations) all possible combinations are included in the final event file. More precisely, the combined weight (w full ) of a production event with weight w prod and a decay event with weight w dec is given by where Br H→bb is the branching ratio of the Higgs boson to bottom quarks, which is treated as an input, while Γ H→bb is the H → bb decay width, which is computed from the decay events. Accordingly, the last factor in Eq. (2.5) integrates to one exactly if one is fully inclusive over the Higgs decay products. This implies that regardless of the accuracy of the last factor (MiNLO vs. MiNNLO PS ) the integrated inclusive cross section will be exactly equal to the cross section obtained from the production events times the branching ratio of the Higgs boson to bottom quarks. Similarly, the precise value of the bottom Yukawa coupling used in the Higgs decay cancels out exactly in the ratio of the last factor in eq. (2.5), since we evaluate the Yukawa coupling at a fixed scale m H . In order to match the final events to a parton shower one needs to take into account that the starting scale of the shower for production and decay is different. In our combined event file we store the starting scale (scalup) of the production stage, while, when interfacing to the shower we recompute on-the-fly the veto scale for the decay (scalup dec ) using the decay kinematics, as explained in section 3.2 of ref. [126].
In order to respect the radiation bound coming from the production stage we simply use the scalup variable stored in the combined event files, as usually done when interfacing Powheg with Pythia. In order to constrain shower radiation off the decay products, we implement a vetoed shower, i.e. we first allow Pythia to generate emissions off the Higgs decay products in all available phase space and, once the shower is completed, we check the hardness of the splittings that were generated in the decays. To do so, we go through all splittings that originated from the Higgs decay products generated by Powheg and, for each, we compute the corresponding hardness using where E rad (E em ) and p rad (p em ) are the energy and momentum of the radiating (emitted) particle, respectively. If for each of these splittings the hardness is smaller than scalup dec , we accept the showered event, otherwise we reject it and we attempt to shower the event again, until the above requirement is fulfilled. After 1000 shower attempts, we reject the event.
Finally, we note that when interfacing to the shower, we make sure that the automatic matrix element corrections to the decay of the parton shower program are switched off, since they are already included with better accuracy in our combined events.

Validation against NNLOPS results
V H production at NNLO+PS has already been calculated with the NNLOPS method, which combines the original MiNLO approach [67] with a reweighting procedure fully differential in the Born phase space to include the NNLO corrections. While such reweighting constitutes a substantial technical obstacle and requires an additional handling of the generated MiNLO events to include NNLO accuracy a posteriori, the physical results are sound. Precisely the complications related to the reweighting are solved by the newer MiNNLO PS method, which includes the NNLO corrections directly during the generation of the events.
NNLOPS results are available so far for W H production without [84] and for ZH production with H → bb decay [126]. It is therefore very useful to validate our MiNNLO PS calculation by comparing to those predictions, which will be done in the following. We will show that MiNNLO PS and NNLOPS results are fully compatible, so that our new MiNNLO PS V H generators supersede the previous NNLOPS ones.

W ± H production
We start by considering W ± H production keeping the Higgs boson stable, and we compare our MiNNLO PS predictions against the results of ref. [84]. Since only results for W + H are presented in ref. [84] we consider only this case.
We consider proton-proton collisions at the LHC with a center of mass energy of √ s = 13 TeV and present predictions for the process pp → e + ν e H. The input parameters are identical to the ones in ref. [84]. In particular, we use the NNLO set of the MMHT2014nnlo68cl parton densities corresponding to a value of α s (m Z ) = 0.118, and we set m W = 80.399 GeV, Γ W = 2.085 GeV and m H = 125.0 GeV. Furthermore we set α EM = 1/132.3489 and sin 2 θ W = 0.2226. The Higgs boson is considered to be on-shell and, in order to compare with ref. [84], contributions where the Higgs boson is radiated off a heavy quark loop are neglected. As already mentioned, these contributions have been shown to be of the order of 1% (see e.g. section I.5.2.c of ref. [146]). The central renormalization and factorization scales are set following the MiNNLO PS procedure. The perturbative uncertainties are obtained from a 7-point scale variation, i.e. by varying µ R and µ F around the central scale by a factor of two with the constraint 1/2 ≤ µ R /µ F ≤ 2. The phase-space definition is fully inclusive, without any fiducial cuts applied.
The total cross sections are 96.72(4) +1.9% −0.6% fb for MiNNLO PS and 96.69(3) +1.3% −1.3% fb for NNLOPS. Thus they agree at the sub-permille level, which is even within the numerical errors and far below the scale uncertainties. There is a small difference in scale uncertainties, which are of the same size, but slightly asymmetric in the case of MiNNLO PS . This is a consequence of the different settings of perturbative scales and the treatment of terms beyond accuracy in the two calculations, and it is therefore not unexpected. Figure 2 compares MiNNLO PS and NNLOPS results at the level of Les Houches events (LHEs) for various differential distributions. In particular, the Higgs rapidity (y H ) and transverse momentum (p T,H ), the rapidity (y W + ) and the transverse momentum (p T,W + ) of the W + boson, as well as the rapidity of the positron (η e ) and the transverse momentum of the W + H system (p T,W + H ) are shown. For all rapidity distributions the central predictions are in excellent agreement between MiNNLO PS and NNLOPS at the level of ∼ 1%, up to statistical fluctuations in the tails. As already observed for the total cross section, some slight differences in the scale-uncertainty bands are visible. Also for the transverse-momentum spectra of the Higgs and W + boson we find decent agreement between MiNNLO PS and NNLOPS. The central predictions differ by less than 1-2%, except for some statistical fluctuations. Finally, the differences in the transverse-momentum spectrum of the W + H system are slightly larger, reaching up to 3-4%. However, this distribution is effectively only NLO accurate at p T,W + H 0 and subject to matching ambiguities at intermediate p T,W + H values, which is indicated also by the enlarged uncertainty bands. MiNNLO PS and through a p T -dependent multi-differential rescaling in NNLOPS). In particular, as already noticed in ref. [77], the shape of the MiNNLO PS uncertainty band more closely resembles that of a resummed result matched to a fixed-order prediction.
We have considered many other observables, and we found the same level of agreement for all of them. In conclusion, with agreement typically at the level of 1 − 2% and at most at the level of a few percent, all within scale variations, the MiNNLO PS and NNLOPS calculations provide remarkably compatible results, and both calculations can be considered to be equivalent, up to small subleading corrections.

ZH production with H → bb decay
The situation for ZH production is very similar. First of all, the inclusive cross sections obtained with the new implementation for the production stage (i.e. before including the H → bb decay and without fiducial cuts) have been thoroughly validated against the completely independent fixed-order calculation implemented in MCFM. We find an agreement at the sub-percent level for the inclusive cross section: The central values are 26.567(1) fb for MCFM, 26.51(1) fb for MiNNLO PS using the Powheg-Box-Res implementation, and 26.63(9) fb using the Powheg-Box-V2 one, where the numbers in brackets refer to statistical uncertainties only. Given that scale uncertainties are at the one-percent level, and that the calculations differ by terms beyond accuracy, the agreement is remarkably good.
We continue by comparing MiNNLO PS predictions to the results of ref. [126], which correspond to the combination of pp → e + e − H production and H → bb decay at NNLO+PS. Our comparison is thus performed for e + e − bb events at the LHE level. As before, we consider √ s = 13 TeV proton-proton collisions at the LHC. The input parameters and fiducial cuts are identical to those in ref. [126]. In particular we use the PDF4LHC15_nnlo_mc parton distribution functions, corresponding to a value of α s (m Z ) = 0. 118 Table 1: Fiducial cuts for e + e − bb production from ref. [126]. We would like to stress, however, that given the way production and decay events are combined in eq. (2.5) the precise value of the Yukawa coupling has no impact, as already pointed out in section 2.4. Scale uncertainties are obtained through customary 7-point scale variations with the constraint 1/2 ≤ µ R /µ F ≤ 2, while correlating the scale variation factors in production and decay. The fiducial cuts are are summarized in table 1. If more than two bottom-flavoured jets (b-jets) are present in the final state, the pair with the invariant mass closer to the Higgs-boson mass is chosen.
The fiducial cross sections are 6.261(7) +0.9% −1.8% fb for MiNNLO PS and 6.348(6) +1.2% −1.4% fb for NNLOPS. Because of the more exclusive setup, the difference between the MiNNLO PS and NNLOPS predictions is slightly larger than for the inclusive cross section and for W + H production in the previous section. Still, the fiducial cross sections agree at the level of 1.4%, which is fully covered by the scale uncertainties. We would like to stress that this level of agreement in the fiducial phase-space volume of the e + e − bb final state is quite remarkable, given the very small scale dependence of this process.
As far as differential distributions are concerned, in figure 3 we consider the invariant mass of the ZH system (m ZH ) and observables related to the bottom quarks from the Higgs decay. In particular, the invariant mass of the b-jet pair (m bb ), its transverse momentum (p T,bb ), and the rapidity difference between the two b-jets (|∆y b,b |) are shown. We note that the two b here refer to the two b-jets with invariant mass closest to the Higgs mass, as discussed above. Also here we find a good agreement between MiNNLO PS and NNLOPS predictions for all observables within the scale-variation bands, with differences, by and large, lower than 1-2%. Besides statistical fluctuations, these differences are induced by the 1.4% shift in the fiducial cross section noticed before, while the agreement in terms of shapes is even better.

Phenomenological results
In this section we present phenomenological results for V H production with H → bb decay in proton-proton collisions at the LHC with a center of mass energy of √ s = 13 TeV. More precisely, we consider both MiNLO and MiNNLO PS results for the processes pp → e ± ν e bb and pp → e + e − bb, and compare our state-of-the-art MiNNLO PS predictions to recent ATLAS data [131]. Besides integrated inclusive and fiducial cross sections we also discuss differential distributions in the fiducial phase space, and we analyze the effects of using different jet-clustering algorithms.

Input parameters and settings
We use the NNLO set of the NNPDF31 [148] parton densities with α s (m Z )=0. 118  The bottom Yukawa coupling is set to y b (m H ) = 1.280 ×10 −2 , which, however, as pointed out in section 2.4, cancels out in the ration of eq. (2.5) when combining production and decay events. We employ Pythia8 [132] with the Monash tune [150] for all matched predictions. We do not include any effect from hadronization, underlying event modelling or a QED shower, except for the comparison to data in section 4.5.
Besides a fully inclusive setup, which we refer to as inclusive, we consider two sets of fiducial cuts. The first one is inspired by the CERN Yellow Report [146] and is refereed to as fiducial-YR. The second set of cuts is taken from ref. [131] and is referred to as fiducial-ATLAS. Both sets of fiducial cuts are summarized in table 2.

Inclusive and fiducial cross sections
We start by discussing integrated cross section both in the fully inclusive case and including fiducial cuts. In table 3 we report MiNLO and MiNNLO PS prediction for both the neutral and the charged-current V H production processes with H → bb decay. For ZH production we report separately the numbers with and without the gluon-fusion contribution that is induced by a closed heavy-quark loop (top and bottom). The inclusive cross sections are in good agreement with the NNLO predictions from vh@nnlo refs. [44][45][46]152] that are obtained in the on-shell approximation with µ R = µ F = m V H . Assuming branching ratios consistent with our calculation (Br H→bb = 0.5824, Br W ± →e ± νe = 0.10898 and Br Z→e + e − = 0.033628), they read σ W + H→e + νebb = 58 fb, σ no gg→ZH ZH→e + e − bb = 16.0 +0.2% −0.3% fb, and σ with gg→ZH ZH→e + e − bb = 17.0 +1.5% −1.2% fb. 1 We note that the inclusive on-shell calculation differs from the MiNNLO PS also with respect to the scale settings and the different treatment of terms beyond accuracy.  Table 3: Integrated cross sections for the e + ν e bb, e −ν e bb and e + e − bb final states in the inclusive and the fiducial-YR setup.
Furthermore, we find in table 3 that MiNNLO PS induces a correction of about 5-6% with respect to MiNLO , both in the fully inclusive and in the fiducial phase space volume, disregarding the loop-induced gg contribution in the case of ZH production. While for the inclusive cross sections MiNLO and MiNNLO PS predictions are compatible within scale variations, the MiNNLO PS corrections in the fiducial phase space are generally not covered by plain scale variations when correlating the scales in production and decay. If we vary the scales in production and decay in an uncorrelated manner, the unnaturally small MiNLO uncertainty band in the fiducial setup increases and the predictions become compatible within scale variations. The scale uncertainties of MiNNLO PS are significantly smaller than the ones of MiNLO in either case, with a reduction by more than a factor of two. Considering the loop-induced gg → ZH contribution, its effect is quite significant. It contributes about almost 8% in the inclusive and almost 11% in the fiducial phase space to the MiNNLO PS cross sections. As pointed out before, the contribution from loop-induced gg diagrams dominates the current theoretical uncertainties. Thus, including NLO QCD corrections to the gg → ZH process in the matching to parton showers is a crucial next step. 2

Distributions in the fiducial phase space
We continue our presentation of phenomenological results by considering differential distributions in the fiducial-YR phase space defined in table 2. Figures 4 and 5 show various observables for e + ν e bb (W + H) and e + e − bb (ZH) production, respectively. 3 In particular, the invariant mass (m bb ) and transverse momentum (p T,bb ) of the bottomquark pair, the absolute rapidity difference between the two bottom quarks (|∆y b,b |), the transverse momentum of the hardest bottom quark (p T,b hard ), the transverse momentum of the positron/electron (p T,e + /p T,e − ), and the invariant mass of the colour-singlet system in the final state (m W H /m ZH ) are shown. We recall that with the two bottom quarks (from the Higgs decay) we actually denote the two jets with at least one bottom quark whose invariant mass is closest to the Higgs-boson mass. The plots are organized as follows: The main frame includes the full MiNNLO PS prediction (blue, solid), the full MiNLO prediction (black, dotted), and the MiNNLO PS prediction with LO decay through Pythia8 (red, dashed) as absolute cross sections of the differential distributions. The lower two insets show the ratio to the MiNNLO PS central prediction, where in the first inset scale uncertainties are evaluated by assuming that the 7-point scale dependence in production and decay is fully correlated, while in the second inset uncorrelated scale variations in production and decay are performed, but excluding the extreme variations (one scale up, the other down). When the Higgs decay is simulated by Pythia8, the scale variation is also handled by Pythia8 using its so-called "automated parton-shower variation" facility [153]. In particular, we vary the renormalisation scale for final-state radiation by a factor 2 around the central value, correlating this variation with the one performed in the production stage. The relative behaviour observed in both the W ± H and the ZH results is very similar. Thus, we can discuss the sets of plots simultaneously.
Looking at the m bb distribution, we notice that MiNLO and MiNNLO PS predictions are very similar as far as the shape is concerned, while the MiNNLO PS corrections increase the cross section by roughly 5%, as already observed for the integrated fiducial cross section in the previous section. MiNLO and MiNNLO PS predictions are fully compatible within scale variations, regardless of whether scale variations are assumed to be correlated or uncorrelated. Interestingly, for this observable the two cases lead to very similar uncertainties, which can be seen by comparing the first and the second ratio inset. Moreover, we observe clearly reduced scale uncertainties for MiNNLO PS compared to MiNLO below the Higgs-mass threshold by roughly a factor of two, but similarly large uncertainties above the threshold. Using Pythia8 to include the H → bb decay at LO compares rather well to using the full MiNNLO PS prediction for production and decay.
For the rapidity difference of the two bottom quarks MiNLO and MiNNLO PS predictions are not compatible at central rapidities when assuming correlated uncertainties. In particular, the MiNNLO PS result features a very small band at the ±1% level. By uncorrelating the scale variation in production and decay more realistic scale bands are obtained, leading to fully compatible MiNLO and MiNNLO PS predictions, where the uncertainties of the latter increase only to a very few percent. This behaviour can be observed also for  Figure 4: Differential distributions for e + ν e bb production with fiducial-YR cuts. K R,prod and K R,dec refer to the variation factors of the renormalization scales of production and decay, respectively.  Figure 5: Differential distributions for e + e − bb production with fiducial-YR cuts.K R,prod and K R,dec refer to the variation factors of the renormalization scales of production and decay, respectively. the remaining distributions, where it is visible almost over the entire plotted ranges. This justifies a more conservative uncertainty estimate by taking uncorrelated scale variations in production and decay. We also observe that the LO H → bb decay through Pythia8 is insufficient to describe this observable accurately, especially in the forward ∆y b,b region. Moreover, in the central rapidity region, the uncertainty bands of the Pythia8-decayed results are not compatible with the MiNNLO PS ones. Considering the other differential observables it is clear that similar features appear. For all of them MiNNLO PS induces a mostly flat ∼ 5% correction on top of MiNLO and leads to a significant reduction of the scale uncertainties. Note that for correlated scale uncertainties some accidental cancellations in the scale variations occur, which are particularly visible from the MiNLO scale-uncertainty band around 200 GeV in both the p T,bb and p T,b hard distributions. The red curve that includes the H → bb decay through Pythia8 at LO is in shape rather similar to the full MiNNLO PS result, but it is roughly 5% above the full MiNNLO PS prediction and thereby outside the quoted scale uncertainties in various phase-space regions.

Impact of jet-clustering algorithms
In this section we address the impact of different jet-clustering algorithms in the definition of b-jets on differential distributions. Our calculations are based on matrix elements with massless bottom quarks in the final state. Therefore, performing an event simulation is particularly important as it allows us to include a finite bottom-quark mass in the generation of Les Houches events or through Pythia when showering the events. This makes it possible to perform any standard jet clustering to define b-jets, which in a massless calculation would otherwise lead to infrared (IR) unsafe results. For massless b-jets, for instance in a fixedorder calculation, it is instead possible to use the flavour-k T algorithm of ref. [147] for an IR-safe definition of b-jets. Figure 6 shows the m bb , p T,bb , p T,b hard , and p T,e + distributions for W + H production in the fiducial-YR phase-space defined in table 2. The different curves correspond to using the anti-k T [151] (blue, solid), the flavour-k T [147] (green, long-dashed), and the k T [154] (red, short-dashed) jet-clustering algorithm, all with R = 0.4. In the bottom panels we show the ratio to the anti-k T result. It is easy to see that the anti-k T and k T algorithms lead to very compatible results with small differences only in certain phase-space regimes, in particular for the m bb distribution. There is a 5-10% difference between anti-k T and k T clustering both at small m bb and at the m H threshold. Such differences are not unexpected and they originate from soft gluons being more likely to be clustered in the b-jets with the k T algorithm.
As far as the flavour-k T algorithm is concerned, the differences to the other two approaches are much more prominent in some regions of phase-space. The flavour-k T algorithm assumes b-jets to originate only from jets that contain an uneven number of bottom quarks, since an even number of bottom quarks in one jet typically originates from g → bb splittings, which in the flavour-k T algorithm shall not be considered as b-jets. The reason for adopting such treatment stems mostly from fixed-order calculations in the massless approximation of bottom quarks, i.e. in the five-flavour scheme (5FS), where the stan-  dard clustering algorithms lead to definitions of bottom-quark jets that are not IR safe, as discussed before. Clear differences of the flavour-k T results to the non-flavour specific clustering algorithms are observed at small m bb and at large values of p T,bb , p T,b hard , and p T,e + . Additionally, around the m bb ∼ m H threshold the flavour-k T curve follows the one of the standard k T algorithm and shows a similar difference to the anti-k T one, as both the flavour and the standard k T algorithm use the same measure in the ordering of the jet clustering. In order to explain these differences we focuse on a comparison between k T and flavour-k T algorithms since the underlying distance measure is more closely related. We remind the reader how the flavour-k T algorithm works. The clustering procedure is the same as any sequential algorithm but the distance measure is modified as follows: where softer means with lower p T . The distance with the beam is defined separating explicitly the beam moving towards positive rapidities B (right-moving beam) and the opposite beam towards negative rapiditiesB (left-moving beam), where p T,B (p T,B ) represents the hard transverse scale associated with the right-moving (left-moving) beam, defined as follows: If the smallest distance is a d ij , then the two (pseudo) particles are clustered together into a (pseudo) particle whose flavour is given by the sum of the two individual flavours. In particular, when a bottom-flavoured particle and anti-particle end up in the same jet, they will give rise to a non-flavoured jet. If the smallest distance is a d ib , the (pseudo) particle is considered as a final jet and removed from the clustering sequence. The large discrepancy between flavour-k T and k T results at small m bb , which is dominated by small ∆R bb values, is due to a different treatment of particular kinematic regions in the two algorithms. Specifically, when two b-quarks are separated by ∆R bb > R, the k T algorithm typically reconstructs two b-jets. On the other hand, the flavour-k T algorithm compares the relative distance d ij with the distance of each bottom quark to the beam d . The latter quantity is by constuction larger than the individual transverse momenta of the bottom and anti-bottom quarks. Since d ij can be smaller than d bB ) even for ∆R bb > R, the two bottom quarks can be clustered even if ∆R bb > R. Similar arguments apply also in the region of large transverse momenta, since when the Higgs boson is highly boosted, its decay products tend to be collimated.
Apart from that, we note that even for bottom quarks with ∆R bb < R the output of the two clustering algorithms can be different: in this case, both algorithms tend to cluster the two bottom quarks together, giving rise to a b-jet in the k T algortihm and a gluon-jet in the flavour-k T one. It is thus more likely that the event passes the fiducial cuts in the k T case, because of the presence of one additional secondary bottom quark. In this context a soft b is often sufficient as it can be promoted to a hard b-jet by clustering with a hard gluon. This clustering mechanism is instead strongly suppressed in the flavour-k T one. In this case the anti-k T and k T algorithms do not actually reconstruct the two hard b-jets originating from the Higgs-boson decay.
Overall, it is clear that in various phase-space regions the flavour-k T clustering leads to significant changes. It is therefore not advisable to use this approach for theoretical predictions, unless the same clustering is used also on the experimental side, which is usually not the case. In fact, the application of the flavour-k T algorithm to experimental data is challenging, as discussed in ref. [155]. Alternatively, one should include a theoretical estimate of the size of the discrepancy due to the use of different algorithms in the theory simulations and measurements. While in fixed-order calculations in the 5FS such approach is mandatory to obtain IR-safe results as far as bottom-quark tagging is concerned, in exclusive event simulations matched to parton showers this problem can be circumvented by including the bottom-quark mass through reshuffling the final-state momenta of the bottom quarks, so that flavour-k T clustering can be avoided. This is an example why matching high accuracy calculations at fixed order and parton showers is important and advantageous. One should bear in mind, however, that such reshuffling is not completely unambiguous and comes with some uncertainties on the bottom-quark kinematics. Thus, ideally the bottom-quark kinematics should be described using massive bottom quarks at amplitude level , i.e. a four-flavour scheme (4FS) calculation, but this is not always feasible at high accuracy with current technology and also comes with other shortcomings [156]. Furthermore, it is important to be aware that in certain configurations a hard reconstructed b-jet can come from a soft bottom quark and that the two selected b-jets may actually not originate from the Higgs-boson decay. These occurrences are less likely when using the flavour-k T algorithm.  [131]. The respective cross sections are reported in table 4. The results correspond to W ± H and ZH production with all leptonic final states, i.e. ± ν bb, + − bb and ν ν bb with ∈ {e, µ, τ }. For this comparison we have included hadronization effects through Pythia8, which however have a negligible impact. It is clear that the measured V H cross sections are fully compatible with our predictions within uncertainties. However, one must bear in mind that this measurement requires relatively large lower cuts on the transverse momentum of the vector boson. Therefore, the experimental error is quite large, being dominated mostly by the limited statistics. In fact, within the current experimental uncertainties the measured cross sections are compatible with zero within one standard deviation. It will be interesting to see how this comparison will develop when the experiments have collected more data. Moreover, the original publications [3,4], where the H → bb decay was observed, involved a cut of only 150 GeV on the transverse momentum of the vector boson, but no unfolding to fiducial cross sections was performed. Therefore, repeating the exclusive data-theory comparison at the fiducial level including a cut of only 150 GeV on vector-boson transverse momentum would be desirable in the future.

Conclusions
We have presented novel predictions for V H production with H → bb decay within the MiNNLO PS framework. NNLO QCD accuracy is retained for both production and decay, while the matching to the Pythia8 parton shower ensures a fully realistic exclusive description of the process at the level of hadronic events. For the W ± H process NNLO+PS accuracy is achieved for the first time for production and decay. Our results have been validated against earlier calculations using the NNLOPS approach, which exploits a numerically intensive multidifferential reweighting to upgrade MiNLO events to NNLO QCD accuracy a posteriori. We find fully compatible results within scale uncertainties when comparing our MiNNLO PS predictions to earlier NNLOPS results for W ± H production with an on-shell Higgs boson and for ZH production including the H → bb at NNLO+PS. The agreement is particularly remarkable considering the fact that the scale uncertainties of both predictions are at the level of only 1-2% in most phasespace regions. Given that the reweighting procedure of the NNLOPS approach induces some (technical) limitations and for V H production it relies on certain approximations in describing the vector-boson decay, our new MiNNLO PS V H generator, which includes NNLO corrections on-the-fly directly during the usual Powheg event generation, lifts those shortcomings and can be considered to supersede the previous NNLOPS implementation. We stress, however, that the NNLOPS predictions can still be useful to assess residual uncertainties.
Phenomenological results have been discussed in detail for both the pp → e ± ν e bb and pp → e + e − bb processes at the level of inclusive and fiducial cross section as well as for differential distributions in the fiducial phase space. We find that the NNLO corrections induced by the MiNNLO PS procedure with respect to MiNLO increase both the inclusive and fiducial cross sections by roughly 5−6%, while the corrections are relatively flat in phase space for the observables under consideration. Moreover, due to the relatively small scale dependence of the V H processes, uncorrelating scale variations in production and decay to estimate the perturbative uncertainties is crucial to obtain MiNLO and MiNNLO PS results that are compatible within their scale bands. Regardless of this, we observe a clear reduction of scale uncertainties, by roughly a factor of two, when including the MiNNLO PS corrections. As far as the H → bb decay is concerned, we have shown that including the full NNLO QCD corrections induces important effects with respect to a LO description through Pythia8 for various observables: For the cuts under consideration, they shift the fiducial cross section by about −5% and induce relevant shape effects. In some cases, these effects are not captured by the scale-variation uncertainty from Pythia8. Finally, our results are in agreement with current ATLAS data, keeping in mind that with the current experimental accuracy the measured V H cross section is still compatible with zero.
While at the level of the matrix elements our calculations involve massless bottom quarks, kinematic bottom-mass effects are included approximately through a reshuffling of the bottom quarks at LHE level. In this setup, we have compared the anti-k T , k T and flavour-k T algorithms. The anti-k T and k T algorithms give rather similar results, and in the bulk region of the phase space this is also the case for the flavour-k T algorthim. However, in some phase-space regions, such as at small m bb and at large transverse momenta, we observe substantial differences in the flavour-k T case. Since experimental analyses usually employ the standard (anti-)k T algorithm, it is important to be aware of these differences and features, especially when using different jet-clustering algorithms in the theoretical calculations and in experimental measurements.
We have performed largely independent implementations in both Powheg-Box-V2 and Powheg-Box-Res, whose results are fully consistent. We reckon that our MiNNLO PS V H generators will be very useful for future analyses of these processes in Run III and, in particular, for the High Luminosity LHC phase. Therefore, we will make all relevant codes to simulate NNLO+PS events for the e ± ν e bb and e + e − bb final states publicly available within the Powheg-Box-Res framework.   Figure 8: Differential distributions for e −ν bb production with fiducial-YR cuts. K R,prod and K R,dec refer to the variation factors of the renormalization scales of production and decay, respectively.