Hadronic decays of Higgs boson at NNLO matched with parton shower

We present predictions for hadronic decays of the Higgs boson at next-to-next-to-leading order (NNLO) in QCD matched with parton shower based on the POWHEG framework. Those include decays into bottom quarks with full bottom-quark mass dependence, light quarks, and gluons in the heavy top quark effective theory. Our calculations describe exclusive decays of the Higgs boson with leading logarithmic accuracy in the Sudakov region and next-to-leading order (NLO) accuracy matched with parton shower in the three-jet region, with normalizations fixed to the partial width at NNLO. We estimated remaining perturbative uncertainties taking typical event shape variables as an example and demonstrated the need of future improvements on both parton shower and matrix element calculations. The calculations can be used immediately in evaluations of the physics performances of detector designs for future Higgs factories.


Introduction
The successful operation of the LHC and the ATLAS and CMS experiments have led to the discovery of the Higgs boson and completion of the standard model (SM) of particle physics [1,2]. Precision test on properties of the Higgs boson including all its couplings with standard model particles becomes one primary task of particle physics at the high energy frontier. There have been proposals for possible Higgs factories to measure the Higgs properties with higher accuracy and to probe possible new physics beyond the standard model. That includes electron-positron colliders like ILC [3], CEPC [4], CLIC [5] and FCC-ee [6], or even a muon collider with the technology advances [7].
In the SM the Higgs boson decays dominantly to hadronic final states which are hard to access at hadron colliders due to huge QCD backgrounds. That includes decay channels of a pair of bottom quarks or charm quarks, a pair of gluons via heavy-quark loops, and four quarks via electroweak gauge bosons, adding to a total decay branching fraction of about 80% [8]. The small backgrounds and low hadron multiplicities at lepton collisions make them an ideal environment to study hadronic decays of the Higgs boson. For a Higgs factory like CEPC, with the designed total luminosity we expect an experimental precision of about one percent for the hadronic decay channels [9]. On the theory side the partial widths have already been calculated to a very high accuracy with the intrinsic uncertainties being at percent level or better [10]. The partial width for H → bb is known up to the next-to-nextto-next-to-next-to-leading order (N 4 LO) in QCD, in the limit where the mass of the bottom quark is neglected [11]. The partial width for H → gg has been calculated to the N 3 LO [12] and N 4 LO [13] in QCD in the heavy top-quark limit. We refer the readers to [14,15] for a complete list of relevant calculations.
On the other hand modelings on hadronic decays at the exclusive level are equally important. For instance the full kinematic and flavor information in hadronic decays will be crucial for correcting for experimental efficiencies [16], especially in the case where the accompanied Z boson also decays hadronically [9,17]. Besides, there are also strong motivations for study of exclusive hadronic decays to test Yukawa coupling of light quarks [18], and to search for possible CP violations in Yukawa couplings [19] or exotic states beyond the SM [20][21][22][23][24]. The fully differential cross sections for H → bb have been calculated to next-to-next-to-leading order (NNLO) in [25,26] and N 3 LO in [27] for massless bottom quarks, and to NNLO in [28][29][30][31] with massive bottom quarks. There have been recent predictions for the hadronic event shapes in the decay of the Higgs boson at next-to-leading order (NLO) [32][33][34] and approximated NNLO [35]. There also exist a NNLO calculation for the Higgs boson decaying into a pair of bottom quarks plus an additional jet for massless bottom quarks [36]. To achieve a fully exclusive description at hadron level, the fixed-order calculations need to be further matched with parton shower for example as implemented in the general purpose MC event generators [37][38][39]. Matching schemes at NLO like MC@NLO [40] and POWHEG [41] have been widely used. There have also been developments on matching NNLO matrix elements with parton shower for specific processes [42][43][44][45][46][47][48][49][50][51][52][53][54][55].
Very recently there have been two implementations towards matching hadronic decays of the Higgs boson at NNLO with parton shower. Ref. [56] presents the matched results for the Higgs boson decaying into bottom quarks with the MiNLO method [57] within POWHEG framework [58] neglecting masses of the bottom quark in the matrix elements. As noticed there neglecting masses of the bottom quark leads to certain mismatches with parton shower. Ref. [59] calculates decays to massless bottom quarks as well as to gluons within the GENEVA framework [60]. In this work we have presented a matched calculation for the Higgs boson decaying into bottom quarks with full bottom-quark mass dependence within the POWHEG framework. Besides the mass effects our calculation also differs with Ref. [56] in the way of merging samples with different multiplicities. In addition we present matched results for the Higgs boson decaying into light quarks and gluons with similar approaches.
The rest of our paper is organized as follows. In Sec. 2, we present the theoretical framework used for our fixed-order calculation as well as the matched calculation. Sec. 3 provides numerical results on the partial width and distributions at various levels. Finally our summary and conclusions are presented in Sec. 4.

