Resummed predictions for hadronic Higgs boson decays

We present the NNLL′ resummed 2-jettiness distribution for decays of the Standard Model Higgs boson to a bb¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{b} $$\end{document}-quark pair and to gluons. The calculation exploits a factorisation formula derived using Soft-Collinear Effective Theory, in which large logarithms of the 2-jettiness are resummed by renormalisation group evolution of the hard, soft and jet contributions to the differential decay rate. We match the resummed predictions to the fixed-order NNLO result using the Geneva framework, extending the validity of the results to all values of the resolution variable and providing a fully exclusive NNLO event generator matched to the Pythia8 parton shower.


Introduction
The lack of any definitive signal of New Physics at the Large Hadron Collider (LHC) suggests that the high-energy physics community must be open to alternative ways to probe Beyond the Standard Model (BSM) effects at collider experiments. In particular, precision measurements of the Standard Model (SM) Higgs sector may be a way to indirectly constrain BSM theories which live at scales beyond our reach, both at the LHC and at future lepton colliders or 'Higgs factories'. It is therefore crucial that theoretical predictions for processes involving the production or decay of the Higgs boson have a precision which matches that of experiment.

JHEP04(2021)254
In this work we consider the hadronic decays of a Higgs boson to a bb-quark pair and to gluons. 1 Although these modes are difficult to observe at a hadron collider due to the large QCD background, they may well prove to be useful probes of Higgs physics at a future lepton collider since they contribute significantly to the total width of the Higgs boson.
The dominant decay mode of the Higgs boson is to a bb-quark pair, with a branching ratio of about 58% [2]. An accurate measurement of this channel would allow an extraction of the Yukawa coupling y b , which is an important input to Higgs studies. The hadronic environment at the LHC makes this a challenging prospect -nevertheless, the decay has been observed recently by both ATLAS and CMS in the V H, or Higgsstrahlung, production channel [3,4]. QCD corrections to the partial width are known up to N 4 LO [5][6][7][8][9][10][11][12], and fully differential NNLO calculations have also been available for quite some time [13,14]. Recently, the fully differential N 3 LO calculation has also been completed [15].
The decay channel to gluons, on the other hand, proceeds via a top-quark loop and contributes around 8% to the total width. Since QCD corrections to this channel are indistinguishable from the bb-quark case at higher orders in perturbation theory, one should consider the classes of processes together to obtain a total hadronic width, as performed in ref. [1]. Nevertheless, at NNLO and in the case of kinematically massless b-quarks the processes can be fully separated, since interference terms between the diagrams vanish. Related complications which arise at N 3 LO have been studied in ref. [16]. In the case of massive b-quarks, these interference terms can no longer be neglected -their impact has been studied in ref. [17] and a full NNLO calculation in the massive case has been carried out in refs. [18,19].
In the limit that M H 2m t , the top-quark loop which couples the Higgs boson to gluons can be integrated out to obtain an effective theory with five light active flavours in which the interaction is local. This simplifies the inclusion of QCD corrections and has allowed calculations to be performed at NNLO [20,21] and, for the total width, at N 3 LO [22] and N 4 LO [23]. The effect of including a finite top-quark mass on the total width has also been studied in ref. [24].
In light of the importance of Higgs physics, several other predictions at various accuracies and using different approximations are available beyond those listed here. A complete review of the state of the theoretical calculations for Higgs boson production and decay processes, including the calculation of electroweak corrections, can be found in ref. [25].
In a recent publication [26], the distributions of the thrust variable in these decay processes were considered and fixed-order computations up to approximate NNLO (which contribute at O(α 3 s ) relative to a Born H → bb/gg process) were performed. In that work, the authors noted the poor convergence of the perturbative series for both processes and were able to show that the approximate NNLO corrections obtained from the singular terms of a SCET-derived factorisation formula could ameliorate the scale dependence of the calculations. They also acknowledged several shortcomings of their calculation, one of which related to the size of the logarithms log n τ /τ which are not resummed in a purely fixed-JHEP04(2021)254 order computation and spoil predictivity in the small τ region. Here, we provide resummed predictions at NNLL accuracy which complement the results of ref. [26]. Compared to NNLL, the resummation at NNLL accuracy incorporates the complete O(α 2 s ) singular structure for T 2 → 0, i.e. all 2-loop virtual and corresponding real corrections, allowing us to consistently match to NNLO.
Using the Geneva formalism developed in refs. [27][28][29], we are also able to construct IR-finite events which combine the advantages of the resummed and fixed-order calculations and are matched to a parton shower. 2 Having a Geneva implementation of the H → bb process will also allow us to produce an NNLOPS generator for the signal process pp → V (H → bb). This can be achieved by combining the Higgs boson decay presented here with our previous calculation of the V H production process [31] in the narrow width approximation. Fixed-order calculations for the full pp → V bb process were performed in refs. [32][33][34] in the massless approximation -more recently, a calculation with massive b-quarks also appeared [19]. An NNLOPS generator for W (H → bb) production via the MiNLO method [35][36][37][38] was presented in [39], while a separate NNLOPS H → bb generator was also made available in ref. [40]. Nonetheless, we believe an independent implementation of the combined corrections to both production and decay in the Geneva framework will provide a useful cross-check of previous results. We leave this development to a future publication. This paper is organised as follows. In section 2, we briefly explain how resummed predictions are obtained from a factorisation formula derived in soft-collinear effective theory and provide numerical results for the resummed T 2 distribution in H → bb and H → gg. In section 3, we briefly recap the main features of the Geneva method relevant for the processes at hand. In particular, we discuss various implementation details, as well as how the matching to the parton shower is achieved. We present our Geneva results in section 4. Finally, we report our conclusions and directions for future work in section 5, while we detail the construction of the phase space mappings used and the analytical NNLO decay rates in appendices A and B respectively.

