Analytic Next-To-Leading Order Calculation of Energy-Energy Correlation in Gluon-Initiated Higgs Decays

TL;DR

This paper presents an analytic NLO calculation of the energy-energy correlation in gluon-initiated Higgs decays within the Higgs Effective Field Theory, providing results useful for future collider measurements and strong coupling determinations.

Contribution

It provides the first full analytic NLO calculation of the EEC in gluon-initiated Higgs decays in the HEFT framework, including asymptotic behaviors.

Findings

Analytic NLO results expressed with classical polylogarithms.

Asymptotic behavior analyzed in collinear and back-to-back limits.

Potential for future collider measurements to probe strong and Higgs sectors.

Abstract

The energy-energy correlation (EEC) function in $e^+e^-$ annihilation is currently the only QCD event shape observable for which we know the full analytic result at the next-to-leading order (NLO). In this work we calculate the EEC observable for gluon initiated Higgs decay analytically at NLO in the Higgs Effective Field Theory (HEFT) framework and provide the full results expressed in terms of classical polylogarithms, including the asymptotic behavior in the collinear and back-to-back limits. This observable can be, in principle, measured at the future $e^+e^-$ colliders such as CEPC, ILC, FCC-ee or CLIC. It provides an interesting opportunity to simultaneously probe our understanding of the strong and Higgs sectors and can be used for the determinations of the strong coupling.

Figures

Tables

Table 1. Table 1 : Fit of eq.  40 to the Pythia data using ROOT with Minuit / Migrad as the minimizer.

Fit region	Fitted ${\alpha }_s$	Fitted $C_1$ in GeV	Minuit / Migrad ${\chi }^2/\text{NDF}$
$-0.40\le \cos \chi \le 0.40$	$0.115\pm 0.033$	$3.36\pm 5.56$	$24.6/38$
$-0.45\le \cos \chi \le 0.45$	$0.114\pm 0.023$	$3.41\pm 3.95$	$30.7/42$
$-0.50\le \cos \chi \le 0.50$	$0.115\pm 0.019$	$3.28\pm 3.26$	$43.1/48$
$-0.55\le \cos \chi \le 0.50$	$0.116\pm 0.012$	$3.28\pm 2.14$	$67.6/52$

Equations

$$e^{+}(k_1)e^{-}(k_2)\to {\gamma }^{\ast }/Z^0\to a(p_a)b(p_b)+X.$$

$$\frac{1}{{\sigma }_{\text{tot}}}\frac{d\Sigma (\chi )}{d\cos \chi }=\sum _{a,b}\int \frac{2E_a{E}_b}{Q^2}\delta (\cos {\theta }_{ab}-\cos \chi )d{\sigma }_{a+b+X},$$

$$e^{+}(k_1)e^{-}(k_2)\to {\gamma }^{\ast }/Z^0\to q(p_1)\bar{q}(p_2)g(p_3).$$

$${\mathcal{L}}_{\text{HEFT}}=-\frac{1}{4}\lambda H\text{Tr}\left(G^{\mu \nu }{G}_{\mu \nu }\right).$$

$$\lambda =\frac{{\alpha }_s}{3\pi }\sqrt{\sqrt{2}{G}_F},$$

$$H\to q(p_1)\bar{q}(p_2)g(p_3)$$

$$H\to g(p_1)g(p_2)g(p_3).$$

$$\frac{1}{{\Gamma }_{\text{tot}}}\frac{d{\Sigma }_H(\chi )}{d\cos \chi }=\sum _{a,b}\int \frac{2E_a{E}_b}{m_H^2}\delta (\cos {\theta }_{ab}-\cos \chi )d{\Gamma }_{a+b+X},$$

$${\Gamma }_{\text{tot}}=\frac{{\lambda }^2{m}_H^3}{8\pi }K(\mu )$$

$$K(\mu )=1+\frac{{\alpha }_s}{2\pi }\left[\frac{73}{2}+11\log \frac{{\mu }^2}{m_H^2}-N_f\left(\frac{7}{3}+\frac{2}{3}\log \frac{{\mu }^2}{m_H^2}\right)\right]+\mathcal{O}({\alpha }_s^2),$$

$$E_a{E}_b\delta (\cos {\theta }_{ab}-\cos \chi )$$

$$\cos {\theta }_{ab}=1-\frac{p_a\cdot p_b}{p_a\cdot Q{p}_b\cdot Q},$$

$$\frac{d^{d}p}{(2\pi {)}^d}2\pi {\delta }_{+}(p^2-m^2)\to \frac{1}{i}\frac{d^{d}p}{(2\pi {)}^d}\left(\frac{1}{p^2-m^2-i\epsilon }-\frac{1}{p^2-m^2+i\epsilon }\right)$$

$$2z\equiv 1-\cos \chi$$

$$E_a{E}_b\delta (\cos {\theta }_{ab}-\cos \chi )=(p_a\cdot Q{)}^2(p_b\cdot Q{)}^2\delta \left(2z{p}_a\cdot Q{p}_b\cdot Q-p_a\cdot p_b\right).$$

$$\frac{1}{2z\,p_{a}\cdot Q\,p_{b}\cdot Q-p_{a}\cdot p_{b}}\biggl{|}_{\textrm{cut}}$$

$$H(Q)$$

$$\sum _{\lambda =1}^2{\epsilon }^{\mu }({\mathbf{p}}_i,\lambda ){\epsilon }^{\ast \nu }({\mathbf{p}}_i,\lambda )=-g^{\mu \nu }+\frac{(p_i^{\mu }{n}^{\nu }+p_i^{\nu }{n}^{\mu })}{p_i\cdot n}-\frac{n^2{p}_i^{\mu }{p}_i^{\nu }}{(p_i\cdot n{)}^2},$$

$$\int \left(\prod _k\frac{d^d{p}_k}{(2\pi {)}^{d-1}}{\delta }_{+}(p_k^2)\right)|\mathcal{M}(H\to gg+X){|}^2\sum _{a<b}2E_a{E}_b\delta (\cos {\theta }_{ab}-\cos \chi )$$

$$\left(\prod _k\frac{d^d{p}_k}{(2\pi {)}^{d-1}}{\delta }_{+}(p_k^2)\right)|\mathcal{M}(H\to gg+X){|}^{2}2E_a{E}_b.$$

$$\int$$

$$=6\int \left(\prod _{k=1}^4\frac{d^d{p}_k}{(2\pi {)}^{d-1}}{\delta }_{+}(p_k^2)\right)|\mathcal{M}(H\to gggg){|}^{2}2E_1{E}_2\delta (\cos {\theta }_{12}-\cos \chi ),$$

$$\left(2z{p}_a\cdot Q{p}_b\cdot Q-p_a\cdot p_b\right){\left[\delta ({\mathcal{K}}_{ab}(z))\right]}^j={\left[\delta ({\mathcal{K}}_{ab}(z))\right]}^{j-1},$$