Effective Lagrangian
Since the top quarks only appear as internal states, one can adopt an effective theory by integrating out the top quark, the interactions can be expressed as to the leading power of inverse of the top-quark mass. We have neglected the masses of light quarks except for the bottom quark. We use the on-shell scheme for renormalization of the gluon field, the bottom quark field and mass, except for the MS running mass in the Yukawa coupling. The renormalization of QCD coupling is carried out with a MS scheme with n l = 5 light flavors. Moreover, it further requires renormalization of the Wilson coefficients including mixing effects, with the one-loop renormalization constants in the MS scheme given by [61], with S = (4π) /Γ(1− ), and β (n l ) 0 = (11C A −4T R n l )/6. The renormalized Wilson coefficients C 1 and C 2 [61][62][63][64][65][66][67][68] carry further logarithmic dependence on the top-quark mass. The results at NNLO are given by . For hadronic decays of the Higgs boson, one usually distinguishes between decaying into bottom quarks and into gluons. Such a separation is only apparent at leading order. In the following when we refer to the bottom-quark channel or gluon channel we mean decays initiated by the coupling C 2 or C 1 respectively. We have not included the cross terms of C 1 and C 2 which are formally of order α 2 S comparing with the decay at leading order [11,13]. In the calculation of decaying into bottom quarks we keep full dependence on the bottom-quark mass in the matrix elements. For the calculation of the gluon channel we set the bottomquark mass to zero. Alternatively for decaying into bottom quarks one may also set the bottom-quark mass in the matrix elements to zero thus neglect associated power corrections. We refer such a calculation as for decays to massless or light quarks in the sense that it can be applied directly to light-quark decay channels induced by various new physics beyond the standard model.

