Thrust distribution in Higgs decays at the next-to-leading order and beyond

TL;DR

This paper provides advanced predictions for the thrust distribution in Higgs boson decays, incorporating next-to-leading and approximate next-to-next-to-leading order corrections, crucial for reducing theoretical uncertainties.

Contribution

It introduces approximate NNLO corrections derived from a factorization formula, enhancing the precision of Higgs decay predictions beyond previous methods.

Findings

Large corrections observed, especially in the gluon channel

Scale variations at lowest orders underestimate higher order effects

Results are essential for controlling perturbative uncertainties

Abstract

We present predictions for the thrust distribution in hadronic decays of the Higgs boson at the next-to-leading order and the approximate next-to-next-to-leading order. The approximate NNLO corrections are derived from a factorization formula in the soft/collinear phase-space regions. We find large corrections, especially for the gluon channel. The scale variations at the lowest orders tend to underestimate the genuine higher order contributions. The results of this paper is therefore necessary to control the perturbative uncertainties of the theoretical predictions. We also discuss on possible improvements to our results, such as a soft-gluon resummation for the 2-jets limit, and an exact next-to-next-to-leading order calculation for the multi-jets region.

Equations

$$T\equiv \underset{\vec{n}}{\max }\frac{{\sum }_i|\vec{n}\cdot {\vec{p}}_i|}{{\sum }_i|{\vec{p}}_i|},$$

$${\mathcal{L}}_{\text{eff}}$$

$$\equiv \frac{{\alpha }_s(\mu )C_t(m_t,\mu )}{12\pi v}H{G}^{\mu \nu ,a}{G}_{\mu \nu }^a+\sum _q\frac{y_q(\mu )}{\sqrt{2}}H{\bar{\psi }}_q{\psi }_q,$$

$$C_t(m_t,\mu )=1+\sum _{n=1}^{\infty }{\left(\frac{{\alpha }_s(\mu )}{4\pi }\right)}^n{C}_t^{(n)}(m_t,\mu ).$$

$$C_t(m_t,\mu )$$

$$\displaystyle+\left(\frac{\alpha_{s}}{4\pi}\right)^{3}\bigg{[}L_{t}^{2}\left(209+46n_{f}-\frac{32}{9}n_{f}^{2}\right)+L_{t}\left(\frac{4834}{9}+\frac{2912}{27}n_{f}+\frac{77}{27}n_{f}^{2}\right)$$

$$\displaystyle\hskip 40.00006pt-\frac{2761331}{648}+\frac{897943\zeta_{3}}{144}+\left(\frac{58723}{324}-\frac{110779\zeta_{3}}{216}\right)n_{f}-\frac{6865}{486}n_{f}^{2}\bigg{]}\,,$$

$$⟨0⟩G^{\mu \nu ,a}{G}_{\mu \nu }^{a}X⟨X⟩{\bar{\psi }}_q{\psi }_{q}0$$

$$\frac{d}{d\ln \mu }{C}_t(m_t,\mu )={\gamma }^t({\alpha }_s(\mu ))C_t(m_t,\mu ),\frac{d}{d\ln \mu }{y}_q(\mu )={\gamma }^y({\alpha }_s(\mu ))y_q(\mu ).$$

$$\frac{1}{{\Gamma }_0^q}\frac{d{\Gamma }_{\text{LO}}^q}{d\tau }$$

$$\frac{1}{{\Gamma }_0^g}\frac{d{\Gamma }_{\text{LO}}^g}{d\tau }$$

$$\displaystyle\hskip 150.00023pt-12\big{(}1-\tau+\tau^{2}\big{)}^{2}\,\ln\frac{1-2\tau}{\tau}\bigg{]}$$

$$\displaystyle+T_{F}n_{f}\,\frac{2}{3\tau}\left[(1-3\tau)(2-15\tau+15\tau^{2})+6\tau\left(1-2\tau+2\tau^{2}\right)\,\ln\frac{1-2\tau}{\tau}\right]\Bigg{\}}\,,$$

$${\Gamma }_0^q(\mu )=\frac{y_q^2(\mu )m_H{C}_A}{16\pi },{\Gamma }_0^g(\mu )=\frac{{\alpha }_s^2(\mu )m_H^3}{72{\pi }^3{v}^2}.$$

$$Z_{{\alpha }_s}=1-\frac{{\alpha }_s}{4\pi }{e}^{-ϵ{\gamma }_E}(4\pi {)}^{ϵ}\frac{{\beta }_0}{ϵ}+\mathcal{O}({\alpha }_s^2),Z_y=1-\frac{{\alpha }_s}{4\pi }{e}^{-ϵ{\gamma }_E}(4\pi {)}^{ϵ}\frac{{\gamma }_0^y}{2ϵ}+\mathcal{O}({\alpha }_s^2),$$