$$\delta ({\mathcal{K}}_{ab}(z))\equiv \frac{1}{2\pi i}\left(\frac{1}{2z{p}_a\cdot Q{p}_b\cdot Q-p_a\cdot p_b-i\epsilon }-\frac{1}{2z{p}_a\cdot Q{p}_b\cdot Q-p_a\cdot p_b+i\epsilon }\right).$$

$$z\to y^2,z\to 1/(1-y^2)$$

$${\hat{z}}_{ab}^n(1-{\hat{z}}_{ab}{)}^m={\int }_0^{1}dz{z}^n(1-z{)}^{m}2p_a\cdot Q{p}_b\cdot Q\delta (2z{p}_a\cdot Q{p}_b\cdot Q-p_a\cdot p_b),$$

$$\int d{\Phi }^{(4)}\mathcal{I}(\{p\})\delta (2z{p}_a\cdot Q{p}_b\cdot Q-p_a\cdot p_b),$$

$$\int d{\Phi }^{(4)}{\hat{z}}_{ab}^n(1-{\hat{z}}_{ab}{)}^m(p_a\cdot Q{p}_b\cdot Q{)}^{m+n}\mathcal{I}(\{p\})$$

Full text

1 Introduction

Many aspirations of the particle physics community for measuring the properties of the recently discovered Higgs-like boson with unprecedented precision and possibly finding new physics beyond the Standard Model are closely linked to the prospects of having a new high-energy $e^{+}e^{-}$ collider in the near future. Precision QCD studies constitute an integral part of the CEPC CEPCStudyGroup:2018rmc ; CEPCStudyGroup:2018ghi , ILC Behnke:2013xla ; Baer:2013cma , FCC-ee Gomez-Ceballos:2013zzn and CLIC Aicheler:2012bya ; deBlas:2018mhx physics programs and it is obvious that our understanding of the strong sector could substantially benefit from new experimental data obtained in the clean environment of a lepton collider.

Over the decades, infrared- and collinear-safe event shape variables turned out to be very useful to confront theoretical calculations in perturbative QCD (pQCD) with the experimental observations. Being functions of the reconstructed 4-momenta of the final state particles, they can be easily extracted from the experimental data. For instance, there are six well-known event shape variables that were studied by the Large Electron Positron (LEP) collider experiments ALEPH Heister:2003aj , DELPHI Abdallah:2004xe , L3 Achard:2004sv and OPAL Abbiendi:2004qz in great details. These are thrust Brandt:1964sa ; Farhi:1977sg , heavy jet mass Clavelli:1981yh , wide and total jet broadening Rakow:1981qn ; Ellis:1986ig ; Catani:1992jc , $C$ parameter Parisi:1978eg ; Donoghue:1979vi and the jet transition variable $Y_{23}$ Catani:1991hj . Although event shape observables are usually discussed in the context of $e^{+}e^{-}$ annihilation, it is also possible to define them in other processes, such as electron-proton Aktas:2005tz , proton-antiproton Aaltonen:2011et or proton-proton Banfi:2010xy collisions.

On the theory side, fixed-order numerical predictions at next-to-next-to-leading order (NNLO) accuracy are available for a wide range of event shape observables Gehrmann-DeRidder:2007nzq ; GehrmannDeRidder:2007hr ; Weinzierl:2009ms ; Ridder:2014wza ; DelDuca:2016ily . Regarding the resummation of endpoint divergences, many results at next-to-next-to-leading logarithmic (NNLL) and even NNNLL accuracy can be found in the literature deFlorian:2004mp ; Becher:2008cf ; Chien:2010kc ; Abbate:2010xh ; Monni:2011gb ; Becher:2012qc ; Banfi:2014sua ; Hoang:2014wka ; Banfi:2016zlc ; Tulipant:2017ybb ; Moult:2018jzp ; Kardos:2018kqj ; Banfi:2018mcq ; Bell:2018gce ; Verbytskyi:2019zhh . Publicly available codes such as Eerad3 Ridder:2014wza or NLOJet++ Nagy:2001fj ; Nagy:2003tz and Event2 Catani:1996jh ; Catani:1996vz implement parton level predictions for event shape variables at NNLO or NLO accuracy respectively, while the effects of parton shower and hadronization can be simulated with dedicated software tools.

One of the tasks of a new $e^{+}e^{-}$ collider will be to measure various event shape variables with even higher precision as compared to what was possible at LEP. However, one should also take into account that the new machine, being a Higgs factory, will have much higher center of mass energy as compared to LEP. This opens up many exciting possibilities to calculate and measure event shape variables that would simultaneously probe our understanding of the strong and the Higgs sectors Gao:2016jcm . One of such observables is the Energy-Energy correlation (EEC) function in gluon-initiated Higgs decays, which is the main subject of this paper.

The original definition of the Energy-Energy Correlation goes back to the late 70s of the last century, when the authors of Basham:1978bw suggested to employ two calorimeters separated by an angle $\chi$ to measure the energies of two hadrons $a$ and $b$ produced in $e^{+}e^{-}$ annihilation

[TABLE]

The EEC is a differential angular distribution that arises from measuring the cosine of the angle $\theta_{ab}$ between all particle pairs $(a,b)$ in the event and weighting each contribution by the particle energies $E_{a}$ and $E_{b}$. The resulting histogram is normalized to unit area and provides a useful way to visualize the energy flow through the calorimeters. The EEC is defined as

[TABLE]

where $d\sigma_{a+b+X}$ is the differential production cross section for the process in eq. (1), $Q$ is the total center of mass energy and $\cos\theta_{ab}=\hat{\bm{p}}_{a}\cdot\hat{\bm{p}}_{b}$. Here $\sigma_{\textrm{tot}}$ denotes the total cross section for $e^{+}e^{-}\to\textrm{hardons}$. The summation is over all $(a,b)$-combinations from the final state hadrons. Notice that we do not include the contributions from self-correlations in our fixed order calculation such that in our results the summation over $a$ and $b$ was replaced by $\sum_{a\neq b}$. The fixed order pQCD predictions for EEC are obtained by calculating eq. (2) at the parton level, with the LO contribution arising from the tree level hard process

[TABLE]

At a future Higgs factory, it appears very natural to measure the EEC in hadronic decays of the Higgs boson. Higgs can decay to a pair of gluons through a top quark loop, or to $b\bar{b}/c\bar{c}$ directly via Yukawa couplings. In this paper we are mainly concerned with the gluon initiated decay of the Higgs boson within the framework of the Higgs Effective Field Theory (HEFT) Wilczek:1977zn ; Shifman:1978zn ; Inami:1982xt ; Kniehl:1995tn . The $b$ or $c$ quark Yukawa initiated decay will be considered in a future work. In the HEFT, the top quark can be integrated out to obtain operators that contain the Higgs field $H$ and two QCD field-strength tensors $G^{\mu\nu}$. The interacting part of the HEFT Lagrangian is then given by

[TABLE]

The corresponding Wilson coefficient $\lambda$ is determined from matching the amplitude for $H\to gg$ in SM to the one in HEFT order by order in the strong coupling $\alpha_{s}$ and the inverse of the top quark mass $1/m_{t}$. To lowest order in $\alpha_{s}$ and $1/m_{t}$ we have

[TABLE]