QCD factorization and fixed-order calculation
The starting point of our calculation is to reproduce the partial width at NNLO in QCD. We rely on the QCD factorization formulas as derived in either heavy-quark effective theory (HQEF) or soft-collinear effective theory (SCET). For instance, in the case where the Higgs boson decays into massive bottom quarks, we choose a principal variable that is the total radiation-energy, E X , in the rest frame of the Higgs boson. In the soft limit, E X m b , the partial width can be factorized as [69,70], where µ is the renormalization scale, m b is the bottom quark mass, and Q = m H is the typical hard scale of the process. The soft function for radiation off massive quarks has been calculated to two-loop level [71]. The hard function at NNLO can be extracted from the two-loop form factors calculated in Refs. [72] and [73]. We have used results of both groups and find full agreement on the hard function. For completeness we include the analytic expressions of the hard function and the soft function at one-loop level in Appendix A. The two-loop results are lengthy and are not shown for simplicity.
In the cases where the Higgs boson decaying into massless quarks or gluons we use τ variable, ≡ 1 − T , with the event shape variable thrust defined as and P i are three-momentum of final state partons, the unit vector n is selected to maximize the projections. Based on SCET, a factorization formula is derived for τ in the two-jet region [74], where H(Q 2 , µ) is the hard function, J(p 2 , µ) is the jet function and S T (k, µ) is the soft function for thrust. Q is the center-of-mass energy, µ is the renormalization scale and Γ 0 is the partial width at tree-level. The NNLO soft function for thrust is calculated in [75][76][77]. The N 3 LO jet functions are calculated in [78][79][80][81]. The hard function can be extracted from the calculations in [82][83][84][85][86]. All these perturbative ingredients have been summarized in Ref. [35]. From the factorization formula and the ingredients up to NNLO one can derive the cumulant of the singular distributions at NNLO as where x ≡ τ or E X /2m b are the principal variables defined above, and Γ (i) s denotes perturbative expansions of the partial width to the i-th power of the QCD coupling α S . Γ (i) s are polynomials of ln x up to order of i for massive bottom quarks and order of 2i for the massless cases. We have exchanged the total radiation energies with a dimensionless variable by taking ratio to twice of the bottom quark mass. Details on derivations of the singular distributions can be found in Ref. [87]. On the other hand from conventional fixed-order calculations we can derive the exact distribution in the three-jet region as In a phase space slicing method, one can obtain the NNLO partial width as given that the cutoff parameter δ is sufficiently small and the power corrections not accounted for can be safely neglected. In calculating NLO corrections to the three-jet rate, we use oneloop results of three-body decays in Ref. [88] for the Higgs boson decaying into massless bottom quarks, and results in Ref. [89] for decaying into gluons. In the case of decaying into massive bottom quarks, we use GoSam 2.0 [90] to generate the one-loop virtual corrections for the three-body decay. Reduction of loop integrals is performed with Ninja [91,92] and scalar integrals are calculated with OneLOop [93,94]. Further we can define the following damping factor, (2.11) which vanishes as x goes to zero and can be expanded perturbatively in the three-jet region. Based on the damping factor we can construct the following damped predictions by reweighting the fixed-order distributions in Eq. (2.9), It is not difficult to show that above integration can reproduce the partial width up to NNLO. We point out that even though we write Eq. (2.12) as an explicit integration over x the reweighting can be applied at the exclusive level for the full phase space of three-jet production, where x can be reconstructed on the event-by-event basis. We use a slightly modified version of the damped three-jet distributions as the matrix element inputs for POWHEG matching. As shown in Appendix B such a distribution agrees with the conventional NLO result in the resolved three-jet region while being normalized to the exact NNLO partial width upon integration over the full three-jet phase space. The damping factors used are inspred by the Sudakov factors in the MiNLO method [56] though they are essentially different. In D(x) we always fix the QCD scale to the hard scale of ∼ m H rather than using a dynamic scale varying with x. Also for the massive case in the exponent there only exist single logarithms of x accompanied by logarithms of the mass of the bottom quark. As will be explained in the following section, introduction of the damping factor is to ensure the decay rate integrable over the full born phase space, and realize an unitarized seperation of the two and three-jet samples in the POWHEG matching. The Sudakov resummation in our approach is completely from the shower MCs, unlike in the MiNLO method.