Resummation from Soft-Collinear Effective Theory
In this section we present, for the first time, the NNLL resummation of the 2-jettiness observable, T 2 , for the decay of a Higgs boson into a pair of b-quarks. We also provide results at the same accuracy for the Higgs boson decay into gluons, which were first presented in ref. [30]. We present numerical results for the dimensionless τ ≡ T 2 /(2 M H ) distribution, where M H is the mass of the Higgs boson.

Formulation
Our basic resolution parameter for the hadronic decays of the Higgs boson is the 2-jettiness, defined as

JHEP04(2021)254
where k runs over all final state particles with momenta p k = (E k , p k ) and n is the unit 3vector resulting from the minimisation procedure. In the case where all final-state particles are massless, it is related to the more familiar thrust T , which was widely studied for e + e − collisions [41,42] and extended to hadronic collisions in ref. [43], by the relation . For the decays that we consider, we work in the rest frame of the Higgs boson and always have that E cm = M H . Exactly like the thrust variable, the 2-jettiness is constrained kinematically (0 ≤ T 2 ≤ M H ) and its value is related to the spatial distribution of the radiation: in the limit T 2 → 0 the final state consists of two pencil-like jets, while for T 2 ∼ M H there are three or more jets distributed in a more spherical configuration. We consider the decay rate differential in T 2 of a Higgs boson to either a bb-quark pair or to a pair of gluons. We consider massless b-quarks with a finite Yukawa coupling to the Higgs y b . In the gluon case, we work in an effective theory in which the top-quark loop that couples the Higgs boson to gluons has been integrated out, leaving an effective local operator Hgg.
The Born level decay rates for the two processes considered are given by It has been shown, both in QCD and SCET, that the differential decay rate factorises in the small T 2 limit [42,[44][45][46] as where the index i = b, g indicates the process in question.
The decay rate has been factorised into a hard contribution H i (M H , µ), a soft function S i (k, µ) and two jet functions J i n (p 2 n , µ) and J ī n (p 2 n , µ). The hard function is defined as the square of the Wilson coefficients which match the full theory (the SM) onto SCET. In the gluon case, an additional matching is required from the heavy-top limit effective theory we are working in onto the SM -thus, the hard functions can be written as The jet functions describe collinear radiation from the Born-level partons along the jet directions n andn, which can be chosen without loss of generality to be orientated alongẑ, viz. n = (1, 0, 0, 1) andn = (1, 0, 0, −1). The soft function accounts for all soft radiation.
Each component in the factorisation theorem must be evaluated at its own characteristic scale in order to prevent the appearance of large logarithms, viz.
However, since the decay rate must be evaluated at a single scale, we evolve the separate functions to a common scale µ via renormalisation group (RG) evolution and in so doing resum the large logarithms of ratios of scales which appear. The resummed spectrum, differential in the Born kinematics, can then be written as

JHEP04(2021)254
where dΓ i B (µ)/dΦ 2 indicates the Born decay rate differential in the two-body phase space and we have used ⊗ to denote the convolutions, dropping the explicit dependence on the convolution variables. The ingredients necessary for NNLL accuracy are all available in the literature and many have been compiled in ref. [26]. For the H → bb case, we take the hard function from ref. [26]. Since we are considering only massless b-quarks, we can use the soft and jet functions as implemented in ref. [27] (and first calculated at NNLO in refs. [47][48][49]) for e + e − → jj, also recycling the evolution kernels from that work. For the H → gg case, the fully expanded hard function, including contributions to the Hgg effective vertex from both SCET and the effective theory where the top-quark is integrated out, appears in ref. [30]. The NNLO jet function [50] and the evolution kernels are also given therein, while we obtain the soft function via a Casimir rescaling of the H → bb case.

Numerical results
We have implemented the resummed calculation, eq.

JHEP04(2021)254
Because of the resummation of large logarithms of τ , our results provide a physical description at small τ , as can be seen in figure 1, improving on the resummed-expanded results presented in ref. [26]. However, the range of validity of the resummed calculation is still limited to low τ values, since the factorisation formula we rely on (eq. (2.3)) is only valid there. Comparing the two decay channels, we see that the H → bb process presents a higher peak, located at a lower value of τ , with a narrower width. The peak of the τ distribution in the H → gg case is instead lower and shifted to larger values of τ , with a broader width. This behaviour is in line with expectations based on a naive analysis of the Casimir scaling of the two processes.
In order to extend the validity of the calculation to all values of τ , one needs to match the resummed result to the fixed-order calculation, which provides the physical behaviour at large τ . Indeed, in this region, the nonsingular contribution becomes sizeable and exponentiating the singular contribution is no longer the correct approach. The matching of fixed-order calculations to resummation has long been established at the level of the resummed observable: the most straightforward approach simply adds the results for the resummed and fixed-order distributions in the τ variable and then subtracts the expansion of the resummed result up to the same order included in the fixed-order result. In this way the calculation is free from doubly counted contributions up to the given perturbative order and includes all the higher-order terms properly resummed.
While the approach just outlined works flawlessly for the τ distribution we are resumming, it is not directly applicable to the construction of a fully exclusive event generator. In the next section, section 3, we show how this can be achieved by means of the Geneva method, allowing us to perform the matching at the fully differential level. 3 3 Implementation in the Geneva framework

Geneva in a nutshell
The Geneva framework allows the matching of a resummed to a fixed-order calculation and thence to parton shower programs such as Pythia [51]. In so doing, it provides theoretical predictions which are accurate over the whole phase space and which describe realistic events of high multiplicity. These can then be hadronised and fed into the analysis routines used by the experimental collaborations. The method for separating events into different multiplicity bins and for performing the matching has already been described thoroughly in refs. [27,29] and related references. Therefore, in this context, we limit ourselves to restating the primary outcomes as applicable to the case of Higgs boson hadronic decays up to NNLL T +NNLO 2 accuracy. We remind the reader that in order to achieve a sensible separation between the exclusive 3-jet and the inclusive 4-jet decay rates, we must also perform the resummation of T cut

JHEP04(2021)254
Dropping the process label i for ease of notation, the Geneva Monte Carlo expressions for the exclusive 2-jet, 3-jet and the inclusive 4-jet rates are given by 4 where the B j , V j and W j are the 0-, 1-and 2-loop matrix elements for j partons in the final state. In the equations above, we have introduced the shorthand notation to indicate that the integration over a region of the M -body phase space is done keeping the N -body phase space and the value of some specific observable O fixed, with N ≤ M . 4 We make a slight abuse of notation in order to highlight the dependence of the dΓ mc i decay rate on the resolution parameters. When an argument contains a single term, e.g. T cut N , it means that the corresponding quantity has been integrated over up to the value of the argument. An argument TN > T cut N implies instead that the corresponding decay rate remains differential in the relevant resolution variable for values larger than the cutoff.

JHEP04(2021)254
The Θ O (Φ N ) term in the previous equation limits the integration to the phase space points included in the singular contribution for the given observable O. For example, when generating 3-body events we use where the map used by the 3 → 4 splitting has been constructed to preserve T 2 , i.e.
and Θ T (Φ 4 ) defines the projectable region of Φ 4 which can be reached starting from a point in Φ 3 with a specific value of T 2 . The usage of a T 2 -preserving mapping is necessary to ensure that the pointwise singular T 2 dependence is alike among all terms in eqs. (3.2) and (3.4) and that the cancellation of said singular terms is guaranteed on an event-by-event basis. The expressions in eqs. (3.3) and (3.5) encode the nonsingular contributions to the 3and 4-jet rates which arise from non-projectable configurations below the corresponding cut. This is highlighted by the appearance of the complementary Θ functions,Θ O , which account for any configuration which is not projectable either because it would result in an invalid underlying-Born flavour structure or because it does not satisfy the T 2 -preserving mapping (see also ref. [31]).
The term V C 3 denotes the soft-virtual contribution of a standard NLO local subtraction (in our implementation, we follow the FKS subtraction as detailed in ref. [52]). We have that with C 4 a singular approximation of B 4 : in practice we use the subtraction counterterms which we integrate over the radiation variables dΦ 4 /dΦ C 3 using the singular limit C of the phase space mapping. U 3 is a LL Sudakov which resums large logarithms of T 3 and U 3 its derivative with respect to T 3 . Its exact form is given by where the Casimir factor C k depends on the flavour content of the 3-jet event, C qqg = 2C F + C A or C ggg = 3C A , and we run the coupling at NNLL order. The term P(Φ N +1 ) represents a normalised splitting probability which serves to extend the differential dependence of the resummed terms from the N -jet to the (N +1)-jet phase space. For example, in eq. (3.2), the term P(Φ 3 ) makes the resummed spectrum in the first term (which is naturally differential in the Φ 2 variables and T 2 ) differential also in the additional two variables needed to cover the full Φ 3 phase space. These splitting probabilities are normalised, i.e. they satisfy The two extra variables are chosen to be an energy ratio z and an azimuthal angle φ. In the soft and collinear limit, z = E sister /(E sister + E daughter ) where the daughter and the sister

JHEP04(2021)254
are assigned to be the pair of particles that are closest according to the N -jettiness metric and which therefore set the value of T N , i.e. which minimise the quantity The daughter particle is defined to be the gluon for q → qg splittings, the quark for g → qq and the softer gluon for g → gg. These definitions in hand, the normalised splitting probability is given by where AP sp (z, φ) is the unregularised Altarelli-Parisi splitting function. 5 The implementation of the splitting probability requires us to construct the full Φ N +1 phase space from Φ N and a value of T N . Similarly, we mentioned above that the real integration in the fixed-order part of the calculation requires us to project from Φ N +1 configurations onto Φ N while preserving the value of T 2 . Both of these tasks demand a map that satisfies eq. (3.8) -the construction of such a map is detailed in appendix A.

Implementation details
In this section we discuss the particulars of the implementation of the Higgs boson decay processes in Geneva. Throughout this section we use the same settings and values for SM parameters as in section 2.2. In the H → bb case, we implement the analytic matrix elements found in ref. [14], while in the H → gg case we interface to the OpenLoops package [53][54][55].

Profile scales
The resummation provided by the RGE of the functions in eq. (2.3) correctly accounts for logarithms of the form log(T 2 /2M H ) which become large in size for small values of T 2 . In the fixed-order region, however, where T 2 is larger, such logarithms are more modest in size and continuing to resum them would introduce undesirable higher-order contributions.
We must therefore switch off the resummation before this happens. This can be achieved by setting all scales to a common nonsingular scale in the fixed-order region, µ NS = µ S = µ J = µ H , which stops the evolution ensuring that the resummed contribution is cancelled out exactly by the resummed-expanded. In order to achieve a smooth transition between the resummation and the fixed-order (FO) regimes, we make use of profile scales µ J (T 2 ) and µ S (T 2 ) which interpolate between the characteristic scales and µ NS [56,57]. Specifically, we have that (3.14) JHEP04(2021)254 where the common profile function f run (x) is given by This form has strict canonical scaling below x 1 and switches off the resummation above x 3 ; for a = 1 it matches the form of the profiles used in e.g. ref. [58]. In order to determine the choice of parameters a, x i it is instructive to examine the relative sizes of the singular and nonsingular contributions as a function of T 2 to determine where the resummation should be switched off. This is done for the two decay channels as shown in figure 2. We see that the singular and nonsingular pieces become similar in size at around τ ≡ T 2 /(2M H ) ≈ 0.3, and therefore set for both processes We notice that in the limit τ → 0 the singular contribution becomes an increasingly good approximation to the fixed-order result, reflecting the proper cancellation of the singular terms between the fixed-order and resummed-expanded parts of the calculation. We set the remaining parameters a = 1, ax 0 = 3 GeV/M H for the gg channel following ref. [30], while for the bb channel we set a = 1/2, ax 0 = 2.5 GeV/M H .

JHEP04(2021)254
The uncertainties associated with the resummed and fixed-order calculations are estimated by varying the profile scales. For the uncertainty arising from the FO part, we adopt the usual prescription of varying µ FO up and down by a factor of 2 and taking the maximal absolute deviation from the central value as a measure of the uncertainty. This preserves everywhere the ratios between the various scales µ H , µ J and µ S and so the arguments of the logarithms which are resummed by the RGE factors are unaffected. In the resummed case, we vary the profile scales for µ J and µ S about their central profiles while keeping µ H = µ FO fixed. Specifically, defining a variation function (see e.g. ref. [59]) we vary the soft-and jet-function scales such that where η = 1/6. In this way the arguments of the resummed logarithms are varied in order to estimate the size of higher-order corrections in the resummed series while maintaining the scale hierarchy More details on the specifics of this prescription may be found in ref. [59]. In addition, we include two more profiles where we vary all x i transition points by ±0.025 simultaneously. We thus obtain 6 profile variations in total and take the maximal absolute deviation in the result from the central value as the resummation uncertainty. The total uncertainty is then obtained as the quadrature sum of the resummation and fixed-order uncertainties. The profiles and their variations are shown in figure 3 for both the bb and gg cases.

Comparing the spectrum and the derivative of the cumulant
Since the profile scales which we just discussed have themselves a functional dependence on T 2 , the integral of the spectrum that one obtains from eq. (2.5) is not exactly equal to the cumulant in eq. (3.1) evaluated at the highest scale.
Choosing canonical scaling, i.e. µ ∝ T cut 2 , we have where T max 2 is the upper kinematical limit. By integrating the spectrum we therefore obtain not only the cumulant but also unwanted additional terms of higher order. Depending on the convergence properties of the perturbative series, these additional terms can be numerically relevant and cause a sizeable difference between the inclusive decay rate obtained from the resummed calculation and the FO result. To obviate this problem, we supplement the JHEP04(2021)254 spectrum in eqs. (3.2) and (3.4) with an additional higher order term. The contribution of this term is restricted to the region of T 2 where the spurious N 3 LL terms are sizeable, and vanishes in the FO region; crucially, upon integration it ensures that the FO rate is recovered. It takes the form: where κ(T 2 ) and µ h (T 2 ) are smooth functions. It is clear that this vanishes in the FO region where µ h (T 2 ) ∼ M H as required -in order to restrict its contribution further, we also choose κ(T 2 ) to tend to zero in this region to minimise its size before exact cancellation is reached and choose the profile scale µ h (T 2 ) to reach M H at a lower value of T 2 than the rest of the calculation. This prevents the accuracy of the tail of the spectrum from being spoiled, while keeping the resulting changes in the peak region contained within its scale uncertainty band. We tune κ(T 2 ) to recover the correct inclusive rate, both for the central FO scale and also for its variations such that the result of integration is identical to a FO calculation for inclusive quantities.

Power-suppressed corrections to the nonsingular cumulant
The integration of the differential decay rate in eq. (3.1) over the Φ 2 phase space produces an NNLO accurate total width. For differential quantities, however, the O(α 2 s ) terms in eq. (3.1) are guaranteed to be NNLO accurate only up to power corrections in T cut 2 since any projective map one could devise could not preserve all Φ 2 quantities simultaneously. This fundamental limit on the accuracy of event generators actually allows us to sidestep the problem of implementing a full NNLO subtraction -since the total width is the JHEP04(2021)254 only quantity that is certain to be NNLO accurate, we can drop all the O(α 2 s ) terms in the cumulant and achieve the correct NNLO width by reweighting. That is, rather than implementing the full form of eq. (3.1), we instead use which requires only a local NLO subtraction. The remaining nonsingular terms take the form where the functions f i (T cut 2 , Φ 2 ) are at worst logarithmically divergent in the small T cut 2 limit. We include the NLO term proportional to f 1 (T cut 2 , Φ 2 ) in eq. (3.21) via an on-the-fly NLO 2 calculation, but neglect the f 2 (T cut 2 , Φ 2 ) piece. The size of this neglected term as a function of the cut is shown in figure 4 for both processes. We see that at our default value of T cut 2 = 1 GeV the missing O(α 2 s ) terms are of a size ∼ 10 −5 GeV in both cases. This amounts to a relative correction of O(0.4%) for the bb channel and of O(1%) for the gg. Smaller power corrections could naturally also be obtained by modifying the factorisation formula eq. (2.3) to include subleading power contributions [60,61] or by lowering further the value of T cut 2 . In this limit, however, the calculation suffers from numerical problems originating from the stability of the matrix elements and of the NLO subtraction procedure close to extreme soft or collinear configurations, which motivates our default choice.
In order to correct for this discrepancy and obtain the correct NNLO inclusive decay width, we may simply rescale the weights of the Φ 2 events in such a way that we match the known analytic result at NNLO. We are thus able to include the effects of the f 2 term JHEP04(2021)254 in eq. (3.22) on the total cross section that would have been present had we implemented eq. (3.1) literally. Since neither eq. (3.1) nor our approach in eq. (3.21) achieves the exact O(α 2 s ) Φ 2 dependence of all observables, our approximation does not inherently limit the accuracy of our predictions.

Interface to the parton shower
We briefly recap the main features of the parton shower interface in Geneva here and refer the interested reader to section 3 of ref. [29] for a more detailed discussion.
The partonic jet decay rates dΓ mc 2 , dΓ mc 3 and dΓ mc ≥4 each include contributions from higher multiplicity phase space points, but only in those cases where T N (Φ M ) < T cut N . In order to make the calculation fully differential in the higher multiplicities, a parton shower is interfaced which adds radiation to each jet decay rate in a unitary and recursive manner. Ideally, the shower should leave the values of the jet rates and their accuracy unaffected, restoring the emissions in dΓ mc 2 and dΓ mc 3 which were integrated over when the jet rates were constructed and also adding extra final-state partons to the inclusive dΓ mc ≥4 . For illustrative purposes, we consider a shower strongly ordered in T N , such that . . . A shower history of this kind could be constructed by taking the output of a shower ordered in a more conventional variable and reclustering the partons using the N -jettiness metric T N .
In general, the requirement of the preservation of the accuracy of the jet rates after applying the shower on a phase space point Φ N sets constraints on the point Φ M reached after the shower. For the cases in which the showered events originate from Φ 2 events, the main constraint is that the integral of the decay rate below the T cut 2 (which is NNLL +NNLO accurate) must not be modified. The emissions generated by the shower must in this case satisfy T 2 (Φ N ) < T cut 2 , so that they recover the events which were integrated over in the construction of the 2-jet exclusive decay rate and add events with more emissions below the cut. In case of a single shower emission we require also that the resulting Φ 3 point is projectable onto Φ 2 , as these are the only configurations at this order which are included in eq. (3.21). Both of these conditions can be implemented with a careful choice of the starting scale of the shower. The preservation of the decay rate below the cut is then ensured by the unitarity of the shower evolution. In practice, we allow for a tiny spillover up to 5% above T cut 2 in order to smoothen the transition. The showering of Φ 3 and Φ 4 events must be treated more carefully in order to preserve the NNLL +NNLO accuracy of the T 2 spectrum. Crucially, we must ensure that the Φ 4 points produced after the first emission are projectable onto Φ 3 using the T 2 -preserving map discussed in appendix A. Since the shower cannot guarantee this, we instead perform the first emission in Geneva (using the analytic form of the LL Sudakov factor and phase space maps) and only thereafter allow the shower to act as usual, subject to the restriction T 4 (Φ N ) ≤ T 3 (Φ 4 ). We apply this procedure only to the Φ 3 events and find that ) .