$${\Gamma }_{\text{V+R}}^i={\int }_{n+1}d{\Gamma }_{\text{real}}^i+{\int }_{n}d{\Gamma }_{\text{virt}}^i={\int }_{n+1}\left(d{\Gamma }_{\text{real}}^i-d{\Gamma }_A^i\right)+{\int }_n\left(d{\Gamma }_{\text{virt}}^i+{\int }_{1}d{\Gamma }_A^i\right),$$

$${\Gamma }_{\text{NLO}}^q$$

$${\Gamma }_{\text{NLO}}^g$$

$$\frac{1}{{\Gamma }_0^q}\frac{d{\Gamma }_{\text{LO,sing}}^q}{d\tau }$$

$$\frac{1}{{\Gamma }_0^g}\frac{d{\Gamma }_{\text{LO,sing}}^g}{d\tau }$$

$$\frac{d\Gamma^{i}}{d\tau}=\Gamma_{0}^{i}(\mu)\,|C_{t}^{i}(m_{t},\mu)|^{2}\,|C_{S}^{i}(m_{H},\mu)|^{2}\int dp_{n}^{2}\,dp_{\bar{n}}^{2}\,dk\,\delta\bigg{(}\tau-\frac{p_{n}^{2}+p_{\bar{n}}^{2}}{m_{H}^{2}}-\frac{k}{m_{H}}\bigg{)}\\ \times J_{n}^{i}(p_{n}^{2},\mu)\,J_{\bar{n}}^{i}(p_{\bar{n}}^{2},\mu)\,S^{i}(k,\mu)\,,$$

$$\frac{d{\Gamma }_{\text{sing}}^i}{d\tau }={\Gamma }_0^i(\mu )\left[\frac{{\alpha }_s(\mu )}{4\pi }{\Delta }_i^{(1)}(\tau ,\mu )+{\left(\frac{{\alpha }_s(\mu )}{4\pi }\right)}^2{\Delta }_i^{(2)}(\tau ,\mu )+{\left(\frac{{\alpha }_s(\mu )}{4\pi }\right)}^3{\Delta }_i^{(3)}(\tau ,\mu )\right],$$

$$\frac{d{\Gamma }_{\text{NNLO,A}}^i}{d\tau }=\frac{d{\Gamma }_{\text{NLO}}^i}{d\tau }+{\Gamma }_0^i(\mu ){\left(\frac{{\alpha }_s(\mu )}{4\pi }\right)}^3{\Delta }_i^{(3)}(\tau ,\mu ).$$

$$\frac{1}{{\Gamma }_0^g}\frac{d{\Gamma }_{\text{LO}}^g}{d\tau }$$

$$\displaystyle+\frac{\alpha_{s}^{2}(\mu)}{\alpha_{s}^{2}(m_{H})}\frac{\alpha_{s}(\mu)}{2\pi}\,\Big{[}\underbrace{C_{A}\,\big{(}4\ln\tau+11\big{)}-T_{F}n_{f}\,\big{(}4\ln\tau+14\big{)}}_{\text{next-to-leading power}}+\mathcal{O}(\tau)\Big{]}\,.$$

$$\frac{d{\alpha }_s(\mu )}{d\ln \mu }=\beta ({\alpha }_s)=-2{\alpha }_s\sum _{n=0}{\left(\frac{{\alpha }_s}{4\pi }\right)}^{n+1}{\beta }_n,$$

$${\beta }_0$$

$${\beta }_1$$

$${\beta }_2$$

$${\beta }_3$$

$${\gamma }^y({\alpha }_s(\mu ))=-\sum _{n=0}{\left(\frac{{\alpha }_s}{4\pi }\right)}^{n+1}{\gamma }_n^y,$$

$${\gamma }_0^y$$

$${\gamma }_1^y$$

$${\gamma }_2^y$$

$$-\left(\frac{556}{27}+48{\zeta }_3\right)C_A{C}_F{n}_f-\frac{70}{27}{C}_F{n}_f^2.$$

$${\gamma }_0^t$$

$${\gamma }_1^t$$

$${\gamma }_2^t$$

$$C_S^i(m_H,\mu )=1+\sum _{n=1}^{\infty }{\left(\frac{{\alpha }_s(\mu )}{4\pi }\right)}^n{C}_S^{i(n)}(L_H),$$

Full text

1 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 Aad:2012tfa ; Chatrchyan:2012xdj . 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. Continuous operation of LHC has shown great success on refined study of the Higgs boson, for example, the recent discovery of the Higgs couplings with top quarks Sirunyan:2018hoz ; Aaboud:2018urx and bottom quarks Aaboud:2018zhk ; Sirunyan:2018kst . On the other hand, the ability of the LHC or high luminosity (HL) LHC are limited on several aspects in the study of Higgs couplings. Due to the huge SM backgrounds, the accuracy of measurements on the Higgs signal strength cannot go below the order of 5% CMS:2018qgz . It is also very difficult to probe Yukawa couplings of the fermions of first two generations Gao:2013nga ; Soreq:2016rae ; Bishara:2016jga ; Bodwin:2013gca ; Kagan:2014ila ; Zhou:2015wra ; Koenig:2015pha ; Perez:2015lra ; Chisholm:2016fzg , as well as possible invisible decay channels present in new physics models. Besides, the sensitivity to Higgs self-interactions are rather weak Goertz:2013kp ; Sirunyan:2018two ; CMS:2018ccd ; Aaboud:2018ftw .