with $G_{F}$ being the Fermi constant. Currently, $\lambda$ is known at the $\textrm{N}^{4}\textrm{LO}$ Baikov:2016tgj accuracy. The effective operator in eq. (4) gives rise to tree-level coupling of the Higgs with $2$, $3$, and $4$ gluons. In this work, we would like to consider the correlations between energies of the partons that arise from a gluon-initiated decay of the Higgs particle at rest, via

[TABLE]

and

[TABLE]

To avoid possible confusion between the EEC from eq. (2) and the one considered in this paper, in the following we will denote them as “standard EEC” and “Higgs EEC” respectively. Then, in analogy to eq. (2) we can define the Higgs EEC as

[TABLE]

where $m_{H}$ is the Higgs boson mass, $\Gamma_{\textrm{tot}}$ is the total decay width for $H\to gg$ and $d\Gamma_{a+b+X}$ denotes the differential decay rate for Higgs decaying into gluons plus anything else. Notice that here the center-of-mass energy squared in $e^{+}e^{-}$ annihilation $Q^{2}$ was replaced with $m_{H}^{2}$, since we are considering the decay of the Higgs particle at rest, so that $Q=(m_{H},0,0,0)$. The normalization of the Higgs EEC with respect to

[TABLE]

ensures the cancellation of $\lambda$ in the final result for the Higgs EEC. The factor $K$ accounts for the corrections to the total decay width $H\to gg$ within the HEFT. For our purposes it is sufficient to use its NLO value Chetyrkin:1997iv

[TABLE]

where $N_{f}$ is the number of light quark flavors.

The derivation of fixed-order pQCD predictions for event shape variables is usually done using numerical methods. While for some observables such as EEC or thrust the LO predictions can be obtained analytically by directly calculating the 3-particle phase space integrals, this “brute-force” approach is clearly not feasible at NLO and beyond. At NLO the main complication arises from the evaluation of the double real emissions. The corresponding 4-particle phase space integrals are very difficult to calculate analytically due to the complexity of the integrand and its dependence on a (often nontrivial) measurement function. Even though the final result for the given infrared- and collinear-safe event shape variable must be finite, the single integrals will suffer from severe soft and collinear divergences. The cancellation of the divergences occurs only at the very end, when real-virtual and double-real contributions are added together.

The availability of reliable numerical NNLO predictions for event shape observables in $e^{+}e^{-}$ annihilation is undoubtedly useful, but the corresponding analytic results (even at NLO) are still interesting and desirable. Currently, the standard EEC is the only QCD event shape observable known analytically at NLO Dixon:2018qgp . There has also been remarkable progress in computing EEC analytically in ${\cal N}=4$ supersymmetric Yang-Mills theory Belitsky:2013xxa ; Belitsky:2013bja . Analytic results are now available not just at NLO Belitsky:2013ofa , but also at NNLO very recently Henn:2019gkr . The goal of this paper is to obtain the analytic NLO result also for the Higgs EEC using the same methods as in Dixon:2018qgp . To the best of our knowledge, no previous numerical or analytic studies of this observable exist in the literature, even at LO. However, we would like to point out that numerical NLO and approximate NNLO results for thrust in hadronic Higgs decays were recently obtained in Gao:2019mlt .

This work is organized in the following way. In section 2 we describe our framework used to calculate EEC-like observables such as the standard EEC or the Higgs EEC analytically at LO and NLO. In particular, we explain the nontrivial procedure of deriving IBP equations for nonlinear propagators and fixing the boundary conditions in the differential equations for the master integrals. Our analytic results for the Higgs EEC are presented in section 3, while section 4 is devoted to the discussion of the asymptotics in the collinear and back-to-back limits. Section 5 contains a comparison of the NLO result to the predictions of Pythia as a way to estimate the importance of the nonperturbative corrections for our observable. Additionally, we employ the Pythia simulation as a toy model for the determination of the strong coupling, which should emphasize the importance of the Higgs EEC for the future, when real data from a new lepton collider should become available. Finally, we summarize the obtained results and present an outlook for future work in this direction in section 6.

For convenience, in the course of our calculation we set $m_{H}$ to unity. This is allowed, as the dependence of the final result on the Higgs mass can be easily restored by dimensional analysis. In order to avoid cluttering the notation, we prefer to use $m_{H}\equiv 1$ in all relations shown in section 2. Notice that in this case the scalar products $p_{a}\cdot Q$ and $p_{b}\cdot Q$ have mass dimension one. The formulas provided in section 3 and all subsequent sections do not make use of this simplficiation and display the full dependence on $m_{H}$.

2 Calculational setup

2.1 Outline

To approach the task outlined in section 1 we directly calculate the matrix elements squared of all the relevant sub-processes and convert the arising 4-particle phase space integrals into loop integrals using the method of reverse unitarity Anastasiou:2002yz ; Anastasiou:2003yy . Those integrals are then reduced to a smaller set of master integrals via Integration-By-Parts reduction Chetyrkin:1981qh ; Tkachov:1981wb , while for the evaluation of the master integrals we employ the method of differential equations Kotikov:1991pm ; Kotikov:1990kg ; Kotikov:1991hm ; Bern:1993kr ; Remiddi:1997ny ; Gehrmann:1999as . Each system of equations is calculated by first turning it to the canonical form Henn:2013pwa and then fixing the boundary constants.

Among all QCD event shape observables, EEC is arguably the simplest object that can be computed in this framework. Unlike thrust, jet masses, jet broadening and the jet transition variable $Y_{23}$, the measurement function of EEC

[TABLE]

does not require us to minimize or maximize a particular kinematic quantity. Furthermore, the measurement function in the contribution from the particle pair $(a,b)$ depends only on the energies and momenta of these particles, which greatly facilitates the calculation of the individual contributions. However, since

[TABLE]

applying the reverse unitarity Anastasiou:2002yz ; Anastasiou:2003yy

[TABLE]

to eq. (11) yields a propagator that is not linear in the loop momentum dependent scalar products. At this point it is convenient to introduce a new dimensionless variable $z$ so that

[TABLE]

and

[TABLE]

This property of the measurement function is the reason why the analytic NLO calculation of an EEC event shape variable is nontrivial. At first glance, the complications related to the appearance of the nonlinear propagator

[TABLE]

make it very challenging to carry out the calculation using the standard techniques. Yet, as will be explained in the following sections, all these difficulties can be efficiently mitigated with a minimal amount of out of the box thinking.

The main feature of our approach is the direct IBP reduction of loop integrals with nonlinear propagators and the subsequent calculation of the corresponding master integrals. A different route was taken in Gituliar:2017umx , where the nonlinear propagator was converted into a product of two linear propagators at the price of introducing an auxiliary parameter $x$. While this linearization significantly facilitates the IBP reduction of the phase space integrals, the evaluation of the auxiliary master integrals using differential equations seems to become more involved. Nevertheless, it would be very interesting to compare the results in our approach Dixon:2018qgp with those using the method of Ref. Gituliar:2017umx , once they become available.

2.2 Amplitudes