JHEP04(2021)254
By choosing Λ 3 ∼ Λ QCD , the Sudakov factor U 3 (T cut 3 , Λ 3 ) becomes vanishingly small and we can relax the shower conditions on the 3-jet contributions. The showered events therefore originate almost exclusively from either dΓ mc 2 or dΓ mc ≥4 . We choose starting scales of T cut 2 and T cut 3 for the Φ 2 and Φ 3 events respectively. For the Φ 4 events, the starting scale t needs to be a measure of the hardness of the splitting, for example the 3-jettiness value T 3 . Here we follow the choice made in ref. [40] and set where the energies are defined in the Higgs boson rest frame. After interfacing to the Pythia8 parton shower, we expect the accuracy for observables other than T 2 to be no worse than that of the standalone Pythia8 shower.

Nonperturbative power corrections and hadronisation
The approach described up to now does not take into account nonperturbative power corrections, which can significantly affect the partonic predictions. The framework of SCET allows these nonperturbative effects to be systematically included via the introduction of a shape function f (k, µ) modifying the soft component as [44,62,63] where S pert 2 is the perturbative soft function. At small T 2 ∼ Λ QCD the shape function gives an O(1) contribution to the cross section, while for larger T 2 values one can show that the dominant contribution is of O(Λ QCD /T 2 ) and results in a overall shift of the T 2 spectrum [56]. The same conclusions can be reached using a dispersive model and an effective value for the strong coupling constant in the nonperturbative regime [64][65][66][67].
The resummed predictions obtained by Geneva at the partonic level only include the perturbative soft function, and we delegate the provision of nonperturbative ingredients to the hadronisation models used in Pythia8. Therefore, after the showering stage, the events are interfaced to the phenomenological hadronisation model in Pythia8 without further constraints on the kinematics of the hadronised event. This means that the hadronisation can potentially cause significant shifts of the T 2 spectrum.
It is known that the 2-jettiness and the thrust observables receive different hadronisation corrections, due to the different treatment of the hadron masses in their definitions [68,69]. Since there are currently no experimental data with which we can compare for these decay channels, in this work we consistently use the definition of T 2 in eq. (2.1) even for hadronised events, despite the larger power corrections compared to schemes with a different mass treatment. This is different from the approach taken for the e + e − → jets study in ref. [27], where the definition based on thrust was used to compare to LEP data.
It is important to notice that we do not include uncertainties from these nonperturbative contributions in the results presented in the next section. In our approach, a crude estimate of their size could in principle be obtained by varying the tune parameters of the  Table 1. Comparison of Higgs boson partial widths obtained from NNLO analytic expressions and at the partonic level from Geneva. Note that, due to the presence of the power corrections displayed in figure 4, the values do not agree exactly within the statistical error and therefore a reweighting must be performed.
of this work. It is worth noting, however, that in any calculation obtained by matching higher-order calculations with parton shower one has to carefully evaluate which parameters are truly encoding nonperturbative effects and should therefore be tuned.