To measure the Higgs properties with higher accuracy and to probe rare decay modes of the Higgs boson, there have been a few proposals to build a future lepton collider that can serve as a Higgs factory. These include ILC Behnke:2013xla , CEPC CEPCStudyGroup:2018ghi , CLIC Lebrun:2012hj and FCC-$ee$ Gomez-Ceballos:2013zzn . At a lepton collider, e.g., the CEPC CEPCStudyGroup:2018ghi , all decay channels of the Higgs boson can be measured in a model-independent way including possible invisible channels, and the total width can be reconstructed. The projected precision on most Higgs couplings are at the percent level thanks to the clean environment An:2018dwb . This is an order of magnitude improvement over the ability of the (HL-)LHC.

Precision experiments require equally precision theoretical predictions. To further scrutinize the SM and to look for possible new physics beyond, it is necessary to calculate higher-order corrections to the production and decay of the Higgs boson. In this respect, there have been enormous advances in recent years. For example, the next-to-next-to-next-to-leading order (N3LO) quantum chromodymamics (QCD) corrections to Higgs boson production via gluon fusion in the heavy top-quark limit Anastasiou:2015ema ; Mistlberger:2018etf and to Higgs boson production via vector boson fusion within the structure function approach Dreyer:2016oyx , the next-to-next-to-leading order (NNLO) corrections to Higgs boson production in association with a jet in the heavy top-quark limit Boughezal:2013uia ; Chen:2014gva ; Boughezal:2015dra ; Boughezal:2015aha , and the next-to-leading order (NLO) corrections to Higgs boson pair production with full top-quark mass dependence Borowka:2016ypz have been known for some time. The two-loop mixed QCD and electroweak corrections have also been calculated recently for the associated production of Higgs boson and a $Z$ boson at electron-positron colliders Gong:2016jys ; Sun:2016bel ; Chen:2018xau .

In this work, we are concerned with the hadronic decays of the Higgs boson. Namely, the final-state consists hadrons initiated by quarks and gluons. This channel is particularly interesting for a future lepton collider, since it is rather difficult to be detected at hadron colliders. This channel also provides a unique place to cleanly study non-perturbative aspects of QCD related to gluon jets. Due to the hadronic nature of this channel, the cross sections receive sizeable QCD corrections. As a result, higher order calculations for various observables in this process are highly demanded. The partial width for $H\to b\bar{b}$ is known up to the next-to-next-to-next-to-next-to-leading order (N4LO), in the limit where the mass of the bottom quark is neglected Baikov:2005rw . The partial width for $H\to gg$ has been calculated to the N3LO in the heavy top-quark limit Baikov:2006ch . We refer the readers to Denner:2011mq ; Spira:2016ztx for a complete list of relevant calculations. At a more exclusive level, the fully differential cross sections for $H\to b\bar{b}$ have been calculated to NNLO in Anastasiou:2011qx ; DelDuca:2015zqa with massless $b$-quarks, and in Bernreuther:2018ynm with massive $b$-quarks.

For hadronic decays, event shapes are a class of good observables. On one hand, they are infrared safe observables which can be theoretically calculated order-by-order in perturbation theory. On the other hand, they can be experimentally constructed from the hadron momenta without the need to specify a jet algorithm. For Higgs boson decay, in particular, one of the authors has proposed to use event shapes such as thrust, hemisphere mass and $C$ parameter to distinguish final states induced by the $Hgg$ coupling and the $Hq\bar{q}$ coupling Gao:2016jcm . This may help to probe possible new physics effects which modifies the light-quark Yukawa couplings. It is also suggested in Li:2018qiy to use jet energy profile to improve the measurement of the $Hgg$ coupling.

In this work, we investigate the thrust distribution in the hadronic decays of the Higgs boson. Such decays can be induced by the effective coupling between the Higgs boson and gluons, and can also be induced by the Yukawa coupling between the Higgs boson and quarks. We discuss these couplings in Section 2. We then calculate the leading order (LO) and the NLO contributions to the thrust distribution in Section 3. We find that the NLO corrections are rather large, and proceed to construct an approximate NNLO prediction in Section 4. We conclude in Section 5.

2 Formalism

In this work, we study the thrust distribution in hadronic decays of the Higgs boson. The thrust $T$ is defined by

[TABLE]

