Four-jet event shapes in hadronic Higgs decays

We present next-to-leading order perturbative QCD predictions for four-jet-like event-shape observables in hadronic Higgs decays. To this end, we take into account two Higgs-decay categories: involving either the Yukawa-induced decay to a $b\bar{b}$ pair or the loop-induced decay to two gluons via an effective Higgs-gluon-gluon coupling. We present results for distributions related to the event-shape variables thrust minor, light-hemisphere mass, narrow jet broadening, $D$-parameter, and Durham four-to-three-jet transition variable. For each of these observables we study the impact of higher-order corrections and compare their size and shape in the two Higgs-decay categories. We find large NLO corrections with a visible shape difference between the two decay modes, leading to a significant shift of the peak in distributions related to the $H\to gg$ decay mode.


Introduction
Future lepton colliders, such as the FCC-ee [1], the CEPC [2], or the ILC [3], are projected to operate as so-called "Higgs factories", producing an unprecedented amount of events which contain a Higgs boson in the final state.In leptonic collisions, Higgs bosons are predominantly produced via the Higgsstrahlung process e + e − → ZH at low centre-of-mass energies (≲ 450 GeV) and via the vector-boson-fusion process e + e − → ℓ lH at high centreof-mass energies (≳ 450 GeV).The clean experimental environment of leptonic collisions, in which contamination from QCD backgrounds such as initial-state radiation and multiparton interactions is absent, will allow for precise measurements of Higgs-boson properties, such as its branching ratios and total width.In particular, it will become possible to have access to so-far unobserved subleading hadronic decay channels such as the decay to gluons or c-quark pairs, which are currently inaccessible in hadron-collider environments, where only the dominant H → b b decay was observed to date [4,5].
In hadronic Higgs decays, three-jet-like event shapes such as the thrust and energyenergy correlators have been discussed as discriminators of fermionic and gluonic decay channels [47][48][49][50].In [51], the full set of the six "classical" three-jet-like event-shape observables related to the event-shape variables thrust, C-parameter, heavy-hemisphere mass, narrow and wide jet broadening, and the Durham three-to-two-jet transition variable (called y 23 ) have been calculated at NLO QCD for hadronic Higgs decays.The results obtained for the two Higgs-decay categories were compared.It was shown that for all event shapes the size of NLO corrections were larger in the Hgg category.The distributions related to the three-jet event-shape variable thrust and y 23 showed the most striking shape differences between Higgs decay categories.It was therefore argued that these event shapes could act as discriminators between the two Higgs-decay categories highlighting where the Hgg category could be enhanced.Recently, a novel method to determine branching ratios in hadronic Higgs decays via fractional energy correlators has been proposed in [52].To the best of our knowledge, four-jet event shapes have not yet been computed in hadronic Higgs decay processes.
With this paper, we aim at delivering, for the first time, theoretical predictions for four-jet-like event shapes in hadronic Higgs decays including perturbative QCD corrections up to the NLO level.The event-shape distributions computed here are related to the eventshape variables narrow jet broadening, D-parameter, light-hemisphere mass, thrust minor, and the four-to-three-jet transition variable in the Durham algorithm (called y 34 ).It is the purpose of this work to compare the impact of higher-order corrections in terms of size, shape, and perturbative stability of these four-jet-like event-shape distributions.
Furthermore, besides its direct applicability in Higgs phenomenology at lepton colliders, our calculation provides necessary ingredients needed for the NNLO and N 3 LO calculations of jet observables in three-jet and two-jet final states of hadronic Higgs decays, which are currently only known for decays to b-quarks [53,54].
The paper is structured as follows.In section 2, we discuss the ingredients of our computation for both Higgs-decay categories up to NLO QCD.Our predictions for the five event-shape distributions are presented in section 3. We summarise our findings and give an outlook on future work in section 4.
2 Hadronic Higgs decays up to NLO QCD Hadronic Higgs decays proceed via two classes of processes; either involving the decay to a quark-antiquark pair, H → qq, or the decay to two gluons, H → gg.Within the scope of this work, we consider a five-flavour scheme and assume all light quarks, including the b-quark, to be massless.In order to allow the Higgs boson to decay to a b b pair, we keep a non-vanishing Yukawa coupling for the b-quark.The coupling between gluons and the Higgs is considered in an effective theory in the limit of an infinitely heavy top-quark.
Based on this, we will classify parton-level processes induced by the two Born-level decay modes into two categories.
In the first class of processes, the Higgs decays to a b-quark pair, mediated by the Yukawa coupling y b between the Higgs and the bottom quark, and the decay is computed g q H q q q q Figure 3: Tree-level Feynman diagrams of four-parton processes in the H → gg decay category contributing to our calculation at LO. Here, q and q ′ denote any flavour u, d, s, c, b, possibly with q ′ = q but q ′ ̸ = q in general.In the second class of processes, the decay proceeds via a top-quark loop.In the limit of an infinitely large top-quark mass, where the top quarks decouple, we compute these decays in an effective theory with a direct coupling of the Higgs field to the gluon field-strength tensor.In this second category, which we shall refer to as the Hgg category, the interaction is mediated by an effective Hgg vertex, which is represented as a crossed dot in the two-parton decay diagram shown in the right-hand part of fig. 1.A summary of the four-and five-parton processes contributing to this decay category, entering our computation at LO and NLO, are presented in the second column of table 1. Feynman diagrams for the tree-level four-parton processes contributing at LO are shown in fig. 3.
Following the nomenclature presented in [51], the preceding discussion can be cast into an effective Lagrangian containing both decay categories as, Here, the effective Higgs-gluon-gluon coupling proportional to α s is given by where we define the Higgs vacuum expectation value v and the top-quark Wilson coefficient The b-quark Yukawa coupling on the other hand is given by which is independent of the top-quark mass M t .Further, the Yukawa mass mb runs with µ R , with the the running of the Yukawa coupling y b with µ R taken into account via the results of [62].
Finally, to conclude some general remarks, we wish to highlight the importance for the present study that the kinematical mass of the b-quark be vanishing.For kinematically massless b-quarks, the operators in the effective Lagrangian given in eq.(2.1) do not interfere or mix under renormalisation, cf.e.g.[48].In particular, this means that the product of Yukawa-induced amplitudes with HEFT-induced amplitudes vanishes to all orders in α s when computing squared matrix elements.Moreover, it implies that the Wilson coefficient C(M t , µ R ) in eq. ( 2.4) and the Yukawa coupling y b (µ R ) in eq. ( 2.3) evolve independently under the renormalisation group.As a result, it is possible to consider two distinct decay categories, H → b b and H → gg, and calculate theoretical predictions and higher-order corrections for each of them separately, as done previously in [51] for the case of three-jet-like event shapes.The effect of interferences between the Yukawa-induced decay H → b b and the HEFT-induced H → gg to order O α 3 s in a framework with massive b-quarks and c-quarks has been studied in [63].
In the remainder of this section we discuss in detail the ingredients of our computation.
We split the discussion in the following parts.In section 2.1, we present the general framework enabling us to compute Higgs decay observables up to NLO, while in section 2.2 and section 2.3 we discuss details of the ingredients and their numerical implementation in the two Higgs-decay categories.In section 2.4, we provide a summary of the checks performed to validate our results.