Results
In this section we present the full Geneva results obtained by matching the resummed calculation to the fixed order. We adopt the same values of SM parameters as in section 2.2 and set T cut 2 = T cut 3 = 1 GeV. We interface to the Pythia8 generator which showers our events 6 and use the e + e − tune 3, turn off QED effects and prevent the decay of b-hadrons. We set the strong coupling used by Pythia8 to α s = 0.118, although ideally, one should perform a dedicated tune to accommodate for this change.
With the setup as described, we verified that we obtain the correct NNLO decay rate up to the power corrections shown in figure 4. The partonic results are presented for each channel in table 1, where the analytic values have been obtained using the formulae appearing in appendix B. In general, the Geneva method also guarantees NNLO accuracy for distributions differential in the Born variables of the process (see for example ref. [31]). In the case of a spin-0 boson decaying into two particles, however, the Born phase space is parameterised by only two angles and is flat in both -there is therefore no non-trivial shape information which can be compared to a fixed-order calculation. We have, however, validated our NLO calculations of H → bbg and H → ggg/H → qqg against aMC@NLO [70] and found perfect agreement. We checked that by increasing the T cut 2 to ∼ 5 GeV we obtain smaller power corrections (see figure 4) -however, since this would limit our higher-order resummed predictions for the shape of the spectrum to T 2 > 5 GeV, in the following we continue to use T cut 2 = 1 GeV and accordingly reweight our events in order to obtain the correct total decay width. 6 The publicly available Pythia8.235 version we used has difficulty parsing events read from an LHEF file in which only one particle appears in the initial state -we therefore add dummy neutrino beams using code provided by S. Prestel to mimic a collider process.  A comparison of the NNLO+NNLL results at the partonic and showered levels is presented in figure 5 for the H → bb process and in figure 7 for the H → gg process, while the corresponding comparisons of the showered and hadronised events are shown in figures 6 and 8. The panels in the plots show three different regions of the 2-jettiness spectrum: the peak (leftmost panels), where resummation effects are expected to be dominant; the transition (centre panels), where the resummed and fixed-order calculations compete for importance; and the tail (rightmost panels), where the resummation is switched off and the fixed-order calculation provides the correct physical description. We observe that in the bb channel the T 2 is well preserved by the shower, while hadronisation effects shift the distribution to higher values of T 2 across all regions. This can be compared to the results obtained in ref. [27], keeping in mind the aforementioned difference between the 2-jettiness definitions used at hadron level and the different energy scale which result in competing contributions to the shift.