where $\vec{p}_{i}$ runs over the 3-momenta of the final state particles, and $\vec{n}$ is a 3-vector with unit norm. It is conventional to introduce the variable $\tau\equiv 1-T$, which we will use extensively later. The limit $\tau\to 0$ corresponds to the final-state configuration of two back-to-back jets, and the limit $\tau\to 1/2$ corresponds to a nearly isotropic event.

Our calculations are based on the effective Lagrangian

[TABLE]

where $\mu$ is the renormalization scale, $v$ is the vacuum expectation value of the Higgs field, $H$ represents the physical Higgs boson after electroweak symmetry breaking, and $G_{\mu\nu}^{a}$ is the field strength tensor of the gluon field. $\psi_{q}$ is the light quark fields namely excluding top quark. The strong coupling $\alpha_{s}(\mu)$ and the Yukawa coupling $y_{q}(\mu)$ are renormalized in the $\overline{\text{MS}}$ scheme with $n_{f}=5$ active flavors, i.e., with the top quark integrated out. The Wilson coefficient $C_{t}(m_{t},\mu)$ comes from integrating out the top quark, whose perturbative expansion can be written as

[TABLE]

The coefficients $C_{t}^{(n)}(m_{t},\mu)$ have been calculated up to N4LO Inami:1982xt ; Djouadi:1991tk ; Chetyrkin:1997iv ; Chetyrkin:1997un ; Chetyrkin:2005ia ; Schroder:2005hy ; Baikov:2016tgj . For our purpose, we need the results up to N3LO, which are given by

[TABLE]

where $L_{t}=\ln(\mu^{2}/m_{t}^{2})$, and we have set explicitly the number of colors $N_{c}=3$ to shorten the expression.

We work in the limit of vanishing light quark masses, $m_{q}=0$, while keeping the Yukawa coupling $y_{q}$ non-zero. This treatment can be justified if new physics beyond the SM leads to a different relation between $y_{q}$ and $m_{q}$ in the low energy effective theory. The zero mass limit is a good approximation as long as $\tau\gg m_{q}^{2}/m_{H}^{2}$.

The massless (chiral) limit brings about a few simplifications to our calculation, which we elaborate in the following. The first immediate effect is that the two operators in eq. (2) do not interfere when computing squared-amplitudes. That is to say, for all final state $X$, the following interference term

[TABLE]

vanishes to all orders in the strong coupling $\alpha_{s}$. This can be easily seen since the QCD interactions preserve chirality in the massless limit, while the quark operator $O_{q}$ couples two quark fields with opposite chirality. Therefore, irrelevant of the final states, it is guaranteed that one of the two matrix elements in the above interference term vanishes.

The second simplification resides in the fact that the two operators in eq. (2) do not mix with each other under renormalization. To see this, it is sufficient to show that the two matrix elements $\braket{q_{L}\bar{q}_{R}}{G^{\mu\nu,a}G_{\mu\nu}^{a}}{0}$ and $\braket{gg}{\bar{\psi}_{q}\psi_{q}}{0}$ are zero. The vanishing of both matrix elements follows from the same argument on chirality in the above. As a result of this observation, the two coefficients $C_{t}(m_{t},\mu)$ and $y_{q}(\mu)$ evolve independently under the renormalization group (RG). We have

[TABLE]

The explicit expressions for the anomalous dimensions $\gamma^{t}$ and $\gamma^{y}$ are known to third order in $\alpha_{s}$, and are collected in Appendix A.

Finally, we note that in the massless limit, the impact of integrating out the top quark on the quark operator is fully absorbed by the Yukawa coupling $y_{q}(\mu)$ defined in the 5-flavor scheme. This is slightly different from the massive case Chetyrkin:1997un , where in addition to the flavor-decoupling in $y_{q}(\mu)$, there is an extra Wilson coefficient $C_{2}(m_{t},\mu)$ coming into play. However, this coefficient arises purely from a similar effect as the operator mixing between $O_{g}$ and $O_{q}$. Since we have shown above that such a mixing is absent when $m_{q}=0$, we can conclude that $C_{2}(m_{t},\mu)$ equals unity to all orders in $\alpha_{s}$.

It is easy to demonstrate the above fact at the two-loop order (where the effect first appears). Consider the matching procedure for the $Hq_{L}\bar{q}_{R}$ amplitude. The matching coefficient comes from 3 contributions in the full theory with a closed top-quark loop: 1) diagrams where the external Higgs field is attached to the top-quark loop, e.g., the first diagram in fig. 1; 2) diagrams where the Higgs filed is attached to the light quark propagator, e.g., the second diagram in fig. 1; and 3) top-quark loop contributions to the renormalization of $y_{q}$ and $\psi_{q}$. The second and third contributions cancel each other if the renormalization constants for Yukawa coupling and quark field, $Z_{y}$ and $Z_{\psi}$, are chosen in the 5-flavor scheme. This cancellation is in fact the very definition of the “5-flavor scheme”, which is obvious if we perform the matching with the external quarks on-shell and the Higgs momentum set to zero. As for the first contribution, it can be immediately seen that the first diagram in fig. 1 vanishes in the massless limit. The absence of the first contribution can be formally proven to all orders, since it is related to the on-shell matrix element $\braket{q_{L}\bar{q}_{R}}{\bar{\psi}_{t}\psi_{t}}{0}$. Such an amplitude must have the form $F^{\mu}(p_{1},p_{2})\,\bar{u}(p_{1})\gamma_{\mu}v(p_{2})$ which is zero due to the equation-of-motion.