The main ingredient of this calculation are the squared matrix elements $|\mathcal{M}(H\to gg+X)|^{2}$ for the real emission processes

[TABLE]

with 3-parton final states and

[TABLE]

with 4-parton final states. To ensure the correctness of the expressions, we generate the corresponding amplitudes in two different ways, using QGRAF Nogueira:1991ex and FeynArts Hahn:2000kx , cf. figure 1.

The FeynArts model for HEFT was created with FeynRules Alloul:2013bka and the Feynman rules for the Higgs-gluon vertices were additionally rederived using the FeynRule function of FeynCalc Mertig:1990an ; Shtabovenko:2016sxi . The same vertices are also employed when obtaining amplitudes from the output of QGRAF.

Custom FORM Vermaseren:2000nd code is used to evaluate and simplify the squared matrix elements in $d$-dimensions, where the color algebra is handled by the Color vanRitbergen:1998pn package. Summations over the polarizations of the gluons are done using the axial gauge

[TABLE]

where for a particular gluon with the 4-momentum $p_{i}$ we use the 4-momentum of another gluon $p_{j}$ as the auxiliary vector $n$.

Some of the terms in $|\mathcal{M}(H\to gg+X)|^{2}$ contain linearly dependent propagators that require partial fractioning. This is also handled by our FORM code, using rules obtained with the ApartFF function of FeynCalc. ApartFF is based on the multiloop partial fractioning algorithm from Feng:2012iq , which is also implemented in the standalone $Apart package.

2.3 Topology identification

As far as the topology identification is concerned, the presence of the nonlinear propagator eq. (16) prevents us from using the established methods Pak:2011xt that rely on the $UF$-representation of the loop integrals. Instead, we identify the identical topologies by trying different combinations of loop momentum renamings $p_{a}\leftrightarrow p_{b}$ and shifts $p_{a}\to Q-\sum_{b\neq a}p_{b}$. Since the sum of the measurement functions is manifestly invariant under such transformations, the topology identification proceeds by considering the full expression

[TABLE]

and applying all the allowed shifts and renamings to each distinct denominator and its loop momentum dependent coefficient. This allows us to identify all subtopologies contained in

[TABLE]

A subtopology is not the final topology of an integral family, since it lacks the nonlinear propagator. Therefore, each subtopology gives rise to $\left(\begin{subarray}{c}k\\ 2\end{subarray}\right)$ integral families, where $k$ is the number of the final state partons. For example, at LO the single subtopology yields $1\times 3$ integral families.

In the case of the subprocess $H\to gggg$ we slightly modify this procedure to maximally exploit the symmetry of the corresponding matrix element squared under $p_{a}\leftrightarrow p_{b}$. This symmetry allows us to write

[TABLE]

so that in this case the number of the final topologies equals the number of the subtopologies. Of course, in this case the set of possible shifts and replacements becomes more restricted. For example, transformations such as $p_{1}\leftrightarrow p_{3}$ or $p_{2}\leftrightarrow p_{4}$ would introduce additional measurements functions in eq. (22), which we would like to avoid. On the other hand, the renamings $p_{1}\leftrightarrow p_{2}$ or $p_{3}\leftrightarrow p_{4}$ do not alter the measurement function in eq. (22) and are, therefore, allowed. This fully symmetric topology identification gives us only 10 integral families in the subprocess $H\to gggg$ as compared to 18 families when using the other method. Unfortunately, other subprocesses of the Higgs decay do not possess such a high degree of symmetry on the level of the matrix element squared so that the fully symmetric topology identification is not helpful there.

All the manipulations related to the topology identification are handled by an in-house Mathematica code. Owing to the small number of the loop momenta (at most 4) and the fact that all quarks are taken massless, on a modern laptop this procedure requires only few minutes per subprocess. For the sake of completeness, all identified subtopologies at LO and NLO are listed in appendix A.

2.4 IBP reduction

There are many publicly available codes that can automatize the procedure of the IBP reduction, but none of those tools can deal with integrals containing nonlinear propagators such as eq. (16) out-of-the box. Nevertheless, the general method of the IBP reduction is perfectly applicable to those integrals. In this sense, we are facing a technical, rather than a conceptual limitation. To overcome this difficulty we split the reduction into two steps, which can be performed using different tools. In the first step we derive the IBP equations for each topology in the usual way by differentiating with respect to the loop momenta. This also applies to the nonlinear propagator. However, the so obtained system of equations turns out to be incomplete, since we miss a relation that would allow us to lower the integer power of the nonlinear propagator. Therefore, we need to augment the system by adding the relation

[TABLE]

with

[TABLE]

Once a solvable system of IBP equations has been obtained for each topology, we can proceed to the second step, where the reduction is performed by explicitly solving these equations for the relevant loop integrals using the Laporta algorithm Laporta:2001dd . Notice that this step can be performed without any knowledge about the explicit form of the propagators in the original topologies. Apart from the IBP equations and the list of the loop integrals we merely need to specify the positions of the cut propagators.

For the technical realization of the procedure described above, the publicly available tools LiteRed Lee:2012cn and FIRE 5 Smirnov:2014hma were used to perform the steps one and two respectively. Although LiteRed actually does not support bases with nonlinear propagators, with a simple trick it can be nonetheless used to derive the corresponding IBP equations. The basic idea is to first create a new basis that omits the nonlinear propagator, so that the NewBasis command can be evaluated without errors. The nonlinear propagator can be then added by directly modifying the list of the propagators Ds and the scalar product replacement rules Toj. These two modifications are sufficient to successfully run the GenerateIBP command of LiteRed and therefore to obtain the IBP equations for the given topology. The extra relation from eq. (23) can be easily added to the existing set of equations by applying the command Toj to a product of the inverse nonlinear propagator and a j-integral with symbolic indices. Once we have the full set of valid IBP equations, it is a trivial exercise to convert them into the FIRE notation and save the result to a text file that will be later used by our FIRE code.

As far as FIRE is concerned, no additional tricks or modifications are required to use as it is a pure IBP solver. In order to run the C++ version of FIRE on the given set of loop integrals we need to generate the so-called start-files that contain all the relevant information about the current topology. The standard way to generate a start-file is to enter the basis by specifying the propagators, the loop and the external momenta. However, it is also possible to bypass this step and to enter only the corresponding IBP equations. This is done by assigning the set of the IBP equations to the FIRE variable startinglist. We would like to stress that this working mode of FIRE is well documented (cf. section 4.1 of Smirnov:2014hma ) and is supported since very early versions of the program. This feature obviously makes FIRE an ideal tool for our purposes.111The newest version of KIRA Maierhoefer:2017hyi also features the ability to solve custom systems of IBP equations. This functionality became publicly available only after this calculation has already been completed. After having marked the cut propagators via the RESTRICTIONS variable, we can proceed with the commands Prepare and SaveStart that generate the corresponding start-file and save it to the disk respectively. In the following we can run the C++ version of FIRE with the obtained start-file in the usual way, just as one would do it for a standard set of propagators.

2.5 Calculation of the master integrals

