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

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.


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. 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 [3,4] and bottom quarks [5,6]. 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% [7]. It is also very difficult to probe Yukawa couplings of the fermions of first two generations [8][9][10][11][12][13][14][15][16], as well as possible invisible decay channels present in new physics models. Besides, the sensitivity to Higgs self-interactions are rather weak [17][18][19][20].
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 [21], CEPC [22], CLIC [23] and FCC-ee [24]. At a lepton collider, e.g., the CEPC [22], 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 [25]. 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-nextto-leading order (N 3 LO) quantum chromodymamics (QCD) corrections to Higgs boson production via gluon fusion in the heavy top-quark limit [26,27] and to Higgs boson production via vector boson fusion within the structure function approach [28], the next-tonext-to-leading order (NNLO) corrections to Higgs boson production in association with a jet in the heavy top-quark limit [29][30][31][32], and the next-to-leading order (NLO) corrections to Higgs boson pair production with full top-quark mass dependence [33] 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 [34][35][36].
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 → bb is known up to the next-to-next-to-next-to-next-to-leading order (N 4 LO), in the limit where the mass of the bottom quark is neglected [37]. The partial width for H → gg has been calculated to the N 3 LO in the heavy top-quark limit [38]. We refer the readers to [39,40] for a complete list of relevant calculations. At a more exclusive level, the fully differential cross sections for H → bb have been calculated to NNLO in [41,42] with massless b-quarks, and in [43] 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 Hqq coupling [44]. This may help to probe possible new physics effects which modifies the light-quark Yukawa couplings. It is also suggested in [45] 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.