In summary, in the limit $m_{q}\to 0$, the hadronic decay of the Higgs boson can be classified at the parton level into two categories, induced by the gluon operator and the quark operator in eq. (2), respectively. These two operators do not mix under renormalization. In the following, we will denote the partonic processes induced by the gluon operator as the $Hgg$ channel, and those induced by the quark operator as the $Hq\bar{q}$ channel. The names might sometimes be misleading, since the two channels can have the same final state particles. For example, the two operators can both induce the $H\to q\bar{q}g$ process. However, according to the discussions around eq. (5), these two amplitudes do not interfere with each other. As a result, from the computational point of view, we can strictly separate the $Hgg$ channel and the $Hq\bar{q}$ channel, and calculate higher order QCD corrections for them independently.

3 The leading order and next-to-leading order results

For the thrust distribution, at LO in $\alpha_{s}$, the $Hgg$ channel contains two partonic subprocesses $H\to ggg$ and $H\to q\bar{q}g$, while the $Hq\bar{q}$ channel has only one subprocess $H\to q\bar{q}g$. The representative Feynman diagrams are depicted in figure 2. The LO result for the $Hgg$ channel has been calculated in Mo:2017gzp . We calculate the LO result for $Hq\bar{q}$ channel and also reproduce the LO result for $Hgg$ channel. The expressions of normalized thrust distribution are given by

[TABLE]

where $\tau\in(0,1/3]$, $\mu$ is the renormalization scale, $\Gamma_{0}^{q}\equiv\Gamma_{0}^{q}(m_{H})$ and $\Gamma_{0}^{g}\equiv\Gamma_{0}^{g}(m_{H})$ are LO partial decay widths at the scale $\mu=m_{H}$, with decay width at a scale of Higgs mass, with

[TABLE]

The NLO corrections to the thrust distribution involve both virtual gluon exchanges and real gluon emissions. The representative Feynman diagrams are show in figure 3. The virtual diagrams contain ultraviolet (UV) divergences which are removed by renormalization of the couplings $\alpha_{s}$ and $y_{q}$. The renormalization constants are given by

[TABLE]

where $\beta_{0}$ and $\gamma^{y}_{0}$ are given in Appendix A; $\epsilon=(4-d)/2$ is the dimensional regulator; and $\gamma_{E}$ is the Euler constant. After renormalization, both real and virtual corrections are separately infrared (IR) divergent, while their sum is finite. In order to implement the cancellation in a Monte-Carlo generator, we adopt the dipole-subtraction method Catani:1996vz . This amounts to introducing an auxiliary function $d\Gamma_{A}$ which has the same singular behaviors in the soft and/or collinear limits. The sum of the virtual and real corrections then be written in the form

[TABLE]

where the integral symbol with subscript $n$ denotes an $n$-body phase-space integration, and $i=q,g$ represent the $Hq\bar{q}$ and $Hgg$ channels, respectively. The two terms in the above formula are both finite, and the integration can be performed numerically. For the $Hgg$ channel, there is an extra contribution from the $C_{t}$ coefficient at NLO. Combining everything, we have the NLO decay rates as

[TABLE]

Based on the above formulas, we construct an in-house Fortran program to compute the differential decay rates. We use the real-emission matrix elements from OpenLoops Cascioli:2011va and the one-loop matrix elements from Refs. Schmidt:1997wr ; DelDuca:2015zqa . The Monte-Carlo integrations are performed with the Cuba library Hahn:2004fe . For the input parameters, we use $\alpha_{s}(m_{Z})=0.1181$, $m_{H}={125.09}$ GeV and $m_{t}={173.5}$ GeV.

In figure 4, we show the LO and NLO thrust distributions in the $Hgg$ and $Hq\bar{q}$ channels, respectively. The error bands reflect the variations of the results when the renormalization scale $\mu$ is varied up and down by a factor of 2 from the nominal choice of $m_{H}$. Note that the LO distributions approach zero when $\tau\to 1/3$, due to phase space constraints. At NLO, with an additional parton emitted, the region $1/3<\tau<(1-1/\sqrt{3})$ opens up. From this figure, one can see that the NLO corrections are rather large for both channels, indicating the bad convergence of the perturbative series. Especially for the $Hgg$ channel, the NLO differential cross section is twice the LO one at $\tau\sim 0.05$. The correction is even more pronounced for larger $\tau$. We also find that the scale uncertainties of the LO results do not overlap with the NLO ones. This indicates that the scale variation of the LO differential cross sections underestimate the theoretical uncertainties. We also show in figure 5 differential cross sections normalized to their central values. At LO the scale variations arise entirely from running of the couplings and show no dependence on kinematics. The scale variations are reduced at NLO for $\tau$ below the kinematic endpoint at LO.