Parton shower matching
The POWHEG method has been widely used for matching NLO computations with parton showers. The basic idea is to single out all different singular regions in the real emissions and to parametrize the phase space with the so-called underlying Born phase space and additional radiation variables. A probability distribution fully differential in the Born phase space is calculated. Meanwhile a Sudakov factor that describes the no radiation probability as a function of radiation k T at each configuration of the Born phase space is also constructed. The Sudakov factor is calculated based on the exact matrix elements at NLO. With the two successive probability distributions, a first hard emission is generated by POWHEG from the Born phase space, and later emissions as handled by parton shower can only be allowed for k T lower than that of the first emission. Here we reproduce the POWHEG formula for the case where only final state emissions are presented [41], where Φ n and Φ rad are the Born phase space and the radiation phase space respectively, p min T is a small cutoff for the transverse momentum of the POWHEG hard emission. TheB(Φ n ) can be thought as projection of the NLO distributions on to the Born phase space and is given byB for a typical subtraction method at NLO like FKS method [95], with R and C being the QCD corrections from real emissions and the counterterms respectively. Precise definition of the Sudakov form factor is given by Further details on the POWHEG framework can be found in Ref. [41]. In standard MC programs, like PYTHIA8 [38], there are similar approaches implementing the so-called matrix element reweighting [96]. The difference is thatB(Φ n ) is replaced by the tree-level matrix element B(Φ n ). Our matched calculations start with the NLO calculations for three-body decays of the Higgs boson, namely with tree-level processes being H → bbg, qqg for decays into bottom quarks or light quarks, and H → ggg(qq) for decays into gluons. One can not simply integrate over the full Born phase space Φ 3 since the tree-level matrix elements have further soft or collinear singularities. We apply the reweighting procedure designed in Eq. (2.12) to suppress those singularities and render theB(Φ 3 ) kernel integrable in the full phase space, denoted as B * (Φ 3 ), with the precise definition given bȳ (2.16) Note the principal variable used in the damping factorD(x) is constructed on the basis of infrared and collinear safe observables, E X for the massive case and τ for the massless cases. We further split the Born phase space and the contributions toB * into the resolved three-jet region and the unresolved one using the τ variable constructed from Φ 3 for both the massive and massless cases, and a merging parameter τ m , (2.17) We choose τ m to be larger than the position where the Sudakov peak locates. Based on above distinctions we feed the POWHEG events with τ (Φ 3 ) > τ m to PYTHIA8 for parton shower as in the nominal matching calculations at NLO. That consists of our three-jet samples. For events with τ (Φ 3 ) < τ m we replace them with randomly generated events from PYTHIA8 for the same hadronic decay channel of Higgs boson, based on parton shower with native matrix elements reweigthing. 1 In this category after parton shower we calculate τ values for splittings from the two daughters of the Higgs boson based on the MC record. We require the larger one of the two τ values is smaller than τ m . Otherwise we repeat the shower. That generates our two-jet samples. With the shower veto we can minimize the overlapping of phase space probed in the three-jet samples and the two-jet samples. In our merging scheme the phase-space overlapping can not be avoided since neither the parton shower nor POWHEG emissions are ordered in τ . In the end we will vary the merging scale τ m and count the variations as part of the perturbative uncertainties. The total normalization of our matched results is fixed to the exact NNLO partial width by construction. Our results have the leading logarithmic accuracy same as PYTHIA8 in the Sudakov region, and maintain the NLO matched with parton shower accuracy in region with three hard jets. We note in the case of massless quarks, the MiNLO method in Ref. [56] provides a smooth matching of the two-jet and three-jet samples without introducing an explicit merging scale. Distributions in the Sudakov region have been reweighted according to results from analytical resummations. Our calculations for the Higgs boson decaying into massless quarks and gluons can be compared with those in Ref. [59]. First in the GENEVA method the decays are seperated into 2, 3 and 4-jet contributions with cuts on the 2 and 3-jettiness, T cut 2 and T cut 3 respectively. Note in Higgs decays the 2-jettiness is equivalent to the thrust variable used here, τ = T 2 /(2m H ). In the 2-jet contributions, the total rate is calculated at NNLO, and is further improved with analytic resummations on T 2 . The 3-jet contributions are based on predictions of analytic resummations on T 2 matched with NLO fixed-order calculations. The 3-jet and 4-jet contributions also include resummations on T 3 . Once matched with parton showers, by choosing a small T cut 2 (∼ 1 GeV) and a dedicated treatment on the first emission, the precisions of analytical resummations on T 2 are maintained. In our approach we choose a much larger value of τ m that is beyond the Sudakov peak. Thus the two-jet samples are dominant. Below τ m the resummations are realized by the Sudakov factors in parton showers.

Numerical results
In the numerical calculations we set the mass of the Higgs boson m H = 125.09 GeV, the vacuum expectation value v = 246.22 GeV, and α S (M Z ) = 0.1181 [97]. The pole masses of the top quark and bottom quark are set to 172.5 GeV and 4.78 GeV respectively. We use the 3-loop running of α S and of the bottom quark mass in a MS scheme with 5 light flavors [98]. We set the renormalization scale to the mass of the Higgs boson unless specified. We use POWHEG-Box-v2 [58] for matching of matrix elements and PYTHIA8.2 [38] for final state QCD showers and hadronizations with the Monash tune [99]. We turn off multiple particle interactions and decays of B hadrons for simplicity. The generated MC event samples are analysed with Rivet program [100].

