Measuring Hadronic Higgs Boson Branching Ratios at Future Lepton Colliders

We present a novel strategy for the simultaneous measurement of Higgs-boson branching ratios into gluons and light quarks at a future lepton collider operating in the Higgs-factory mode. Our method is based on template fits to global event-shape observables, and in particular fractional energy correlations, thereby exploiting differences in the QCD radiation patterns of quarks and gluons. In a constrained fit of the deviations of the light-flavour hadronic Higgs-boson branching ratios from their Standard Model expectations, based on an integrated luminosity of $5\,\text{ab}^{-1}$, we obtain $68\%$ confidence level limits of $\mu_{gg}=1 \pm 0.05$ and $\mu_{q\bar{q}}<21$.


Introduction
Electron-positron colliders operating as Higgs-boson factories constitute one of the main options for future accelerator-based high-energy experiments [1,2].Similar to the LEP experiments of the 1990's and their precision determination of the Z-boson properties, the primary target of such a facility will be the profound analysis of the Higgs-boson properties.This encompasses in particular measuring its width, branching ratios, and kinematic distributions of its decay products with unrivalled precision, searching for minute deviations from Standard Model predictions, usually formulated in the language of Effective Field Theories, see, e.g.[3][4][5].One of the observables that attracted considerable attention is the branching ratio of the Higgs boson, and hence its coupling, to gluons, which at a hadron collider such as the LHC can be accessed only through its (dominant) gluon-fusion production mode.Due to the much cleaner final states produced in lepton-lepton annihilation, at a future e + e − collider this decay can be studied directly in events where the Higgs boson decays into jets.Ignoring for the moment decays into intermediate W and Z-boson pairs, the gluon decay mode can be extracted by vetoing QCD events with displaced vertices which emerge in the weak decays of b and c hadrons [6], a technique also underpinning measurements of branching ratios to heavy quarks at the LHC [7][8][9][10].Having thus eliminated the two main competitor decay modes with similar characteristics, the remaining QCD events are due to the Higgs boson coupling to the light degrees of freedom.Realising that the Yukawa coupling to the u, d, and s quarks is negligible in the Standard Model yields the coupling of the Higgs boson to gluons.Results based on this strategy indicate branching ratio measurements at relative per-cent accuracy [11,12].In this study we will advocate an alternative approach to this measurement, that is agnostic towards the presence of displaced vertices and is underpinned by fits to event-shape observables alone.This is in line with previous ideas to use event-or jet shapes as theoretically well understood taggers [13,14], but goes arXiv:2306.03682v2[hep-ph] 29 Nov 2023 a step further by never explicitly assigning a flavour to a given event.Such a strategy relies on two wellestablished properties of QCD radiation patterns, namely firstly that gluons carry two colours resulting in a ratio C A /C F = 9/4 of colour charges with respect to quarks and hence about twice as many emissions, and secondly that the finite masses of the heavy quarks shield the collinear divergence of gluon emission, thereby depleting their QCD radiation in this region, a phenomenon also known as the "dead cone" effect [15].Combining both effects allows the placement of direct constraints to the sum of the light-quark Yukawa couplings.
In particular we will use fractional energy correlations [16] which are geared towards a systematic study of the collinear regions of the radiation pattern.Our studies here supplement the existing strategy for the measurement of the Higgs-boson branching ratio to gluons, see also Ref. [17] for a recent study using jet charge as a discriminating variable.They furthermore provide an alternative to first attempts to measure the Yukawa coupling to light quarks through rare decays such as H → ϕγ at the LHC [18,19].In relying on event-shape variables alone, hadronisation effects constitute the dominant systematic uncertainty when fitting Monte Carlo results to (synthetic) data.To account for this we re-tune the new cluster-hadronisation model [20] of the SHERPA event generator [21] to LEP data, and quantify the resulting uncertainties through repeated tunes with varying input data.We refer to the resulting sets of alternative hadronisation parameter values as replica tunes.Our discussion is structured as follows: In Section 2 we detail the setup of our analysis, and in particular the event-shape observables we use.This is followed by Section 3 where we describe our simulations with SHERPA and we put special emphasis on the re-tuning of its fragmentation model.In Section 4 we discuss the results emerging from fits to various event-shape distributions, with and without soft-drop grooming the hadronic final state, which we present as allowed two-dimensional regions of values for the deviations µ gg and µ qq of the Higgs boson branching ratios into gluons and light quarks.We conclude and summarise our study in Section 5.