To summarize, our NLO calculation reveals a few unsatisfactory features which make us believe that even higher order corrections are phenomenologically important. To obtain the full NNLO thrust distribution for, e.g., the $Hq\bar{q}$ channel, one needs to calculate, among others, the two-loop virtual corrections to the $H\to q\bar{q}g$ process, the one-loop virtual corrections to the $H\to q\bar{q}gg$ process, and the tree-level $H\to q\bar{q}ggg$ process. One also needs to combine these contributions, either analytically or numerically, in order to cancel the infrared divergences. Before get into such an involved computation, it is useful to estimate the size of the NNLO corrections. The rest of this paper will be devoted to the calculation of thrust distributions at approximate NNLO based on a factorization formula in small-$\tau$ limit. The factorization formula can also be used to resum large logarithms appearing in small-$\tau$ region, where the perturbative expansion is doomed to fail. This will be left to a future work in preparation.

The factorization formula deals with singular terms of the form $\ln^{n}\tau/\tau$ in the thrust distributions. Before going into the NNLO corrections, we can extract such singular terms in the LO results from eq. (7). The results are given by

[TABLE]

In figure 6, we compare numerically the singular terms at LO against the exact results by plotting their ratios. From there one can see the singular terms dominate at small-$\tau$ region. They remain as the leading contributions up to $\tau\sim 0.25$, where the non-singular terms contribute about 30% and 20% for $Hgg$ and $Hq\bar{q}$ respectively.

4 Factorization at small $\tau$ and approximate NNLO

In this section, we briefly introduce the factorization formula at small $\tau$, and use it to derive an approximate NNLO formula for the thrust distribution. In the $\tau\to 0$ limit, the final state hadrons form two nearly back-to-back jets in the rest frame of the Higgs boson. In this reference frame, it is convenient to choose two light-like vectors $n=(1,0,0,1)$ and $\bar{n}=(1,0,0,-1)$ to represent the directions of the two jets. The momenta of the two jets are then labeled by $p_{n}$ and $p_{\bar{n}}$. The factorization formula can be obtained using the language of soft-collinear effective theory (SCET) Bauer:2000ew ; Bauer:2000yr ; Bauer:2001yt ; Beneke:2002ph ; Beneke:2002ni ; Becher:2014oda , following the derivations for the $e^{+}e^{-}\to q\bar{q}$ process Schwartz:2007ib ; Becher:2008cf ; Bauer:2008dt . The factorized form is given by

[TABLE]

where $i=q,g$ denote the $Hq\bar{q}$ and $Hgg$ channels, respectively. We have defined $C_{t}^{g}(m_{t},\mu)\equiv C_{t}(m_{t},\mu)$ and $C_{t}^{q}(m_{t},\mu)\equiv 1$, corresponding to the matching coefficients discussed in section 2.

The formula eq. (13) involves several ingredients, which we introduce in the following. The hard Wilson coefficients $C_{S}^{i}(m_{H},\mu)$ comes from integrating out the hard fluctuations at the scale $\mu\sim m_{H}$. They are defined as the matching coefficient from the full theory eq. (2) to SCET. They can be obtained from the $Hq\bar{q}$ and $Hgg$ form factors, which are know up to the 3-loop order Harlander:2003ai ; Gehrmann:2005pd ; Moch:2005tm ; Gehrmann:2010ue ; Gehrmann:2014vha . From these results, the Wilson coefficients $C_{S}^{q}$ and $C_{S}^{g}$ can be extracted up to the next-to-next-to-next-to-leading order (N3LO). The jet functions $J_{n}^{i}(p_{n}^{2},\mu)$ and $J_{\bar{n}}^{i}(p_{\bar{n}}^{2},\mu)$ describe collinear emissions along the directions of the two jets. The typical jet masses are given by $p_{n}^{2}\sim p_{\bar{n}}^{2}\sim\tau m_{H}^{2}$. Both the quark jet function and the gluon jet function have been calculated to the N3LO Becher:2006qw ; Becher:2010pd ; Bruser:2018rad ; Banerjee:2018ozf . The soft functions $S^{i}(k,\mu)$, on the other hand, describe soft emissions with typical momenta $k\sim\tau m_{H}$. The quark soft function has been known analytically up to the NNLO Schwartz:2007ib ; Fleming:2007xt ; Kelley:2011ng ; Monni:2011gb . For our purpose, we also need the scale-dependent part of the N3LO soft function, which can be obtained through its RG equation. Note that the scale-independent part of the N3LO soft function was also extracted numerically, albeit with large uncertainty Bruser:2018rad . Up to the N3LO, the gluon soft function can be obtained from the quark one by a Casimir scaling $C_{A}/C_{F}$. The explicit expressions for the above ingredients are collected in Appendix B.