Partial width at fixed order
We start with results on the total partial decay width of the Higgs boson decaying into bottom quarks and gluons up to NNLO in QCD. For bottom quarks we show separately results with full bottom quark mass dependence and results neglecting bottom quark mass except in the Yukawa coupling. In Fig. 1 we plot dependence of the NLO and NNLO width on the cutoff of the phase space slicing variable, which is chosen to be the total radiation energies normalized to twice of the Higgs boson mass, for the case of massive bottom quarks. All predictions have been normalized to the LO partial width for simplicity while the absolute values will be summarized in later section. The scattering points with error bar represent our calculations with MC statistical uncertainties. They clearly approach the genuine NLO/NNLO predictions taken from Ref. [30] that are represented by the horizontal lines. In general the differences due to the neglected power corrections are within the MC errors which are about one per mille, if the cutoff is below 10 −3 .  Figure 1: Dependence of the total partial width calculated at NLO and NNLO on the phase space slicing parameter for the Higgs boson decaying into bottom quarks with full bottom-quark mass dependence. The phase space slicing variable is chosen to be the total radiation energies normalized to twice of the Higgs boson mass. The horizontal lines represent corresponding reference predictions from the literature, see text for details. All predictions are normalized to the LO partial width.
In Fig. 2 we present similar results for the Higgs boson decaying into massless bottom quarks and gluons. In this case the phase space slicing variable is the τ variable constructed  [101] for bottom quarks and gluons. Again we see excellent agreements between our numerical results and the analytical results in the literature. We find a cutoff parameter of the order ∼ 10 −2 is good enough to ensure the power corrections being negligible especially for the gluon channel. We can further compare the K-factors of the bottom-quark channel in the calculations with full mass dependence and using massless approximation. At both NNLO and NLO the K-factors with respect to LO are larger by about half percent in the calculation with full mass dependence. We summarize the partial decay width at NNLO for the bottom-quark channel and gluon channel in Table. 1. The results are calculated with our numerical program and further cross checked with available calculations in the literatures [30,101]. We also include values of the running QCD coupling and the bottom quark mass at the chosen renormalization scales. With the full bottom-quark mass dependence the NNLO width is reduced by about half percent due to suppressions of the phase space. Here and below we have not included interference contributions between the two operators in Eq. (2.1) which can be calculated separately.
We move to the damped predictions from integrating over the three-body phase space.   Differences between the damped NNLO predictions and the exact NNLO predictions have been normalized to the leading order width times α 2 S . The damped calculations reproduce the exact results to a precision of one per mille in general. In Fig. 4 we further compare the exact and damped predictions on the distribution of radiation energies E X (τ variable) for decaying into bottom quarks (gluons). We note that for a NNLO calculation of the total partial width it only predicts above distributions at NLO. The distributions from exact calculations diverge when E X or τ goes to zero due to the soft and collinear singularities. The damped predictions render a strong suppression in the soft and collinear region and ensure the distributions being smoothly integrable in the full phase space.
In the following calculations when matching the NNLO calculations with parton shower we use a damping parameter of x 0 = 0.5(0.05) for the massive (massless) case. That ensures the total partial width normalized automatically to the NNLO fixed-order predictions and also in resolved three-jet region the underlying parton-level results agree with the NLO fixed-order predictions.