Analysis of event shapes in e + e − → ZH
We here propose an analysis at a future lepton collider operating at the working point for Higgs-strahlung production, i.e. √ s ≳ m H +m Z .Our goal is to consider hadronic Higgs-boson decays, where we separate the branchings to gluons, light (up, down and strange) quarks, charm quarks and bottom quarks.Experimentally, those will at first all be seen as hadronic channels.However, due to the difference in the QCD radiation pattern between quarks and gluons, and the imprints of finite quark masses, one can expect differences in observables such as the well-studied event shapes.Following the selection cuts of Ref. [6], we identify Z-boson candidates as pairs of opposite-sign leptons within ±5 GeV of the nominal Z mass.The reconstructed Z-boson is required to have at least a transverse momentum of p T,Z > 10 GeV and a longitudinal momentum of at most 50 GeV.To suppress irreducible backgrounds from ZZ events, we require for the opening angle between the two leptons θ l + l − < 100 • .We additionally ask for a total hadronic mass of all other particles to be at least m had > 75 GeV.In order to select events where the hadronic final state is likely to originate from a decaying Higgs-boson, we constrain the recoil mass of the lepton pair, defined as see also [2], to be similar to the Higgs-boson mass.In practice we use 120 GeV < m recoil < 130 GeV.We base the calculation of event-shape observables on charged-particle tracks1 and consider the family of fractional energy correlations [16] FC with x = 0.5, 1, 1.5.The sums run over all charged tracks i and j with respective energies E i,j and threemomenta ⃗ q i,j .All energies and angles are evaluated in the Higgs-boson rest frame, which we reconstruct as the full charged final state, excluding the two leptons from the Z-boson decay.We can analyse the behaviour of this class of observables in response to a single soft-gluon emission off a hard parton from the hadronic decay as commonly done in the context of resummation calculations [16].In terms of the soft-gluon transverse momentum k t and rapidity η relative to the hard parton, the fractional energy observables scale like Hence x = 1 corresponds to a purely transverse-momentum like scaling; larger (smaller) values of x will give a higher weight for pairs of particles with smaller (larger) opening angles.The fractional energy correlations are very similar to jet angularities that are studied at the LHC in the context of quark and gluon tagging, see for example [22][23][24][25][26].The choice x = 1.5 corresponds to the Les-Houches angularity [27].By the Heaviside function in Eq. (2.2), the fractional energy correlations are implemented as sum over the contributions from two hemispheres defined by the axis ⃗ n T .We follow the original definition and use the thrust variable to define the reference axis, and hence the hemispheres.There are further standard observables that are written in this way, for example total hemisphere broadening, or total mass We also considered properties of individual hemispheres, for example the mass of the heavier hemisphere and the broadening of the wider one, but did not observe any noteworthy increase in performance and hence focus on the ones described above.In terms of Eq. (2.3), the broadening scales with b = 0 and behaves similar to FC 1 while the mass would correspond to FC 0 , i.e. b = 1, which we do not analyse here.
As an additional handle, we employ soft-drop grooming [28] with the goal to reduce hadronisation corrections.While the soft-drop grooming algorithm was originally developed to mitigate the contamination of jets from effects that are typically simulated as underlying event or multiple parton interactions, it has been shown to be effective in mitigating non-perturbative corrections to event-shape observables in leptonic and hadronic collisions as well [29][30][31].We apply the algorithm individually to the two event hemispheres.To this end, we recluster their respective constituents using the Cambridge/Aachen jet algorithm [32,33], then undoing the last clustering step between the subjets i, j and checking for the soft-drop condition If this condition is satisfied, the procedure terminates.Otherwise, the softer of the subjets, with smaller energy, is discarded and the procedure is repeated for the harder one.This continues until Eq.(2.6) is true, or the remaining subjet consists of only one track.Note that other references include an angular dependence in the soft-drop condition, whereas we here only consider the β = 0 case of [28], which is equivalent to the modified mass drop tagger [34,35].We also restrict our study to the conventional case z cut = 0.1.
The observables are then calculated on the remaining particles after grooming, however, normalised by the hadronic energy before grooming.This treatment is necessary in order to define collinear safe observables [30].
For our final results, histograms for the differential distribution of event shapes v are constructed as sums over the individual decay channels where the sum runs over the hadronic decay modes of the Higgs-boson, into light (q q), charm (cc) and bottom (b b) quarks, gluons (gg), as well as two hadronically decaying vector bosons.In the last term we add the irreducible background from ZZ production.The factors µ i parametrise deviations from the Standard Model (SM) partial Higgs-boson decay widths, with the SM corresponding to µ i = 1 ∀ i.Ultimately, we aim for an experimental determination of the coefficients µ i .In the following, we explore the possibility to set limits on simultaneous deviations of µ gg , µ q q and µ b b from 1, while leaving the total cross section unchanged.
To achieve this, we scan different points in µ gg and µ q q , fixing µ cc = µ W W = µ ZZ = 1, and imposing the constraint (2.8) 3 Simulation with SHERPA To simulate particle-level events we use the SHERPA event generator.The main physics aspects of the framework are documented in [21], while we here work with a pre-release version 3.0β [36].We use SHERPA's default dipole shower based on Catani-Seymour factorisation [37], that supports finite parton masses in the splitting kernels and branching kinematics.Parton-to-hadron transitions are described by SHERPA's built-in cluster-hadronisation model [20] and hadron decays are treated by its internal decay package [21,38].We will comment below on a dedicated hadronisation-parameter tune based on sensitive measurements from LEP experiments.We analyse our simulated data with the RIVET package [39] and use the CONTUR tool [40,41] for statistical analyses and the calculation of exclusion limits.

