Probing light quark Yukawa couplings through angularity distributions in Higgs boson decay

We propose to utilize angularity distributions in Higgs boson decay to probe light quark Yukawa couplings at $e^+e^-$ colliders. Angularities $\tau_a$ are a class of 2-jet event shapes with variable and tunable sensitivity to the distribution of radiation in hadronic jets in the final state. Using soft-collinear effective theory (SCET), we present a prediction of angularity distributions from Higgs decaying to quark and gluon states at $e^+e^-$ colliders to ${\rm NNLL}+\mathcal{O}(\alpha_s)$ accuracy. Due to the different color structures in quark and gluon jets, the angularity distributions from $H\to q\bar{q}$ and $H\to gg$ show different behaviors and can be used to constrain the light quark Yukawa couplings. We show that the upper limit of light quark Yukawa couplings could be probed to the level of $\sim15\%$ of the bottom quark Yukawa coupling in the Standard Model in a conservative analysis window far away from nonperturbative effects and other uncertainties; the limit can be pushed to $\lesssim 7-9\%$ with better control of the nonperturbative effects especially on gluon angularity distributions and/or with multiple angularities.


Introduction
After the discovery of the Higgs boson at the Large Hadron Collider (LHC) [1,2], precision measurements of the properties of the Higgs and proving that the Higgs boson is indeed responsible for electroweak symmetry breaking and mass generation have become a forefront goal of high energy physics.Determining the Yukawa couplings of the Higgs boson to fermions is one of the avenues to verify the Standard Model (SM) of particle physics.Since the Yukawa couplings are completely determined by the fermion mass in the SM, i.e. y q = √ 2m q /v with v = 246 GeV, it is virtually impossible to probe the light quark Yukawa couplings directly due to the smallness of their mass.But these light quark Yukawa couplings may receive large modifications from new physics (NP) beyond the SM [3], and thus these NP effects could be probed by careful measurements of the Yukawa couplings.
Due to the large QCD backgrounds for the hadronic decay of Higgs boson at the LHC, a direct measurement of the light quark Yukawa couplings is challenging.Several approaches have been proposed to constrain the light quark Yukawa couplings indirectly.For example, one can measure the light quark Yukawa coupling through i) rare modes of Higgs boson decays, e.g.h → J/Ψγ (ϕγ, ργ, ωγ) [4][5][6]; ii) Higgs production in association with a charmtagged jet [7]; iii) global analysis of the Higgs data [8][9][10][11]; iv) the transverse momentum (p T ) distribution of Higgs boson or jet in Higgs production processes [12][13][14][15][16][17]; v) the kinematic shapes of Higgs pair production [18], off-shell Higgs production [19,20] and triple heavy vector boson production [21]; etc.The above proposals demand accurate calculations of light quark meson formation, charm tagging efficiency and faking rates of light quarks, or precise knowledge of the p T spectrum of Higgs boson.Compared to hadron colliders, e + e − colliders provide direct access to all possible decay channels of the Higgs boson due to the clean environment.Several plans for future lepton colliders have been proposed, including the CEPC [22], ILC [23], CLIC [24] and FCC-ee [25].Through the Higgs and Z boson associated production, the inclusive cross section could be measured to 0.5% accuracy at √ s = 250 GeV with integrated luminosity of 5 ab −1 [22].Such high accuracy of the total cross section offers a possibility to measure the light Yukawa couplings directly at the e + e − colliders.
The main difficulty of measuring the light quark Yukawa couplings at e + e − colliders is due to contamination from Higgs boson decays to gluons, obscuring jets initiated by light quarks.To suppress the gluon background, it is necessary to use some form of quark and gluon jet discrimination.A wide variety of quark and gluon discriminants have been proposed, e.g., in Refs.[26][27][28][29][30][31][32][33][34][35].To first approximation, the main underlying feature is that the initiating energetic quark radiates soft or collinear gluons in proportion to its color factor C F = 4/3, while initiating hard gluons radiate additional gluons proportional to the factor C A = 3, making gluon-initiated jets "broader" or more "diffuse" than quark-initated jets.Good discriminants tease out these and other more subtle differences between gluon and quark jets.Ideally, the discriminants have properties that can also be predicted reliably from first principles in QCD, though powerful methods also exist that are totally data-driven and require no input from QCD at all, e.g.[36].
One class of observables that can discriminate broad features of quark and gluon jets and can also be predicted to high accuracy in QCD are hadronic event shapes (e.g.[37][38][39][40][41]).Thus we explore in this paper their potential to constrain the light quark Yukawa couplings at e + e − colliders.Similar ideas have been discussed in Ref. [42], e.g.event shapes including thrust [43], heavy hemisphere mass [44,45], C parameter [46,47], broadening [48] and Durham 2to-3-jet transition parameter [49], and jet energy profile [50][51][52][53].Ref. [42] showed that these event shapes can provide a much stronger sensitivity for the light quark Yukawa couplings compared to methods proposed for the LHC, e.g.y q < 0.091y b by event shapes for q = u, d, s at 95% confidence level at CEPC [42], while it only can be constrained to 0.4 ∼ 0.5y b by Higgs p T spectrum at the LHC [12,13], where y b is the bottom quark Yukawa coupling in the SM.Recently, Ref. [20] studied the potential of off-shell Higgs production to further improve the Yukawa coupling limit, projecting improvements of 20% or more compared to Higgs+jet production depending on assumptions about uncertainties on background.
Event shapes in e + e − collisions to hadrons have already been computed to very high accuracy, up to N 3 LL ′ resummed accuracy matched to NNLO fixed-order calculations, e.g.[38,39,41,54].Recently, theoretical predictions of event shape observables in Higgs boson decays have also been computed to high accuracy, e.g. the thrust distribution has been calculated to NLO accuracy plus NNLO singular terms [55], the energy-energy correlation from Higgs decaying to gluon mode has been computed to NLO accuracy [56].Ref. [57] also computed the 2-jettiness distribution from the decay of the Higgs boson to a b b quark pair and to gluons, to NNLL ′ +NNLO accuracy and even approximate N 4 LL accuracy in [58].In this work, we propose to use a class of event shapes, angularities [59], to separate the q q channel from the gg channel in Higgs decays and to improve measurements of the light quark Yukawa couplings at lepton colliders.(The thrust and total jet broadening used in Ref. [55] are special cases of the angularities.)We take advantage of recent results that make possible the prediction of angularity distributions to NNLL ′ accuracy in resummed perturbation theory [60].We also newly compute the LO fixed-order O(α s ) corrections to general angularities in Higgs decay.We find that using any single angularity distribution allows determination of the light quark Yukawa couplings to similar or somewhat better precision as other individual event shapes studied in Ref. [55], depending on whether we use a conservative analysis window at larger τ a away from large nonperturbative effects and other corrections, yielding a potential bound of y q ≲ 0.15y b , or a more aggressive window that enters further into the nonperturbative region but thus takes advantage of the very different peaks for the quark and gluon angularity distributions, yielding y q ≲ 0.7 − 0.9y b .We also perform a very preliminary study using Pythia of the ability of double differential distributions of angularities to improve this reach further, attempting to utilize their additional power over single angularity distributions to distinguish quark and gluon jets.We find that, at the present time, realizing large additional improvements is difficult due to the backgrounds from H → b b, cc and Zq q.However, if those backgrounds could be further suppressed, say by a factor 10, both single-and multi-angularity distributions could yield even better limits on y q .
The paper is organized as follows: In Section 2, we obtain the factorization of angularity distributions from Higgs boson decay in the formalism of SCET, allowing large logs in the distributions to be resummed through renormalization group evolution.Both the decay processes H → q q and H → gg are calculated to NNLL+O(α s ) accuracy.Although higher orders are now possible (e.g.[61]), for our illustrative study here, we do not include higher-order corrections.The numerical predictions are given in Section 3. We then show the precision with which light quark Yukawa couplings could be measured using angularity distributions in Section 4. Finally, we conclude in Section 5.The one-loop angularity distributions in Higgs boson decay and some technical details of our analysis are discussed in the Appendix.
2 Factorization and resummation of event shapes

