Hadronic decays of Higgs boson at NNLO matched with parton shower

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

$$\mathcal{L}_{eff}=\mathcal{L}_{QCD}-\frac{H}{v}C_{1}G_{\mu\nu}^{a}G^{\mu\nu,a}-\frac{H}{v}C_{2}m_{b}\bar{\psi}_{b}\psi_{b},$$

(1)

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 $\overline{\rm MS}$ running mass in the Yukawa coupling. The renormalization of QCD coupling is carried out with a $\overline{\rm MS}$ scheme with $n_{l}=5$ light flavors. Moreover, it further requires renormalization of the Wilson coefficients including mixing effects,

$$C_{1}^{0}=Z_{11}C_{1},\quad\quad C_{2}^{0}=Z_{21}C_{1}+C_{2},$$

(2)

with the one-loop renormalization constants in the $\overline{\rm MS}$ scheme given by hep-ph/9708255 ,

	$\displaystyle Z_{11}$	$\displaystyle=$	$\displaystyle 1-\frac{\alpha_{S}^{(n_{l})}(\mu)}{4\pi}S_{\epsilon}\frac{2\beta_{0}^{(n_{l})}}{\epsilon},$
	$\displaystyle Z_{21}$	$\displaystyle=$	$\displaystyle-\frac{\alpha_{S}^{(n_{l})}(\mu)}{4\pi}S_{\epsilon}\frac{12C_{F}}{\epsilon},$		(3)

with $S_{\epsilon}=(4\pi)^{\epsilon}/\Gamma(1-\epsilon)$, and $\beta_{0}^{(n_{l})}=(11C_{A}-4T_{R}n_{l})/6$. The renormalized Wilson coefficients $C_{1}$ and $C_{2}$ Kataev:1981gr ; Inami:1982xt ; Dawson:1990zj ; Djouadi:1991tka ; Kataev:1993be ; hep-ph/9405325 ; Spira:1995rr ; hep-ph/9708255 carry further logarithmic dependence on the top-quark mass. The results at NNLO are given by

	$\displaystyle C_{1}$	$\displaystyle=$	$\displaystyle-\frac{\alpha_{S}^{(n_{l})}(\mu)}{12\pi}\Big{\{}1+\frac{\alpha_{S}^{(n_{l})}(\mu)}{4\pi}11+\left(\frac{\alpha_{S}^{(n_{l})}(\mu)}{4\pi}\right)^{2}\Big{[}L_{t}(19+\frac{16}{3}n_{l})$
			$\displaystyle+\frac{2777}{18}-\frac{67}{6}n_{l}\Big{]}+\mathcal{O}(\alpha_{S}^{3})\Big{\}},$
	$\displaystyle C_{2}$	$\displaystyle=$	$\displaystyle 1+\left(\frac{\alpha_{S}^{(n_{l})}(\mu)}{4\pi}\right)^{2}\Big{[}{40\over 9}-{16\over 3}L_{t}\Big{]}+\mathcal{O}(\alpha_{S}^{3}),$		(4)

with $L_{t}=\ln(\mu^{2}/m_{t}^{2})$.

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 $\alpha_{S}^{2}$ comparing with the decay at leading order Davies:2017xsp ; Herzog:2017dtz . 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 bottom-quark 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.

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}\ll m_{b}$, the partial width can be factorized as 1408.5150 ; 1410.3165 ,

$$\frac{1}{\Gamma_{0}}\frac{d\Gamma}{dE_{X}}=H(Q^{2},m_{b}^{2},\mu)\int\mathrm{d}kS(k,\mu)\delta(E_{X}-k),$$

(5)

where $\mu$ 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 1408.5134 . The hard function at NNLO can be extracted from the two-loop form factors calculated in Refs. hep-ph/0508254 and 1712.09889 . 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 $\tau$ variable, $\equiv 1-\mathcal{T}$, with the event shape variable thrust defined as

$$\mathcal{T}=\max_{{n}}\frac{\sum_{i}|{P}_{i}\cdot{{n}}|}{\sum_{i}|{P}_{i}|},$$

(6)

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 $\tau$ in the two-jet region 0803.0342 ,

$$\frac{1}{\Gamma_{0}}\frac{d\Gamma}{d\tau}=H(Q^{2},\mu)\int\mathrm{d}p_{L}^{2}\mathrm{d}p_{R}^{2}\mathrm{d}kJ(p_{L}^{2},\mu)J(p_{R}^{2},\mu)S_{T}(k,\mu)\delta(\tau-\frac{p_{L}^{2}+p_{R}^{2}}{Q^{2}}-\frac{k}{Q}),$$