With the aid of FIRE we can successfully reduce the loop integrals in each topology to a small set of master integrals. For this particular calculation no further steps are required, as, upon some loop momentum shifts and renamings, all of the resulting master integrals can be related to the masters that were already computed in Dixon:2018qgp . However, since the explicit evaluation of these integrals requires some labor and due to the fact that Dixon:2018qgp contains almost no technical details regarding this step, we would like to provide more explicit explanations in this work.

In section 2.1 we have already voiced our preference for solving the master integrals using the method of differential equations. We can regard the application of this technique as a two step process. First, we need to turn each system of differential equations into a canonical form Henn:2013pwa , provided that such a form exists. The algorithms of Lee Lee:2014ioa or Meyer Meyer:2016slj ; Meyer:2018feh can be used to construct the corresponding transformation matrix iteratively. Once the canonical form has been obtained, we can trivially determine the loop integrals at arbitrary order in $\varepsilon$ up to a set of unknown boundary constants. To determine those constants we need to find suitable boundary conditions, which, depending on the process of interest, may turn out to be very challenging.

To implement the first step of the method, we can again use the tricks from section 2.4 and employ LiteRed to differentiate with respect to $z$ via the command Dinv. The resulting loop integrals are then passed to FIRE. After another IBP reduction we obtain a closed system of differential equations in $z$. The search for the canonical form can be automatized using Fuchsia Gituliar:2017vzm , although other publicly available codes such as Epsilon Prausa:2017ltv or CANONICA Meyer:2017joq are also suitable for our purposes.

Some of our topologies can be converted into the canonical form only after performing a nonrational transformation of $z$, which is not automatically determined by the public Python version of Fuchsia222 It is worth noting that in the meantime the development C++ version of Fuchsia offers a helper tool that can suggest a suitable nonrational transformation, based on Lee:2017oca .. Luckily, in our case the required transformations can be directly inferred from the analytic result for the standard EEC in $\mathcal{N}=4$ SYM Belitsky:2013ofa . The functional dependence of the classical polylogarihtms appearing in the final result on $\sqrt{z}$ readily suggests what kind of transformations are required here. Consequently, we find that using one of the two transformations

[TABLE]

all systems of equations can be straightforwardly converted into the canonical form. The solutions of canonical basis can be easily written down iteratively order by order in $\epsilon$, all in terms of harmonic polylogarithms (HPLs) Remiddi:1999ew with argument $z$ or $y$, which can be conveniently manipulated using the Mathematica package HPL Maitre:2005uu .

Then we can proceed to the second step, i. e. the determination of the boundary constants. The difficulty to find suitable boundary conditions for the EEC observable makes this step the most time consuming part of the whole calculation. This is mainly due to the fact that the fixed order result diverges for $z\to 0$ and $z\to 1$ so that we cannot impose any regularity conditions in these limits.

Therefore, we have to use several different methods to determine the values of the integration constants up to $\mathcal{O}(\varepsilon^{0})$ which is required for the NLO result. First of all, in the collinear limit $z\to 0$ we can predict the power of the leading $1/z$ singularity by applying collinear power counting to each of the loop integrals. Working with light-cone coordinate $p=(p^{0}+p^{3},p^{0}-p^{3},p_{\perp})$, when two of the four partons are collinear ($p_{a}\sim p_{b}\sim(\lambda^{2},1,\lambda)Q$ with $\lambda\ll 1$), the other two can be both anticollinear $p_{c}\sim p_{d}\sim(1,\lambda^{2},\lambda)Q$, both hard $p_{c}\sim p_{d}\sim(1,1,1)Q$, anticollinear and collinear, or anticollinear and hard. Therefore, for each loop integral and each region we can assign a definite scaling in $\lambda$ to $z$, the measure and the occurring scalar products. For practical purposes, it is most convenient to work with the symmetric parametrization of the 4-particle phase space from Gehrmann-DeRidder:2003pne . When the integral expanded around $z=0$ contains stronger poles than predicted by the power counting, the integration constants of those terms must be fixed in such a way, that these contributions vanish. Another requirement we impose is that the master integrals which are pure functions of uniform transcendental weight before converting to canonical basis must vanish in the unphysical limit $z\to\infty$. While the physical values of $z$ lie within the region $0\leq z\leq 1$, we can always perform an analytic continuation to the whole $z$ plane.

The most powerful but also most time-consuming procedure to determine the integration constants involves matching to the inclusive 4-particle phase-space master integrals from Gehrmann-DeRidder:2003pne . Our starting point is the obvious identity

[TABLE]

with $\hat{z}_{ab}=p_{a}\cdot p_{b}\,/(2\,p_{a}\cdot Q\,p_{b}\cdot Q)$, where $m$ and $n$ are some nonnegative integers and the Dirac delta corresponds to the nonlinear propagator from eq. (16) that appears in our master integrals. These integrals can be represented as

[TABLE]

where $d\Phi^{(4)}$ denotes the 4-particle Lorentz-invariant phase space and $\mathcal{I}(\{p\})$ is the part of the integrand that does not depend on $z$. Multiplying the integrand with $z^{n}(1-z)^{m}\,2p_{a}\cdot Q\,p_{b}\cdot Q$ from eq. (25) and with $(p_{a}\cdot Q\,p_{b}\cdot Q)^{m+n}$ and integrating over $z$ from [math] to $1$ we obtain the relation

[TABLE]

where $m$ and $n$ must be chosen such, that the integral over $z$ is convergent. The right-hand side of eq. (27) is evaluated as follows. First of all, we multiply the original integrand from eq. (26) by $2(p_{a}\cdot Q\,p_{b}\cdot Q)^{m+n+1}$ and reduce it to master integrals. All the resulting master integrals are already present in the corresponding system of differential equations, so that we know their solutions up to the unknown boundary constants. Then, we multiply the resulting linear combination of the master integrals by $z^{m}(1-z)^{n}$ and employ HyperInt Panzer:2014caa and HPL Maitre:2005uu to perform the integration over $z$. Now let us turn to the left-hand side of eq. (27), where the loop integrals do not depend on $z$ and contain no nonlinear propagators. Hence, it is trivial to reduce those integrals to the master integrals of the inclusive 4-particle phase space, which are known analytically since long time Gehrmann-DeRidder:2003pne . Equating both sides of eq. (27) provides us with additional relations between the integration constants. Since each evaluation of the right-hand side of eq. (27) requires an IBP reduction of integrals containing the nonlinear propagator, the practical application of this method to fix the boundary constants is very time consuming. Moreover, to avoid spending too much on time on each single reduction, the integer values of $m$ and $n$ should be chosen sufficiently small, but still large enough to cure the divergences in $z\to 0$ and $z\to 1$ on the right-hand side of eq. (27). Luckily, the values with $m,n\leq 1$ are enough for this calculation. We would also like to remark that this method bears similarity with the procedure employed in Gituliar:2017umx . However, in our case it was much easier to apply, as our master integrals depend only on $z$, while the auxiliary master integrals from Gituliar:2017umx are functions of $z$ and the linearization parameter $x$.