Predictions at parton level
We proceed with predictions of matching the NNLO fixed-order calculations with parton shower using the POWHEG method and the merging scheme outlined in earlier sections. We have checked the predictions for a variety of event shape variables but only show here results for τ variable due to limited space. In Fig. 5 we show our nominal predictions for the Higgs boson decaying into massive bottom quarks with the default choice of the merging scale τ m = 0.075, the renormalization scale µ R = m h , and the scale in parton shower µ P S = p T . From left to right, we show separately the variations due to changes of above three scales in turn. In each plot the upper panel shows the normalized distribution and the lower panel shows ratios of different predictions. We notice that even there are massive particles in the final state, the thrust or τ variable can be constructed in a similar way as in the massless case. The normalized τ distribution shows a clear dependence on the merging scale in the designed merging region around τ ∼ 0.075. For alternative merging scales of τ m = 0.05, 0.10, the normalized distribution can change by 15% at most because of the differences of the three-jet and two-jet samples in the merging region. However, when going away from the merging region, the distributions are less affected. The change of renormalization scale shows its impact in the region of τ larger than the merging scale. Variation of the renormalization scale by a factor of two leads to a change of the distribution by about 15% in that region. In the two-jet region the normalized distributions vary by only a few percents in order to maintain the total partial width to the NNLO values at the associated scales.
The variation of shower scale can shift the height and position of the peak of the distribution significantly that are determined by the two-jet sample. Changing the square of the shower scale by a factor of two serve as an estimate of the uncertainties due to QCD resummation in the region of τ below the merging scale. We also notice a non-negligible change of the distribution for τ well above the merging scale. That is because in the three-jet sample after the first hard emission from POWHEG the subsequent radiations are handled by parton shower and are thus affected by the choice of shower scale. The impact of these subsequent emissions is especially pronounced in the tail region of τ where fixed-order predictions are limited due to phase space constraints.
For our NNLO calculation of the Higgs boson decaying into massless bottom quarks, in principle we can also match it with parton shower if neglecting the mismatch of the bottom quark mass used in matrix element calculations and that used in parton shower. Here we instead apply our massless calculations to the Higgs boson decaying into light quarks, for example, strange quarks or even up and down quarks. In such cases both quark masses in the matrix element calculations and in the parton showers can be safely neglected. We show results for the Higgs boson decaying into up quarks in below since the matched predictions will be independent of the flavor of quarks at the parton level. In Fig. 6 we plot the τ distribution and its dependence on the three scales in a same format as at most 10% in the merging region and are negligible for τ above the merging scales. The distribution also show a slightly smaller variation on the renormalization scale possibly due to mass dependence in the POWHEG matching method.
Lastly we show our matched predictions for the Higgs boson decaying into gluons in Fig. 7. The τ distribution peaks at a larger value and has much broader shape due to stronger radiation strength of gluons. Thus the normalization of the three-jet sample is relatively larger than in the case of quarks since we use the same merging scales in the gluon channel and the quark channel by default. The size of variations due to the merging scale is within 10% but is extended to a broader region comparing to the quark channels. The variations on renormalization scale are about 10∼20% for large τ values. The region dominated by the two-jet sample also show a larger dependence on the renormalization scale comparing to the quark channels in order to balance the variation of normalizations of the three-jet sample. The distribution shows a much larger dependence on the scale in parton shower as expected. The variations can reach 10% even for τ around 0.2 indicating possible large contributions from missing higher-order matrix elements.
To conclude we have used variations of three scales, namely the merging scale, renormalization scale, and shower scale to estimate the perturbative uncertainties of our matched predictions. We find the variations well cover the full kinematic region with each scale dominates in its designed region, i.e., the shower scale for small-τ region, the merging scale for transition region, and the renormalization scale for large-τ region. In the following we will design our matched calculation with the default choices of the three scales as the nominal pre- diction and take the quadratic sum of variations induced by the three scales as the estimate of full perturbative uncertainties. As mentioned earlier our matched calculations are fully exclusive though we have chosen τ as the principal variable for merging. We have checked explicitly our predictions for distributions of the total hemisphere broadening or the heavy hemisphere mass and found similar behavior for dependence of the distributions on the three scales. Our prediction has the same logarithmic accuracy as the usual parton shower MC in the Sudakov region while it maintains the NLO matched with parton shower accuracy in the region with three hard jets and the NNLO accuracy for the total partial width. Comparisons of our parton-level predictions with results of analytic resummations are included in Appendix C that indicates the logarithmic accuracies of shower MCs are well preserved.