(7)

where $H(Q^{2},\mu)$ is the hard function, $J(p^{2},\mu)$ is the jet function and $S_{T}(k,\mu)$ is the soft function for thrust. $Q$ is the center-of-mass energy, $\mu$ is the renormalization scale and $\Gamma_{0}$ is the partial width at tree-level. The NNLO soft function for thrust is calculated in Schwartz:2007ib ; Fleming:2007xt ; Kelley:2011ng . The N³LO jet functions are calculated in Becher:2006qw ; Becher:2010pd ; Bruser:2018rad ; Banerjee:2018ozf . The hard function can be extracted from the calculations in Harlander:2003ai ; Gehrmann:2005pd ; Moch:2005tm ; Gehrmann:2010ue ; Gehrmann:2014vha . All these perturbative ingredients have been summarized in Ref. 1901.02253 .

From the factorization formula and the ingredients up to NNLO one can derive the cumulant of the singular distributions at NNLO as

$$\Gamma_{s}(x)\equiv\int_{0}^{x}dx\frac{d\Gamma_{s}}{dx}=\Gamma_{0}(1+\Gamma_{s}^{(1)}(x)+\Gamma_{s}^{(2)}(x)),$$

(8)

where $x\equiv\tau$ or $E_{X}/2m_{b}$ are the principal variables defined above, and $\Gamma_{s}^{(i)}$ denotes perturbative expansions of the partial width to the $i$-th power of the QCD coupling $\alpha_{S}$. $\Gamma_{s}^{(i)}$ 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. 1210.2808 . On the other hand from conventional fixed-order calculations we can derive the exact distribution in the three-jet region as

$$\frac{d\Gamma_{3j}(x)}{dx}=\Gamma_{0}\left(\frac{d\Gamma_{3j}^{(1)}(x)}{dx}+\frac{d\Gamma_{3j}^{(2)}(x)}{dx}\right).$$

(9)

In a phase space slicing method, one can obtain the NNLO partial width as

$$\Gamma^{NNLO}=\Gamma_{s}(\delta)+\int_{\delta}dx\frac{d\Gamma_{3j}(x)}{dx},$$

(10)

given that the cutoff parameter $\delta$ 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 one-loop results of three-body decays in Ref. 1501.07226 for the Higgs boson decaying into massless bottom quarks, and results in Ref. hep-ph/9707448 for decaying into gluons. In the case of decaying into massive bottom quarks, we use GoSam 2.0 1404.7096 to generate the one-loop virtual corrections for the three-body decay. Reduction of loop integrals is performed with Ninja 1203.0291 ; 1403.1229 and scalar integrals are calculated with OneLOop 0903.4665 ; 1007.4716 .