FCC-ee setup
We assume the operating conditions for a Future Circular Collider, FCC-ee, running at a centre-of-mass energy of √ s = 240 GeV [1,2].We simulate the processes e + e − → Z(→ µ + µ − )H(→ q q) at order α 3 EW y 2 q separately for q = u, d, s, collectively referred to as light-quark decays, q = c and q = b, and e + e − → Z(→ µ + µ − )H(→ gg).We here assume the one-loop decay H → gg in the heavy top-quark limit, treating it through an effective ggH vertex [42][43][44].We take into account the Higgs-boson decays to W W * and ZZ * by generating a sample of e + e − → Z(→ µ + µ − )H events where the Higgs is forced to decay into W q q′ or Zq q and the on-shell vector bosons likewise decay into quarks.We rescale our leading-order results to the branching ratios from [45], corresponding to by appropriately adjusting the event weights in the final samples.While we neglect contributions from q = u, d, we estimate the H → ss branching ratio by scaling the H → cc result [46], i.e. q=u,d,s, Note that we handle c and b quarks as massive in the parton-shower evolution [37].Finally, we simulate a sample for resonant and non-resonant di-boson production, i.e. e + e − → µ + µ − q q at order α 4 EW , which we refer to as ZZ background.We do not include any higher-order corrections at this stage of the simulation.The computed cross sections are scaled to an integrated luminosity of L = 5 ab −1 to account for the full projected statistics accumulated by the FCC-ee at this centre-of-mass energy, and we include another factor of two to emulate the use of both electron and muon channels.We construct histograms for the considered set of observables according to Eq. (2.7).The statistical errors are scaled accordingly with the estimated number of entries for a given bin i, N i , as √ N i .This is combined with a covariance matrix for the systematic variations of hadronisation-model parameters, see below, derived as where the average is taken over runs with different tuning parameters.
In Fig. 1 we show example predictions for the fractional energy correlations F C 1.5 and F C 0.5 .We illustrate the SM distribution obtained from our simulations with SHERPA as well as the variations corresponding to two representative sampling points in the (µ gg , µ q q ) plane, i.e. {µ gg = 1, µ q q = 4} and {µ gg = 1.18, µ q q = 1}.