Given the factorization formula eq. (13), it is straightforward to obtain the leading singular terms for the thrust distribution by expanding the formula in terms of $\alpha_{s}(\mu)$. Up to the NNLO, the singular part of the thrust distribution can be formally written as

[TABLE]

with $i=q,g$. The explicit expressions of the coefficients $\Delta_{i}^{(n)}(\tau,\mu)$ can be found in Appendix C.

With the above formula, we can now perform a comparison similar to fig. 6 for the NLO corrections. This is shown in fig. 7. Again we see that the $\Delta_{i}^{(2)}$ term serves as a very good approximation of the exact NLO correction up to $\tau\sim 0.2$. This leads us to believe that the $\Delta_{i}^{(3)}$ term should also provide a good description of the NNLO correction in this region. Therefore, we define our Approximate-NNLO (NNLO${}_{\text{A}}$) thrust distribution as

[TABLE]

Namely, we add the NNLO singular contribution from $\Delta_{i}^{(3)}$ to the exact NLO result calculated in the previous section.

In fig. 8, we show the approximate NNLO results for the $Hgg$ and $Hq\bar{q}$ channels in the region $0.05\leq\tau\leq 0.25$. In the upper plots we show the absolute distributions, while in the lower plots we show the ratios of the differential cross sections to the LO central values. We see that the NNLO corrections are still quite large. Especially for the $Hgg$ channel, the NNLO correction can reach about 50% of the NLO differential cross section. Nevertheless, the NNLO band now marginally overlaps with the NLO one, indicating that the perturbative series starts to converge. We can therefore expect that the scale variations of the NNLO results provide a relatively honest estimate of the perturbative uncertainties due to missing higher order corrections.

To see more clearly the relative scale variations at each order, we show in fig. 9 the ratios of the integrated cross sections in the bin $\tau\in[0.1,0.2]$ to their central values at $\mu=m_{H}$. The slopes of the curves indicate how strong the predictions depend on the unphysical renormalization scale $\mu$. We observe that the scale dependence consistently decreases as we go to higher orders in perturbation theory. However, for the $Hgg$ channel, the variation of the cross section is still at the level of $\pm 10\%$ when $\mu$ is varied in the range $[m_{H}/2,2m_{H}]$, which calls for further improvement to match the precision of future $e^{+}e^{-}$ colliders.

Finally, it should be noted that the factorization formula (13), and hence the leading singular term in Eq. (14), captures only the leading power (LP) contribution enhanced by $1/\tau$. Recently, there have been a lot of efforts to calculate the next-to-leading power (NLP) corrections for various processes. In particular for thrust distribution, this has been considered in Moult:2018jjd . It will be interesting to include such higher power contributions in the approximate NNLO formula. This will improve the accuracy of the approximate formula for moderate $\tau$, and will also extend its range of validity to larger values of $\tau$. While this is beyond the scope of the current work, it is straightforward to perform a power expansion in $\tau$ for the LO distribution using the analytical expressions (7). For example, in the $Hgg$ channel, the result is given by

[TABLE]

In Fig. 10, we study the convergence of the power expansion for the LO distributions with the central scale choice $\mu=m_{H}$. We show the ratios of the first 4 orders in the power expansion to the exact LO result. It can be seen that in the $Hgg$ channel, the NLP contribution brings the approximate result much closer to the exact one. On the other hand, in the $Hq\bar{q}$ channel, the NLP result accidentally behaves worses than the LP one for $\tau>0.15$. Only by including even higher power corrections can one obtain a reliable approximation to the exact LO result. It would be interesting to see in the future whether the same conclusions can be drawn for the NLO and NNLO results.

5 Conclusion and Outlook

In this paper, we have presented predictions for the thrust distribution in hadronic decays of the Higgs boson to quarks and gluons. Our calculation is based on a low energy effective theory with $Hgg$ effective coupling and $Hq\bar{q}$ Yukawa couplings by integrating out the top quark. We have calculated the NLO QCD corrections to both channels and find large impacts on the differential cross sections. Especially for the di-gluon case, the NLO corrections can be as large as the LO results (corresponding to a $K$-factor $\sim 200\%$). The scale variations of the LO fail to predict the genuine perturbative uncertainties, and are barely reduced by the inclusion of the NLO corrections. Besides, the NLO calculation provides a new leading contribution to the large $\tau$ region $1/3<\tau<(1-1/\sqrt{3})$, in which the LO distribution vanishes.