JHEP04(2021)254
In the gg channel, the effect of the shower on the T 2 spectrum is greater, especially at the lowest values of T 2 , but it preserves the shape of the distribution to within the scale variation bands closer to the peak and in the transition region. We also notice that the partonic and showered predictions give a negative cross section for very small values of T 2 , below the JHEP04(2021)254  nonperturbative freeze-out of the profile scales. This behaviour should not be concerning as it happens in a region where the perturbative resummed results are already questionable and, as mentioned before, we do not include any nonperturbative uncertainty. We have verified that the size of the negative value is augmented by the nonsingular corrections which we include via an additive approach. Indeed, when examining the resummed T 2 distribution alone, the behaviour at small T 2 remains negative but is compatible with a value of zero to within the quoted uncertainties.
A peculiar feature is observed in the first bin of figure 7, which contains the cross section below T cut 2 and is positive. This is again a consequence of the missing nonsingular corrections in eq. (3.21), which are included by the reweighting procedure. Since these are particularly large for this process, see figure 4, their effect is to change the sign of the cumulant below T cut 2 . We also observe somewhat larger effects on the spectrum due to hadronisation compared to the bb case, particularly in the peak region. The seemingly unusual behaviour at small T 2 is a consequence of the already discussed shift of the spectrum after hadronisation resulting in a smearing of the first bin. We stress again that the small error bands reported are due to the lack of nonperturbative uncertainties. Figure 9. Jet broadening and the JADE two-to-three differential jet rate at the partonic, showered and hadronised levels for H → bb. Figure 10. Jet broadening and the JADE two-to-three differential jet rate at the partonic, showered and hadronised levels for H → gg.