LEP1 setup and Tuning
For the present study we perform dedicated tunes of SHERPA's new cluster-fragmentation model AHADIC++ [20], focusing on event-shape observables.We provide non-perturbative (tuning) uncertainties through replica tunes.
Similar to previous tunes, we concentrate on observables for hadronic final states in electron-positron annihilation accurately measured at LEP1.Our simulations rely on next-to-leading order (NLO) QCD matrix elements for e + e − → q q + {0, 1}j dressed by the parton shower, using the SHERPA implementation of the MEPS@NLO formalism [47].The contributing tree-level amplitudes are obtained from the built-in matrix element generators COMIX [48] and AMEGIC++ [49], while the required one-loop amplitudes are obtained from OPENLOOPS [50].For the hard-scattering amplitudes we consider q = u, d, s, c massless and only take into account mass-effects for q = b, therefore adding explicitly the tree-level contribution for the 4b final state.As jet-separation parameter in the merging prescription we use log 10 Q 2 cut /s = −2 [51].We employ the APPRENTICE tuning tool [52] in combination with analyses provided by RIVET [39].We simultaneously tune all parameters listed in Table 2 in Appendix B, starting from rather wide parameter ranges, decreasing the intervals in a sequence of tunes.APPRENTICE uses actual generator runs with varied hadronisation parameters to construct a bin-wise polynomial surrogate of the Monte-Carlo response for observables of interest.To narrow down the tuning ranges in subsequent iterations, we construct multiple such surrogates each with a different set of generator runs and in this way find equivalent tunes, with results of similar quality.The outcome of such set of tunes is then used to shrink the parameter ranges for the next iteration.After the final iteration, the obtained family of equivalent tunes is used to estimate the remaining non-perturbative uncertainties by interpreting its members as replica tunes, i.e. by re-running the Monte-Carlo simulation with the alternative parameter values.Our observable selection for the new generator tunes is similar to the ones used for the initial AHADIC++ tunes presented in [20], as well as its predecessor [53].It consists mostly of mean and differential chargedparticle multiplicities [54], event shapes like thrust and its minor and major variants [55,56], as well as the b-quark fragmentation function [57,58].Furthermore, we consider jet-rates for the Durham algorithm [59], and a selection of multiplicities of identified hadrons [60], thereby aiming for a general tune suitable for event-shape and jet observables.The complete list of analyses and differential distributions can be found in Table 1 of Appendix B. Figure 2 shows exemplary results for observables used in the tuning, including the final SHERPA prediction and the corresponding non-perturbative tune uncertainty indicated by the light-blue band.To allow for sufficient freedom in the variations for all 16 parameters considered in the tuning, we provide 50 replica   tunes.They represent our tuning uncertainties by having each replica tuned with a different, random subset of Monte-Carlo runs.We extract the uncertainty bands in Fig. 2 by re-running the simulation for each of the replica tunes, and plot the envelope of the resulting deviations.We find good agreement between our SHERPA predictions and data, with deviations of the central tune to the data being largely covered by our estimated non-perturbative uncertainties.