Formalism
In this work, we study the thrust distribution in hadronic decays of the Higgs boson. The thrust T is defined by where p i runs over the 3-momenta of the final state particles, and n is a 3-vector with unit norm. It is conventional to introduce the variable τ ≡ 1 − T , which we will use extensively later. The limit τ → 0 corresponds to the final-state configuration of two back-to-back jets, and the limit τ → 1/2 corresponds to a nearly isotropic event.
Our calculations are based on the effective Lagrangian where µ 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 a µν is the field strength tensor of the gluon field. ψ q is the light quark fields namely excluding top quark. The strong coupling α s (µ) and the Yukawa coupling y q (µ) are renormalized in the MS scheme with n f = 5 active flavors, i.e., with the top quark integrated out. The Wilson coefficient C t (m t , µ) comes from integrating out the top quark, whose perturbative expansion can be written as The coefficients C (n) t (m t , µ) have been calculated up to N 4 LO [46][47][48][49][50][51][52]. For our purpose, we need the results up to N 3 LO, which are given by where L t = ln(µ 2 /m 2 t ), 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 τ m 2 q /m 2 H . 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.2) do not interfere when computing squared-amplitudes. That is to say, for all final state X, the following interference term 0|G µν,a G a µν |X X|ψ q ψ q |0 (2.5) vanishes to all orders in the strong coupling α 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.2) do not mix with each other under renormalization. To see this, it is sufficient to show that the two matrix elements q LqR |G µν,a G a µν |0 and gg|ψ q ψ 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 , µ) and y q (µ) evolve independently under the renormalization group (RG). We have The explicit expressions for the anomalous dimensions γ t and γ y are known to third order in α 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 (µ) defined in the 5flavor scheme. This is slightly different from the massive case [49], where in addition to the flavor-decoupling in y q (µ), there is an extra Wilson coefficient C 2 (m t , µ) 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 , µ) equals unity to all orders in α 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 LqR 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 ψ q . The second and third contributions cancel each other if the renormalization constants for Yukawa coupling and quark field, Z y and Z ψ , 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 q LqR |ψ t ψ t |0 . Such an amplitude must have the form F µ (p 1 , p 2 )ū(p 1 )γ µ v(p 2 ) which is zero due to the equation-of-motion.
In summary, in the limit m q → 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.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 Hqq 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 → qqg process. However, according to the discussions around eq. (2.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 Hqq channel, and calculate higher order QCD corrections for them 1 1 Figure 2: Representative Feynman diagrams for the Hgg channel (left) and the Hqq channel (right) for the thrust distribution at LO. independently.

The leading order and next-to-leading order results
For the thrust distribution, at LO in α s , the Hgg channel contains two partonic subprocesses H → ggg and H → qqg, while the Hqq channel has only one subprocess H → qqg. The representative Feynman diagrams are depicted in figure 2. The LO result for the Hgg channel has been calculated in [53]. We calculate the LO result for Hqq channel and also reproduce the LO result for Hgg channel. The expressions of normalized thrust distribution are given by where τ ∈ (0, 1/3], µ is the renormalization scale, Γ q 0 ≡ Γ q 0 (m H ) and Γ g 0 ≡ Γ g 0 (m H ) are LO partial decay widths at the scale µ = m H , with decay width at a scale of Higgs mass, with 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 α s and y q . The renormalization constants are given by where β 0 and γ y 0 are given in Appendix A; = (4 − d)/2 is the dimensional regulator; and γ 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 [54]. This amounts to introducing an auxiliary function dΓ 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 where the integral symbol with subscript n denotes an n-body phase-space integration, and i = q, g represent the Hqq 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 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 [55] and the one-loop matrix elements from Refs. [42,56]. The Monte-Carlo integrations are performed with the Cuba library [57]. For the input parameters, we use α 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 Hqq channels, respectively. The error bands reflect the variations of the results when the renormalization  scale µ 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 τ → 1/3, due to phase space constraints. At NLO, with an additional parton emitted, the region 1/3 < τ < (1 − 1/ √ 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 τ ∼ 0.05. The correction is even more pronounced for larger τ . 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 Figure 6: Comparison between the exact results and the singular terms at LO. others, the two-loop virtual corrections to the H → qqg process, the one-loop virtual corrections to the H → qqgg process, and the tree-level H → qqggg 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-τ limit. The factorization formula can also be used to resum large logarithms appearing in small-τ 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 τ /τ in the thrust distributions. Before going into the NNLO corrections, we can extract such singular terms in the LO results from eq. (3.1). The results are given by 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-τ region. They remain as the leading contributions up to τ ∼ 0.25, where the non-singular terms contribute about 30% and 20% for Hgg and Hqq respectively.

Factorization at small τ and approximate NNLO
In this section, we briefly introduce the factorization formula at small τ , and use it to derive an approximate NNLO formula for the thrust distribution. In the τ → 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 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 pn. The factorization formula can be obtained using the language of soft-collinear effective theory (SCET) [58][59][60][61][62][63], following the derivations for the e + e − → qq process [64][65][66]. The factorized form is given by where i = q, g denote the Hqq and Hgg channels, respectively. We have defined C g t (m t , µ) ≡ C t (m t , µ) and C q t (m t , µ) ≡ 1, corresponding to the matching coefficients discussed in section 2.
The formula eq. (4.1) involves several ingredients, which we introduce in the following. The hard Wilson coefficients C i S (m H , µ) comes from integrating out the hard fluctuations at the scale µ ∼ m H . They are defined as the matching coefficient from the full theory eq. (2.2) to SCET. They can be obtained from the Hqq and Hgg form factors, which are know up to the 3-loop order [67][68][69][70][71]. From these results, the Wilson coefficients C q S and C g S can be extracted up to the next-to-next-to-next-to-leading order (N 3 LO). The jet functions J i n (p 2 n , µ) and J ī n (p 2 n , µ) describe collinear emissions along the directions of the two jets. The typical jet masses are given by p 2 n ∼ p 2 n ∼ τ m 2 H . Both the quark jet function and the gluon jet function have been calculated to the N 3 LO [72][73][74][75]. The soft functions S i (k, µ), on the other hand, describe soft emissions with typical momenta k ∼ τ m H . The quark soft function has been known analytically up to the NNLO [64,[76][77][78]. For our purpose, we also need the scale-dependent part of the N 3 LO soft function, which can be obtained through its RG equation. Note that the scale-independent part of the N 3 LO soft function was also extracted numerically, albeit with large uncertainty [74]. Up to the N 3 LO, 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. (4.1), it is straightforward to obtain the leading singular terms for the thrust distribution by expanding the formula in terms of α s (µ). Up to the NNLO, the singular part of the thrust distribution can be formally written as with i = q, g. The explicit expressions of the coefficients ∆ (n) i (τ, µ) 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 ∆ (2) i term serves as a very good approximation of the exact NLO correction up to τ ∼ 0.2. This leads us to believe that the ∆ (3) i term should also provide a good description of the NNLO correction in this region. Therefore, we define our Approximate-NNLO (NNLO A ) thrust distribution as i (τ, µ) .  In fig. 8, we show the approximate NNLO results for the Hgg and Hqq channels in the region 0.05 ≤ τ ≤ 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 τ ∈ [0.1, 0.2] to their central values at µ = m H . The slopes of the curves indicate how strong the predictions depend on the unphysical renormalization scale µ. 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 ±10% when µ 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 (4.1), and hence the leading singular term in Eq. (4.2), captures only the leading power (LP) contribution enhanced by 1/τ . 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 [81]. 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 τ , and will also extend its range of validity to larger values of τ . While this is beyond the scope of the current work, it is straightforward to perform a power expansion in τ for the LO distribution using the analytical expressions (3.1). For example, in the Hgg channel, the result is given by next-to-leading power +O(τ ) .  Figure 10: The ratios of the LP, NLP, next-to-next-to-leading power (NNLP), and nextto-next-to-next-to-leading power (NNNLP) results to the exact LO.
In Fig. 10, we study the convergence of the power expansion for the LO distributions with the central scale choice µ = 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 Hqq channel, the NLP result accidentally behaves worses than the LP one for τ > 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.

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 Hqq 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 ∼ 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 τ region 1/3 < τ < (1 − 1/ √ 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 τ 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 τ up to ∼ 0.25 at LO and NLO. We then use the formula to give an approximate NNLO prediction for the thrust distribution in the range τ ∈ [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 τ . In this region, the singular terms ln n τ /τ in Eq. (4.2) are too large at each order in α 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 τ 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.

A Ingredients relevant for LO and NLO calculations
The β-function is defined as where the coefficients are given by [79] Here the color factors are C A = N c , C F = (N 2 c − 1)/(2N c ), T F = 1/2 and n f = 5 is the number of light quarks. For β 2 and β 3 we have substituted N c = 3 to shorten the expressions.
The anomalous dimension of the Yukawa coupling y q (µ) is the same as the anomalous dimension of quark masses. It is given by with the coefficients given by [71] γ y 0 = 6C F , The anomalous dimension of C t (m t , µ) is actually not used in our calculation, since we always evaluate the coefficient at the renormalization scale µ as in eq. (3.5). We nevertheless give it here [50] γ t 0 = 0 ,

B Ingredients relevant for the leading singular terms
We expand the hard Wilson coefficients C i S in eq. (4.1) as The NLO and NNLO coefficients are given by [70,71] C g(1) We now turn to the jet function J i (s, µ) in eq. (4.1). In practice, it is more convenient to work with its Laplace transform with γ E the Euler constant. The transformed jet function can be expanded as For our purpose, we need the NLO and NNLO coefficients, as well as the L J -dependent part of the N 3 LO coefficients. They are given by [65,74,75,80] and 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, µ) is similar. We define its Laplace transform as Again, we need the expansion coefficients ofs i (L S , µ) up to the NNLO and the L Sdependent terms at N 3 LO. They can be written as [65,77] s q(1) (L S ) = C F −8L 2 S − π 2 , s q(2) (L S ) = C F n f − 32 9 L 3 S +