General framework
Given an infrared-safe observable O, the differential four-jet decay rate of a colour-singlet resonance of mass M , like the Standard-Model scalar Higgs boson, can be written at each perturbative order up to NLO level (i.e., including corrections up to the third order in the strong coupling α s ) as where d B and d C denote the differential LO and NLO coefficients, respectively.Subject to the order of the calculation, n, the differential decay rate in eq.(2.5) is normalised to the LO (n = 0) or NLO (n = 1) partial width, Γ (1) or Γ (0) , respectively.
Schematically, the LO coefficient B can be determined as while the NLO coefficient C is obtained as Here, dΓ R and dΓ V denote the real and virtual (one-loop) correction differential in the four-parton and five-parton phase space, respectively.The real and virtual subtraction terms dΓ S NLO and dΓ T NLO , on the other hand, ensure that the real and virtual corrections are separately infrared finite and make them amenable to numerical integration.The notation Φ 4 (Φ 5 ) represents a kinematic mapping from the five-parton to the four-parton phase space, specific to the subtraction term dΓ S NLO .The choice of subtraction terms is in principle arbitrary and subject only to the requirements that the virtual subtraction term dΓ T NLO cancels all explicit poles in dΓ V , the real subtraction term dΓ S NLO cancels all implicit singularities in dΓ R , and the net contribution of the subtraction terms to the decay (2.8) The last requirement is equivalent to demanding that dΓ T NLO be the integral of dΓ S NLO over the respective one-particle branching phase space.To compute our predictions for the four-jet like hadronic event shapes in Higgs decays, we here rely on the antenna-subtraction framework [40,64,65] to construct all subtraction terms and implement our calculation in the publicly available2 EERAD3 framework [14].This code has previously been used to study event shapes [12,13] and jet distributions [66] in e + e − → 3j at NNLO and was recently extended to include hadronic Higgs decays to three jets at NLO [51].Our new implementation builds upon the latter and is done in a flexible manner, utilising the existing infrastructure, such as phase-space generators of the e + e − → 3j calculation at NNLO.In particular, the real-radiation contributions proportional to α 2 s of the three-jet computation of Higgs-decay observables enter our current calculation of four-jet-like observables at the Born level.We implement all matrix elements and subtraction terms in analytic form, enabling a fast and numerically stable evaluation of the perturbative coefficients up to NLO level.
Four-particle tree-level matrix elements in the Hb b category are taken from the NLO H → b bj process in [51], which have been calculated explicitly using FORM [67].Treelevel five-parton matrix elements are taken from the N 3 LO H → b b and NNLO H → b bj calculation in [53,54], in turn calculated using BCFW recursion relations [68].Similarly, one-loop four-parton amplitudes are taken from the same calculation [53,54].These have been derived analytically by use of the generalised unitarity approach [69], using quadruple cuts for box coefficients [70], triple cuts for triangle coefficients [71], double cuts for bubble coefficients [72], and d-dimensional unitarity techniques for the rational pieces [73].
As per eq.(2.5), the differential four-jet rate depends on the partial two-parton width.
At NLO, the inclusive width of the H → b b decay reads where the LO partial width is given by (2.10) The running of the Yukawa coupling y b with µ R present in the partial width is taken into account via the results of [62].