Results
Let us now turn to the statistical analysis of event shapes, measured as described in Sec. 2. We add the histograms corresponding to our analysis to the CONTUR framework [40], and use its statistical analysis modules to compute confidence levels for the exclusion of different points in the (µ gg , µ q q ) plane.The two dimensional exclusion limits for µ gg versus µ q q based on the three fractional energy correlations F C 1.5 , F C 1 and F C 0.5 are shown in Fig. 3, respectively.As already indicated in the introduction, QCD radiation tends to mainly populate the soft and collinear regions of phase space, where the dead-cone effect associated to the finite and relatively large masses of the c and b quarks most visibly manifest themselves, and where differences due to different colour charges (the C F of quarks versus the C A of the gluons) lead to directly observable differences in the numbers of particles emitted.Accordingly, the Les-Houches angularity F C 1.5 tends to be the most sensitive observable, since it gives the largest weight to collinear emissions.Nevertheless, all three choices x = 1.5, 1, 0.5 are able to limit µ gg to be within 1 ± 0.10 based on a 68% confidence limit.However, based on F C 1.5 one should be able to set a stronger limit on µ gg to be within 1 ± 0.07 and additionally limit µ q q < 33, while we can only exclude µ q q values larger than 45 based on F C 1 .Finally, F C 0.5 appears to not be sensitive to µ q q within the range we consider.Including soft-drop grooming of the hemispheres does not result in any significant improvements, as shown in Fig. 4, the equivalent of Fig. 3 but with grooming included.In fact, the limits worsen slightly, which could to some extend have been anticipated, since there are competing effects at work.Grooming will remove some information from the radiation pattern, but on the other hand has the potential to reduce the impact of hadronisation corrections and hence the associated systematic uncertainty.Apparently, this reduction is not sufficient to compensate for the loss in information, at least with the grooming parameters we have considered here.One could imagine that an optimisation of z cut and the inclusion of angular dependence in the soft-drop condition could lead to more competitive results.In addition it is certainly worth stressing that the combination of probable future refinements of the hadronisation models and the drastically increased data set of a potential FCC-ee (10 12 events vs 10 7 at LEP-I) will most likely significantly reduce the uncertainty related to the modelling of the parton-to-hadron transition.To illustrate this, we present in Appendix A Fig. 7 selected exclusion-limit plots for the scenario of negligible non-perturbative uncertainties.As anticipated, the limits improve, resulting in µ gg = 1 ± 0.05 and µ q q < 25 for plain F C 1.5 , and µ gg = 1 ± 0.06 and µ q q < 28 for its soft-drop groomed variant.One of the major differences of our procedures so far, compared to traditional tagging methods, is that we effectively tag any event as a whole.When individually tagging jets, or hemispheres for that matter, one would want to include a requirement that both tags are compatible with the desired final state.To mimic this, we consider a measurement of the fractional energy correlations but on each hemisphere separately.
We then derive exclusion limits based on the corresponding two-dimensional histograms.While we hope to expose additional information in this way, it should be clear that this is a more involved observable definition.In particular, joint resummed calculations of several observables are far less advanced than what would be available for the distributions considered above.The resulting confidence levels can be found in Fig. 5.They are somewhat improved compared to the baseline in Fig. 3.In particular, for F C 1.5 we obtain µ gg = 1 ± 0.05 and µ q q < 21.This suggests that combining individual results from the two hemispheres into  one measurement reduces the impact of systematic uncertainties due to non-perturbative corrections, i.e. hadronization.This is further highlighted in the appendix, Fig. 7 where we elucidate the impact of vanishing hadronization uncertainties due to the improved understanding and tuning of the corresponding models.Lastly, we analyse some examples of more traditional event shapes that were measured by the LEP experiments and enter our tuning.We focus on the sum of the hemisphere masses and the total broadening.The hemisphere masses would scale similar to F C 0 , which we do not show, in the infrared limit, while total broadening is similar to F C 1 .Correspondingly, the limits derived from the masses, shown in the leftmost subplot of Fig. 6, is even weaker than that obtained from F C 0.5 , consistent with the decline depending on x.
In the middle of Fig. 6 we show the same plot but for total broadening, where we find a similar behaviour to F C 1 , as expected.Finally, in the rightmost plot we show the limits obtained from considering the twodimensional distribution of the broadening of both hemispheres.Again, we observe slightly improved limits compared to the one-dimensional distribution, in line with what was observed for F C 1 .

