Energy-energy correlation in hadronic Higgs decays: analytic results and phenomenology at NLO

In this work we complete the investigation of the recently introduced energy-energy correlation (EEC) function in hadronic Higgs decays at next-to-leading order (NLO) in fixed-order perturbation theory in the limit of vanishing light quark masses. The full analytic NLO result for the previously unknown EEC in the $H \to q \bar{q} + X$ channel is given in terms of classical polylogarithms and cross-checked against a numerical calculation. In addition to that, we discuss further corrections to predictions of the Higgs EEC event shape variable, including quark mass corrections, effects of parton shower and hadronization. We also estimate the statistical error on the measurements of the Higgs EEC at future Higgs factories and compare with the current perturbative uncertainty.

worth noting that by making use of the AdS/CFT duality in N = 4 SYM one can also obtain a strong-coupling limit result for the EEC [26,27]. In QCD, the collinear limit of EEC can be understood by using the recently available factorization theorem [28], which also improves the resummation beyond the leading logarithmic (LL) accuracy [29,30]. The back-to-back limit features an all-order factorization formula [31] that makes use of the transverse-momentum dependent (TMD) factorization [32][33][34]. In this limit, resummed predictions are currently known at the N 3 LL accuracy [31,[35][36][37].
The answer to the question raised in the above paragraph lies in the very definition of the energy-energy correlator. As we will see below, the Dirac delta that introduces correlations between energies of partons or final state hadrons can be straightforwardly converted to a loop-momentum dependent (albeit nonlinear) propagator and subjected to the standard methods of computing higher-order corrections, such as integration-by-parts (IBP) reduction [38,39] and differential equations [40][41][42][43][44][45]. Moreover, these steps can be carried out using off-the-shelf software packages for loop computations: The specifics of our observable (e.g. custom IBP equations for loop integrals with nonlinear propagators) can be encoded in the Mathematica scripts used to invoke the existing tools, so that the tools themselves do not require any modifications. The existence of numerical NNLO results [36,46] (making use of the CoLoRFulNNLO method [47][48][49]) as well as the availability of public codes (e.g. Event 2 [50,51], NLOJet++ [52,53], Eerad3 [54]) capable of evaluating the EEC numerically greatly facilitate the cross-checks of new analytic results.
It is important to stress that when speaking of "EEC" we do not limit ourselves to the original definition of this event shape variable for electron-positron annihilation to partons via the reaction e + e − → qq + X. For example, Transverse-Energy-Energy Correlations (TEEC) [55] have already been studied in the context of proton-proton [56] and electronproton [57] collisions. The back-to-back limit of TEEC can be investigated using recently obtained factorization theorems for hadron-hadron [58] and electron-hadron [59] colliders.
Recent considerations of the three-point [60,61], four-point [62] and multi-point energy correlators [63,64] as well as two-point gravitational energy correlators [65] represent further exciting extensions of the original EEC concept and signalize an increased interest of the theorist community in such novel event shape variables.
For phenomenological purposes, EEC can be employed as a tool to determine the value of the strong coupling constant (cf. e.g. [66] for a recent study) by comparing the available theoretical predictions to the existing electron-positron collider measurements. In [18] it was suggested that a new event shape variable, denoted as the Higgs EEC, could provide an intriguing connection between the strong and the Higgs sectors by defining an observable equally accessible to experimentalists analyzing the data from a future Higgs factory and to theorists calculating the corresponding predictions. Furthermore, this observable could be potentially used for the purpose of α s determinations from hadronic Higgs decays. A high-energy lepton collider, be it CEPC [67,68], ILC [69,70], FCC-ee [71] or CLIC [72,73], would be capable of copious production of Higgs bosons in the clean environment of e + e −annihilations. It is, therefore, not unreasonable to expect that in the future we might witness a high precision measurement of the Higgs EEC using data collected at a leptonic Higgs factory.
The analytic NLO results presented in [18] concerned only the H → gg + X channel calculated in the Higgs Effective Theory (HEFT) [74][75][76][77] with massless quarks. The goal of this work is to present analytic results also for the channel H → qq + X, thus completing the fixed-order investigation of the Higgs EEC at NLO. Being the largest Higgs decay branching ratio, Higgs decaying into bottom quarks has received much attention from the theory community. For example, the partial decay width of H → qq has been calculated to N 4 LO [78][79][80], and the fully differential decay width for the same process is known to N 3 LO [81][82][83] for massless quarks and to NNLO [84] for massive quarks. Some interesting results obtained very recently are the calculation of Higgs decaying into two bottom quarks and an additional jet at NNLO [85], the study of the Higgs decay into four bottom quarks at NLO [86] and the investigation of the thrust distribution for Higgs going into a pair of bottom quarks or gluons plus an additional jet at NLO and approximate NNLO [87].
For the sake of clarity, in the following we will denote the H → gg +X and H → qq +X contributions as Hgg EEC and Hqq EEC respectively. The Higgs EEC is then understood to contain both channels. The original observable from [6] will be referred to as the standard EEC.
Following [18], we define the Higgs EEC as with Γ tot being the total decay width for H → hardons, whereas dΓ a+b+X describes the differential decay rate of a Higgs decaying into two hadrons plus anything else. Furthermore, we have cos θ ab =p a ·p b , where (E a , p a ) T and (E b , p b ) T denote the 4-vectors of the hadrons a and b respectively. Finally, χ is the angle between two calorimeters measuring the energies of a and b, while m H stands for the Higgs boson mass. By summing over all available final state hadron pairs (a, b) and weighting their contributions to the energy flow by the product of their energies divided by the square of the Higgs mass, we obtain a differential angular distribution normalized to unit area.
To calculate the Higgs EEC in perturbation theory we replace the hadrons by partons and exclude self-correlations, so that the contributions with a = b are removed from the summation in eq. (1.1). The interacting part of the relevant Lagrangian reads where the first term stems from the HEFT with λ being the corresponding Wilson coefficient (known up to N 4 LO [88]). The second term is the Standard Model Yukawa interaction for quarks, with y q being the Yukawa coupling for the quark flavor q.
To facilitate the analytic calculation we choose to work in the massless quark limit, while keeping nonvanishing Yukawa couplings. The top quark contributions are thus omitted and we have only 5 active quark flavors. As has already been observed in [87], the chiral symmetry of massless QCD ensures that in this approximation there is no interference between the H → gg + X and H → qq + X channels. The respective operators also do not mix under the renormalization so that both pieces can be treated separately. Since the gluonic channel has already been computed in [18], our sole remaining task is to calculate the contribution from Higgs decaying to a quark-antiquark pair and one or two additional partons. The 3-parton final state corresponds to the LO result, while the 4-parton states are needed for the NLO.
We normalize the Hqq EEC contribution with respect to the total decay width for H → qq given by where C A stands for the number of colors and K(µ) encodes higher order corrections in α s . The K-factor for H → bb in the limit where the bottom mass is set to zero is currently known at O(α 4 s ) [78][79][80], and the full scale dependence up to O(α 3 s ) can be found in [89]. This normalization prescription ensures that Hqq EEC does not depend on y q , while the dependence on m H enters only through log(µ/m H ) and vanishes for the renormalization scale choice µ = m H .
Our paper is organized as follows. We describe the technical details of our Higgs EEC calculation for the H → qq + X channel in section 2 and subsequently present the obtained analytic results (including the asymptotic behavior in the collinear and back-toback limits) in section 3. Section 4 explores the phenomenological implications of the Hqq EEC observable. Finally, our conclusions and possible future extensions of this work are summarized in section 5.

Technical framework
Our calculation essentially follows the path that has already been outlined in [17] and explained in details in [18], so that we keep the following description short.
First of all, we need to obtain matrix elements squared |M(H → qq + X)| 2 for real, double-real and real-virtual corrections to the Higgs decaying into a quark-antiquark pair. The real and double-real contributions follow directly from squaring the corresponding tree-level amplitudes with 3-or 4-parton final states respectively A visualization of the double-real contributions using the cut diagram notation is shown in figure 1. Working in the rest frame of the decaying Higgs particle, we have Q = (m H , 0, 0, 0) T . The real-virtual piece follows from the interference of the tree-level and 1-loop 3-parton final states. The Higgs EEC observable without the overall normalization factor is obtained by multiplying |M(H → qq + X)| 2 with the measurement function Since the real-virtual piece involves only a massless 3-particle phase space, it is sufficiently simple to be integrated directly via HyperInt [90]. However, the NLO double-real contribution leaves us with a large number of complicated and badly divergent 1 phase-space integrals. We choose to handle them by employing the method of reverse unitarity [91,92] which effectively trades the measurement function for the following nonlinear cut propaga- The occurring loop integrals can then be reduced using IBP techniques. The resulting master integrals can be solved via differential equations by finding a canonical form [93] for each of the systems and then determining the integration constants using suitable boundary conditions. In practice, we generate the Higgs decay amplitudes using QGRAF [94] and Fey-nArts [95]. FeynCalc [96][97][98], FORM [99] and Color [100] are used to prepare the squared matrix elements, evaluate them in d-dimensions and carry out the color algebra.
We also employ FeynHelpers [101] and Package-X [102,103] for the calculation of the real-virtual matrix element. To avoid dealing with ghost contributions we make use of the axial gauge when summing over the gluon polarizations. The obligatory topology identification step proceeds by considering all possible ways to exchange loop momenta p a ↔ p b or to perform a shift p a → Q − b =a p b . Notice that the invariance of the sum of the measurement functions for different partons under these manipulations lead to a significant simplification of this task. In the first step, instead of looking at the full integrand it is convenient to omit the Dirac delta from the measurement function and enumerate the occurring subtopologies. In the second step we augment each identified subtopology with the corresponding nonlinear cut propagator. In the case of a 4-parton final state, one subtopology gives rise to 6 integral families, stemming from the parton pairs (1, 2), (1, 3), (1, 4), (2, 3), (2,4) and (3,4). The subprocesses with qqg, qqqq and qqq q final states contain only one subtopology each, given by and respectively. The most complicated double-real piece stemming from the qqgg final state involves 3 following subtopologies that lead to 18 integral families. The search for a minimal set of subtopologies as well as the generation of the final integral families is done using in-house Mathematica scripts. Custom codes written on top of FeynCalc and LiteRed [104] are used to handle linearly dependent propagators via partial fraction decomposition and to derive symbolic equations for the IBP reduction. Then, the IBP-reduction is carried out with FIRE [105,106], where we submit our custom IBP equations to the program via the variable startinglist and mark all cut propagators through the RESTRICTIONS setting. Finally, we map the obtained master integrals to the set of integrals that was calculated in [17]. Just as in the case of the Hgg EEC, we find no new master integrals that cannot be expressed as a linear combination of masters from the standard EEC integral basis at NLO.
Upon adding all contributions together and carrying out the UV-renormalization of the real-virtual contribution, we end up with a manifestly finite result, as expected from the IR-safe property of the EEC event shape variables.

Analytic results at NLO
The main result of this work is the analytic expression for the Hqq EEC at O(α 2 s ) given by with N c being the number of colors. The overall prefactor 1/K(µ) stems from the normalization prescription given in eq. (1.3), while A Hqq (z) and B Hqq (z) denote the LO and NLO coefficients respectively. One may wonder why the coefficient of log µ m H in the numerator of eq. (3.1) is proportional to β 0 + 6C F . The origin of this term can be traced back to the usual strong coupling constant renormalization and the additional Yukawa renormalization [87,107], Notice that K(µ) to order O(α s ) is given by Using eq. (3.3), one could expand eq. (3.1) to O(α 2 s ), obtaining a result with the coefficient of log µ m H being exactly proportional to β 0 . The LO piece is directly proportional to C F and can be written as The NLO coefficient B Hqq (z) can be decomposed into where B Hqq,lc (z), B Hqq,nlc (z) and B Hqq,N f (z) stand for the leading color, next-to-leading color and the N f pieces respectively. The color structure of the NLO coefficient is identical to the one observed in the standard EEC. This is not surprising, as both observables are quark-initiated quantities.
The analytic structure of the color coefficients precisely follows the pattern known from the standard EEC and the Hgg EEC. We again find the same set of building block functions g (j) i , were j denotes the pure transcendental weight The coefficients of these functions (except for some coefficients multiplying g 3 ) are rational polynomials of the form Every color component also contains a term proportional to 1/z 7/2 g 3 and those pieces are symmetric under √ z → − √ z, which appears to be a universal feature of EEC observables at NLO [17,18,21]. On the other hand, it is interesting to observe that the highest power of z in the numerators of the rational polynomials is only 7, at variance with 8 in the case of the Hgg EEC and 9 for the standard EEC at NLO. Furthermore, the largest value of the power k being 5 is true also for the standard EEC, while it can go up to 6 for the Hgg EEC. Presumably, these small differences between the Hqq EEC and the standard EEC can be largely attributed to the different vertex structures of L int : the former is initiated through a scalar-fermion coupling, while the latter starts via a vector-fermion interaction.
The analytic results for the separate color components at NLO read as follows 4 .
In the same manner we can also explore the back-to-back limit. Notice that the presence of large logarithms from soft and collinear emissions signals the necessity of a proper resummation using the existing techniques [31,32,58,108]. Expanding around z = 1 we find where the leading power terms can be also obtained using the formalism of [31].

Phenomenological applications
In the following we present a brief discussion on phenomenological applications of the Hqq EEC event shape variable in Higgs boson decays. We verify our analytic formulas by comparing them to a numerical result that was obtained using Monte Carlo (MC) integration. In the numerical calculation we used independent matrix elements that were automatically generated with GoSam 2.0 [109], while the real corrections were treated using the dipole subtraction method [110]. We set the strong coupling constant to α s (M Z ) = 0.1181 in the calculations. Analytic and numerical results for the Hqq EEC at LO and NLO are shown in figure 3, where the underlying process is the decay of the Higgs into massless quarks and all distributions are normalized to the total partial width at LO. To simplify the comparison and improve the visual quality of the plot, we choose slightly different cos χ values for the curves describing the analytical and numerical distributions. As can be inferred from the plot, within the MC errors we find a perfect agreement between our analytic and numerical predictions both at LO and NLO.
A direct comparison to the future experimental data requires additional corrections to the fixed-order theory prediction, which we discuss below.
First of all, the effect of the finite bottom-quark mass m b can be nonnegligible in a fixed-order calculation. Since bottom quarks are treated as massless in our analytic result, it is important to estimate the impact of this simplification. To this end, we performed another NLO numerical calculation of the Hqq EEC, where the bottom-quark was treated as massive with m b = 4.78 GeV. The ratio of the Hqq EEC results with massive and massless quarks is shown in the upper panel of figure 4. We observe that in the range of | cos χ| < 0.95 the bottom-quark mass corrections reduce the distribution by about 4% in the bulk region and can reach 10% to 20% in the back-to-back and collinear regions. The latter is not surprising, as it is well known that collinear radiations are suppressed due to the finite quark mass. Second, for a meaningful Higgs EEC prediction we must also include parton shower and hadronization corrections. To account for that, we match our NLO calculation with massive bottom quarks to parton shower using POWHEG-BOX-V2 [111,112] and PYTHIA 8.2 [113]. In the PYTHIA setup we use the Monash tune [114] and, for the sake of simplicity, force all B hadrons to be stable. Figure 4 shows that the parton shower can substantially enhance the EEC distribution in the whole cos χ range, with the corrections amounting to almost 40% in the collinear region. This observation hints that fixed-order NNLO QCD corrections to the Hqq EEC could be potentially large. Furthermore, the hadronization corrections are equally significant and can reach more than 10%.
Finally, we estimate the perturbative uncertainty of the matched NLO predictions by varying the renormalization scale and the square of the parton shower scale independently by a factor of two around their nominal values, chosen as m H for the renormalization scale and k 2 T for the square of the shower scale. We add the two scale variations in quadrature and plot the uncertainty band in the lower panel of figure 4. The total uncertainty from the scale variations lies between 5% and 10% in the plotted region. Given the existence of NNLO numerical calculations for massless bottom quarks [85], we expect that this uncertainty can be, in principle, significantly reduced in the future.
In addition to that, the lower panel of figure 4 also contains the projected experimental uncertainties. In this case we incorporate only the statistical errors and assume the total number of events being 4 × 10 5 , which corresponds to the number of H → bb decays that CEPC [115] is expected to collect during its first 7 year data taking period. We estimate the statistical errors by first generating 40 ensembles of events and then calculating the standard deviation of the EEC in each bin from the values predicted by all ensembles. This procedure is meant to account for the strong statistical correlations among different bins that typically arise when studying EEC-like observables: a single event generates multiple histogram entries, hence simultaneously contributing to many bins. In our case, for all bins the observed uncertainties constitute at most 0.5%. Experimental systematic errors, that we choose to ignore here, can be attributed to the signal extraction from the SM background as well as event reconstruction and detector resolution. Although they are expected to be dominant over the statistical errors, a thorough estimation of these uncertainties is beyond the scope of the present paper.

Summary
Higgs EEC is a novel event shape variable that can be measured by reconstructing 4vectors of the final state particles originating from hadronic Higgs decays. This observable opens an interesting perspective of α s determinations from Higgs precision measurements at future Higgs factories and is therefore of great relevance for experimentalists interested in exploring Higgs phenomenology at high-energy lepton colliders.
In this work we employed methods pioneered in [17] to calculate Higgs EEC in the H → qq + X channel at NLO in the fixed-order perturbation theory. This result can be combined with the already available computation in the H → gg +X channel [18] to obtain the full Higgs EEC in the limit of vanishing light quark masses. The analytical structure of the Hqq EEC is very similar to that of the Hgg EEC and the standard EEC: All 3 results can be calculated using the same set of master integrals and written in terms of the same building block functions that involve classical polylogarithms up to weight 3.
As far as the phenomenology of the Hqq EEC is concerned, we employed numerical methods to study the importance of the effects missing in the analytic calculation: finite bottom-quark mass, parton shower and hadronization. On the one hand, the corrections due to the finite bottom-quark mass turn out to be numerically rather small, apart from the collinear and back-to-back regions. On the other hand, parton shower and hadronization effects can lead to enhancements of tens of percents beyond the NLO fixed-order predictions. The remaining scale variations are at the level of 5% to 10%. At the same time, the projected statistical uncertainties on the measurements of the Higgs EEC at future Higgs factories are at sub-percent level. Therefore, we conclude that improved perturbative calculations and a more accurate modeling of the hadronization are mandatory in order to match the future experimental precision.
Theoretical investigations of EEC-like observables continue to expand our understanding of the mathematical underpinning of perturbative QCD. The multitude of results made available in the recent years corroborate that the study of EEC has become a very active field of research within the phenomenology of the strong interactions at high energies. Even though every new calculation raises the bar a bit higher, there is obviously still a lot of work left to be done. At NLO one could consider other underlying processes that lead to hadronic decays or try to incorporate effects of massive quarks, while at NNLO we still lack the full fixed-order result even for the standard EEC. Taking a broader view, it would be very rewarding to search for techniques that could enable us to obtain NLO analytic results for event shape variables other than the EEC. Given the amount of progress in the field made in the last few years, we may very well expect to witness even more exciting findings in the years to come.
2 + g  where r = √ z/ √ 1 − z . The function g (2) 5 is real-valued in z = ir,z = −ir and is known as the Bloch-Wigner function. The building block functions in eq. (B.3) and their specific combinations in eq. (B.2) are identical to those appearing in the B qq int term of the standard EEC [17]. This is not surprising, since both the standard EEC and the Hqq EEC are quarkinitiated observables.
As far as the asymptotics is concerned, in the collinear limit B Hqq ,qq int reads We observe that the leading power terms (i.e. the coefficient of 1/z in eq. (B.4)) are the same as the corresponding terms of the standard EEC [17]. These terms can be predicted by the jet calculus approach at the next-to-leading logarithm (NLL) accuracy [116]. Looking at the back-to-back limit by expanding B Hqq ,qqint up to next-to-leading power we again find that the coefficient of 1/(1−z) precisely reproduces the corresponding results for the standard EEC [17].