Finally, we also demand that after substituting all master integrals into the full NLO real correction, in the limit $z\to 0$ it may not diverge stronger than $1/z$ Konishi:1978yx ; Richards:1982te . When combined together, the above constraints are sufficient to determine all the required boundary constants.

2.6 Real-virtual correction to Higgs EEC

In the previous sections we were mainly concerned with the double-real corrections to the Higgs EEC. While this is undoubtedly the most important and also the most complicated ingredient of our analytic calculation, to obtain the full NLO result we must also include the real-virtual corrections.

The evaluation of these corrections does not pose any difficulties. The matrix element squared for the corresponding loop diagrams is currently known analytically at two loops, for both the dimension $5$ effective operator Gehrmann:2011aa and dimension $7$ effective operators Jin:2018fak . As a cross-check, we have explicitly recomputed the 1-loop contributions to the subprocesses $H\to ggg$ and $H\to q\bar{q}g$. The amplitudes were generated with FeynArts and evaluated using FeynCalc. To calculate the 1-loop integrals analytically we employed the FeynHelpers Shtabovenko:2016whf interface between FeynCalc and Package-X Patel:2015tea ; Patel:2016fam .

To carry out the renormalization, we included the counter terms for the dimension 5 effective operator to our FeynRules HEFT model and generated the corresponding counter term diagrams using FeynArts. Using the ability of FeynCalc and Package-X to explicitly distinguish between UV and IR poles at 1-loop, we have verified that the inclusion of the counter term diagrams removes all UV poles in the real-virtual corrections, leaving us with IR poles only.

The obtained results agree with the expressions arising from the analytic continuation and proper UV renormalization of the real-virtual contributions from Kilgore:2013gba . The final integrals over the massless 3-particle phase space can be carried out directly using the HyperInt Panzer:2014caa package. The still present IR poles in real-virtual cancel against the remaining IR poles in the double-real piece. As expected, the final result for the Higgs EEC is manifestly finite.

3 Full analytic results

The final analytic NLO result for the Higgs EEC can be written as

[TABLE]

with $A_{H}(z)$ and $B_{H}(z)$ being the LO and NLO coefficients respectively, while $\beta_{0}=11/3C_{A}-4/3N_{f}T_{f}$. In QCD we have $C_{A}=N_{c}=3$, $C_{F}=(N_{c}^{2}-1)/(2N_{c})=4/3$ and $T_{f}=1/2$, while $N_{c}$ denotes the number of colors. The overall prefactor $1/K$ arises from the normalization with respect to the total decay width for $H\to gg$ in HEFT, cf. eq. (9). Figure 2 shows the Higgs EEC at LO and NLO for $N_{f}=5$ with the corresponding uncertainties from varying the renormalization scale $\mu$ between $2m_{H}$ and $m_{H}/2$. We set $m_{H}=125.0\textrm{ GeV}$ and use the following values of the strong coupling $\alpha_{s}$ obtained with RunDec Chetyrkin:2000yt ; Herren:2017osy at 4-loop accuracy

[TABLE]

Before presenting the explicit result for the LO and NLO contributions let us, for the sake of clarity, decompose the corresponding coefficients $A_{H}(z)$ and $B_{H}(z)$ into different color pieces that appear in the full result. For $A_{H}(z)$ we will use the subscript “lc” to denote the leading color contribution ($\sim N_{c}$). In the case of $B_{H}(z)$ “lc” stands for the component proportional to $N_{c}^{2}$, “nlc” for the next-to-leading color part ($\sim N_{c}$) and “nnlc” for the next-next-to-leading color piece ($\sim 1/N_{c}$). The components proportional to the number of flavors $N_{f}$ will be named accordingly. For LO we find

[TABLE]

with

[TABLE]

The color decomposition of the NLO coefficient $B_{H}(z)$ reads

[TABLE]

It is interesting to observe that the Higgs EEC exhibits a much richer color structure as compared to the standard EEC. While the LO coefficient of the standard EEC is directly proportional to $C_{F}$, eq. (30) depends not only on $C_{A}$ but also on $N_{f}$. Notice that the $N_{f}$-piece in the standard EEC is an NLO effect, while in the Higgs EEC it appears already at LO. Furthermore, unlike Eq (32) the $B(z)$ coefficient of the standard EEC contains no terms proportional to $N_{f}^{2}$. These differences can be largely attributed to the fact that the Higgs EEC is a gluon-initiated observable, while the standard EEC is a quark-initiated quantity. For example, the splitting of a gluon into a quark-antiquark pair is a LO effect in the Higgs EEC, but for the standard EEC it can occur only at NLO and beyond.

The mathematical structure of the explicit expressions for the color coefficients from eq. (32) is very similar to what has already been observed in the analytic NLO result for the standard EEC. Each term is essentially a product of a rational function of polynomials in $z$ multiplying a building block function $g_{i}^{(j)}$ of pure transcendental weight $j$. These functions are

[TABLE]

All these functions have only $z$ and $\sqrt{z}$ dependence, but intermediate results do have explict dependence on $i\sqrt{z}/\sqrt{1-z}$, which cancels out when combining real-virtual contributions with $q\bar{q}gg$ cut or $gggg$ cut. This is a strong check of the correctness of our final results.

The numerator of the rational function is an univariate polynomial with the highest power of $z$ being 8, while the denominator is a function of type $(1-z)^{m}z^{k}$ with an integer $m$ between 0 and 1, $k$ between 0 and 6. The only exception to this rule arises from the single terms present in $B_{H,\text{lc}}(z)$, $B_{H,\text{nlc}}(z)$ and $B_{H,\text{nnlc}}(z)$ which explicitly depend on $\sqrt{z}$ through $g_{3}^{(2)}$ and the corresponding coefficients. The products of those with $g_{3}^{(2)}$ are symmetric with respect to $\sqrt{z}\to-\sqrt{z}$, a property that has already been observed in the standard EEC result Dixon:2018qgp and the result in $\mathcal{N}=4$ SYM Belitsky:2013ofa .

It is worth noting that $B_{H,N_{f}^{2}}(z)$, which is the simplest piece of $B_{H}(z)$, turns out to be much simpler than the $B_{N_{f}}(z)$ part of the QCD result. While the latter contains one weight 3 function and an explicit dependence on $\sqrt{z}$, the former depends only two weight 1 and three weight 2 functions and is free of square roots.

In the following we list the explicit values of the color components of $B_{H}(z)$.

[TABLE]

To assess the relative importance of the different color components for the final result, we plot $A_{H}(z)$ and $B_{H}(z)$ together with their color coefficients in figure 3 and figure 4 respectively.

4 Asymptotics

Let us now explore the asymptotics of the Higgs EEC by expanding the full NLO result in the back-to-back and collinear limits. For the limit $z\to 0$ expanded up to $\mathcal{O}(z)$ we find

[TABLE]

The other limit, $z\to 1$ corresponds to the back-to-back limit. In this limit, the large logarithms originate from soft and collinear radiations. The resummation of large logarithms in this limit is well understood, and is closely related to various transverse-momentum dependent distribution Collins:1981uk ; Dokshitzer:1999sh ; Moult:2018jzp ; Gao:2019ojf . The expansion up to $\mathcal{O}(1-z)$ yields