Conclusions
In this article we have presented a proposal for the extraction of hadronic branching ratios of the Higgs boson at a future lepton collider operated as a Higgs factory.To separate couplings of the Higgs boson to light degrees of freedom, i.e. those to gluons and light quarks from the ones to massive charm and bottom quarks, our ansatz does not require the reconstruction of displaced vertices from weak decays of heavy flavour hadrons.Instead, we propose to instrument event-shape observables that are sensitive to the differences in the radiation patterns of light/heavy quarks and gluons.In particular, we have focused on fractional energy correlations with and without soft-drop grooming of the two event hemispheres with respect to the thrust axis.Based on dedicated simulations with the SHERPA event generator framework we showed that stringent limits on the deviation of the Higgs-boson branching ratios into light-quarks, gluons and bottom-quarks can be obtained.Assuming the full projected luminosity of 5 ab −1 for the FCC-ee collider at √ s = 240 GeV, we were able to derive limits as tight as µ gg = 1 ± 0.05 , µ q q < 21 , while Eq.(2.8) has to hold.
The best sensitivity we obtained for the fractional energy correlation F C 1.5 that is geared to explore the collinear region of the hard shower initiator.To address the systematic uncertainty from the non-perturbative fragmentation of partons to hadrons we performed a dedicated tune to event-shape and jet observables measured by the LEP experiments, thereby deriving, for the first time, a whole family of replica tunes that allow us to estimate the residual model-parameter uncertainties of SHERPA's cluster fragmentation.The obtained results motivate further studies and refinements.In our study, we deliberately ignored any information on potential secondary vertices that could in principle be used for a better separation of the heavy-flavour contributions.Furthermore, additional selection cuts could help to further reduce the contribution from Higgs-boson decays to vector bosons that subsequently decay hadronically, populating the region of rather large event-shape values, thereby diluting the gluon signal.Similarly, the parameters of the grooming algorithm can certainly be further optimised to better balance the reduction of non-perturbative corrections against the imprints of the QCD radiation pattern in the hadronic final states.Concerning the theoretical predictions, for the considered event-shape variables all-orders analytical predictions could be derived, for example at next-leading-logarithmic accuracy based on the CAESAR formalism [16,30,61,62], or at next-to-next-to-leading-logarithmic level through ARES [63,64].Alternatively, accurate resummed predictions for event shapes have been obtained using effective field theory techniques, see for example [65,66], which have also been applied to jet observables closely related to the energy correlations used here [67][68][69][70][71]. Furthermore, progress has recently been made on including finite quark masses in resummed calculations [72][73][74].To further reduce non-perturbative uncertainties we envisage dedicated analyses of hadronic final states at the FCC-ee prior to attempts to measure Higgs-boson couplings that should enter the tuning and help to further constrain the model-parameter uncertainties.

A Exclusion limits assuming vanishing non-perturbative uncertainties
We highlight the effect of vanishing non-perturbative uncertainties in Fig. 7, where we show exclusion plots for the three procedural variants of the F C 1.5 observable considered before.

B Central tune parameters and uncertainty ranges
The re-tuning of SHERPA's hadronisation model was based on RIVET analyses of LEP 1 measurements and observables detailed in Tab. 1.We list the the initially considered parameter range, the central tune value and corresponding uncertainty intervals, determined from 50 replica tunes, in Tab. 2. For the full list of AHADIC++ model parameters and their physical interpretation we refer the reader to Ref. [20].

Figure 1 :
Figure 1: Predictions for fractional energy correlations F C 1.5 (left) and F C 0.5 (right) at the FCC-ee for the SM and two hypotheses for modified couplings of the Higgs-boson to QCD partons.The bands indicate the combined statistical and systematic errors of the SM prediction.

Figure 3 :
Figure 3: Exclusion limits based on fractional energy correlations (from left to right) F C 1.5 , F C 1 , F C 0.5 .

Figure 5 :
Figure 5: Exclusion limits based on fractional energy correlations (from left to right) F C 1.5 , F C 1 , F C 0.5 , measured individually on the two hemispheres.

Figure 6 :
Figure 6: Exclusion limits based on (from left to right) the sum of the hemisphere masses, the total broadening, and the 2-dimensional distribution of broadenings in the two hemispheres.

Figure 7 :
Figure 7: Examples of limits obtained from variants of the F C 1.5 observable when neglecting nonperturbative uncertainties, without and with soft-drop grooming (left, middle) and measured individually on the two event hemispheres (right).The transparent bounds are identical to the main text, while the solid ones are obtained when including statistical uncertainties only.
[58] Shown are the charged-particle multiplicity n ch measured by ALEPH[54](left), thrust T as measured by DELPHI[56](center), and the B-hadron energy fraction x B measured by OPAL[58](right).Each of the SHERPA predictions corresponds to 10 7 events, ensuring that the statistical errors are negligible and the depicted uncertainties are dominated by the variations of non-perturbative model parameters.
Figure 4: Exclusion limits based on soft-drop groomed fractional energy correlations (from left to right) F C 1.5 , F C 1 , F C 0.5 .