Further we can define the following damping factor,

	$\displaystyle\mathcal{D}(x)=$	$\displaystyle\exp{[\Gamma^{(1)}_{s}(x)+\Gamma^{(2)}_{s}(x)-(\Gamma_{s}^{(1)}(x))^{2}/2}]$
	$\displaystyle\equiv$	$\displaystyle\exp{[D^{(1)}(x)+D^{(2)}(x)]},$		(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. (9),

$$\Gamma^{damp.}=\int dx\mathcal{D}(x)\left(\frac{d\Gamma_{3j}(x)}{dx}-D^{(1)}(x)\Gamma_{0}\frac{d\Gamma^{(1)}_{3j}}{dx}\right).$$

(12)

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. (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 1912.09982 though they are essentially different. In $\mathcal{D}(x)$ we always fix the QCD scale to the hard scale of $\sim 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.

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 0709.2092 ,

$$d\Gamma=\bar{B}({\Phi}_{n})d{\Phi}_{n}\Big{(}\Delta({\Phi}_{n},p_{T}^{min})+\Delta({\Phi}_{n},k_{T}({\Phi}_{n+1}))\frac{R({\Phi}_{n+1})}{B({\Phi}_{n})}\mathrm{d}{\Phi}_{rad}\Big{)},$$

(13)

where $\Phi_{n}$ and $\Phi_{rad}$ are the Born phase space and the radiation phase space respectively, $p_{T}^{min}$ is a small cutoff for the transverse momentum of the POWHEG hard emission. The $\bar{B}({\Phi}_{n})$ can be thought as projection of the NLO distributions on to the Born phase space and is given by

	$\displaystyle\bar{B}({\Phi}_{n})$	$\displaystyle=[B({\Phi}_{n})+V({\Phi}_{n})]$
		$\displaystyle+\int\mathrm{d}{\Phi}_{rad}[R({\Phi}_{n+1})-C({\Phi}_{n+1})],$		(14)

for a typical subtraction method at NLO like FKS method hep-ph/9512328 , 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

$$\Delta({\Phi}_{n},p_{T})=\exp\Big{(}-\int\frac{[\mathrm{d}{\Phi}_{rad}R({\Phi}_{n+1})\theta(k_{T}({\Phi}_{n+1})-p_{T})]^{\bar{{\Phi}}_{n}={\Phi}_{n}}}{B({\Phi}_{n})}\Big{)}.$$

(15)

Further details on the POWHEG framework can be found in Ref. 0709.2092 . In standard MC programs, like PYTHIA8 Sjostrand:2014zea , there are similar approaches implementing the so-called matrix element reweighting hep-ph/0611247 . The difference is that $\bar{B}(\Phi_{n})$ is replaced by the tree-level matrix element $B(\Phi_{n})$.

Our matched calculations start with the NLO calculations for three-body decays of the Higgs boson, namely with tree-level processes being $H\rightarrow b\bar{b}g$, $q\bar{q}g$ for decays into bottom quarks or light quarks, and $H\rightarrow ggg(q\bar{q})$ for decays into gluons. One can not simply integrate over the full Born phase space $\Phi_{3}$ since the tree-level matrix elements have further soft or collinear singularities. We apply the reweighting procedure designed in Eq. (12) to suppress those singularities and render the $\bar{B}(\Phi_{3})$ kernel integrable in the full phase space, denoted as $\bar{B}^{*}(\Phi_{3})$, with the precise definition given by

	$\displaystyle\bar{B}^{*}({\Phi}_{3})$	$\displaystyle=\mathcal{\tilde{D}}(x(\Phi_{3}))[B({\Phi}_{3})(1-\tilde{D}^{(1)}(x(\Phi_{3})))+V({\Phi}_{3})]$
		$\displaystyle+\int\mathrm{d}{\Phi}_{rad}[\mathcal{\tilde{D}}(x(\Phi_{4}))R({\Phi}_{4})-\mathcal{\tilde{D}}(x(\Phi_{3}))C({\Phi}_{4})].$		(16)

Note the principal variable used in the damping factor $\mathcal{\tilde{D}}(x)$ is constructed on the basis of infrared and collinear safe observables, $E_{X}$ for the massive case and $\tau$ for the massless cases. We further split the Born phase space and the contributions to $\bar{B}^{*}$ into the resolved three-jet region and the unresolved one using the $\tau$ variable constructed from $\Phi_{3}$ for both the massive and massless cases, and a merging parameter $\tau_{m}$,

	$\displaystyle d\Gamma$	$\displaystyle=\bar{B}^{*}({\Phi}_{3})\theta(\tau(\Phi_{3})-\tau_{m})d{\Phi}_{3}\Big{(}\Delta({\Phi}_{3},p_{T}^{min})+\Delta({\Phi}_{3},k_{T}({\Phi}_{4}))\frac{R({\Phi}_{4})}{B({\Phi}_{3})}\mathrm{d}{\Phi}_{rad}\Big{)}$
		$\displaystyle+\bar{B}^{*}({\Phi}_{3})\theta(\tau_{m}-\tau(\Phi_{3}))d{\Phi}_{3}.$		(17)

We choose $\tau_{m}$ to be larger than the position where the Sudakov peak locates. Based on above distinctions we feed the POWHEG events with $\tau(\Phi_{3})>\tau_{m}$ to PYTHIA8 for parton shower as in the nominal matching calculations at NLO. That consists of our three-jet samples. For events with $\tau(\Phi_{3})<\tau_{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.¹¹1The first-order matrix elements reweighting was implemented in PYTHIA8 for the Higgs boson decaying into quarks but not for decaying into gluons. In this category after parton shower we calculate $\tau$ values for splittings from the two daughters of the Higgs boson based on the MC record. We require the larger one of the two $\tau$ values is smaller than $\tau_{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 $\tau$. In the end we will vary the merging scale $\tau_{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. 1912.09982 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. 2009.13533 . First in the GENEVA method the decays are seperated into 2, 3 and 4-jet contributions with cuts on the 2 and 3-jettiness, ${\mathcal{T}}^{{\rm cut}}_{2}$ and ${\mathcal{T}}^{{\rm cut}}_{3}$ respectively. Note in Higgs decays the 2-jettiness is equivalent to the thrust variable used here, $\tau={\mathcal{T}}_{2}/(2m_{H})$. In the 2-jet contributions, the total rate is calculated at NNLO, and is further improved with analytic resummations on ${\mathcal{T}}_{2}$. The 3-jet contributions are based on predictions of analytic resummations on ${\mathcal{T}}_{2}$ matched with NLO fixed-order calculations. The 3-jet and 4-jet contributions also include resummations on ${\mathcal{T}}_{3}$. Once matched with parton showers, by choosing a small ${\mathcal{T}}^{{\rm cut}}_{2}\,(\sim 1\,{\rm GeV})$ and a dedicated treatment on the first emission, the precisions of analytical resummations on ${\mathcal{T}}_{2}$ are maintained. In our approach we choose a much larger value of $\tau_{m}$ that is beyond the Sudakov peak. Thus the two-jet samples are dominant. Below $\tau_{m}$ the resummations are realized by the Sudakov factors in parton showers.

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. 1911.11524 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}$.

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 $\tau$ variable constructed from all final state partons. The reference values of NLO/NNLO width are taken from Ref. 1707.01044 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 $\sim 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 1707.01044 ; 1911.11524 . 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. (1) which can be calculated separately.

We move to the damped predictions from integrating over the three-body phase space. We demonstrate the effectiveness of the procedure by comparing the damped predictions with the exact NNLO predictions for different choices of the renormalization scale, $\mu=\{m_{h}/2,\,m_{h},\,2m_{h}\}$, and the QCD coupling constant, $\alpha_{S}(m_{Z})=\{0.110,\,0.1181,\,0.130\}$. Fig. 3 shows the comparison for the total partial width of the Higgs boson decaying into bottom quarks with full mass dependence, and of the Higgs boson decaying into gluons. We include results for three different choices of the damping parameter, $x_{0}=\{0.5,\,1,\,2\}$ for the case of massive bottom quarks with $x=E_{X}/2m_{b}$, and $x_{0}=\{0.05,\,0.1,\,0.2\}$ for the case of gluons with $x=\tau$. Precise definition of the damping parameter is given in Appendix B. Differences between the damped NNLO predictions and the exact NNLO predictions have been normalized to the leading order width times $\alpha_{S}^{2}$. 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}$ ($\tau$ 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 $\tau$ 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.

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 $\tau$ 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 $\tau_{m}=0.075$, the renormalization scale $\mu_{R}=m_{h}$, and the scale in parton shower $\mu_{PS}=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 $\tau$ variable can be constructed in a similar way as in the massless case.

The normalized $\tau$ distribution shows a clear dependence on the merging scale in the designed merging region around $\tau\sim 0.075$. For alternative merging scales of $\tau_{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 $\tau$ 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 $\tau$ below the merging scale. We also notice a non-negligible change of the distribution for $\tau$ 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 $\tau$ 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 $\tau$ distribution and its dependence on the three scales in a same format as Fig. 5. The genuine feature is similar to that of the Higgs boson decaying into bottom quarks and we only highlight a few differences comparing to the case of bottom quarks. The normalized distribution shows less dependence on the merging scale in full range of $\tau$. For instance, the variations are at most 10% in the merging region and are negligible for $\tau$ 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 $\tau$ 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$\sim$20% for large $\tau$ 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 $\tau$ 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-$\tau$ region, the merging scale for transition region, and the renormalization scale for large-$\tau$ region. In the following we will design our matched calculation with the default choices of the three scales as the nominal prediction 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 $\tau$ 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.

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 Sjostrand:2014zea and the Monash tunes Skands:2014pea . We show the normalized $\tau$ 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 $\tau$ values which are well known. In the tail region of $\tau$ where the distribution falls off rapidly, the hadronization and decay of hadrons tend to push the parton-level kinematics to larger $\tau$ 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 $\tau$ 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 $\tau$ 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 $\tau$. The overall picture on distribution of the total hemisphere broadening is very similar to that of the $\tau$ 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 $\tau$ distribution.

In Fig. 10 we present a comparison of the matched predictions on the normalized $\tau$ 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 fixed-order predictions have been truncated at $\tau=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 $\tau$ between 0.1 and 0.2, and 40% for $\tau\sim 0.3$. That is in agreement with the findings in Ref. 1901.02253 that higher-order logarithmic corrections enhance the distributions for massless quarks by a similar amount. The fixed-order predictions fail completely at close to the three-jet kinematic limit of $\tau=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%$\sim$10% in the distribution of bottom quarks for $\tau\gtrsim 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 $\tau$ 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.