Predictions at hadron level
In the final step we further include hadronization for our matched predictions which is mandatory for any phenomenological study. We use the default hadronization model in PYTHIA8.2 [38] and the Monash tunes [99]. We show the normalized τ distributions for the Higgs boson decaying into bottom quarks, up quarks and gluons in Fig. 8. We compare the nominal predictions at parton level and at hadron level with colored bands indicate the full perturbative uncertainties. Predictions from PYTHIA8 at hadron level are also included for comparison. In each plot the upper panel shows the normalized distributions and the lower panel shows ratios with respect to the nominal prediction at hadron level.
We found the hadronization corrections can be as large as tens percents when close to the  peak region and shift the Sudakov peak to larger τ values which are well known. In the tail region of τ where the distribution falls off rapidly, the hadronization and decay of hadrons tend to push the parton-level kinematics to larger τ values and thus increase the distribution significantly. In between the peak and tail region the hadronization corrections are at the level of 1% for the quark channels and 5% for the gluon channel. The hadronization corrections are almost the same for the decays into bottom quarks and into up quarks with corrections to the former being slightly smaller due to the presence of the bottom quark mass. The size of the perturbative uncertainties are not affected by hadronization. They range between 10% to 30% from the peak region to the tail region, and are even larger when τ is to the left of the peak region. It is interesting that the native hadron-level predictions from PYTHIA8 lie within our uncertainty bands for the full kinematic region considered and for all three decay channels. However they show a harder spectrum in general comparing to our nominal predictions.
In Fig. 9 we show similar results for distributions of the total hemisphere broadening B T . The distribution peaks at a value of B T that is almost twice of the value of the τ distribution for both the quark channel and the gluon channel due to different dependence on soft and collinear emissions. In the tail region the B T distribution falls even more rapidly than τ . The overall picture on distribution of the total hemisphere broadening is very similar to that of the τ distribution including for the size of perturbative uncertainties and the size of the hadronization corrections. One difference is that the hadronization corrections tend to decrease the distribution in the tail region which is opposite to the case of the τ distribution.
In Fig. 10 we present a comparison of the matched predictions on the normalized τ distribution to the fixed-order predictions (NNLO for the total partial width or NLO for the distribution) for the decays to bottom quarks and to up quarks respectively. The fixedorder predictions have been truncated at τ = 0.1. By checking on the ratios of the fixed-order predictions to the matched predictions at parton level it is evident that by matching to parton shower the distributions are enhanced by about 20% for τ between 0.1 and 0.2, and 40% for τ ∼ 0.3. That is in agreement with the findings in Ref. [35] that higher-order logarithmic corrections enhance the distributions for massless quarks by a similar amount. The fixedorder predictions fail completely at close to the three-jet kinematic limit of τ = 1/3 where multi-parton emissions play an essential role. In the lower panel of Fig. 10 we plot ratios of the distributions for up quarks to that for bottom quarks with different calculations. We observe a suppression of 5%∼10% in the distribution of bottom quarks for τ 0.1 for both the fixed-order calculation and the two matched calculations at parton and hadron level. On the other hand the distribution of bottom quarks peaks at a slightly smaller τ value than that of up quarks. In the merging region it shows the distribution of bottom quarks is 10% higher than that of up quarks which could be an artifact due to merging uncertainties.