The above observations indicate that higher order corrections beyond NLO are needed to reduce the perturbative uncertainties of theoretical predictions, in order to match the experimental precision at a future Higgs factory. As a first step, we have derived an approximate formula based on a factorization theorem valid in the small $\tau$ limit. The formula captures the leading singular terms arising from soft and collinear emissions. We show that the formula provides a reasonable approximation to the exact result for $\tau$ up to $\sim 0.25$ at LO and NLO. We then use the formula to give an approximate NNLO prediction for the thrust distribution in the range $\tau\in[0.05,0.25]$. We find that the NNLO corrections are still quite sizable and important. They also reduce the scale uncertainties significantly. Therefore, the NNLO results must be taken into account for future experiments.

A couple of improvements over the results in this work are ongoing. First of all, the fixed-order predictions presented in this work cease to be valid in the region of very small $\tau$. In this region, the singular terms $\ln^{n}\tau/\tau$ in Eq. (14) are too large at each order in $\alpha_{s}$, such that the perturbative convergence is spoiled. An all-order resummation of these singular contributions is mandatory to arrive at reliable predictions. The ingredients for such a resummation at the next-to-next-to-next-to-leading logarithmic accuracy are available, and can be readily applied. The second improvement concerns the large $\tau$ region. The approximate NNLO formula obtained in this work is not valid there. An exact NNLO calculation would be necessary to correctly describe the tail of the thrust distribution. These improvements will be presented in our forthcoming articles.

Acknowledgments

This work was supported in part by the National Natural Science Foundation of China under Grant No. 11575004 and 11635001. The work of J. Gao was sponsored by CEPC theory program and by the National Natural Science Foundation of China under the Grant No. 11875189 and No.11835005. The work of W. Ju was supported in part by the China Postdoctoral Science Foundation under Grant No. 2017M610685. Y. Gong and W. Ju would like to thank SJTU for hospitality during their stay.

Appendix A Ingredients relevant for LO and NLO calculations

The $\beta$-function is defined as

[TABLE]

where the coefficients are given by vanRitbergen:1997va

[TABLE]

Here the color factors are $C_{A}=N_{c}$, $C_{F}=(N_{c}^{2}-1)/(2N_{c})$, $T_{F}=1/2$ and $n_{f}=5$ is the number of light quarks. For $\beta_{2}$ and $\beta_{3}$ we have substituted $N_{c}=3$ to shorten the expressions.

The anomalous dimension of the Yukawa coupling $y_{q}(\mu)$ is the same as the anomalous dimension of quark masses. It is given by

[TABLE]

with the coefficients given by Gehrmann:2014vha

[TABLE]

The anomalous dimension of $C_{t}(m_{t},\mu)$ is actually not used in our calculation, since we always evaluate the coefficient at the renormalization scale $\mu$ as in eq. (11). We nevertheless give it here Chetyrkin:2005ia

[TABLE]

Appendix B Ingredients relevant for the leading singular terms

We expand the hard Wilson coefficients $C^{i}_{S}$ in eq. (13) as

[TABLE]

where

[TABLE]

The NLO and NNLO coefficients are given by Gehrmann:2010ue ; Gehrmann:2014vha

[TABLE]

and

[TABLE]

We now turn to the jet function $J^{i}(s,\mu)$ in eq. (13). In practice, it is more convenient to work with its Laplace transform

[TABLE]

where

[TABLE]

with $\gamma_{E}$ the Euler constant. The transformed jet function can be expanded as

[TABLE]

For our purpose, we need the NLO and NNLO coefficients, as well as the $L_{J}$-dependent part of the N3LO coefficients. They are given by Becher:2008cf ; Becher:2009th ; Bruser:2018rad ; Banerjee:2018ozf

[TABLE]

and

[TABLE]

The $L_{J}$-independent terms $c^{J}_{3q}$ and $c^{J}_{3g}$ are known, but are not relevant to the calculations in this work.

The case for the soft function $S^{i}(k,\mu)$ is similar. We define its Laplace transform as

[TABLE]

where

[TABLE]

Again, we need the expansion coefficients of $\tilde{s}^{i}(L_{S},\mu)$ up to the NNLO and the $L_{S}$-dependent terms at N3LO. They can be written as Becher:2008cf ; Kelley:2011ng

[TABLE]

where again the constant term $c_{3q}^{S}$ is not relevant for this work. The expression for the gluon soft function can be obtained from the quark one by a Casimir scaling.

Appendix C Leading singular terms up to NNLO

In Eq. (14), the singular parts of thrust distributions are expressed in terms of the coefficients $\Delta_{i}^{(n)}(\tau,\mu)$. Here we give their explicit expressions, where we set the number of colors $N_{c}=3$ for simplicity. We also set $\mu=m_{H}$ to get rid of the scale-dependent logarithms, and one can easily recover them through the RG equation. For the $Hgg$ channel, the results are given by

[TABLE]

where $L_{HT}=\ln(m_{H}/m_{t})$. For the $Hq\bar{q}$ process, we have

[TABLE]