JHEP04(2021)254
Finally, in figures 9 and 10 we show the results for distributions other than the 2-jettiness that we use as input to our Geneva implementation, for the bb and gg cases respectively. We consider the JADE clustering metric y 23 for separating two exclusive jets from three or more [71,72] and the jet broadening (B T ) [73,74] event shape defined as follows where the sum runs over all final state particles andn T is the thrust axis.

JHEP04(2021)254
It is important to remark that we do not expect the Geneva method to provide a higher formal accuracy for these observables, but it is nonetheless interesting to observe the effects of our predictions at the various stages. In general the bb decay channel is better behaved after showering, providing results that are compatible with the predictions at the partonic level over the majority of the phase space.
We notice that, at small values of the jet broadening, some unappealing artefacts appear. We have verified that these features are a consequence of the additive matching to the fixed-order calculation, used to include the nonsingular corrections. They are not present, for example, when examining the jet broadening distribution of events obtained from the T 2 resummed calculation alone. We therefore conclude that the secondary peak which appears constitutes an effect beyond the perturbative accuracy of our formulae, despite being numerically large. In the future, one could explore whether by exploiting the freedom in the handling of the recoil by the 3 → 4 mapping one could ameliorate this undesirable effect -however, we do not pursue this further here. We also see deviations for particular values of the observables after hadronisation and hadron decays are included. In particular we notice a significant shift in the JADE y 23 observable for both decay channels, which is not unexpected for this specific jet-clustering measure.