Summary
In summary we have presented calculations of hadronic decays of the Higgs boson at NNLO matched with parton shower, including the decays to bottom quarks, light quarks and gluons. We use phase space slicing methods to reproduce the NNLO partial widths. The matched calculations are based on the POWHEG framework together with PYTHIA8. We keep the full dependence on the bottom quark mass in the matrix elements thus ensuring a consistent matching with parton shower. The calculations provide descriptions on exclusive decays of the Higgs boson at leading logarithmic accuracy in the Sudakov region and NLO matched with parton shower accuracy in the three-jet region, with normalizations fixed to the NNLO partial width. As examples we show detailed results on distributions of event shape variables like thrust and total jet broadening. We found parton shower plays important roles even in the three-jet region where it increases the NLO distributions by 20% in general and even more when close to the boundary of the phase space of fixed orders. By comparing distributions for bottom quarks and light quarks we found the normalized distributions are suppressed by 5 to 10% in the three-jet region for bottom quarks. Effects of hadronizations can be significant in the Sudakov region or far tail region, and are slightly different for bottom quarks and light quarks. We estimate the total perturbative uncertainties on the distributions and found they are more than 10% in the bulk region of the distributions especially for the decay to gluons. Predictions at such an accuracy will not be sufficient for measurements at future Higgs factories. However, they can certainly help with evaluations of the physics performances of the detector designs which are under investigation. On another hand, with recent progresses on developments of parton shower MCs [102][103][104][105][106][107][108] and given the available matrix elements at higher orders [36], we expect further improvements on modeling of hadronic decays of the Higgs boson in the near future.

A Soft and hard functions at one-loop
For QCD radiations off two massive quarks the soft function can be calculated using the eikonal approximation and can be expanded order-by-order in α S . The Laplace transformed soft function at one-loop can be expressed as [71] S(z, L) where L = ln(κ/µ), and γ E is the Euler-Mascheroni constant. We have replaced the bottom quark mass with the dimensionless variable, The cusp anomalous dimension is given by [109][110][111], and the constant piece of the soft function is The hard function can be obtained from the form factor of the Hbb vertexF renormalized in the MS scheme where Z H removes fully the remaining infrared divergences in the form factor. At one-loop the hard function can be expressed as where the renormalization scale dependence in the second line cancels with that from the soft function and the third line is due to conversion from the on-shell scheme to the MS scheme for the bottom-quark Yukawa coupling. The two-loop hard function can be obtained in a similar way with the two-loop form factor calculated in [72,73].

B Damping factor
In this appendix we include further details on the damped predictions used. We can rewrite it as where we have simply added and subtracted a same term from the original definition in Eq. (2.12). Now the integrand of the first line becomes a total derivative. In the second line each piece in the brackets is now regularized for x goes to zero and the damping factor in front can be expanded order-by-order. To be specific, we find In deriving the last step one can imagine adding an arbitrarily small cutoff as the lower limit of the integration, and the formula simply goes back to Eq. (2.10) for the phase-space slicing method. In all our derivations we have fixed the renormalization scale rather than using a dynamic scale varying with x. Before passing the reweighted three-jet matrix elements to POWHEG we would like to further eliminate those spurious O(α 3 S ) terms in Eq. (B.2). In addition we require distributions in the resolved three-jet region to be less modified by the reweighting. These can be achieved by using a slightly modified damping factor as below, where x 0 is a variable to control the functional range of the reweighting, and x is defined as is a constant term tuned numerically to eliminate spurious terms of O(α 3 S ). It depends on the preselected value of x 0 as well as on the renormalization scale.

C Comparison to analytic resummations
We compare our predictions for the Higgs boson decaying into massless quarks or gluons at the parton level to results from analytic resummations. The factorization formulas and renormalization group evolutions can be found in Ref. [74] for thrust and [112] for total hemisphere broadening. We do not reproduce details of the analytic resummations but only show numerical results at leading logarithmic (LL) and next-to-leading logarithmic (NLL) accuracy with a choice of α S (M Z ) = 0.118. In Fig. 11 we show the cumulate width as a function of τ normalized to its value at the 3-jet kinematic limit, τ 0 = 1/3, for the Higgs boson decaying into massless quarks and gluons respectively. The colored bands indicate the variations when changing the square of the shower scale by a factor of two. We find that when close to the Sudakov region, the MC predictions in general lie between the LL and NLL results, and are more closer to the NLL ones. Fig. 12 shows similar comparisons for B T . We again find good agreements between the MC predictions and the analytic resummations. We conclude that the usual logarithmic accuracy in parton shower MCs are well preserved in our matched predictions.