Factorization of cross section in SCET
The angularities are defined as [59,62], where, for us, we will take Q = m H = 125 GeV, the sum is over all final state particles i.The pseudorapidity η i and transverse momentum p i T of each particle i is measured with respect to the thrust axis [43,63] in the rest frame of the Higgs boson.The parameter a defining each angularity τ a is a continuous parameter a < 2 for infrared safety, though we will focus on a < 0.5 in this work, since soft recoil effects which complicate the resummation become important as a → 1 [59,[64][65][66][67][68].The angularities have found a wide range of applications in e + e − collisions producing hadrons (e.g.Refs.[60,[69][70][71][72][73]), and their distributions have been calculated to NNLL ′ resummed and O(α 2 s ) matched accuracy [60].To describe the decay of Higgs boson, we will use the following effective Lagrangian, where µ is the renormalization scale, b is the color index (summed over) of the gluon field with field strength tensor G b µν .In our calculation, we ignore the masses of the light quarks, but keep the Yukawa coupling y q itself.The coupling of the Higgs to the gluon fields comes from integrating out the top quark.
For small values of τ a ≪ 1, the degrees of freedom in the final state at leading power are collinear quarks and gluons in the two back-to-back directions in the Higgs rest frame, and (ultra)soft gluons radiated between them.We can predict the cross section in this regime by matching the operators in Eq. (2.2) onto operators in SCET: Here χ n,n and B n⊥,n⊥ are the SCET gauge invariant collinear quark and gluon fields respectively [74][75][76][77].Computing the cross section in this regime, the event shape distributions in e + e − collisions can be factorized into hard, jet and soft functions [59,60,69,70,78], where i = q, g corresponds to H → q q, gg, t n,n a are the natural arguments of the jet functions of dimension 2 − a [70,78], and k s is the natural dimension-1 argument of the soft function.The decay spectra in Eq. (2.4) are normalized to the leading-order partial decay widths, (2.5) is the hard coefficient of Higgs boson decaying to quarks or gluons, and is given by the square of the matching coefficients in Eq. (2.3), (2.6) These encode virtual fluctuations at scales µ ∼ m H that give quantum corrections to the decay operators on the left-hand sides of Eq. (2.3) that are integrated out of the lower-scale effective theory of collinear and soft excitations.At leading order in QCD factorization, soft wide-angle radiation can be shown to factor for the dynamics inside collinear jets [79], such that the sum of all soft emissions couple only to Wilson lines encoding the light-cone directions n, n of jets and their color representations.In SCET this is encapsulated in a field redefinition of collinear fields with soft Wilson lines so that soft gluons no longer couple directly to any collinear particles [77].All collinear radiation and splitting inside jets are described by the jet functions, J i n,n (t n,n a , µ), defined here by [71,78] where the traces are over color and Dirac (Lorentz) indices, the "label" momentum operators P µ fix the large label momentum components of the SCET collinear quark (gluon) fields χ n (B n⊥ ) [76] (here, Q = m H ), and τ n a is an operator that fixes the angularity of final states produced by the collinear quark (gluon) fields [69].Meanwhile, the soft functions S i (k s , µ) are defined by where Y n (Y n ) are soft Wilson lines in the fundamental (adjoint) representation along the direction of n µ , with A s in the fundamental representation, and Y n is similar but in the adjoint representation.In Eqs.(2.7) and (2.8), the operators τ n,s a acting on a collinear or soft final state returns the contribution to the angularity τ a of that state, defined by its action on the collinear (soft) states X n,s , and can also be constructed in terms of the energy-momentum tensor in QCD [69,[80][81][82].
In the next subsection we will review perturbative calculations of the above functions appearing in the factorized decay spectra and use them to evolve each piece and sum large logarithms in the full decay spectra.

RG Evolution and Resummation
The prediction for the decay spectrum Eq. (2.4) in fixed order QCD perturbation theory contains logs of τ a at every order in α s , which become large for τ a ≪ 1 and need to be resummed.This can be achieved by RG evolution of each piece of the factorized cross section (see, e.g.[78,83]).
The Yukawa coupling y q (µ) and strong coupling α s (µ) obey the following renormalizationgroup (RG) equations, (2.11) The anomalous dimension γ y and the β function have expansions in α s that we express: where the coefficients γ n y , β n are given in Eqs.(B.1) and (B.2).The one-loop hard functions are [84], The soft functions for quark and gluon jets are known to one [85] and two loops [40,86] for a = 0.For generic values of a, the soft function was computed to one loop order in [70]: The one-loop result in Laplace space is, Here the color factor C i = C F,A for quark and gluon respectively and ν a is the Laplaceconjugate variable to the angularity τ a .The two-loop soft functions for a ̸ = 0 have recently become computable using the program SoftSERVE [87][88][89], which were used for predictions of e + e − angularities to NNLL ′ accuracy in [60].To at least two-loop order, whether the directions of the n, n Wilson lines in Eq. (2.8) are incoming or outgoing does not affect the perturbative results for the soft functions [90], so we can also use the results of [88,89] here.We give the results for the non-cusp anomalous dimensions in Eq. (B.5).The 2-loop constant terms for the quark channel can be found in [60,89], though we will not use them in this paper, where we go only to NNLL+O(α s ) accuracy.
The jet functions for massless quark and gluon jets are also known to one [91,92], two [93,94] and even three loops [95,96] for a = 0.For generic values of a, the one-loop jet functions for quarks were computed in [70], and for gluons in [71,97].The result can be expressed in Laplace space, where the non-cusp anomalous dimension coefficient γ i is given by, The coefficient of the constant term is a function of a, f i (a), which is defined as, Here T F = 1 2 , and n f = 5 is the number of active quark flavors.The remaining integrals in the definitions of f q,g are easily evaluated numerically.
The hard, soft and jet functions H, S, J obey the renormalization group equations (RGEs), where F = H, S, J.The anomalous dimension γ i F is given by, Here κ H = 4, κ S = 4 2−a , κ J = − 2 1−a , j H,S = 1 and j J = 2 − a.Both the cusp and non-cusp anomalous dimensions can be expanded as, The coefficients Γ i n and γ i F n are given in App.B. For the two-loop jet functions for generic a, the cusp parts of the anomalous dimensions are known, while the non-cusp anomalous dimension for quark jets can be obtained from the hard and soft anomalous dimensions, by RG consistency, (2.21) The two-loop constant terms for a ̸ = 0 for quark jets were obtained in [60] from numerical computations of the QCD singular cross section using EVENT2 [98,99] together with the numerical results for soft function constants in [89].Newer, highly precise results for γ J (a) and 2-loop constants for a ̸ = 0 for quark jet functions have been presented in [100].Similar computations could be done for the gluon jet function for arbitrary a but lie outside the scope of this paper.
It is convenient to present the resummed results for the cumulative distribution, whose resummed form in momentum (τ a ) space is conveniently expressed in terms of the Laplace-transformed jet and soft functions acting with derivative operators on a resummation kernel [78,101,102]: Note that we use a form of the cusp evolution kernel proposed in [60] that keeps the resummed cross section explicitly independent of the factorization scale µ appearing in the original factorized cross section Eq. (2.4), at every order in resummed perturbation theory, While the other evolution kernels are defined as follows, Explicit expansions order by order for these kernels are given in [60].The sums of individual hard, jet, and soft evolution kernels Ω i , K i and K i γ used in Eq. (2.23) are defined, (2.26)

Fixed-order matching
In order to obtain a reliable prediction for large values of τ a , we also need to match our calculation to the full QCD fixed-order distribution.To O(α s ), the full QCD distribution is, We follow the method in Ref. [70] to calculate the coefficient A i a (τ a ) numerically.The detail of the calculation can be found in the appendix.Analytical results are only known for a = 0, found, e.g. in Ref. [37,55,103 The fixed-order angularity in SCET at O(α s ) is given by, (2.28) For quark final state, For gluon final state, where the function f g was given in Eq. (2.17).We define the difference away from τ a = 0 between full QCD and SCET as, Note that the coefficient α s /(2π) is not included in the above definition.We need to add this remainder function to the resummed distribution in order to obtain the NLL ′ + O(α s ) and NNLL + O(α s ) accuracy.The numerical results of r i a (τ a ) for values of a equal to {−1.0, −0.5, 0.0, 0.25, 0.5} are shown in Fig. 1.
The kink in Fig. 1 at τ max a = 1 3 1−a/2 is due to the fact that the full QCD distribution vanishes above the maximum kinematic value of τ a and only the singular part will give a contribution above this value.As an important consistency check of this matching technique, we integrate over τ a to get the fixed-order correction for the Higgs partial decay width.The total partial decay width can be written as, (2.32) The numerical results of Γ i H1 are summarized in Table .1.It shows our results are in agreement with Ref. [104] for all the values of a of we are interested in.
Figure 1.The O(α s ) nonsingular remainder functions for Higgs decaying to quarks (orange) and gluons (purple).Their relative size compared to the singular and total fixed-order cross sections is shown in Fig. 2.
H→qq, a=-1 The singular and nonsingular distributions for a = (−1.0,0.0, 0.5).The value of τ a where the singular (black) and nonsingular (blue) cross determines the transition point τ a = t 2 (a) between resummation and fixed-order regions that will be used in the scale profile functions later.

Nonperturbative shape function
The soft function in the factorization theorem will also receive corrections from nonperturbative hadronization effects.It can be parameterized into a soft shape function f mod (k) as [105][106][107], Here S i PT is the perturbative soft function and ∆i a = ∆i 1−a is a gap parameter with ∆i ∼ Λ QCD and ∆i is an a-independent parameter.This scaling of ∆i a tracks the known scaling 1/(1 − a) of the first moment of the shape function [62,108].The shape function can be expanded in a complete set of orthonormal basis functions [109], where (2.35) Here P n are Legendre polynomials.In our calculation, we will only keep b 0 = 1 and set b n = 0 for n > 0. The parameter λ i corresponds to the first moment of the f mod .However, it has been shown in Ref. [107] that the gap parameter ∆a has a renormalon ambiguity.To ensure good perturbative convergence of the soft function Eq. (2.33) and resulting cross section, as well as a stable definition of ∆i a , it is necessary to subtract/add a series removing this ambiguity from both S PT and ∆i a .In order to cancel the ambiguity, we can split the gap parameter as, ∆i where δ i a can be perturbative expansion with a same renormalon ambiguity as S i PT and ∆ i a is a renormalon free parameter.The R is the subtraction scale, which is defined through, where ν = ν a /m H , the same scheme adopted in, e.g.[39,60,110].There are multiple other schemes to define the series δ i a , see e.g.[111], and one could study the quantitative effects of varying schemes (cf.[112]), though this lies outside the scope of this paper.The shift Eq.(2.36) results in the renormalon-free soft function, The non-perturbative effects will shift the perturbative cross section.The shift parameter Ω i (µ, R) is given by,

.39)
In our calculation, we choose Ω with the reference scale R ∆ = 1.5 GeV.We also assume the Casimir scaling for the parameter This is purely an assumption, though probably fairly good, that we make to simplify our illustrative analysis below.A definitive study should measure Ω q,g separately.The gap parameter ∆ i a can be evolved to any other subtraction scale R and soft scale µ S , using the formalism in [113,114], formulas also summarized in Ref. [60].

Scales in resummation
From the arguments of the logarithms in the hard, jet and soft functions, the canonical scales should be, (2.40) However, the canonical scales do not properly take into account the transition from the resummation region into the fixed-order region where τ a is not small or into the nonperturbative region for τ a ≤ Λ QCD /m H .The use of profile scales has been proposed to smooth the transition between those different scale regions [39,109] and various specific forms for these profile functions are possible, e.g.[39,60,97,115].We use: where e H simply controls the normalization of the hard scale, while e J,S control variations of the shape of the jet and soft scales as a function of τ a above and below their default, canonical shapes.The ranges over which we vary them are given in Table 2. (Actually we keep e S = 0 as variations of other parameters below will make varying e S redundant.)The running scale is defined as, (2.42) The function ζ smoothly transitions between regions, and is defined as where the coefficients are . (2.44) The parameters t i are used to control the transition of different regions and we parameterized ), determined by the crossing of singular and nonsingular parts of the cross section in Fig. 2. The blue triangle is from H → q q, while the red box is for H → gg.The black line is the fitting function t 2 = 0.295 1−0.637a .them as, Here −1 is the angularity of the spherically symmetric configuration.The parameters t 0 and t 1 control the transition between the non-perturbative and resummation regions.In the nonperturbative region, we freeze the scale at µ run ∼ 1 − 3 GeV to allow the a shape function to smoothly describe this region and avoid the Landau pole in α s (µ).The parameter t 2 was determined by the point where singular and nonsingular contribution become comparable.As an example, we show the value of t 2 for a = (−1.0,0.0, 0.5) in Fig. 2, and it is determined by the crossing points of blue and black lines.It shows the t 2 is almost same for H → q q and H → gg, see Fig. 3.The particular functional form in Eq. (2.45) is the same as that found for 1-jettiness in DIS in [115], and the same form happens to fit the transition points plotted in Fig. 3 as well.Therefore, we will use the same t 2 formula for both Higgs decaying to quark and gluon states.Now, the exact values of t 0,1,2,3 are somewhat arbitrary, thus we vary them using the parameters n 0,1,2,3 to estimate theoretical uncertainties, across ranges given in Table 2.
For the renormalon subtraction scale, we choose R(τ a ) = µ S (τ a ) with µ 0 → R 0 and initial value of R 0 < µ 0 for τ a < t 0 .The remainder function depends on another scale µ ns which is used to estimate the higher-order effects from fixed order calculation [41], The values we pick for R 0 for quark and gluon τ a distributions are also given in Table 2.

Numerical results
Below, we present the numerical results of angularity distributions from Higgs boson decaying to quarks and gluons at NLL ′ + O(α s ) and NNLL + O(α s ) accuracy.In order to obtain a comprehensive theory uncertainty, we should consider all the scales and parameters from the profile function in our calculation.For a conservative estimation, we use the scan method to calculate the uncertainty band [39], and take 64 random selections of profile and scale parameters within the ranges shown in Table 2.Note that the soft scale parameter e S has been fixed to zero and the variation of soft scale µ S is from other parameters, e.g.n 0 , n 1 , µ 0 , r.In Fig. 4, we show the integrated distribution Γ i Hc (τ a ) for five values of angularity parameter a = (−1.0,−0.5, 0.0, 0.25, 0.5) to NLL ′ + O(α s ) and NNLL + O(α s ) accuracy.Note that the total partial decay width Γ i Htot is defined with the renormalization scale µ = m H . From NLL ′ + O(α s ) to NNLL + O(α s ), the scale uncertainties can reduce a little and the correction e H e J e S n 0 (q) n 0 (g) 0.5 ↔ 2 −0.5 ↔ 0.5 0 1 ↔ 2 GeV 2.8 ↔ 3.5 GeV n 1 (q) n 1 (g) n 2 n 3 µ 0 (q) 8.5 ↔ 11.5 GeV 25 ↔ 28 GeV 0.9 ↔ 1.1 0.8 ↔ 0.9 1.0 ↔ 1.2 GeV µ 0 (g) R 0 (q) R 0 (g) r 2.2 ↔ 3.0 GeV µ 0 (q) − 0.4 GeV µ 0 (g) − 1.8 GeV 0.75 ↔ 1.33 Table 2.The parameter ranges for our theoretical uncertainty estimation [60].Note that we use the different parameters n 0,1 , µ 0 and R 0 for quark and gluon final states.and NNLL + O(α s ) accuracy (green and orange for H → q q and light blue and purple for H → gg), including a nonperturbative shape function, compared to LO (red), NLO (gray) and approximate NNLO (blue).The fixed order prediction is from [55].Note that the normalization factor Γ Htot is defined in NLO.
is sizable for small a, e.g. a = −1.It is also clear that the quark and gluon decay modes from Higgs boson show different shapes in the small τ a region since the different color structure of quark and gluon.Therefore, the precise study of the angularity event shapes offers the possibility to distinguish the quark and gluon jets from Higgs boson decay and further can be used to constrain the light quark Yukawa couplings.
The differential distribution of angularities can be obtained by taking the derivative of Eq. (2.23); see Fig. 5 for the numerical results.We note that the gluon distributions peak at much larger values compared to the quark case.It could be understood from the Casimir scaling C A /C F in Sudakov factor.The distributions from gluon are also much broader compared to the quark cases due to the stronger QCD radiation.We also find that the larger a would shift the peak to a larger value.This is because varying a will change the proportions of two-jet-like events and three-or-more-jet-like events and further to change the peak position of the τ a distributions.In Fig. 7, we compare our theoretical calculation of Higgs decaying to quarks and gluons to PYTHIA [116] at parton and hadron levels.It shows that PYTHIA predictions at hadron level agree with our theoretical calculation very well, but the parton level results do not in the nonperturbative region at very small τ a .We should note that the design of profile functions is somewhat of an art.The choices of the parameters are based on obtaining properties of the theoretical predictions such as smoothness and convergence of uncertainty bands.
Since the thrust from Higgs boson decay has been calculated up to approximate NNLO (i.e.full NLO and singular NNLO) accuracy, it is also useful to compare our resummed results with the fixed order prediction in Ref. [55]; see Fig. 6.It shows that that our resummed prediction (which includes a shape function) agrees with the NNLO prediction very well in Higgs decaying to quark or gluon states at sufficiently large τ 0 , but there are deviations in the small τ 0 region, where especially for gluons, where the effects of resummation and of the nonperturbative shape function cause the distribution to peak at some value of τ 0 .However, if we focus on the region τ 0 ∈ [0.1, 0.20], our results agree with the approximate NNLO very well for both quark and gluon final states.
Our theoretical predictions show in Figs.4-7 include a soft shape function Eq. (2.38), with parameters chosen simply as described there.No attempt has been made to tune this for Higgs decay, they were guided simply by similar values they take in event shapes in e + e − → hadrons (e.g.[39,41]), for purposes of producing illustrative numerical predictions.When a serious comparison to data is performed, the data can be used to constrain these soft shape function model parameters and test properties such as the universal scaling properties of the leading nonperturbative shift governed by Ω q,g 1 [62,108].

Probing light quark Yukawa couplings
Next, we will use angularity distributions to probe light quark Yukawa couplings at CEPC and our results are easy to generalize to other e + e − colliders.At CEPC, the Higgs boson is dominantly produced through the Higgs and Z boson associated production, i.e. e + e − → HZ.
The signal we are interested in is Higgs decaying into q q with q = u, d, s and Z boson decaying into lepton pairs.The major SM backgrounds are H(→ gg, b b, cc)Z, Zq q and H(→ V V * → 4j)Z, with V = W ± , Z1 .To suppress the backgrounds from heavy quarks, we could use flavor tagging techniques.Ref. [42] has shown in this case that the heavy flavor backgrounds can be removed mostly if we require two non-b and c jets in the final state.The background Zq q can be highly suppressed after we include the kinematical cuts, e.g.recoil mass [22,117].As shown in Ref. [42], after the kinematical cuts and requiring two light jets in the final state, we could get the number of background events for the H → gg down to The background of H → V V * will contribute to the tail region of the angularity distributions and the number of events after including above analysis is N 4q b = 0.6 N g b [42].It is clear that the gluon background is the major obstacle for probing light quark Yukawa couplings.Therefore, we propose to use the hadronic angularity distributions of Higgs boson to separate the gluon background from the signal.We should note that the kinematical cuts in this section just modify the normalization but not the shape of the τ a distributions.Thus, we could use the normalized angularity distributions to study the Yukawa couplings without input any kinematical cuts for the signal and backgrounds.
In this work, we assume NP effects only change the branching ratio (BR) of Higgs decaying to light quarks, then we define the ratio, The 1-loop corrections Γ q,g H1 to the Higgs decay widths were given in Eq. (2.32) and Table 1.The total number of signal events is N s = R q N g b .In the SM, R q ≃ 0 due to the smallness of the quark mass.We divide the angularity into N bin = 10 bins and use the binned likelihood function to estimate the sensitivity for the Quark NLL ' +(αS) Gluon NLL ' +(αS) Quark NNLL+(αS) Gluon NNLL+(αS) . The normalized distributions of angularity from PYTHIA at hadron level for H → V V * → 4j with a = (−1.0,0.0, 0.5) (red lines).The orange and purple lines are from our theoretical prediction for H → q q and H → gg.The perpendicular dashed lines denote the analysis regions [τ min (a), τ max (a)] for probing light quark Yukawa couplings.
hypothesis with R q against the hypothesis with R q = 0 [118], where b i and n i are the number of the backgrounds and observed event in the ith bin, respectively, and s i (R q ) is the number of signal events in the ith bin for the parameter R q .We set n i = s i (0) + b i .The number of signal events in each bin is determined by our theoretical calculation at NNLL + O(α s ) accuracy.The number signal events in the ith bin is, The normalized angularity distributions from H → b b, cc are almost same with H → q q except for very small τ a region due to the mass effect of heavy quark.In order to avoid the possible heavy quark mass effect, the impact of non-perturbative function and the contamination of high tail backgrounds, we can very conservatively require τ min (a) < τ a < τ max (a) with τ min (a) = 15  125 × 3 a and τ max (a) = 0.295 1−0.637a (see the vertical dashed lines in Fig. 5).Therefore, we expect the angularity distributions from heavy quarks should be same with the light quarks in this region.The distributions from Zq q should have the same shape as those predicted in Ref. [60] at NNLL ′ + O(α 2 s ) accuracy.The background H → V V * can be estimated by event generator Madgraph [119] and PYTHIA [116] at the leading order; see Fig. 8.It is clear that the four-quark background from gauge boson pair is only sensitive to the tail region and far away from the peak of the signal.The total backgrounds in ith bin is, Here j = g, HF, ZZ, 4q denotes the different background processes.We define the test ratio of likelihood function, The parameter q yields the exclusion of the hypothesis 1 with R q ̸ = 0 versus the hypothesis 0 with R q = 0 at the qσ confidence level.Thus, where n i (R q ) = s i (R q ) + b i .For simplicity, we normalize the light quark Yukawa couplings to the SM bottom quark y b ≡ y b (m H ). Figure 9 displays the expected 1σ (green) and 2σ (orange) confidence level exclusion limit on Yukawa coupling y q ≡ y q (m H ) from the angularity distributions.It shows that the angularity event shapes, in the conservative τ a window we considered above, could give a fairly strong constraint for the light quark Yukawa couplings, i.e. y q ≲ 0.15y b (for a = −1) to y q ≲ 0.22y b (for a = 0.5), at the 2σ confidence level.
The theoretical uncertainties will change the upper limit of light quark Yukawa couplings from 5% to 14%.We note that the limit for y q obtained in this way is considerably larger than the results in Ref. [42].It could be understood from our analysis strategy, i.e. in this window, we have τ a > τ min (a), which is far away from the peak region of the signal, while it is not in Ref. [42].From the PYTHIA prediction (see Fig. 7), we know the typical nonperturbative hadronization effects will shift the peak of angularities by about ∆τ a ∼ 0.01.Furthermore, the nonperturbative corrections to our predictions for light quark angularities remain small (or within the universal shift region) to considerably smaller values of τ a .If we push our analysis region a bit more aggressively into the peak region, e.g. to a lower limit of τ a > t g 0 (a) ∼ 3.5/125×3 a , then we obtain the upper limit y q ≲ 0.09y b at 2σ confidence level for all our choices of a.We can try to be even more aggressive than this, and push the lower limit for τ a left of the peak of the quark angularity distributions, even though nonperturbative and heavy-quark mass corrections are larger there-this is for illustrative, motivational purposes only.Choosing the analysis region τ a > t q 0 (a) = 1/125 × 3 a , we obtain the stronger limit y q ≲ 0.07y b .We summarize the various limits we obtain on y q for the different analysis windows in Table 3.We stress that in this work we have used a simple educated guesses for the quark and gluon shape function parameters, to which this analysis region will be sensitive.If our guesses are close to accurate, though, the above limit is indicative of the power of a single angularity distribution to put a limit on y q .A definitive limit would require further control of the shape function parameters through universality/scaling arguments and/or extraction from an independent dataset.Alternatively, to further improve the measurements and reduce the impact of the non-perturbative effects, we could use the soft-drop grooming technique [120][121][122], and this could be considered in a future work.Nevertheless, the potential limits indicated in Table 3 show the promise of angularities in separating quark signal from gluon background to obtain strong limits on y q .
Since the angularities depend on the continuous parameter a, it is also useful to combine the multiple angularity distributions in the analysis.Here we build upon ideas previously developed in [123].As an example, we show the normalized double differential distributions of angularities with a = −1.0 and a = 0.5 from PYTHIA at hadron level for H → q q and 1σ CL exclusion 1σ CL exclusion+th.unc.2σ CL exclusion 2σ CL exclusion+th.unc.3.
H → gg in Fig. 10.Performing a likelihood analysis in the regions τ min (a) < τ a < τ max (a) with the same min and max limits for each a as in the single-angularity analysis above, we find that limits on the Yukawa couplings could be further improved about 2% ∼ 3% compared to the single angularity analysis at 2σ confidence level.Although the two different angularities will improve the power of distinguishing the quark and gluon jets, the backgrounds from H → b b, cc and Zq q will become important since both of them have a similar angularity distributions as the signal.Therefore, the next step to further constraint the light quark Yukawa couplings could be from improving the b-tagging efficiency, detector resolution and so on.If both H → b b, cc and Zq q could be further reduced about 10 times than current analysis, we estimate that the results from two different angularities could be improved about 12% compared to the single angularity analysis.(Note that these estimates are obtained in our most conservative analysis windows.)Methods for resummation of two angularities have been developed and applied to predictions in [72,73].Further development of such predictions for double differential distributions of angularities in Higgs decays and their applications to phenomenology would seem worthwhile.) from PYTHIA at hadron level for H → q q and H → gg.On the two side axes we show the corresponding single differential spectra.We analyse regions with the same conservative cuts τ min (a) < τ a < τ max (a) for each τ a as in the one-dimensional analyses used in Table 3 to estimate the sensitivity for the y q .

Summary
In this paper, we studied a class of event shape variables angularities from Higgs boson decay at the e + e − colliders based on soft-collinear effective theory.Both the quark and gluon final states are calculated to NLL ′ + O(α s ) and NNLL + O(α s ) accuracy.From NLL ′ + O(α s ) to NNLL + O(α s ) accuracy, we found a sizable correction for small a angularties.The differential angularity distributions also show a Casimir scaling C A /C F from quark to gluon.We compared the predictions resulting from this analysis to PYTHIA at both parton and hadron level, and find good agreement for hadron level results.
Based on the difference of angularity distributions between quark and gluon final state, we proposed to test the light quark Yukawa couplings through angularity distributions at lepton colliders.As an example, we show that the CEPC with √ s = 250 GeV and an integrated luminosity of 5 ab −1 could give a constraint y q ≲ 0.15 − 0.22y b for q = u, d, s, for different values of the angularity parameter a, at 2σ confidence level using a conservative analysis region [τ min (a), τ max (a)] far away from small τ a where hadronization and b-mass effects are larger.The theoretical uncertainty for this upper limit is around 10%.In order to further improve the results, it is important to extend the analysis region to small τ a .It shows that the upper limit could be reduced to y q ≲ 0.09y b , similar to [42], if we push the analysis region down to [t g 0 (a), τ max (a)], or even smaller y q ≲ 0.07y b in the region [t q 0 (a), τ max (a)].However, the theoretical prediction in the small τ a region is strongly dependent on the nonperturbative model assumptions and this issue could be overcome by gaining better knowledge of the gluon soft shape function, in particular, or by utilizing soft-drop grooming techniques.Utilizing multiple angularities at once also shows promise in improving the potential limit on y q further.Note: As this paper was being finalized, Ref. [61] appeared, computing Higgs angularity distributions to NNLL ′ resummed and NLO fixed-order accuracy.We have not yet compared to the technical results of this paper; here we have focused more extensively on the phenomenological application of the results to the determination of Yukawa couplings.
There are two subprocesses H → gg(→ q q) and H → ggg that can contribute to A g a (τ a ).For H → gg(→ q q) channel we have,

)). (A.4)
The coefficient A g a (τ a ) = A gq q a (τ a ) + A ggg a (τ a ).The thrust axis in these 3-body final states is always along the direction of the particle with the largest energy.Therefore, the phase space in the x 1 , x 2 plane can divided into three regions, as shown in Fig. 11.The value of the angularity τ a (x 1 , x 2 ) is dependent on which parton has the largest energy.In a region where x i > x j,k , the angularity τ a (x 1 , x 2 ) is given by, τ a (x 1 , x 2 ) (A.5) In the following, we will use phase space region C where x 3 > x 1,2 as an example to discuss the calculation of A i a (τ a ).The integration in phase space A and B can be related to the integration over region C by the appropriate change of variables.The angularity in phase space C can be written as, (A.6) It turns out to be convenient to change the integration variables from x 1 , x 2 to τ a and w, where w = 1 − x 1 2 − x 1 − x 2 . (A.7) Therefore, x 1 and x 2 are given by, x 1 (w, τ a ) = 1 − w + w τ a w 1−a/2 (1 − w) a/2 + w a/2 (1 − w) 1−a/2 The integration region of w is determined by solving these inequalities, which are equivalent to the conditions, τ a = F a (w min ), τ a = F a (1 − w max ), (A.10) whose solutions give the lower and upper limits of the integral over w, w min and w max , respectively, where the function F a is given by (1 + w) 1−a/2 (w 1−a + (1 − w) 1−a ).(A.11) These limits are themselves functions of τ a , and cease to have a solution above the maximally allowed value of τ a = τ max a , which is τ max a = (1/3) 1−a/2 at O(α s ). 2 Then, for example, the integral for A q a in Eq. (A.2) is expressed: w min (τa) dw J(w, τ a ) (1 − x 1 )(1 − x 3 ) x 1,2 =x 1,2 (w,τa) , (A.12) where x 1,2 are expressed in the form Eq. (A.8), and x 3 = 2 − x 1 − x 2 .The factor J is the Jacobian associated with the variable transformation Eq. (A.8), J(w, τ a ) = (A.13) The endpoints w min,max in Eq. (A.10), the Jacobian Eq. (A.13), and the integral Eq. (A.12) are all easily evaluated numerically, which is how we have computed the A i a 's and the resulting remainder functions illustrated in, e.g.Fig. 1.The gluon channel contributions A g a given by Eqs.(A.3) and (A.4) are expressed and evaluated similarly to Eq. (A.12).