Conclusions
In this work we have resummed the 2-jettiness at NNLL for hadronic Higgs boson decays in the bb and gg channels via a SCET approach. Compared to previous fixed-order results, we observe the expected improved behaviour in the small T 2 region, where the physical Sudakov peak is now described correctly. We have also implemented these processes in the Geneva framework, which has allowed us to match the resummed calculations with NNLO fixed-order predictions and to a parton shower. This has required an examination of the interplay of the singular and nonsingular contributions, in order to determine the region in which resummation effects are dominant and hence design profile scales which provide a smooth transition between the resummed and fixed-order regimes. As a result we have produced NNLO accurate event generators interfaced to the Pythia8 parton shower for the two processes, which provide accurate predictions in all regions of phase space.
We compared predictions at the partonic, showered and hadronised levels, finding the expected good agreement for the total decay rates and for the T 2 distribution up to the showered level. We observed larger differences due to the hadronisation, especially in the gg channel.
The completion of this work will eventually allow us to combine our H → bb result with the Geneva V H production generator in the narrow width approximation, yielding a full NNLOPS generator for the signal channel of the l + l − bb final state. Given the recent observation [3,4] of the Higgsstrahlung process by the ATLAS & CMS experiments at the LHC, this will constitute an important phenomenological result. It will also allow a direct comparison with the only other existing NNLOPS generator for this process [40]. In light of the findings in ref. [26] regarding the convergence of the perturbative series and the N 3 LO results at fixed order which are also available for this decay channel, it might also be JHEP04(2021)254 interesting to consider building an event generator at N 3 LOPS level. Another avenue for development might be the inclusion of a finite b-quark mass in the generator, given recent work on fixed-order calculations [19]. We leave this to future consideration.

A Constructing a 2-jettiness-preserving map
The map used for 3 → 4-body splittings and 4 → 3-body projections presented in this section was first developed and applied to the process e + e − → jj in ref. [27]. Here, we detail the construction of the map as used in that work and in addition provide the translation to the splitting variables T 2 , z and φ needed for the Higgs boson decay case.
We start by considering the case of a splitting, which takes as input N -body phase space points Φ N and generates from them (N+1)-body phase-space points Φ N +1 . Since we wish to calculate the NLO distribution in 2-jettiness, T 2 , while still generating exclusive Φ 3 points, we must use a map that produces Φ 4 points with the same value of T 2 as the Φ 3 points with which we started. Unfortunately, the construction of such a map is challenging since T 2 is a global variable. A more manageable approach is to seek a map which preserves not the exact 2-jettiness, T 2 , but instead a related quantity, the fully recursive 2-jettiness, T FR 2 , defined by the following procedure: 1. Recursively cluster the starting phase space point Φ M down to a Φ 3 point using the N -jettiness metric for final state particles The quantity we obtain through this procedure has the same singular structure as the exact T 2 , with any differences being captured by the nonsingular contributions. Starting from the 3-parton phase space point Φ 3 , which is the input of the splitting map, we label its momenta as

JHEP04(2021)254
The thrust axis will lie along the direction of the hardest parton (i.e. along p 1 ), and we have that When we split to a 4-parton event and then cluster back to a massive 3-parton event, the thrust axis is still determined by the most energetic of the three partons and we have that If we are to preserve T FR 2 , clearly eqs. (A.3) and (A.4) must be equal and so the hardest parton in the massive 3-parton event obtained after reclustering must be parton 1. We can then split the massive leg to produce a 4-parton point. The emitter may or may not be the hardest parton -these two cases must be treated separately.
We will now detail how the Φ 4 point is obtained from the Φ 3 point in the two separate cases while preserving T FR 2 . In addition, we will show in each case that taking the singular limits of the Jacobian of the transformation reproduces the limits of the FKS Jacobian and that our fixed-order subtractions in eq. (3.9) therefore survive unaltered.

A.1 Case 1: the emitter as hardest parton
We deal first with the case in which the emitter is the hardest parton, which we call the FR primary (FRp) map. In this case the emitter is p 1 and, denoting the sum of the momenta of the split pair by k, we must have that k p 1 in order to keep the thrust axis in the same direction. We must also have | k| = | p 1 |. These conditions therefore fix the sum of the three-momenta of the split pair.
In order to proceed with the actual construction of the split configuration we use the same choice of variables as in the FKS approach, and therefore adopt a similar notation: we label the momenta in Φ 3 byk 1 ,k 2 ,k 3 , with the emitter chosen ask 3 . For the Φ 4 momenta we use k 1 , . . . , k 4 with the split pair k 3 , k 4 and k = k 3 + k 4 . The recoil momenta are defined as k rec =k 1 +k 2 , k rec = k 1 + k 2 . (A.5) As discussed above, the splitting preserves the three-momentum of the emitter which constrains the momentum of the split pair and the recoil: We must now determine k 0 and define the recoil constituents such that they remain massless and sum to k rec . Since we have that | k| = | k 3 | =k 0 3 , we may obtain an expression for k 0 : Recalling the definitions of the FKS variables Φ FKS rad ≡ {ξ, y, φ}