[TABLE]

Note that the appearance of $\zeta_{2}\log(2)$ in the constant term at next-to-leading power, which comes from both $B_{H,\text{lc}}$ and $B_{H,\text{nnlc}}$, is different from standard EEC Dixon:2018qgp , where it originates solely from the next-to-leading color part. We note that the leading power terms in the $z\to 1$ limit are in full agreement with the factorization prediction based on the formalism in ref. Moult:2018jzp .333Private communication with Anjie Gao. The subleading power corrections of the Higgs EEC are closely related to the perturbative power corrections for Drell-Yan/Higgs $p_{T}$ distribution at hadron collider. In the latter case, perturbative power corrections in the presence of rapidity divergence Ebert:2018gsn has been studied very recently. We expect the same method can be used to compute the perturbative power corrections for EEC.

5 Estimates of nonperturbative corrections

As every QCD observable, Higgs EEC necessarily contains nonperturbative corrections to the result obtained in pQCD. Even though they are considered to be suppressed as the inverse of the relevant energy scale (in our case $m_{H}$), these contributions are necessary for a meaningful comparison between theoretical predictions and experimental data. In practice, the size of the nonperturbative effects can be estimated using a suitable model (e. g. DWM Dokshitzer:1999sh ) or by simulating parton shower and hadronization with software tools such as Pythia Sjostrand:2014zea . The latter approach is also what we choose to assess the influence of the nonperturbative corrections on the Higgs EEC.

In our Pythia setup we consider the process $e^{+}e^{-}\to H\to gg$ at $\sqrt{s}=m_{H}$ and generate 5000 events, which, according to An:2018dwb roughly corresponds to what CEPC expects to collect in the $H\to gg$ decay channel over a data taking period of 7 years. The tiny size of the Higgs to electrons coupling is irrelevant here, as the Higgs EEC is obtained from the decay of the Higgs particle at rest and does not depend on its production mechanism.

Pythia itself implements only the LO hard matrix element for $H\to gg$ and employs parton shower to approximate higher order corrections. This is sufficient to simulate gluon-initiated $H\to q\bar{q}g$ final states, but turns out to be problematic in the case of $H\to ggg$. The reason for this is that the gluon splitting is not the only contribution to the 3 gluon final state. The other contribution arises directly from the 3-gluon-Higgs interaction vertex and is not negligible. However, the latter is not implemented in the current version of Pythia, so that the simulation of the $H\to ggg$ final state is obviously incomplete.

In principle, one could overcome this issue by matching the current Pythia two jet + parton shower simulation to the full LO matrix element that incorporates both 2-gluon and 3-gluon vertices. However, since we are mainly interested in a qualitative comparison between Pythia and our analytic result, we choose not to go this route and content ourselves with what is already implemented in the code of Pythia.

Despite of these shortcomings, we believe that the comparison of our NLO results to the Pythia predictions is still interesting, especially in view of the lack of any experimental data for the Higgs EEC event shape observable. A more careful and meaningful comparison to a proper numerical simulation that uses the NLO hard matrix element shall be addressed in a future work.

The Higgs EEC distribution is calculated according to Abreu:1996na

[TABLE]

where $\chi_{ab}$ is the angle between the directions of the 3-momenta of the final state particles $a$ and $b$, $\Delta\cos\chi$ is the histogram bin width and $\cos\chi$ is the lower edge of the bin. The overall normalization by the total number of events $N_{\textrm{events}}$ and the bin width ensures that the area under the histogram is unity.

We employ Pythia 8.2 and ROOT 6.14 Brun:1997pa to simulate gluon-initiated Higgs decays both to hadrons and to partons. In the latter case the hadronization effects are disabled via the master switch HadronLevel:Hadronize = 0. Only statistical uncertainties are included. The corresponding plot is shown in figure 5. The discrepancy between the Pythia data and the NLO result can be attributed to the missing $H\to ggg$ vertex in Pythia hard matrix element, sizeable NLO corrections and nonperturbative effects.

For simplicity, we choose to ignore the systematic errors of this simulation, so that the displayed error bars stem only from statistical uncertainties. Since a single event with more than 2 finite state particles generates multiple histogram entries, the errors in different angular bins are statistically correlated. Thus, it does not seem appropriate to calculate the errors assuming that the simulated data obeys an uncorrelated Poisson distribution. Instead, we adopt the approach that was used by the TOPAZ experiment when measuring EEC at the TRISTAN collider Adachi:1989ej . To be more specific, for each Pythia configuration we generate 50 additional samples with the same settings and the same number of events but different random seeds. Then, for each bin $i$ we can calculate the standard deviation $\sigma_{i}$ via

[TABLE]

where $x_{j}(i)$ is the content of the $i^{\textrm{th}}$ bin in the $j^{\textrm{th}}$ sample, $n=50$ and the mean for the $i^{\textrm{th}}$ bin, $\mu_{i}$ is given by

[TABLE]

Comparing the size of the errors in our simulations we observe the unphysical effect that the curve with hadronization seems to have smaller errors than the one without. This suggests that the systematic errors (especially for hadronization) might be much larger than the statistical uncertainties. Looking only at the central values, it is interesting to observe that in the $\cos\chi$ region between $-0.5$ and $0.5$, both curves seem to lie very close to each other. However, given the possibility of underestimated uncertainties, we must abstain from deriving any conclusions on the size of the nonperturbative corrections.

To check whether our observation is a genuine Pythia effect or a possible artifact of an accidental fine-tuning, we varied several important parameters responsible for simulating the hadronization effects according to Skands:2014pea . In particular, we set the effective value of $\alpha_{s}(M_{Z})$ in the final-state radiation to the physical value $0.118$ (TimeShower:alphaSvalue = 0.118), activated the CWM Catani:1990rr resummation of subleading terms using the 2-loop running of $\alpha_{s}$ (TimeShower:alphaSuseCMW = on, TimeShower:alphaSorder = 2) and varied the value of the IR shower cutoff (TimeShower:pTmin) by setting it to $1.0\textrm{ GeV}$ and $2.0\textrm{ GeV}$. As can be inferred from figure 6, the effects of these changes in the parameter sets are rather mild and they seem to support the initial observation that, judging from the central values only, Pythia possibly indicates comparably small influence of the hadronization for $|\cos\chi|\leq 0.5$. Those effects, however, become more prominent in the collinear and back-to-back regions, $z\to 0$ and $z\to 1$ respectively. Of course, these naive comparisons should not be interpreted as a definitive statement on the size of the hadronization effects in the Higgs EEC. To make such a statement one would need to perform a much more careful numerical simulation and error analysis. Here we merely describe our observations and speculate about their possible interpretation.

Finally, it is tempting to attempt a fit of the analytic NLO result to the Pythia prediction in order to perform a toy determination of $\alpha_{s}$. Our motivation is clear: It is well known that the standard EEC measurements at LEP and other experiments were often used to determine the value of the strong coupling. Therefore, should a future lepton collider measure the $H\to gg$ channel with sufficient precision, it will not take long until experimentalists and theorists will attempt to extract the value of $\alpha_{s}$ from the Higgs EEC measurements. The purpose of our toy fit is to attract the attention of the high energy physics community to this scenario and to ensure that when the real data becomes available, all other ingredients for a rigorous extraction of $\alpha_{s}$ will be already there.