B Ingredients for NNLL resummation
In this Appendix we collect results needed to compute angularity distributions to NNLL accuracy.The coefficients of β function in the MS scheme are given by [124][125][126], ) and the coefficients of the anomalous dimension γ y of the Yukawa coupling y q in Eq. (2.12) are given by [127]: The cusp anomalous dimensions up to 3-loop order [128,129] Γ i 0 = 4C i , The non-cusp anomalous dimension for the hard function up to 2-loop level [84,101,[130][131][132][133], γ q H0 = −12C F , γ q H1 = −2C F  which vanish for a = 0.The integral representations can easily be evaluated numerically to high accuracy for any value of a, and the relevant values for our work are given in Table 4.

Figure 3 .
Figure 3.The value of t 2 (a) with a = (−1.0,−0.75, −0.5, −0.25, 0.0, 0.25, 0.5), determined by the crossing of singular and nonsingular parts of the cross section in Fig.2.The blue triangle is from H → q q, while the red box is for H → gg.The black line is the fitting function t 2 = 0.295 1−0.637a .

Figure 4 .
Figure 4. Integrated angularity distributions for values of a = (−1.0,−0.5, 0.0, 0.25, 0.5) at NLL ′ + O(α s ) (green for quark and blue for gluon) and NNLL + O(α s ) (orange for quark and purple for gluon), including shape function effects.The theoretical uncertainties have been estimated with the scan method.

a=- 1 Figure 5 .Figure 6 .
Figure 5. Differential angularity distributions for values of a = (−1.0,−0.5, 0.0, 0.25, 0.5) at NLL ′ + O(α s ) (green for quark and blue for gluon) and NNLL + O(α s ) (orange for quark and purple for gluon) for Higgs decay, including shape function effects.The theoretical uncertainties have been estimated with the scan method.The perpendicular dashed lines denote the analysis regions [τ min (a), τ max (a)] for probing light quark Yukawa couplings.
GeV with an integrated luminosity of 5 ab −1 .The background of Higgs decaying to heavy flavor quarks (H → b b, cc) is about N HF b = 0.1N g b , and the number of Zq q events is N ZZ b

Figure 10 .
Figure 10.The normalized double differential distributions of two jet angularities (τ −1.0 , τ 0.5 ) from PYTHIA at hadron level for H → q q and H → gg.On the two side axes we show the corresponding single differential spectra.We analyse regions with the same conservative cuts τ min (a) < τ a < τ max (a) for each τ a as in the one-dimensional analyses used in Table3to estimate the sensitivity for the y q .

Figure 11 .
Figure 11.The three body phase space.The parameter x i is defined asx i = 2E i /m H and x 1 + x 2 + x 3 = 2.

Table 1 .
]. NLO QCD correction for the Higgs partial decay width for selected values of a.