JHEP04(2021)254
we may substitute in and solve the quadratic equation; one obtains Having determined k 0 in terms of ξ and y, we must carefully examine which solutions are kinematically allowed. The specific Φ FKS rad variables determine which (if either) of these roots are permitted. In addition to ensuring that the solutions are real, we must also have that k 0 >k 0 3 and that k 0 3 > 0. The reality constraint gives The effect of the remaining two constraints is determined by the sign of y. For y > 0, only the positive root is a valid solution (since k 0 3 < 0 for the negative root), and we have a stronger constraint on ξ: For y < 0, the positive root is valid over the range in ξ set by the real constraint, and the negative root is valid for ξ > 2k 0 3 /E cm : It remains for us to construct the four momenta of the Φ 4 event. We define 16) and the parameter We assign the recoil by defining a boost B t along k rec with magnitude β t and a constant scaling of momenta α, so that Boosting along the recoil direction and then rescaling momenta allows us to keep the recoil three-momentum fixed. We can solve for the parameters α and β t using

JHEP04(2021)254
which gives two constraints: These can be solved in terms of δ andk 0 3 to obtain For the splitting to exist, we must also ensure that 0 < β t < 1. Rewriting β t as , (A. 24) we see that for k 0 , k 0 rec >k 0 3 (i.e. for timelike k, k rec ), the condition on β t is satisfied. Specifically, we requirek The upper bound on k 0 implies the additional constraint where x = k 2 /E 2 cm . This can be translated into a constraint on ξ and y. We may then split k into k 3 and k 4 using the FKS variables in the same way as for the FKS splitting.
We can also invert the procedure and construct the projective map from Φ 4 to Φ 3 . Again we must preserve the three-momentum of the split pair, so that We need only now define the individual partons in the recoil, which we achieve by using the same boost technique as before. Defininḡ where the inverse boost is now along − k rec , as before we can obtain two constraints: Solving, we naturally recover the same α and β t as in eq. (A.22). In this case, however, the constraints 0 < β t < 1 and 0 < α < 1 are automatically satisfied so that the projection from any Φ 4 point onto a Φ 3 point is well-defined.
Substituting eq. (A.10) into eq. (A.34), one can verify that taking the limit α → 1 and expanding about ξ = 0 or y = 1 one obtains the soft or collinear limits of the FKS Jacobian, see e.g. section 5 of [52]. This means that one can use the same counterterms as those of the FKS subtraction to obtain a local cancellation of the infrared divergences.

A.2 Case 2: the emitter as a softer parton
In the case where the emitter is not the most energetic particle, the FRp map is no longer appropriate because we no longer need to keep the thrust-axis aligned with the emitter.
In this case, we can use instead the Catani-Seymour (CS) map [75]. For example, if we assume the emitter is p 2 and perform the splitting considering p 3 as the spectator parton, the hardest parton p 1 is then unchanged by the splitting and the quantity (p 2 + p 3 ) 2 is preserved. 7 This means that the thrust axis remains along p 1 and that the value of T FR 2 is also unchanged. It is left for us to show that the singular limit of the Jacobian when using the CS map with FKS variables is the same as in the FKS case up to an overall rescaling.
To describe the splitting in this case we adopt the CS notation, wherep ij is the emitter andp k is the recoil in the Φ 3 phase space. The daughters of the splitting are labelled p i and p j , while p k is the recoil in the Φ 4 phase space: p ij +p k → p i + p j + p k = p ij + p k .
(A. 35) For the case at hand, we begin by factorising Φ 4 into the 4-parton CS phase space and a radiation part where y ij,k = p i · p j p i · p j + p i · p k + p j · p k = (p i + p j ) 2 (p i + p j + p k ) 2 (A.37) and x, Ω 2 are a set of variables which parameterise the splitting p ij → p i + p j . We now wish to express the {x, cos θ} in terms of the FKS variables {ξ, y} which we do using the 7 When the emitter is instead p3 the rôles of the emitter and spectator are interchanged but the same quantities are preserved.

JHEP04(2021)254
defining relations of the CS and FKS variables: Solving, we find that The Jacobian of this transformation is given by which are exactly the soft and collinear limits of the Jacobian in the usual FKS map. Once again, the consequence is that the subtractions are precisely the same as in the FKS case and we are therefore able to use the CS mapping consistently with the FKS subtractions.

A.3 Recasting the mapping for use in the splitting functions
The mapping which we have constructed in this appendix is used not only in the fixed-order pieces of eqs. (3.2) and (3.4) but also to make the resummed spectrum fully differential in Φ 4 via the splitting function defined in eq. (3.13). We must therefore be able to construct JHEP04(2021)254 a Φ 4 phase space point using the mapping given a Φ 3 point and values of three splitting variables. 8 This can be achieved à la FKS, but in order to do so we must express our splitting variables in eq. (3.13) in terms of the FKS variables. We consider four momenta p 12 , p 3 , p 4 before the splitting producing a configuration p 1 , p 2 , p 3 , p 4 afterwards and assume the hierarchy E 1 < E 2 < E 3 < E 4 . Our splitting variables are defined to be the azimuthal angle φ, the 3-jettiness T 3 and an energy ratio z: where | p 12 | = | p 1 + p 2 |. Rewriting, we have that It remains for us to determine the quantity | p 12 | in terms of the Φ 3 momenta. This depends on whether the FRp or CS map is being used. In the FRp case, this is rather straightforward -the FRp map preserves the value of | p 12 | by construction and so we have that | p 12 | = | p 12 |. Thus, the last term in the numerator of eq. (A.56) disappears and the expression simplifies. In the CS case, matters are slightly more complicated. The map preserves the four-momentum of the most energetic particle p 4 while p 3 is the spectator; from the definition of the CS variables, we have that