For simplicity, we model the nonperturbative correction to the analytic NLO result by employing the “simple power correction” ansatz from Abdallah:2003xz which amounts to fitting

[TABLE]

to the Pythia data, where $\Sigma_{H,\textrm{pert}}(\chi,\alpha_{s})$ corresponds to eq. (28) with $N_{f}=5,N_{c}=3$ and $\alpha_{s}$ being a fit parameter. Another fit parameter is $C_{1}$ which is assumed to be a constant that accounts for the nonperturbative corrections. Regarding the fit region, it is important to ensure that the values of $\cos\chi$ are sufficiently far away from the collinear and back-to-back regions, where the pure fixed-order result ceases to be valid.

The results of our binned maximum likelihood fit to the 5000 events simulated with Pythia are summarized in table 1. Despite this being a toy fit, the fitted values of $\alpha_{s}$ for different $\cos\chi$-regions are consistent with $\alpha_{s}(m_{H})=0.113$ within the (admittedly very large) error bounds. The fitted curve for $|\cos\chi|\leq 0.5$ is shown in figure 7.

6 Summary

In this paper we have initiated the study of the Energy-Energy Correlation in gluon-initiated Higgs decays, that can simultaneously probe our understanding of the strong and Higgs sectors. Using the techniques developed in our previous work on the standard EEC in $e^{+}e^{-}$ annihilation Dixon:2018qgp and working in the framework of the Higgs EFT, we obtained the fully analytic result for Higgs EEC at NLO in fixed-order perturbation theory. This result bears strong similarities to the analytic NLO result for the standard EEC, in the sense that both observables are described by the same set of master integrals and can be written using the same basis of building block functions, made of classical polylogarithms up to weight 3.

Regarding the phenomenological relevance of the new observable, large QCD background accompanying Higgs decays at hadron colliders makes it very unlikely that the decay channel $H\to gg$ can be measured by the LHC experiments in the near future (if at all). However, such a measurement will be certainly possible in the clean environment of a future lepton collider, such as ILC, CEPC, FCC-ee or CLIC. Given sufficiently high statistics, the Higgs EEC could be also used for the determination of the strong coupling constant, as it was the case in the studies of the standard EEC conducted by the LEP experiments.

In order to attract more interest to the Higgs EEC from the experimental and theoretical side, we provided a very rough estimate the size of the hadronization effects using Pythia 8. The simulated data hints that the nonperturbative corrections might be comparably small in the central region. Since the simulation is done using only the LO hard matrix element for $H\to gg$ implemented in Pythia and suffers from the absence of the $H\to ggg$ vertex and the error bars in our plots stem only from statistical uncertainties these results should be obviously taken with a grain of salt. Nonetheless, by comparing our analytic result to the Pythia prediction and using a naive model for nonperturbative effects, we also performed a toy determination of $\alpha_{s}$, obtaining a value that is (within the very large error bounds) compatible with the current PDG world average. The purpose of this exercise is to motivate the high energy physics community to explore the phenomenology of the Higgs EEC and other related event shape observables in greater details.

For the future, it would be very useful to obtain the NNLO correction to the Higgs EEC, at least numerically. This is important for the phenomenological studies, especially the determination of $\alpha_{s}$. A fully analytic NNLO result is also very desirable, but as of now the complexity of the corresponding expressions poses a major obstacle to this endeavor. Although EEC currently appears to be the most convenient event shape variable for fixed-order analytic calculations, one may still wonder about the feasibility of analytic NLO results for other related observables such as thrust or the $C$ parameter. Once such results are obtained for parton production in $e^{+}e^{-}$ annihilation, it would be a simple exercise to obtain the corresponding predictions also for the Higgs decays.

Acknowledgements.

We thank Peter Skands for his detailed explanation of the Pythia 8 options related to the hadronization and parton showers. We thank Anjie Gao for checking the results in the back-to-back limit, and Han-tian Zhang for useful discussions. The work of M.X.L., V.S., T.Z.Y., and H.X.Z. was supported in part by the National Science Foundation of China (11135006, 11275168, 11422544, 11375151, 11535002) and the Zhejiang University Fundamental Research Funds for the Central Universities (2017QNA3007).

Appendix A Subtopologies

A.1 LO

Single subtopology for $H\to ggg$

[TABLE]

Single subtopology for $H\to q\bar{q}g$

[TABLE]

A.2 NLO

Ten (sub)topologies for $H\to gggg$

[TABLE]

Two subtopologies for $H\to q\bar{q}gg$

[TABLE]

Single subtopology for $H\to q\bar{q}q\bar{q}$

[TABLE]

Single subtopology for $H\to q\bar{q}q^{\prime}\bar{q}^{\prime}$

[TABLE]

Appendix B Asymptotics of the color components

B.1 Collinear limit

[TABLE]

B.2 Back-to-back limit

[TABLE]

Appendix C Identical-quark interference contributions

In this appendix, we provide the identical-quark interference terms that correspond to the $q\bar{q}q\bar{q}$ cut diagram from figure 1 $(\text{c})$, which contributes to $B_{H,\textrm{nnlc}}$. In addition, $B_{H,\textrm{nnlc}}$ also receives contributions from the $q\bar{q}gg$ cut diagram from figure 1 $(\text{b})$. We can write the real contributions to $B_{H,\textrm{nnlc}}$ as

[TABLE]

Similarly with $B_{\textrm{qq}_{\textrm{int}}}$ from the standard EEC (Dixon:2018qgp, ), there are no virtual corrections to $B_{H,\textrm{qq}_{\textrm{int}}}$ at NLO. Thus, this contribution is separately IR finite and gauge invariant and can be written as

[TABLE]

where we introduced two additional building block functions of pure transcendental weight 3, in addition to those already presented in eq. (33)

[TABLE]

It is interesting to observe that these two functions always appear together as $(g_{6}^{(3)}-2g_{7}^{(3)})$, which is also true for $B_{\textrm{qq}_{\textrm{int}}}$ of the standard EEC (Dixon:2018qgp, ). However, $B_{H,\textrm{qq}_{\textrm{int}}}$ is much simpler than $B_{\textrm{qq}_{\textrm{int}}}$. While the latter contains 11 building block functions, including those with an explicit dependence on $\sqrt{z}$ and $i\sqrt{z}/\sqrt{1-z}$, the former can be expressed in terms of only 8 functions, all of which depend solely on $z$.

Expanding the result for $B_{H,\textrm{qq}_{\textrm{int}}}$ in the collinear limit we get

[TABLE]

while the expansion the back-to-back limit yields

[TABLE]

We would like to stress the fact that both limits are free of leading-power singularities $1/z$ and $1/(1-z)$ respectively. This is mainly because quark-antiquark pairs can only arise in the splitting of gluons originating directly from the decay of the Higgs and those gluons cannot be soft.