Effective-theory contributions
In the Hgg category, we take the four-parton tree-level matrix elements from the NLO H → 3j calculation in the Hgg category presented in [51].Five-parton tree-level and fourparton one-loop amplitudes are obtained by crossing the ones used in the pp → H + 1j NNLO calculation in NNLOJET [74,75].The five-parton tree-level amplitudes are based on the results presented in [76][77][78][79], while the four-parton one-loop amplitudes are based on the results presented in [76,[79][80][81][82][83][84][85][86].The four-parton decay rate further receives contributions from the O(α s ) expansion of the top-quark Wilson coefficient C(M t , µ R ).We implement these as a finite contribution to the virtual correction, As per eq.(2.5), the differential four-jet rate depends on the partial two-parton width.
At NLO, the inclusive width of the H → gg decay reads Γ where the LO partial width is given by (2.13)

Validation
We have performed numerical checks of all matrix elements on a point-by-point basis against automated tools.Tree-level four-and five-parton matrix elements are tested against results generated with MADGRAPH5 [87,88], whereas four-parton one-loop matrix elements are validated using OPENLOOPS 2 [89].In all cases, we have found excellent agreement between our analytic expressions and the auto-generated results.
Real subtraction terms dΓ S NLO have been tested numerically using so-called "spike tests", first introduced in the context of the antenna-subtraction framework in [90] and applied to the NNLO di-jet calculation in [91].To this end, trajectories into singular limits are generated using Rambo [92] and the ratio dΓ S NLO /dΓ R is evaluated on a point-by-point basis.Where applicable, azimuthal correlations are included via summation over antipodal points.In all relevant single-unresolved limits, we have found very good agreement between the subtraction terms and real-radiation matrix elements.
To confirm the cancellation of explicit poles in the virtual corrections dΓ V by the virtual subtraction terms dΓ T NLO , we have both checked the cancellation analytically and numerically.

Results
In the following, we shall present numerical results for the five different four-jet event shapes: thrust minor, light-hemisphere mass, narrow jet broadening, D-parameter, and the four-to-three-jet transition variable in the Durham algorithm.We discuss our numerical setup and scale-variation prescription in section 3.1, before defining the four-jet event-shape observables in section 3.2.Theoretical predictions are shown and discussed in section 3.4.

Numerical setup and scale-variation prescription
We consider on-shell Higgs decays with a Higgs mass of M H = 125.09GeV.We work in the G µ -scheme and consider electroweak quantities as constant parameters.Specifically, we consider the following electroweak input parameters yielding a corresponding value of α = 1 128 .As alluded to above, we keep a vanishing kinematical mass of the b-quark throughout the calculation, but consider a non-vanishing Yukawa mass.The running of the Yukawa mass with µ R is calculated using the results of [62], corresponding to mb (M H ) close to 2.61 GeV.
We choose the Higgs mass as renormalisation scale, µ R = M H , and apply a scale variation µ R → k µ µ R about this central scale with k µ ∈ 1 2 , 2 .We wish to point out that this prescription also affects the normalisation of our distributions via eq.(2.5).We use one-and two-loop running for the strong coupling α s , at LO and NLO respectively, obtained by solving the renormalisation-group equation at the given order, as detailed in [12,14].For the strong coupling, we choose a nominal value at scale M Z given by corresponding to the current world average [93].

Four-jet event-shape observables
We consider the five different four-jet event-shape observables related to the following event shapes thrust minor, light-hemisphere mass, narrow jet broadening, D-parameter, and the Durham four-to-three-jet transition variable.They are defined as follows [40]: Thrust minor We define thrust minor as where ⃗ n Min is given by ⃗ n Min = ⃗ n T × ⃗ n Maj .Here, ⃗ n T is the thrust axis [94,95], defined as the unit vector which maximises and ⃗ n Maj is a unit vector for which in addition ⃗ n Maj • ⃗ n T = 0 holds.The thrust minor measures the distribution of particles orthogonal to the plane defined by ⃗ n T and ⃗ n Maj .It vanishes for topologies with less than four particles, as three particles always span a plane parallel to the thrust and thrust-major axis and two particles exactly align with the thrust axis.
Light-hemisphere mass Starting from the thrust axis ⃗ n T , events are divided into two hemispheres (jets), H 1 and H 2 , separated by an axis orthogonal to ⃗ n T .For each hemisphere, the hemisphere mass is calculated as with E vis the visible energy in the event.The light hemisphere mass is then given by the smaller of the two hemisphere masses [96], For massless particles, the light-hemisphere mass is non-vanishing only for topologies with at least four particles, as the lighter hemisphere will otherwise always consist of a single particle.
Narrow jet broadening Jet broadenings measure the distribution of the transverse momenta of particles with respect to the thrust axis.The narrow jet broadening is defined as the jet-broadening value of the more narrow hemisphere, with the hemisphere broadenings B i given by [97] for the two hemispheres H 1 and H 2 defined by the thrust axis.Since a single particle has a vanishing hemisphere broadening, the narrow jet broadening is zero for events with less than four particles, for which always one of the hemispheres contains only a single particle.

D-parameter
The D-parameter measures the planarity of an event and vanishes for planar configurations.As such, it is non-zero only for events containing at least four particles, as three (or less) particles always span a plane.It is defined as [35] D = 27λ 1 λ 2 λ 3 , (3.9) given in terms of the three eigenvalues λ 1 , λ 2 , and λ 3 of the linearised momentum tensor, Four-to-three-jet transition variable The four-to-three-jet transition variable denoted by y 34 corresponds to the jet resolution parameter y cut at which an event changes from a four-jet to a three-jet event according a specific clustering algorithm.Here, we specifically consider the Durham algorithm, in which the distance measure reads We use the so-called E-scheme, in which the four-momenta of the two particles are added linearly in each step of the algorithm, All observables introduced in section 3.2 are non-vanishing only for at least four resolved particles.In the region of phase space where only three jets can be resolved, i.e., where the observables of section 3.2 become small, the differential rate becomes enhanced by large logarithms of the form log m (1/O).In this phase-space region, fixed-order calculations become unreliable and a faithful calculation of the differential decay width requires the resummation of the large logarithms.
To avoid large logarithmic contributions spoiling the convergence of our fixed-order calculation in the three-jet region, we restrict the range of validity of our predictions as follows: we impose a small cut-off y 0 = 10 −5 on linear distributions and y 0 = e −7 on logarithmically binned distributions.We consider this minimal value for all observables mentioned above in section 3.2 with the exception of the distribution related to the fourto-three-jet transition variable y 34 , where we require a cut off at y 0 = e −10 .The y 0 cut-offs are imposed on both four-and five-parton configurations and ensure the reliability of our predictions in the whole kinematical region considered.
In addition to the observable cutoff y 0 , we introduce a technical cut-off t cut and require that all dimensionless two-parton invariants in the five-parton states to be above this cut, s ij /s > t cut .This technical cut-off parameter improves the numerical stability of the antenna-subtraction procedure by avoiding large numerical cancellations between the real subtraction term and the real matrix element.By default, we employ a technical cut-off of t cut = 10 −7 .
We note that there is a subtle interplay of the observable cut-off y 0 and the technical cut-off t cut .We have studied this interplay for all observables considered in section 3.2 by varying the technical cut-off between t cut ∈ {10 −8 , 10 −7 , 10 −6 } and verified that this variation leaves all distributions above the observable cut-off y 0 unaffected.We further wish to emphasise that the independence on t cut also validates the correct implementation of our subtraction terms.

Comparison of predictions in both Higgs-decay categories
In this subsection, we present theoretical predictions at LO and NLO QCD for the event shapes defined in section 3.2.For each event shape, we show two binnings in the observable; a linear binning to highlight the general structure and a logarithmic binning to emphasise the behaviour of the distribution in the infrared region, in particular the position of its peak.In all cases, we present results according to eq. (2.5).In particular, we normalise by the LO (NLO) partial two-jet decay width for LO (NLO) distributions. 3umerical predictions for the five event-shape observables are shown in figs.All distributions exhibit the usual characteristics of event-shape observables at fixed order.
The LO distributions diverge towards positive infinity in the infrared limit on the lefthand side of the plots.In these regions of phase-space, one of the four particles in the LO Born configuration becomes unresolved and the four-particle configuration assumes the shape of a planar three-jet event.At NLO, the distributions develop a peak close to zero and diverge towards negative infinity in the infrared limit.As a consequence, all distributions show these characteristic behaviour towards the infrared region.As mentioned before, an accurate description of the observables in this phase-space region requires the resummation of large logarithms, which we do not include in our predictions.Instead as detailed in section 3.3, we restrict ourselves to provide predictions above a minimum value of the observables.Generically, we observe rather large NLO corrections with K-factors between 1.3 and 2.3 and scale uncertainties of similar size in both categories.It may be surprising that these corrections are numerically slightly bigger for the H → b b decay category.The reason for this is an interplay of two effects.On the one hand, predictions in the Hgg category generally receive larger NLO corrections, which results in numerically bigger NLO coefficients dC/dO.On the other hand, the NLO distribution are normalised by 1/Γ (1) , as opposed to 1/Γ (0) in the LO case, cf.eq.(2.5).Because the NLO correction to the partial two-particle decay width is again numerically bigger for H → gg decays than for H → b b decays (∼ 1.6 compared to ∼ 1.2 at µ R = M H ), the NLO distributions are subject to a more sizeable scaling.For the observables considered here, this has the effect that the full NLO correction becomes numerically smaller in the H → gg decay category.
Thrust minor In fig. 4 we present results for the thrust minor at LO and NLO in the two decay categories.In both decay modes, we observe rather large NLO corrections, with Kfactors at the central scale around 1.9 in the Hb b category and 1.7 in the Hgg category.In addition, there is a visible shape difference between the two decay modes, as can be inferred  Light Hemisphere Mass    Four-to-three-jet transition variable Figure 8 contains results for the four-to-threejet transition variable denoted by y 34 , using the Durham jet algorithm to build jets at LO and NLO for the two decay categories.We note that, different to the other event shapes, we plot results for this four-jet resolution distribution only with a logarithmic binning.Nevertheless, we observe similar results as before.Events in the Hgg category a peak close to the limit before diverging towards negative infinity.
We have shown that all of the event shapes receive large NLO corrections with Kfactors at the central scale between 1.7 − 2.3 in the Hb b decay category and 1.3 − 1.8 in the Hgg decay category, depending on the observable.Interestingly, the NLO corrections are slightly smaller in H → gg decays.This can be explained by the normalisation to the NLO two-particle decay width, which receives larger corrections in the H → gg case and as such leads to a stronger scaling in this decay category.If, instead, a normalisation to the LO two-parton decay width was chosen, we would find larger NLO K-factors for distributions in the H → gg category.In both decay categories, the largest NLO corrections can be observed for the D-parameter and thrust minor.The size of the NLO corrections for the narrow jet broadening, light-hemisphere mass, and the four-jet resolution is slightly smaller.
Concerning the shape of the distributions inside a given Higgs-decay category, we have shown that all event-shape distributions display very similar behaviour at LO, whereas the shape visibly changes at NLO level.For all observables considered here, the shape of the NLO corrections is rather flat for Hb b decays, while it is more curved in Hgg decays.
We further find visible shape differences between the Yukawa-induced and gluonic decay categories at NLO.Specifically, we observe a characteristic shift of the peaks of the NLO distributions away from the infrared limit for H → gg decays, accompanied by a similar shift of the intersection of the LO and NLO predictions.
Our calculation provides crucial ingredients for the computation of higher-order QCD corrections for event-shape observables in hadronic Higgs decays.Most imminently, it gives access to precise predictions including next-to-leading order corrections to four-jet like event-shape observables in all hadronic Higgs-decay channels, as needed for phenomenological studies of Higgs decay properties, at a future lepton collider.While we have focussed only on decays to bottom quarks here, it is straightforward to extend our implementation to decays to other quark flavour species.We have implemented our computation in such a way that it facilitates the matching to resummed predictions in the future.Our calculation also provides part of the real-virtual and double-real contributions to the NNLO calculation of three-jet event-shape observables.To this end, the NLO subtraction terms used here need to be suitably extended to NNLO-type subtraction terms.As we have here constructed all subtraction terms within the antenna-subtraction framework, this is a straightforward extension and will be subject of future work.

Figure 1 :
Figure 1: Feynman diagrams of the two decay categories of a Higgs boson decaying into two jets, H → b b (left) and H → gg (right).

Figure 2 :
Figure 2: Tree-level Feynman diagrams of four-parton processes in the H → b b decay category contributing to our calculation at LO. Here, q denote any flavour u, d, s, c, b.
using the Standard-Model Lagrangian.The two-parton decay diagram associated to this decay mode is shown in the left-hand panel of fig. 1.In the remainder of this paper, we shall refer to these decay modes as belonging to the Hb b category.The relevant four-and five-parton processes contributing to this decay category, entering our calculation at LO and NLO, respectively, are shown in the first column of table 1. Representative Feynman diagrams for the tree-level four-parton processes contributing at LO are shown in fig. 2.

4 to 8 .
Each figure contains four plots; we present results with linear binning in the top row and results with logarithmic binning in the bottom row; results in the H → b b decay category are shown in the left-hand column, while results in the H → gg decay category are shown in the right-hand column.Each plot contains LO and NLO predictions, shown with a dashed and solid line, respectively.Predictions obtained by the scale variation are shown with lighter shading.The infrared region is located towards the left-hand side of each plot, whereas the hard multi-particle region is located on the right-hand side.

Figure 4 :
Figure 4: Thrust minor in the H → b b (left) and H → gg (right) decay category.The top row shows linearly binned histograms, the bottom row shows the same histograms with logarithmic binning, see main text.

Figure 5 :
Figure 5: Light hemisphere mass in the H → b b (left) and H → gg (right) decay category.The top row shows linearly binned histograms, the bottom row shows the same histograms with logarithmic binning, see main text.

Figure 7 :
Figure 7: D-parameter in the H → b b (left) and H → gg (right) decay category.The top row shows linearly binned histograms, the bottom row shows the same histograms with logarithmic binning, see main text.

Table 1 :
Partonic channels contributing to the decay H → 4j at LO and NLO.