NNLOPS description of the H$\to\! b\bar{b}$ decay with MINLO

We present an event generator that describes the Higgs boson decay into $b$-quarks at next-to-next-to-leading-order (NNLO) in QCD and allows for a consistent matching to parton shower. For this purpose, we work within the POWHEG framework and employ the MINLO method to produce an NNLO-accurate H$\to\! b\bar{b}$ event sample which can further be interfaced with a generic Higgs boson production mode. The resulting events, that combine Higgs production and decay, can be matched to parton shower by means of a vetoed shower. As an example, we present a short phenomenological study for associated Higgs production (HZ) and vector boson fusion (VBF) processes. We release the source code that was used to obtain these results as an $\text{h_bbg}$ user process in POWHEG-BOX-V2. We also provide tools that enable interfacing the decay events with events for an arbitrary Higgs production mode and subsequent matching of the combined Les Houches (LH) events with a parton shower.


Introduction
Since the discovery of the Higgs boson by ATLAS and CMS in 2012, the study of Higgs properties became a central part of the LHC physics programme.The mass of the Higgs boson is already measured with an accuracy of about 0.2% [1], hence, within the Standard model (SM), all couplings of the Higgs to other particles are predicted.Nevertheless, beyond the Standard Model (BSM) effects can modify these couplings, therefore constraining Higgs properties becomes one of the priorities of the LHC physics programme.For this reason, increasing experimental and theoretical precision of measurements and predictions for production and decay rates is crucial to enhance the sensitivity of searches for New Physics.
The largest branching fraction of the Higgs boson is to b-quarks.Studies of very rare Higgs production modes obviously benefits from the large H → b b branching fraction the most.For instance, in double-Higgs production, which is the crucial process to constrain the Higgs self-coupling directly, the most promising decay modes (4b, 2b + 2τ , 2b + 2γ) involve at least one Higgs boson decaying into b-quarks.This is true both for the High-Luminosity (HL) and High-Energy (HE) LHC [2] programmes.
The H → b b decay mode is notably hard to study experimentally because b-quarks typically give rise to b-jets which are more complicated final states compared to decay modes involving charged leptons.Furthermore, extracting the H → b b signal is challenging because of large QCD backgrounds, coming from gluons splitting to b-quarks and/or from other SM processes.The sensitivity can be enhanced by using the resonant kinematics of the signal, by requiring b-jets with high transverse momentum and by exploiting different properties of jets originating from a boosted colour singlet or a boosted gluon.A lot of progress was achieved in this direction in recent years and refined techniques were suggested and are currently used in experimental measurements, see e.g.refs.[3,4].Both ATLAS and CMS reported evidence for the H → b b decay [5,6] followed by an observation [7,8].Significant increase in statistics of the upcoming Run III and the HL-LHC programme will allow for a comprehensive use of boosted techniques and considerably improve the significance of these measurements.
On the theoretical side, an accurate description of this decay is also challenging since higher-order QCD corrections are very large.The NLO QCD corrections to the Higgs decay to massive b-quarks have been known for a long time [9][10][11][12][13][14].NNLO corrections for massless b-quarks were computed more recently [15,16] with even first N 3 LO results appearing [17,18].Today also NNLO corrections to the decay to massive b-quarks are known [19,20].Ref. [21] also presents the exact computation of the top-Yukawa correction to H → b b decay, which appears first at NNLO.
The main Higgs production processes, gluon-fusion (ggf), Vector boson Fusion (VBF) and associated HW/HZ production are known to NNLO accuracy, or better, see e.g.Ref. [22].Other Higgs production modes are also expected to be known to NNLO accuracy in the upcoming years.In view of this, combining the best accuracy in production and decay seems natural.This was done in specific cases, e.g. for HW/HZ in refs.[21,[23][24][25].It is well-known that corrections to production and decay can be discussed separately to a very good approximation.In fact, since the Higgs boson is a scalar, there are no spin-correlation effects.Furthermore, the fact that it is a colour singlet implies that nonfactorizable corrections between production and decay vanish exactly at NLO.At NNLO, the exchange of two gluons between production and decay is possible, but these contributions are O(α 2 s Γ H /M H ) for sufficiently inclusive observables, i.e. for observables that are not sensitive to soft gluon emissions [26,27].
The description of the H → b b decay at NNLO matched to a parton shower is the main focus of this paper.While NNLO corrections to the decay were already studied, higherorder parton shower effects can still have large effects on kinematic distributions of phenomenological relevance, in particular when exclusive fiducial cuts are applied.Moreover, including NNLO corrections into a fully exclusive event generator is particularly important for analysis performed by the experimental collaborations.
In this work we use the MiNLO method [28,29] to build a generator (henceforth named Hbbg-MiNLO) that produces fully exclusive events for the H → b b decay.Subsequently the accuracy of these events is upgraded to NNLO by means of a reweighting procedure.Due to the considerations made above, it is possible to interface such decay events with any Higgs boson production process and further match them to a parton shower.We explain how we implement this combination of production and decay events, and what restrictions need to be imposed on a parton shower evolution such that the fixed-order accuracy is preserved.
The paper is organised as follows.In Sec. 2 we focus on the decay of the Higgs to b quarks and we first explain how to achieve predictions that are NLO accurate both for H → b b and H → b bg observables using MiNLO.Since the environment is particularly simple in this case, we additionally present a numerical verification of the MiNLO accuracy by performing a check in the α s → 0 limit.We then detail the reweighting procedure to obtain NNLO accuracy for the H → b b decay width.In Sec. 3 we describe the combination of Hbbg-MiNLO generated decay events with a generic Higgs production process and the subsequent combination with a parton shower.In Sec. 4 we show two explicit examples where we combine NNLO H → b b decay events with the HZ Higgs production process, as well as with the VBF Higgs production process.We discuss the advantages of the approach presented in this article, and, for HZ production, we also compare our results with those presented in Ref. [30].We summarise our findings in Sec. 5.

H → b b at NNLOPS accuracy
In this section we describe a construction of an event generator that describes the Higgs decay into b-quarks at NNLO QCD and allows for inclusion of parton shower (PS) effects.This is achieved by first merging the H → b b and H → b bg generators using the MiNLO procedure and then performing a reweighting of the Hbbg-MiNLO events to the exact NNLO H → b b partial width.
We focus here on the steps necessary to reach NNLO QCD accuracy.As a prerequisite we extended POWHEG to incorporate a process without initial-state radiation, where QCD partons are present only in the final state.We leave the discussion of the combination of NNLO-reweighted events with a generic Higgs production mode, as well as interface with PS for later sections.

Merging H → b b and H → b bg generators using MiNLO
We start by considering the decay of the Higgs boson into a pair of massless b-quarks accompanied by an additional gluon, H(p 1 ) −→ b(q 1 ) + b(q 2 ) + g(q 3 ). (2.1) Our first goal is to construct an event generator that treats the H → b bg processes at NLO accuracy and at the same time allows for integrating out the extra gluon, yielding an NLO accurate description of the H → b b decay channel.We follow the usual MiNLO procedure [28,29] that enables us to simulate the H → b b and H → b bg processes with a single event generator without introducing an additional merging scale.The method uses an NLO fixed-order calculation of the H → b bg process together with information encoded in the Sudakov form factor of a resummed prediction for some infrared safe observable that vanishes for Born level H → b b events.
In this case, we use a resummed prediction for the three-jet resolution parameter y 3 .The y 3 observable corresponds to a value of a y cut parameter within a clustering algorithm that separates between two-and three-jet configurations or, equivalently, a maximal value of y cut for which two jets are obtained.For reasons that will be explained in this subsection, we use the Cambridge algorithm [31,32] as a clustering algorithm for the H → b b(g) merging procedure.Its definition relies on two variables where M 2 H is the squared centre-of-mass energy, i.e. the mass of the on-shell Higgs boson, E i , E j and θ ij denote energies and the angle between particles i and j in the Higgs boson rest-frame.During the clustering sequence, pairs of particles (i, j) are sorted according to the ordering variable v ij , and the pair with a minimal value of v ij is selected.The procedure checks whether a test variable y ij is smaller than a parameter y cut , and if this is the case the particles i and j are recombined into a pseudo-particle and the algorithm repeats the procedure with all remaining (pseudo-)particles.On the other hand, if y ij > y cut , the softer of the two pseudo-particles is considered a jet and removed from the list of pseudo-particles.The clustering procedure stops with all remaining (pseudo-)particles declared as jets.
The resolution parameter y 3 classifies the hardness of events, i.e.H → b b events with additional hard, large angle gluons yield a large value of y 3 , close to the kinematical upper bound of 1/3, while events with soft or collinear gluons have y 3 1, which leads to large logarithms L ≡ − ln y 3 in the differential cross section.Events with H → b b Born kinematics have y 3 = 0.
In order to merge the H → b b and H → b bg calculations with the MiNLO method, we need a resummed calculation for y 3 that includes all terms of order α 2 s L [29].A next-tonext-to-leading-logarithmic (NNLL) accurate resummation of the y 3 observable has been performed in Ref. [33] and we use this result to extract the required Sudakov form factor ∆(y 3 ).We decided to use the Cambridge algorithm as it is the one for which the NNLL resummation for y 3 has the simplest structure.To reconstruct the Sudakov form factor, we use the formulae given in App.B of Ref. [34], setting a = 2, b = 0 for the pure radiator function, and use the required NNLL multiple emission corrections outlined in Eq. (4) of Ref. [33].Further details about the implementation of the NNLL ingredients within the MiNLO formalism are given in App. A.
Given the Sudakov form factor, the usual MiNLO formula for the POWHEG B function [35] has the form where Φ b bg is the phase space of the three-body decay, B H→b bg , V H→b bg and R H→b bg denote the Born, virtual and real matrix elements respectively, and ∆ (1) (y 3 ) is the O(α s ) expansion of the Sudakov form factor ∆(y 3 ).R, V and ∆ (1) contain an additional power of α s , which is also evaluated at q 2 t = y 3 M 2 H .1 As shown in Ref. [29], the integration over the radiation phase space leads to an NLO accurate description of H → b b observables.In particular, we obtain an NLO accurate result for the H → b b partial width A proof for the specific case at hand is given in App.B. We stress that this particular process offers a unique opportunity to verify the absence of spurious O(α 3/2 s ) terms in a simple way.Since no initial state partons are involved and we do not need to rely on external input such as parton distribution functions, we can easily consider the limit where α s → 0 and check that the result has the correct scaling, i.e. that the difference between Γ MiNLO H→b b and the nominal NLO result Γ NLO H→b b is of O(α 2 s ).We performed such a numerical check by manually changing the reference value of the strong coupling, α s (M Z ), at the level of POWHEG input card. 2 For the evaluation of the strong coupling at other scales, such as the renormalisation scale µ R or the resolution scale q t , we use an internal POWHEG routine for the running of the strong coupling implemented with the standard β function.
We present our results of the numerical check in Fig. 1.We show the difference between the Hbbg-MiNLO result and the pure NLO Higgs partial width to b-quarks.We also normalise the difference to the LO width and further divide it by a factor of α 2 s .We write where α s is meant to be α s (M H ). We study the behaviour of δ(α s ) as a function of the strong coupling.
If spurious O(α 3/2 s ) terms were present in the MiNLO result, the difference defined in Eq. (2.5) would feature an increasing behaviour when approaching the α s → 0 limit.On the contrary, we see that all curves in Fig. 1 flatten out in the limit of small α s .The three curves correspond to different values of the resummation scale Q introduced in App.A

Reweighting
In order to achieve NNLO accuracy for H → b b-observables, we perform a reweighting of Hbbg-MiNLO events with a reweighting factor W(Φ b b), where Φ b b is the phase space of the underlying H → b b process.In this particular case, because of the scalar nature of the Higgs boson, the decay is isotropic in the Higgs boson rest frame and therefore W(Φ b b) is just a number given by the constant ratio (2.6) In order for the NNLO reweighting not to affect hard events, which are already described at NLO accuracy by MiNLO, we modify the reweighting factor as follows, cf.[29], where and the function h(y 3 ), is used to classify the hardness of an event, i.e. h(y 3 ) approaches one when the b b-pair is accompanied by soft or collinear emissions, while it should satisfy h(y 3 ) 1 for hard events.In order to achieve these constraints, we choose y 3,ref = e −4 and y 3,pow = 2.

Production-decay interface and parton shower matching
In this section we describe how to combine the NNLO-reweighted events produced with the Hbbg-MiNLO code with a generic Higgs production process.We also discuss how to deal with small off-shell effects.Furthermore, we briefly describe how such hybrid events can be consistently interfaced with a parton shower keeping the fixed-order accuracy of the event sample, both in production and decay.Detailed instructions and practical aspects are given in the manual accompanying the public code.
We stress that the reweighting to NNLO, both in production and/or decay, is optional, and if these steps are skipped one still obtains NLOPS accuracy in production and decay.For instance, one might consider including NNLOPS in decay, which is numerically feasible, while keeping only NLOPS accuracy in the production, in case the multi-dimensional reweighting to get NNLOPS accuracy is numerically very intensive or even unavailable, e.g. for VBF or ttH processes.

Combination of a generic Higgs production events with the H → b b events
The starting point is to assume that one has access to a set of Les Houches events produced with POWHEG-BOX-V2 with an undecayed Higgs boson and possibly other final state particles.Furthermore, one also has available a set of events with a Higgs-boson decay into b-quarks plus additional QCD radiation prepared using the Hbbg-MiNLO code, as described in the previous section.In the following we describe the sequence of steps needed to obtain a combined event that involves the production and the decay of a Higgs boson in a melded event.
Momenta: Let us consider a single event record event_Hprod with a Higgs boson in the final state and possibly other final-state particles.We denote the Higgs boson momentum by p H . On the other hand we have an event record event_Hdec that describes a decay of an on-shell Higgs boson, in its rest frame, to b-quarks accompanied by additional radiation.From these two sets of momenta we generate the full event as follows: • We pull out the momenta of the Higgs decay products q i from the event_Hdec record.
Note that momenta q i are specified in the Higgs boson rest-frame.
• From the momentum of the Higgs boson in the production stage, p H , we compute a factor λ = p 2 H /M H .For an on-shell Higgs boson λ = 1, but if the Higgs boson in the event_Hprod was generated with its virtuality distributed according to a Breit-Wigner shape then λ = 1.
• In case λ = 1, we reshuffle the q i momenta such that the energy-momentum conservation is restored, i.e. q i = S(λ)q i , such that i q i = λM H .3 • We boost momenta q i from the Higgs boson rest frame to the laboratory frame, • The final event record event_Hprod_x_dec is obtained by copying momenta and colour connections of all particles from the event_Hprod, changing the status of the Higgs boson to be an intermediate resonance (istup = 2), and appending the p i momenta of the decay products with their colour connections at the end of the list.
Note that such a treatment of small off-shell effect neglects terms suppressed by (Γ H /M H ).
Weights: Each production and decay event comes with its own set of weights due to scale variations, i.e. {w prod,i } and {w dec,j }.The weights of the final event_Hprod_x_dec are obtained by multiplying w prod,i × w dec,j .The most general set of weights is obtained by keeping all possible (i, j) combinations.Nevertheless, simpler options are possible to consider, for instance, by correlating the production and decay scale variations, as done for the results presented in this paper.
PS radiation bound: As usual, POWHEG-BOX events come with an upper bound for additional radiation that is generated during parton shower evolution, the scalup variable.
In the final event_Hprod_x_dec record we keep the original scalup from the production stage, i.e. scalup_prod, while the veto scale for the decay, scalup_dec, can be easily reconstructed using the decay kinematics.We leave the discussion of further details to the next subsection.

Interface to parton shower
A consistent matching of an (N)NLO-accurate event generator with a parton shower requires a special treatment of the parton shower evolution.Within the POWHEG-BOX framework, the hardest emission has already occurred at the level of event generation and each event (event_Hprod as well as event_Hdec, in our case) is equipped with an upper bound for the subsequent radiation generated by a parton shower.This quantity is commonly known as the scalup variable.
As explained in the previous subsection, the full event record, event_Hprod_x_dec, contains only the scalup_prod value.Nonetheless the bound scalup_dec related to the H → b b decay stage can be reconstructed from the decay kinematics using the POWHEG definition of hardness in final-state radiation [36].In particular, we identify the hardest splitting that occurred in the decay using the information contained in the colour connections.We denote by p rad the radiated parton and by p em the parton that has emitted it, and then calculate the hardness of the splitting t as where E rad and E em are the energies of the two partons in the Higgs boson rest frame.We finally set scalup_dec to be t, whereas if the decay event has no additional radiation besides the "Born" gluon in Φ bbg , we set scalup_dec to be an infrared cutoff of O(Λ QCD ).Such a treatment of the full event event_Hprod_x_dec provides a hardest emission for the production stage and one for the decay stage.This setup corresponds to the allrad option described in App.A of Ref. [36], hence, we follow the procedure introduced in Ref. [36] to shower event_Hprod_x_dec events with the Pythia8 parton shower.More specifically, in order to respect the radiation bound coming from the production stage we simply pass the scalup variable, equal to scalup_prod, to the shower program, as is always done when interfacing POWHEG with Pythia.In order to constrain the radiation emitted off the decay products through showering, we implement a vetoed shower, i.e. we let Pythia8 generate emissions off the decay products in all the accessible phase space and, after the shower is completed, we check the hardness of the splittings that were generated.First, we find the splittings that originated from all the Higgs decay products generated by POWHEG, then, for each of them, we calculate the corresponding hardness hardness_dec using Eq.(3.1).If for each of those splittings, the corresponding hardness is smaller than the veto scale scalup_dec, we accept the shower history; otherwise we reject it and we attempt to shower the event again, until the condition hardness_dec < scalup_dec is met.A very similar result can be obtained by using built-in facilities of Pythia8, namely by using the UserHooks class [37,38].

Practical implementation and phenomenological studies
In this section we present two concrete applications of the method described in this article, which have NNLOPS accuracy in the H → b b decay.We first consider the associated Higgs production process pp → HZ(→ ) described here at NNLOPS accuracy both in production and decay, and compare results obtained here with those of Ref. [30], which have the same NNLOPS accuracy for production but only NLOPS accuracy for the decay. 4The two results also differ in the treatment of the interface to the parton shower in the decay, as explained in detail below.Next, we consider vector boson fusion Higgs production, described at NLOPS accuracy for the production stage.We interface the production with NNLOPS accurate H → b b decay and we compare these results with those obtained by letting Pythia8 handle the Higgs decay.

Associated Higgs production
Input parameters and fiducial cuts: As far as the Higgs production process is concerned, we used a setup identical to Ref. [30] which facilitates a direct comparison to those results.In particular, we consider 13 TeV LHC collisions and use PDF4LHC15_nnlo_mc parton distribution functions [40][41][42][43] For the decay, our implementation of Hbbg-MiNLO uses the matrix elements from Ref. [16].The bottom Yukawa coupling in H → b b decay is evaluated in the MS scheme at the decay renormalisation scale µ R,dec = M H .The running bottom Yukawa coupling is computed from the on-shell Yukawa coupling using an O(α 2 s ) conversion that is implemented along the lines of RunDec package [44,45].The numerical value of the bottom quark MS mass at the renormalisation scale is m b (M H ) = 3.152 GeV and the corresponding Yukawa coupling is y b (M H ) = 1.280 • 10 −2 .Note that at variance to the Ref. [30] the b-quarks are treated as massless inside the Hbbg-MiNLO generator.
The production events were reweighted to the NNLO accuracy using NNLO fixed-order predictions from MCFM-8.0 [46] with central factorization and renormalization scales set to the sum of the Higgs boson and the Z boson mass, µ = M H + M Z .The NNLO prediction for the decay was obtained from an analytical result of Ref. [47] and uses the Higgs boson mass as a central renormalization scale.When interfacing the fixed-order predictions with a parton shower we use Pythia8 [48] using the Monash tune [49].
The scale uncertainty is obtained by correlating the scale variation factors (K r , K f ) of the MiNLO and NNLO predictions that enter the reweighting procedure.The theoretical uncertainty is estimated by performing the usual 7 point scale variation, i.e. it corresponds to an envelop of 7 scale combinations, according to 1/2 ≤ K r /K f < 2.
We consider the same fiducial cuts as in Ref. [30].In particular, we require two charged leptons with |y | < 2.5 and p t, > 7 GeV.The harder lepton should satisfy p t, > 27 GeV, the invariant mass of the leptons-system should be in the range 81 GeV < M ¯ < 101 GeV.Additionally we require at least two b-jets with |η b | < 2.5 and p t,b > 20 GeV.Jets are defined using the flavour-k t algorithm [50], where we consider only b-quarks to be flavoured, and all other light quarks to be flavourless.

Results:
We start by showing the invariant mass of the two b-jets in Fig. 2. The jets are reconstructed using the flavour-k t algorithm [50] with R = 0.4 (upper plots) and R = 0.7 (lower plots).If more b-jets are present in the final state, the pair with the invariant mass closes to M H is selected.Plots on the right-hand-side include loop-induced gg → HZ contribution.The cuts applied here are those described earlier in this subsection with an additional cut of p t,Z = 150 GeV to select boosted Higgs-decays.The plots show predictions at fixed-order (green curve) as implemented in MCFM-8.0 [46], which account for NNLO corrections in production and NLO corrections in the decay, the NNLOPS predictions of Ref. [30], which also have NNLO accuracy in production and only NLO corrections to the decay (red) and the predictions of this work (black), which are NNLO accurate both in production and decay.The labels "PS1" and "PS2" make the reader aware of the fact that differences between the two NNLOPS predictions are not only due to the NLO vs NNLO treatment of the decay, but also to the different handling of the matching to the parton shower, which for PS1 uses a single overall scalup bound for a whole event, as described in Sec.3.2 of Ref. [30], while the PS2 method comes with a separate scalup bound for production and decay stages, as outlined in Sec.3.2 of this work.
We first focus on the plots without gluon-induced contributions.The large differences between pure fixed order predictions (green) and NNLOPS approaches, in particular in the low invariant mass region, has already been discussed in detail in Ref. [30].Here, we focus instead on the comparison between NLO-decay-PS1 (red) and NNLO-decay-PS2 (black) predictions.The most striking difference between the two is in the size of the uncertainty bands.As was already discussed in Ref. [30], the uncertainty of the NLO-decay-PS1 result is underestimated due to the well-known feature of the POWHEG simulation, i.e. the scale is varied at the level of the B function, which is inclusive over radiation, while the M b b spectrum is sensitive to secondary radiation.On the other hand the NNLO-decay-PS2 has a much more realistic uncertainty estimate, as the uncertainty is driven by the renormalization We shown results at pure fixed-order, NNLO in production and NLO in decay as implemented in MCFM-8.0 [46] (green), the NNLOPS predictions of Ref. [30], which have only NLO corrections to the decay (red) and the predictions of this work (black).Predictions of Ref. [30] and this work also differ in the handling of the matching to the parton shower, see text for more details.
scale variation in Eq. (2.3).We also note that the size of the relative uncertainty in the fixed-order predictions and the NNLO-decay-PS2 predictions are comparable.
Within their uncertainties, NLO-decay-PS1 and NNLO-decay-PS2 predictions agree very well in the whole M b b spectrum with the exception of the region close to M b b = M H .In particular, we notice that just above (below) the peak the NLO-decay-PS1 prediction lies below (above) the NNLO-decay-PS2 prediction.These differences are not unexpected.Notably, the hardest emission from the decay was treated differently in the two results, i.e. in the NLO-decay-PS1 approach it was generated using the POWHEG Sudakov, while in the NNLO-decay-PS2 simulation it is controlled by the MiNLO Sudakov.Furthermore, the composition of the Les Houches events in the two cases is rather different.For instance, in the NLO-decay-PS1 approach one can start with pure H → b b decay events, not accompanied by additional radiation, where subsequent emissions from the decay generated completely by the parton shower, with a value of scalup determined by the hardest POWHEG emission in production.In NNLO-decay-PS2 events instead one can have events with two emissions from the decay, which can further be combined with a production event involving up to two additional emissions.As a result of the above dissimilarities between the LH events, different values of scalup are used in the two approaches.While both methods are formally correct, numerical differences are not surprising.Altogether, the NNLO-decay-PS2 approach has the advantage that more emissions are generated using exact matrix elements which are properly matched to PS, whereas in the NLO-decay-PS1 approach part of this task was left to the parton shower alone.Similar effects are observed when the gluon-induced contribution is included, as presented in the right-hand plots in Fig. 2. In this case, an additional difference between the two approaches is due to the fact that all the decays in the NNLO-decay-PS2 event-sample are generated using the interface method described in Sec.3.1, whereas in the case of NLO-decay-PS1 the decay for these events, which are formally already O(α 2 s ), was always generated by Pythia.
Next, in Fig. 3 we show the transverse momentum of the two b-jets used for the Higgs reconstruction, p t,b b, without (left) and with (right) a p t,Z cut of 150 GeV.Again the main features of these plots were already described in Ref. [30], see discussion of Fig. 9 therein.Here, we focus on the comparison between the predictions of that paper and this work.We note that, in this case, both the NLO-decay-PS1 and the NNLO-decay-PS2 predictions have realistic scale-uncertainty bands.The uncertainty band of the latter result is slightly narrower.This is expected since, for the new result, NNLO corrections in the decay are included.
In Fig. 4 we show the transverse momentum of the hardest of the two b-jets, p t,b (hard).Both plots are considered without the boosted Higgs-decay selection, i.e. the p t,Z > 150 GeV cut is omitted.The left (right) plot is without (with) gluon-induced terms.We observe that the shape of the two NNLOPS predictions are significantly different from the MCFM-8.0one.This is mostly due to final state radiation generated by the parton shower which falls out of the jet, thereby softening its spectrum.
In the right-hand side plot, we also notice that, when the gluon-induced terms are included, the uncertainty bands of all predictions gets larger, as can be seen clearly above ∼ 100 GeV.As was already discussed in Ref. [30], due to their kinematics, this is the region where these events contribute, and the widening is due to the fact that they have a LO-like uncertainty band.
Finally, in the region up to ∼ 80 GeV, the NLO-decay-PS1 uncertainty band is considerably smaller than the one of the NNLO-decay-PS2 result.This region is populated mostly by events with at least one relatively hard emission from the decay.When such an emission is generated as the POWHEG hardest emission (NLO-decay-PS1), the scale uncertainty is underestimated (see discussion about Fig. 2).The NNLO-decay-PS2 result provides a more accurate uncertainty band.We have also checked that the relative size of this band in this region is slightly smaller than the MCFM-8.0one, as expected given the higher perturbative accuracy of the NNLO-decay-PS2 description of the decay.

Vector boson fusion Higgs production
To demonstrate the flexibility of our program, we consider now the NLO production of a Higgs boson via vector boson fusion (VBF) matched to parton shower, as implemented in POWHEG [51].We consider the Higgs boson decaying to b-quarks, and study the effect of including NNLOPS corrections to the decay.We compared our predictions to a more standard treatment where Pythia8 fully handles the Higgs decay and all subsequent radiation off the decay products.
The scale variation procedure we use for VBF is similar to the one used for HZ events, except that here we further restrict to a 3-points scale variation by correlating renormalisation and factorisation scale variation, i.e.K r = K f .We also fully correlate the renormalization scale variation for the production and the decay.
In order to show results in a fiducial phase space which is similar to the one probed in a real analysis, we apply cuts that enhance the vector boson fusion signal with respect to its typical backgrounds.In particular we require two tagging jets with p t,j > 30 GeV, |y j | < 5, y j1 • y j2 < 0, M jj > 600 GeV and |y j1 − y j2 | > 4.2.Furthermore, we require the presence of two b-jets with p t,b > 20 GeV and |y b | < 2.5.As for the HZ results shown above, jets are defined using the flavour-k t algorithm, where we consider only b-quarks to be flavoured, and all other light quarks to be flavourless.
Results: As for the results shown in the previous section, the more interesting observables to show the effect of the NNLOPS treatment of the H → b b decay are those obtained from the b-jet kinematics, and we focus on the invariant mass of the two b-jets (M b b), shown in Fig. 5, and the transverse momentum of the hardest of the two b-jets (p t,b (hard)) in Fig. 6.
Jets are reconstructed using the flavour-k t algorithm [50], and we show results for R = 0.4 (left plots) and R = 1 (right plots).As done for HZ, if more than two b-jets are present in the final state, the pair with the invariant mass closes to M H is selected.
For VBF, we compare our new result (NLOPS for VBF production, NNLOPS for H → b b) against the result where the H → b b decay is handled by Pythia8.These results are labeled in the plots as "NLO+NNLO-dec" (green) and "NLO+MC-dec" (red), respectively.
The NNLOPS treatment of the decay has its most important impact on the invariant mass of the two b-jets that originate from the Higgs decay, shown in Fig. 5. First we focus on the shapes of the central predictions.First of all, one can notice that, above the peak, the cross section drops very quickly, much more than in the HZ case, see Fig. 2.This is due to the nature of VBF production, in particular when tight VBF cuts are applied.In fact, for M b b, the region above the peak is dominated by radiation from production being clustered together with the b-jets, i.e. a configuration which is clearly suppressed by the VBF cuts, where the light jets are required to be far from the central rapidity region, where the Higgs boson is typically found.
Far from the peak, the two results are in rather good agreement, and we observe a pattern not too different from what was observed in Ref. [30], see Fig. 10 therein.We attribute the shape difference observed just below the peak (100 GeV M b b < M H ) to differences between the Pythia8 and the MiNLO Sudakov.The distribution in the left plot of the Fig. 5 peaks at a slightly smaller value of M b b than it does in the right plot.This behaviour is attributed to the jet radius used in the clustering sequence, i.e. for a bigger jet radius (R = 1.0, right plot) more radiation is clustered within b-jets moving the peak The VBF cuts described in the text are applied.The jets are reconstructed with the flavour algorithm of Ref. [50] for R = 0.4 (left) and R = 1.0 (right).
to higher M b b values than for smaller jet radius (R = 0.4, left plot).The differences close to the peak, driven by soft emissions off the b-quarks, can be associated with the different treatment of their masses in the two results, i.e. the "NLO-NNLO-dec" features massless b-quarks while the "NLO-MC-dec" one treats the Higgs decay to massive b-quarks.Note that the difference seems to be particularly large in the first bin to the right of the peak.This is not surprising since the absolute values of the cross-section in the neighbouring bins differ by almost two orders of magnitude and even a small migration of events from the peak induces large effects in the next bin.
The size of the uncertainty band also requires some explanation.It is about 15-20% in the case of the "NLO+NNLO-dec" prediction below the Higgs mass, and about 5% above.This has to be contrasted with an uncertainty band which is at most a few percent throughout the whole mass spectrum when the decay is handled by Pythia8.In the latter situation, the smallness of the band is due to a well-known POWHEG feature, namely the fact that the only source of scale dependence comes from the POWHEG B function, which is an inclusive quantity with respect to the radiation kinematics.For VBF, the NLO scale dependence is very small, and since we are not performing any scale variation within Pythia8, the uncertainty band for M b b is very small.On the other hand, the scale uncertainty band for the "NLO+NNLO-dec" result is as expected: below the peak the result features a scalevariation uncertainty in line with the fact that this is the kinematic region dominated by radiation off the b-quarks, and thereby the size of the uncertainty is of NLO type.This is the same pattern as observed in the HZ case.
Next, in Fig. 6, we present the distributions of the transverse momentum of the hardest b-jet.Both results are in very good agreement.Nevertheless, we observe a slight distortion of the "NLO+NNLO-dec" result with respect to the "NLO+MC-dec" one for small values of the p t,b (hard) in case of R = 0.4 jet radius (Fig. 6, left plot).Once again, this can be attributed to out-of-cone radiation which is more pronounced in the "NLO-NNLO-dec" result.Such an effect causes events with moderate initial value of the p t,b to fail to fulfill fiducial requirements after the QCD radiation due to the NNLO treatment and PS evolution.We also notice that the scale uncertainty of the "NLO+NNLO-dec" band widens in the region of small transverse momenta consistently with the fact that in this region perturbative uncertainties are bigger than for the hard part of the spectrum.
As a final remark, we point out that all of the results presented in this work do not contain an assessment of PS uncertainties.This is particularly relevant for the "NLO+MCdec" case in the VBF example, where there is no estimate of the uncertainty related to the Higgs boson decay to b-quarks.A comprehensive assessment of this type of uncertainties goes beyond the scope of this paper.Nevertheless, the results of this paper which treat the Higgs boson decay at NNLO accuracy provide a more realistic description of the shape and uncertainties for various observables.

Conclusions
In this work, we presented a simple method to describe the H → b b decay with state-ofthe-art NNLO accuracy consistently matched to a parton shower.We also provide a simple tool that allows to combine decay events with any available (N)NLOPS description of a Higgs production process.
The method exploits the factorization between production and decay of the Higgs boson, and relies on generating Higgs production events at (N)NLO accuracy, where the Higgs boson is undecayed, and combining them with Higgs decay events that describe the decay of the Higgs boson to b-quarks at NNLO accuracy.These combined LH events can then be showered using any transverse momentum ordered parton shower.Since radiation from production and decay factorize, the parton shower requires two separate upper bounds for the radiation, one for production and one for the decay.The former one is stored in the event file, while the latter can be easily computed from the decay kinematics.
As a proof of principle, we studied the effect of including NNLOPS corrections to the decay in the case of HZ Higgs production and VBF Higgs production processes.We find that parton shower effects are moderate once NNLO corrections to the decay are included, nevertheless including higher-order matrix elements provides a more realistic description of the decay and a better estimate of the theoretical uncertainty.
We have released the source code of the program, as an h_bbg user process in the POWHEG-BOX-V2, that enables a user to generate Hbbg-MiNLO events and further reweight them such that the NNLO H → b b decay width is reproduced.Although the generation process is straightforward and extremely fast, a ready-to-use sample of event files is also available upon request.We also include a tool to merge decay events with a Higgs production mode and a parton shower driver that allows for a consistent interface with Pythia.
Currently, the NNLO corrections to the decay are treated in the massless approximation.Nevertheless, NNLO corrections to the H → b b decay with massive b-quarks can also be included in a similar manner.At the moment, the code can be interfaced to any single-Higgs production mode but extension to double-Higgs production, where one or both Higgs bosons decay to b-quarks, is possible.We leave this for future work.from the App.B of Ref. [34], by setting a = 2 and taking the limit b → 0. They read In the expressions above, the constant part of the radiator (i.e.R(y 3 → 0)) has been removed, as these corrections are already reconstructed by the H → 4 partons matrix elements of the real corrections to the H → b bg process.Moreover, we also neglect terms proportional to the A 3 coefficient as they give rise to terms O(α 2 s ) terms after integration over radiation which are beyond our accuracy.
The values of the β-function coefficients read The relevant anomalous dimensions read where with N f being a number of light quark flavours.The B 2 anomalous dimension requires a slightly more detailed treatment.Our starting point is the Drell-Yan value [52,53] We are now ready to assemble the Sudakov form factor used in Eq. (2.3).Furthermore, in order to probe the uncertainty related to missing higher-order terms, the renormalisation scale dependence in the Sudakov radiator R(y 3 ), that is α s (M H ) → α s (µ R ) change, can be implemented along the lines of Eqs.(B.7)-(B.9) of Ref. [33].Similarly, we can introduce the resummation scale dependence that allows to estimate the uncertainty related to missing higher-order logarithmic terms, i.e.
L −→ L = ln(κ 2 res /y 3 ), (A.14) with κ res = Q res /M H .This shift is compensated by the following changes in the g i (x) functions, which guarantee that the logarithmic accuracy is preserved: We set the value of the central resummation scale to be Q res = M H .For the results presented in this work the resummation scale was not varied, i.e. κ res = 1.As discussed in App.A, the radiator functions include all logarithms related to the exponentiation of the one-gluon result, but also all NNLL real-radiation multiple emission effects which are incorporated through a modification of the B 2 coefficient.Σ In order to see which terms can be dropped without spoiling the NLO accuracy of the result, we remind the reader that the power-counting is fixed by the integral

B Analytic proof of NLO accuracy of
It is then easy to see that, provided the expansion is done using α s = α s (y 3 Q 2 ), the largest term beyond O(α 2 s ) comes from the A 3 coefficient in the Sudakov form factor.After taking the derivative, these terms contribute as α 3 s L. Hence, upon integration, according to Eq. (B.10) they give rise to an NNLO correction.This means that all O(α 3 s ) terms in D(α s , y 3 ) can be dropped without spoiling the NLO accuracy of the result after integration.From Eq. (B.8) it is clear that D(α s , y 3 ) also contains a term α 2 s B 2 , which upon integration gives a correction O(α 3/2 s ) and therefore has to be included.On the other hand Eq. (B.7) can be expanded in powers of the coupling constant as dΣ(y 3 ) dL = α s dΣ (1)  By explicitly taking the derivative of Eq. (B.2), the above terms can be writen as dΣ (1) (y 3 ) dL = dR (1) (y 3 ) dL = D (1) (y 3 ) (B.12) and dΣ (2) (y 3 ) dL = R (1) (y 3 )D (1) (y 3 ) + D (2) (y 3 ) .(B.13) Accordingly, dropping terms that are order α 3 s and higher, one can rewrite Eq. (B.7) as dΣ(y 3 ) dL = e R(αs,y 3 ) α s 1 − α s R (1) (y 3 ) dΣ (1)  We now clearly see that Eq. (B.14) agrees with the standard MiNLO formula, which then also must have NLO accuracy provided the Sudakov form factor includes the B 2 resummation coefficient.We stress that it is crucial that the expansion in Eq. (B.9) is performed around α s = α s (y 3 Q 2 ).Would one expand around α s (Q 2 ) then O(α 3 s L 2 ) terms would be present in Eq. (B.9) which would contribute as O(α 3/2 s ).Accordingly, this sets the scale of the strong coupling constant in the MiNLO framework to µ 2 = y 3 Q 2 .
. We set M H = 125.0GeV, Γ H = 4.14 MeV, M Z = 91.1876GeV, Γ Z = 2.4952 GeV, M W = 80.398 GeV and Γ W = 2.141 GeV.Moreover we set G F = 1.166387 • 10 −5 GeV −2 , sin 2 θ W = 0.23102, and α EM (M Z ) = 1/128.39.The H → b b branching ratio is set to Br(H → b b) = 0.5824 [22].For the contributions where the Higgs boson is radiated from a heavy-quark loop we use the pole mass of the heavy quark.We set the pole mass of the bottom quark to m b = 4.92 GeV and the pole mass of the top quark to m t = 173.2GeV.

→→→→RatioFigure 2 .
Figure 2. Invariant mass of the two b-jets reconstructed using with the flavour-k t algorithm with R = 0.4 (upper plots) and R = 0.7 (lower plots).Left (right) plots are without (with) gluoninduced terms.We shown results at pure fixed-order, NNLO in production and NLO in decay as implemented in MCFM-8.0[46](green), the NNLOPS predictions of Ref.[30], which have only NLO corrections to the decay (red) and the predictions of this work (black).Predictions of Ref.[30] and this work also differ in the handling of the matching to the parton shower, see text for more details.

→→RatioFigure 3 .
Figure 3.As Fig. 2 but for the transverse momentum of the two b-jets used for the Higgs reconstruction, p t,b b, without (left) and with (right) a p t,Z cut of 150

RatioFigure 4 .
Figure 4.As Fig. 2 but for the hardest of the two b-jets used for the Higgs reconstruction, p t,b (hard), without (left) and with (right) the gluon-induced terms, and without the cut on p t,Z .

RatioFigure 5 .
Figure 5. Invariant mass of the two b-jets in VBF Higgs events with Higgs decaying to b-quarks.The VBF cuts described in the text are applied.The jets are reconstructed with the flavour algorithm of Ref.[50] for R = 0.4 (left) and R = 1.0 (right).

RatioFigure 6 .
Figure 6.Transverse momentum of the hardest b-jet in VBF Higgs events with Higgs decaying to b-quarks.The VBF cuts described in the text are applied.The jets are reconstructed with the flavour algorithm of Ref. [50] for R = 0.4 (left) and R = 1.0 (right).