Analytic Computation of Three-point Energy Correlator in QCD

The energy correlator measures the energy deposited in multiple detectors as a function of the angles among them. In this paper, an analytic formula is given for the three-point energy correlator with full angle dependence at leading order in electron-positron annihilation. This is the first analytic computation of trijet event shape observables in QCD, which provides valuable data for phenomenological studies. The result is computed with direct integration, where appropriate parameterizations of both phase space and kinematic space are adopted to simplify the calculation. With full shape dependence, our result provides the expansions in various kinematic regions such as equilateral, triple collinear and squeezed limits, which benefit studies on both factorization and large logarithm resummation.


Introduction
One of the event shape observables that attracts lots of recent interest in quantum chromodynamics (QCD) and the collider physics community is energy correlators. Traditionally, energy-energy correlation (EEC) measures the energy deposited in two detectors as a function of the angle between these two detectors [1,2]. As observed in [1], the fact that energy weights suppress the soft divergence makes EEC less sensitive to soft gluon emissions. More recently, EEC is generalized to a broader class of observables called energy correlators. In particular, the three-point energy correlator (EEEC), which depends on the three angles among the detectors, contains the nontrivial shape information of the scattering process [3][4][5]. In perturbative theories, EEEC is defined as where i, j and k run over all final-state particles, Q is the total energy of the electronpositron annihilation, and dσ is the differential cross section. For convenience, we normalize the distribution to the born cross section. EEEC is infrared finite in the tree-level γ * → 4 jets process, which allows us to perform the calculation in d = 4 dimension.  Figure 1: (a) A graph on the three-point energy correlator. The three detectors are separated by finite angles θ 12 , θ 13 and θ 23 , capturing outgoing particles at specific angles from the hard interaction and summing their energies. (b) The "zongzi"-shaped kinematic space {x 1 , x 2 , x 3 }, which is constrained by the four-particle phase space.
Energy correlators are almost the simplest infrared (IR) safe jet observables to compute analytically. The leading order (LO) EEC in QCD is obtained since 1970s [1,2]. Recently, EEC is also computed analytically to next-to-leading order (NLO) in QCD [6][7][8] and NNLO in N = 4 super Yang-Mills (SYM) theory [9,10]. At the same time, the collinear limit of the LO EEEC in both N = 4 SYM and QCD is studied in [3] and the complete LO N = 4 SYM result becomes available very recently [5]. In this paper, we calculate the complete LO EEEC in QCD, which shares a similar function space and analytic structure as in N = 4 SYM.
There is also lots of progress in studying energy correlators with effective field theories (EFTs), such as Soft-Collinear Effective theory (SCET) [11][12][13][14][15], which proves to be essential in jet substructure. As summarized in [6], EEC is both singular in the collinear and backto-back limits, and large logarithms in both limits could possibly spoil the perturbation theory. Regarding the collinear region, the resummation has been achieved to the nextto-next-to-leading logarithm (NNLL) accuracy in QCD [16] and N = 4 SYM [17]. In the back-to-back limit, EEC is resummed to NNLL accuracy and matched to NNLO fixed-order prediction [18][19][20][21][22], while a new factorization formula is also introduced in [23], allowing the resummation to N 3 LL [24]. With the recently derived four loop rapidity anomalous dimension in QCD [25,26], EEC is also resummed to N 4 LL accuracy in the back-to-back limit [26]. EEC can also be studied at a hadron collider, the simplicity of the soft function allows the NNLL resummation in the back-to-back limit [27]. More interestingly, the collinear factorization can also be generalized to EEEC observable, where the distribution is factorized into the convolution of a hard function and a jet function. While the factorization is straightforward in SCET, the resummation becomes subtle due to multiple variables. One way is to project the full kinematic region into a one-dimension space, which is referred to as the projected energy correlators [4]. The projected N -point correlator is defined as and its collinear logarithms can be resummed to NNLL accuracy [28]. It would be also interesting to study EEEC in other kinematic limits. Since the shape dependence of EEEC provides more information on the jet substructure, several singular regions besides collinear remain unexplored: equilateral limit (x 1,2,3 ∼ η), squeezed limit (x 1 ∼ 0, x 2,3 ∼ x), coplanar limit and so on. Our fixed-order calculation allows one to extract both leading power (LP) and next-to-leading power (NLP) expansions, which benefit the large logarithm resummations. In a word, the energy correlator is a bridge to precision standard model tests and new physics searches. The energy correlators attract lots of attention on the phenomenological side these days. In Ref. [29], both the shape dependence and the scaling behavior of EEEC, as well as the ratio of projected energy correlators with respect to EEC are measured with the CMS open data. The close agreement between theoretical prediction and CMS open data proves that energy correlators will play an important role in precision QCD measurement and jet substructure, and it would be interesting to perform the measurement at the Large Hadron Collider (LHC). Besides, energy correlators enable measurements in hadronic environments to be theoretically predicted by means of modern loop computation techniques and track functions [4,30,31]. Traditionally, the calculation of track-based observables (e.g. angularities) requires the full functional form of track functions T (x) [32], of which the renormalization group evolution is described by complicated nonlinear equations. However, it is found recently that energy correlator is advantageous for studying track information since it only needs a finite number of track functions moments, which are just numbers and hence do not take part in the phase space integration. It is also suggested that an energy correlator can be applied to top quark mass measurement at the LHC [33].
It has been observed that N -point energy correlators can be written as (N + 2)-point Wightman correlation function of energy flux operators and source operators that produces the localized excitation [34]. Explicitly, EEEC can be alternatively defined by where the energy flux operator is given by integrated stress-energy tensor T µν along the direction n i [35][36][37][38]: and for electron-positron collision, the source operator O is the electromagnetic current. In conformal field theory (CFT), the light-ray operator product expansion (OPE) [39,40] of the energy flux operators reveals that the collinear behavior of EEC is determined by the spin-3 non-local operators [34]. Recently, the squeezed limit of EEEC has been investigated, where the light-ray OPE is developed at leading twist in QCD, in order to understand the transverse spin structure in the squeezed limit [41]. In fact, this spin structure gives rise to a quantum interference at colliders: when rotating the squeezed detector by an angle φ with respect to the third detector, the interference between the intermediate virtual gluon with different helicity leads to a cos(2φ) dependence [42]. Furthermore, standard CFT tools like conformal blocks and Lorentz inversion formula [43,44] are also developed to organize the power correction of triple-collinear EEEC [45,46], opening a new window to studying jet substructure. More recent progress can be found in [47,48]. An outline of this paper is as follows. In Section 2, we introduce the calculation method for three-point energy correlators at leading order. Briefly speaking, we directly integrate the tree-level matrix elements over the four-particle phase space and express the result in terms of transcendental polylogarithmic functions. With our parameterization, the nonanalytic structure in the phase space factorizes and EEEC is reduced to a two-fold integral that can be calculated directly. We discuss the structure of the analytic expression and the numerical checks in Section 3. In Section 4, we extract the equilateral limit, the triple collinear limit and the squeezed limit contributions. The analytic formula for equilateral EEEC and its endpoint behaviors is given for all partonic channels. For the triple collinear limit, we also present a method that allows us to directly extract the subleading power corrections from expanding the EEEC integrand. We summarize in Section 5.

Calculation setup
The leading order EEEC arises from the tree-level process γ * → 4 partons. Given the appearance of the non-standard measurement function in Eq. (1.1), it is not easy to directly apply the modern loop techniques like Integration-by-parts (IBP) [49] and differential equations [50,51]. While for the cases that only involving one non-standard cut propagator like δ(x 1 − (1 − cos θ ij )/2), a method was proposed in Refs. [6][7][8] to allow for a generalized IBP reduction in LiteRed [52,53] and Fire [54,55]. The appearance of three non-standard cut propagators in Eq. (1.1) makes the application of the method in Refs. [6][7][8] much less efficient. Instead of trying to improve the efficiency of the same method, we take the EEEC definition Eq. (1.1) and calculate the phase space integral directly. The main feature of our method is appropriate parameterizations of the four-particle phase space dPS 4 and the kinematic space {x 1 , x 2 , x 3 }, which makes the direct integration possible. Since we only care about EEEC at LO, it is safe to perform the computation in the d = 4 dimension.

Amplitudes and topology identification
We start by calculating the matrix elements squared |M| 2 for γ * → 4 partons with QGRAF [56] and FORM [57], where the color algebra is handled by the Color package [58]. The calculation includes three subprocesses: where q andq stand for non-identical quarks compared with the quarks q andq. In Fig. 2, we present some typical diagrams for the matrix elements. We also compute the same matrix elements squared in FeynArts [59] and FeynCalc [60,61] as a crosscheck. For both of them, we adopt the axial gauge when summing the gluon polarizations where for a particular parton with the momentum p i , the momentum of another parton p j is used as the auxiliary vectorn. By Lorentz invariance, the matrix elements squared are expressed in terms of the standard Mandelstam variables s ij = (p i + p j ) 2 . Our topology identification is a bit different from standard QCD calculations, where the established methods require the UF-representation of the Feynman integrals [62]. In our calculation, it is enough to permute the final state momenta p 1,2,3,4 or equivalently, permute s ij , and classify terms that are invariant under such transformations. Importantly, we have to carry the energy weights E i E j E k together since they are not invariant under particle renaming. Since we are not going to use the topologies for IBP reduction, the point of topology identification is to reduce the integrand and simplify the phase space integration. After obtaining the reduced matrix elements |M(p 1 , p 2 , p 3 , p 4 )| 2 , we rename Figure 2: Some typical graphs on the matrix elements squared |M| 2 for γ * → 4 partons. The first graph corresponds to the double gluon emissions, while in the second graph, the gray lines represent the non-identical quark pair. The last two graphs show the interface between identical quark pairs. the particles such that the energy weights all become E 1 E 2 E 3 : where Π ijk is a short-hand notation for the measurement function In the second line of Eq. (2.3), we make the summation on 4 explicit and rename the momentum labels such that there is no energy weight E 4 in the expression. Then we only need to calculate the unsymmetrical part in the last line since the result is fully symmetric in x 1 , x 2 , x 3 (x 1,2,3 ).

Phase space parameterization
The most challenging part of calculating EEEC is the calculation of the phase space integral. Recall that the massless four-particle phase space measure in d dimension [63] is given by The main difficulty of the direct integration method then comes from the non-trivial constraint Θ(−∆ 4 ), which corresponds to a complicated region of the four-particle phase space. To resolve this problem, we first introduce the energy fractions of three final state particles in the center of mass frame of γ * , Notice that q = p 1 + p 2 + p 3 + p 4 , the above equation becomes where and in the following we set Q 2 = 1. Although the energy fractions break the symmetry of renaming final state particles, the complicated constraint Θ(−∆ 4 ) decouples from the integrals. For example, where ∆ 4 is factorized into integration variables dependent and non-dependent parts Here ∆ 4 ≤ 0 becomes the constraint for the kinematic space {x 1 , x 2 , x 3 }. Fig. 1b shows the allowed kinematic regions, and as we will see in Section 4, the shape dependence of EEEC is encoded in different limits of this region. Note that in the triple collinear limit, ∆ 4 is further reduced to where √ x 1 , √ x 2 and √ x 3 can be interpreted as the lengths of three sides for a triangle due to Helen's area formula.
In summary, EEEC is simplified to an integral over the energy fraction z 1 , z 2 and z 3 . While z 3 is integrated by the δ function in Eq. (2.11), the remaining two-fold integral can be finished using a package called HyperInt [64].

Direct integration
Without the constraint from the δ function in Eq. (2.11), the ranges for integration variables z 1 and z 2 are both from 0 to 1. With the constraint, the integration regions become nontrivial. Explicitly, the result is found to be where we use f (x 1 , x 2 , x 3 , z 1 , z 2 , z 3 ) to represent the EEEC integrand. In our calculation, the integration of z 2 in Eq. (2.14) can be easily carried out with standard mathematical tools like Mathematica or Maple. From the result, we identify the following two possible square roots that will appear in the final result of the LO EEEC, To perform the computation of the remaining one-fold integral with respect to z 1 , we need to rationalize the square roots in Eq. (2.15) by parameterizing the kinematic space {x 1 , x 2 , x 3 }. Explicitly, we introduce a complex variable z and its congugatez as well as a purely imaginary variable t via , (2.16) such that the two square roots are rationalized Note that as observed in Eq. (2.13) and in Ref. [3], the second square root disappears in the triple collinear limit and we no longer need t variable. While in the triple collinear limit, z turns out to be a nice variable that characterizes the triangle shape dependence of EEEC and manifests the S 3 × Z 2 symmetry, {z, t} are not good variables for the full shape dependence and for phenomenological studies eventually. So we will change back to the angular distances x 1,2,3 after finishing the calculation. Using the z, t parameterization, we partial fraction the integrand and format the denominators to be linear functions in the last integration variable z 1 . Subsequently we can evaluate the final integration in HyperInt [64]. It gives us the result in terms of Goncharov polylogarithms (GPLs) [65][66][67], up to transcendentality-two. The GPL is defined iteratively by G(a 1 , · · · a n ; x) ≡ x 0 dt t − a 1 G(a 2 , · · · a n ; t) , There are lots of analytic calculations at two-loop order found involving GPLs, both loop integrals and phase space integrals. It is conjectured that GPLs up to transcendentalitythree can be expressed in term of logarithms and classical polylogarithms Li n (x) with n ≤ 3 [68]. For transcendentality-four GPLs, one also need the special function Li 2,2 (x, y).
For EEEC at LO, we only need GPLs up to transcendentality-two. The conversion from low transcendental weight GPLs to polylogarithms can be done with public packages like PolyLogTools [69] or gtolrules.m [70]. We use the latter package to achieve the conversion. The same results can be obtained by the direct integration from the definition in Eq. (2.19), and we also modify the arguments to meet Mathematica's branch prescription for polylogarithms. After the conversion, our results are expressed in terms of classical polylogarithms.
To simplify the expression, we first collect the transcendental functions with the same rational coefficients. This constructs a raw transcendental function space in terms of classical polylogarithms. All the rational functions are simplified by the MultivariateApart package [71], which implements the partial fraction algorithms for multiple variables. However, simplifying the raw transcendental function space is in general not easy given the three variables. We start by applying transcendentality-two identities to simplify the individual base. A typical set of dilogarithm identities is as follows: which all comes from the well-known five-term identity [72]: It turns out useful to use the five-term identity to simplify complicated arguments. Then we add back all permutation terms of x 1,2,3 (what we call symmetrization) as specified in Eq. (2.3), and reorganize the result such that all bases and the corresponding coefficients are real. A better transcendental basis was already presented in Ref. [5] for the N = 4 SYM EEEC at LO. So for the last step, we try to project our function basis to the basis in Ref. [5]. Explicitly, we symmetrize the N = 4 function basis and construct a new linear independent transcendental weight-two basis. By evaluating both basis at a same numerical point and applying PSLQ algorithm [73,74], we managed to find the linear relations between the elements of these two function bases and successfully simplify our full result in QCD.
As a crosscheck, we evaluate the original result from HyperInt numerically using public GPL libraries like GiNac [75] and FastGPL [76], and compare with the predictions of our final analytic expression. It is interesting to ask how we can simplify the expression in the first step. On the one hand, this requires one to know the singularities of the result and to rule out all spurious poles and branch cuts. Landau equation [77] or Polynomial reduction [78] provides a sufficient set of possible singularities, but it is challenging to figure out the minimal set, especially in the high-dimensional complex hyperplane. Some of the progress can be found in Refs. [79,80]. On the other hand, given the singularities of the integrals, there are still ambiguities how to choose the arguments of polylogarithms. To our knowledge, there is no public algorithm to search for the best arguments that make the expression shortest. It is possible that this can be done with symbol [81] or even with Machine Learning [82] in the future, but it is out of scope of this paper.

Results
In this section, we present the full result for three-point energy correlator in QCD, in terms of the angular distance variables x 1,2,3 and the short-hand notations for the square roots: At LO, there are three color channels where we normalize the distribution to the born-level cross section, and α s = g 2 s /(4π) with g s being the strong coupling constant. All channels contain functions up to transcendentalitytwo and take the form The EEEC function space is composed of 7 logarithmic bases f 1,··· ,7 and 21 polylogarithmic bases g 1,··· ,21 . The transcendental weight-one bases are with the explicit S 3 permutation symmetry. The transcendental weight-two bases are Re log Re log Re log .
Although we use the short-hand notation Re and Im to make the expressions compact, all the bases are analytic functions themselves. The corresponding coefficients R (1) i and R (2) i are rational functions in terms of x 1,2,3 and s 1 , s 2 . We provide these coefficients in the ancillary file. We emphasize that EEEC encodes both scaling information and non-trival shape dependence since it is a three-parameter jet observable. Unlike collinear EEEC, the longest angular distance x L = max{x 1 , x 2 , x 3 } does not factorize out. Instead, the kinematic space is fully determined by the second square root s 2 2 ≤ 0, i.e.
To verify our analytic result, we consider two special cases {x 1 , x 2 , x 3 } = {3y, 2y, y} as well as {x 1 , x 2 , x 3 } = { 11 4 y, 5 3 y, y} and calculate them in Event2 [83,84] and NLOJet++ [85]. The obtained results are in good agreement with our analytic results (see Fig. 3). We also pick the first configuration and separate different color structures as well as the identical quark pair contribution in Event2. The detailed comparison can be found in Fig. 4.  .7), the kinematic spaces are cutoff at y = 1 3 and y = 959 2640 ≈ 0.36 respectively. Fifty billion points are sampled and the internal cutoff is set to 10 −14 in Event2. Ten billion events are generated in NLOJet++. Our result can be useful for phenomenological studies in precision QCD and jet physics. With a simple function basis as well as the simplified rational coefficients, evaluating its numerical values to high precisions is much faster than the raw GPL expression or a Monte Carlo program. As an example, it is easy to numerically evaluate our analytic expression in Mathematica to 200 digits precision within 4 seconds for a regular point in a single core machine. The simplicity of the result strongly encourages us to compute EEEC in QCD for gluon-initiated or bb-initiated Higgs decays analytically in the future.

Kinematic analysis
Given the complete shape dependence of the three-point energy correlator, it is interesting to investigate its behavior under different kinematic limits. Fig. 5 shows several typical regions that could be useful for understanding the singularities and resummation. They are • Triple collinear limit:

and its permutations
• Back-to-back limit: x 1 ∼ 1 and its permutations • Coplanar limit: s 2 → 0 Alternatively, one can borrow the variables {s, τ 1 , τ 2 } from [5], which is related to the angular distance via Here we put the three points in a circle with radius √ s on the celestial sphere and τ 1,2 corresponds to the angle between two of them (see Fig. 6). One can also extract the kinematic limits using the new coordinate, and particularly, it is more convenient to expand the coplanar limit via s → 1.
Back-to-Back x 1 1 Figure 5: Various kinematic limits in the {x 1 , x 2 , x 3 } "zongzi"-shaped space. We denote the triple-collinear limit, squeezed limits and back-to-back limits using different colors. The coplanar limit corresponds to the boundary of the kinematic space itself. The full 3D dynamic figure can be found in the ancillary file.
From phenomenological perspective, one can also slice the kinematic space and apply a constraint on the angular distance x 1,2,3 . Since the demonstration of the full EEEC requires a 3D density plot, from which it is difficult to read information, applying the kinematic constraints helps to reduce the dimension of the plots. The most simplest case is the equilateral EEEC, where three angular distances are the same Two other typical choices are isosceles configuration x 1 = x 2 and the right configuration Figure 6: (a) A graph on the {s, τ 1 , τ 2 } coordinate on the celestial sphere [5]. (b) The kinematic limits of EEEC under the {s, τ 1 , τ 2 } coordinate. s → 0 and s → 1 lead to the triple collinear and coplanar limit respectively, while τ i → 1 represents squeezed limit.
x 3 +x 2 = x 1 . When dealing with data from experiments or simulation programs like Pythia [86][87][88], it is straightforward to apply these constraints directly in the event selection. There are also overlaps between the expansion of kinematic limits and the configuration constraints. For example, one can study both the triple collinear limit and coplanar limit in the equilateral configuration. Applying both expansion and slicing together gives a clearer picture on specific jet substructure that we want to understand.
In this section, we will focus on the analytic result of equilateral EEEC and triple collinear EEEC at next-to-leading power, as well as the squeezed limit. We present the full expression for the equilateral EEEC, with equilateral function space included. To give a concrete example of applying both the kinematic expansion and configuration constraint together, we discuss the x → 0 (collinear limit) and x → 3 4 (coplanar limit) singular behaviors, which could be interesting to the studies of trijet events at colliders. Regarding the triple collinear limit, we present the NLP correction analytically. It turns out that the collinear function space only contains one transcendental weight-one function and three weight-two functions under the S 3 symmetry. Finally, we extract the LP squeezed limit and discuss the ambiguity of the definition under the triple collinear limit. The geometry reveals that the squeezed limit is actually path-dependent.

Equilateral limit
It is straightforward to extract the equilateral EEEC from our analytic formula. Alternatively, one can take the equilateral limit before performing the integral. Explicitly, Eq. (2.11) becomes which can be evaluated directly. The Heaviside function suggests the equilateral EEEC is cutoff at x = 3 4 . The analytic result is written as where the normalization factor is N = 1 4πx √ 3−4x , and 1 4 and 1 2 are symmetry factors due to identical particles. In Eq. (4.3), we also introduce another variable t = √ 3 − 4x to make the result more compact. The n f contribution is given by where the transcendentality-two part T 2 (t) is 4 , (4.5) with the corresponding function space 3 = 2(tan −1 t) 2 + ζ 2 , g (4.6) As a Single-valued Harmonic Polylogarithm (SVHPL), Bloch-Wigner function satisfies The results for the other two partonic channels are given as follows, 6 , (4.10) 5 , (4.12) 6 , (4.13) where we need five more function bases 8 = π log t 2 + 1 , g 9 = ζ 2 . (4.14) In Fig. 7, we also show the equilateral EEEC result from Event2, which has good agreement with our analytic formula. Even slicing the kinematic space with equilateral constraint gives us an interesting result. There are two singular limits: the collinear limit x → 0 and coplanar limit x → 3 4 . The collinear expansion is given by where the coefficients are linear combinations of Clausen functions and Riemann zeta functions: Here κ ≡ Cl 2 π 3 = Im Li 2 e iπ 3 is Gieseking's constant, with Cl 2 (φ) = − φ 0 log |2 sin x 2 |dx being the Clausen function. This is a transcendentality-two number that is typical in the trijet computation (e.g., the one-loop trijet soft function [89]). Another interesting feature is from the n f color factor, i.e., the coefficient of n f for F 4 in Eq. (4.16) doesn't involve κ while κ still shows up for F 1 , F 2 , F 3 . We will find a similar feature in the triple collinear limit in the next subsection.
Due to the kinematic cut Θ 3 4 − x , there is no back-to-back limit in the equilateral configuration. Instead, the three detectors are separated by an angle 2π 3 on the same plane, which refers to the coplanar limit. The coplanar expansion includes both fractional power divergence and logarithmic divergence: (4.17) and the corresponding coefficients are where we need one more transcendentality-two number Li 2 (−3). The C A + 2C F structure in R 1 implies that the leading logarithm in the cumulant can possibly be predicted by a Sudakov form factor [90]. It would be interesting to study the singularity structure in both limits and resum the large logarithms in the future. In the collinear limit, our result provides the regular terms that complete the two-loop equilateral EEEC jet function. To recover its close form in , one might need to compute equilateral EEEC to higher orders in expansion. The soft gluon enhancement appears in the coplanar limit, and similarly, our fixed-order calculation provides the needed ingredients for its resummation. Some of the similar analysis for another trijet event shape observable D-parameter can be found in Refs. [91][92][93]. Either way, equilateral EEEC contains valuable information on understanding the symmetric trijet events in electron-positron collisions.
Like two-point energy correlator (EEC), equilateral EEEC has nice analytic properties and is free of Sudakov shoulders [94]. For event shape observables like thrust, C parameter and heavy jet mass, the range of the parameter grows order by order in perturbation theory, and the incomplete cancellation between real emissions and virtual corrections leads to divergences or kinks at fixed orders. To obtain a precise measurement of the strong running coupling α s , one will have to resum the Sudakov shoulder logarithms that fall into the relevant regions [89]. However, in equilateral EEEC, since the three particles are separated by the same angle, the maximum angle is 2π 3 when all three particles fall into the same plane. This geometry constraint remains the same in higher-order perturbation theory so that the IR cancellation is guaranteed by Kinoshita-Lee-Nauenberg (KLN) theorem [95,96].

Triple collinear limit at next-to-leading power and beyond
The factorization theoroem at leading power (LP) has been well understood for different observables and different physical processes. It allows for the resummation of large logarithms to very high accuracy. Oppositely, much less is known for the factorization theorem and its violation at next-to-leading power (NLP). The complete factorization framework is still not established for NLP observables. In this subsection, we focus on the NLP contribution from the direct calculation point of view, where only a few cases have been carried out [22,[97][98][99][100][101][102]. More discussions can be found in Ref. [103].
At LP, the triple collinear EEEC is factorized as the hard function H = {H q , H g } and the jet function J = {J q , J g }, which both live in the flavor space. In momentum space, the EEEC jet function also decouples into the triple collinear phase space [104,105] and the 1 → 3 splitting functions [105][106][107]. Benefiting from the decoupling, the calculation was performed in Ref. [3], which becomes the first analytic calculation of a three-parameter jet substructure observable. From our EEEC result with full angle dependence in Section 3, it is straightforward to extract NLP contribution in the triple collinear limit. While the NLP contribution itself can provide comparison data for the study of NLP factorization, it can not provide hints toward NLP factorization. It is more interesting to explore a similar decoupling of the phase space and the integrand to NLP, and extract the NLP corrections from a direct computation.
For the purpose of extracting the triple collinear limit, we perform the following rescaling and expand the corresponding formula in λ order by order. To decouple the phase space measure and the integrand to NLP, we start from Eq. (2.14) and reformulate it as in the following, 1 σ tot where g (x 1 , x 2 , x 3 , z 1 , z 2 ) is used to represent the EEEC integrand. Notice that the upper bound of z 2 depends on x 3 , this makes the decoupling non-trivial at NLP. At LP, it is safe to expand the upper bound and the integrand separately, where the leading terms in λ directly gives us the decoupling at LP, as computed in Ref. [3]. At NLP, one can not just expand the integrand g (x 1 , x 2 , x 3 , z 1 , z 2 ) to the next-to-leading term while only keeping the leading term in the upper bound. One may expand both the upper bound and the integrand to next-to-leading terms, however, it doesn't make the computation simpler and also mixes a part of NNLP and beyond into the NLP contribution.
To extract the exact NLP contribution and make the decoupling explicit, we separate the interval of z 2 integration into the following three intervals, On the right-hand side of Eq. (4.21), the first term corresponds to the triple collinear phase space measure and contributes to LP and beyond, the second term starts to contribute at NLP, and the third term only contributes to NNLP and beyond. The right-hand side of Eq. (4.21) is formulated in a way that the i-th term contributes only at N i−1 LP and beyond. It is similar to Ref. [108] where the operators are organized in a way such that the i-th type operators contribute only at α i s and beyond. Let us focus on the second term, together with the integrand, we have where · · · only contributes to NNLP and beyond. In summary, the contribution up to NLP can be written as The intervals of both integration variables z 1 , z 2 in the above equation don't involve any kinematic variables x 1,2,3 . Therefore, we can safely expand the integrand in the triple collinear limit, which makes the computation of the NLP contribution as easy as LP. We emphasize that the contact term that is proportional to δ z 2 − (1 − z 1 ) is crucial to get correct result at NLP. The above method to extract NLP for triple collinear EEEC may also be useful to compute the NLP contributions for other observables. It is also straightforward to generalize the above method to NNLP and beyond. An alternative method to extract the NLP contribution is to first integrate over z 2 in Eq. (4.20) without performing any expansion. By simplifying the resulting integrand and performing the expansion, one can integrate over z 1 and obtain the NLP contribution. The simplified integrand with z 1 dependence is also useful as a cross-check of our final analytic formula, such that we also provide it in the ancillary file.
We use both methods to extract the triple collinear limit to NLP and find the same result. The validity of both methods is verified by the fact that no poles in the integration variables are generated when performing the expansion. We also compare the result with the final full analytic formula by setting x 1,2,3 to very small numbers, and we find the difference is indeed an NNLP contribution.
The result up to NLP in the triple collinear limit can be written as where λ is just used to track the expansion order and should be set to 1 at the end. The leading contribution C LP agrees with the result in Ref. [3]. The subleading term C NLP is new and also contains three color channels: The n f contribution is given below: Similarly, for other color factor contributions, we have and is composed of one logarithm log(x 1 ) and three transcendental weight-two functions: with Bloch-Wigner function D − 2 (z) defined in Eq. (4.7) and z introduced in Eq. (2.16). The simplicity of the above NLP results encourages us to go to NNLP and beyond. Using the second method, we expand the integrand to N 10 LP. Interestingly, we find that for n f contribution all polylogarithmic functions disappear at N 3 LP and beyond, which correspond to positive powers of λ in Eq. (4.24). In particular, we present the first few terms: Only π 2 remains at higher powers of the n f channel, while in other color channels, the polylogarithmic functions still show up. The simplicity may imply that there are some hidden symmetries, we leave it to future study.

Squeezed limit
Another interesting kinematic region is the squeezed limit: x 1 → 0, x 2 ∼ x 3 ∼ η and its permutations. Using our analytic formula for EEEC, it is straightforward to extract this limit. At LP, we find, with (4.32) The functional form is similar to N = 4 SYM, made of a single logarithm log(1 − η).
In fact, there are ambiguities in the definition of the squeezed limit. In the above expansion, we apply the isosceles constraint first and take the zero limit of the third angular distance. However, this is not a unique choice since we can start with other configuration constraints. The ambiguity obtains a geometry interpretation if we study the squeezed limit under the triple collinear limit. As observed in Refs. [41,42], the squeezed limit is accompanied by an angular dependence. If we adopt the z variable defined in Eq. (2.16), one of the squeezed limits is z → 0, and the expansion reads with r and t = e iθ introduced in Fig. 8a. Our definition x 1 → 0, x 2 ∼ x 3 ∼ η then corresponds to approaching z = 0 via the path θ = π 2 . In other words, the z point is forced to fall into a circle located in (1, 0) with the radius 1 (See Fig. 8b). This gives 1 σ tot (4.34) which agrees with the η → 0 limit of B(η). Another choice of the path is through θ = π 4 , as discussed in Ref. [3], with which the expansion is different: (4.35) Interestingly, we get identical results for Eq. (4.34) and Eq. (4.35) if we take the N = 1 SYM limit by setting T F = 1/2, n f = N c , C F = N c , C A = N c . It can be explained by looking at Eq. (4.33), there the expression becomes t-independent for the squeezed limit at LP, i.e., the coefficient of 1/r 2 .
While extracting the higher power corrections in the squeezed limit is pretty straightforward once a unique definition is given, the geometry interpretation is invalid beyond LP. Nevertheless, our result indicates studying the overlap among kinematic limits themselves (triple collinear limit and squeezed limit in this case) is also theoretically important. The structure of singularities becomes clear in such a joint kinematic limit and they will be useful in investigating jet substructure. It is also interesting to ask how we can organize the power corrections under joint kinematic limits in general. We leave these possible directions for future studies. Figure 8: (a) The triangle formed by the three angular distance √ x 1 , √ x 2 and √ x 3 under the triple collinear limit. We introduce the distance r from the origin to the top point and the angle θ between this side and the x-axis. (b). The squeezed limit with isosceles constraints x 1 → 0, x 2,3 ∼ η. In the triple collinear limit, this path becomes a circle |z − 1| = 1.

Summary
The energy correlator is one of the most important event shape observables widely used in both precision QCD and collider physics. Proposed in the 1970s, EEC has been playing an important role in various aspects of QCD measurements and jet physics studies, such as the precise measurement of strong coupling α s . Three-point energy correlator, which captures more information about the scattering events, can be more powerful for probing jet substructure. In this paper, we calculate the three-point energy correlator at leading order in electronpositron collisions and initialize the studies of EEEC kinematics. Instead of rewriting measurement functions as cut propagators and using IBP reduction, we approach the calculation with direct phase space integration. With appropriate parameterization of the phase space dP S 4 and kinematic space x 1,2,3 , we factorize out the Heaviside Θ function from the integral and rationalize all square roots, which allows performing the remaining integration. The QCD result is very similar to N = 4 SYM EEEC, in the sense that they share the same function space that is composed of polylogarithmic functions up to transcendental weight-two.
The simplification of the leading order EEEC involves two steps. Since HyperInt expresses the result in terms of GPLs, we need to convert it into polylogarithms and modify their arguments with dilogarithm identities, in order to meet Mathematica's branch prescription. Then we construct the raw function space by collecting rational coefficients and map it to N = 4 SYM EEEC function space in Ref. [5]. The linear relations between these two function spaces allow us to reduce the leading order EEEC in terms of the latter function space. With the simple function space as well as simplified rational coefficients, the file size of our analytic formula is small and the numerical evaluation is very fast. The simplicity strongly encourages us to analytically compute EEEC for gluon-initiated or bb-initiated Higgs decays in the near future.
Given the multiple angular distance dependence, the EEEC kinematic space becomes more interesting. Various kinematic limits remain unexplored. In Section 4, we discuss the equilateral limit x 1 = x 2 = x 3 = x, triple collinear limit x 1 ∼ x 2 ∼ x 3 → 0 and squeezed limit x 1 ∼ 0, x 2 ∼ x 3 ∼ η, and the analytic results in all limits become very simple. Under equilateral limit, the angular distance x is cutoff at x = 3 4 , which corresponds to the coplanar configuration. Regarding the triple collinear limit, we present a method that allows us to directly compute the next-to-leading power corrections from expanding the EEEC integrand. The NLP result is simple and shares the same collinear function space as the LP. We also discuss the overlap between the triple collienar limit and the squeezed limit, where the ambiguity of the squeezed limit definition receives a geometry interpretation.
In fact, EEEC provides a large playground for studying factorization as well as its violation for different configurations. Using our analytic result, it is straightforward to extract the needed ingredients like jet functions. Deriving the factorization theorem for specific limits and performing resummation for EEEC could be theoretically interesting and phenomenologically important. The simple mathematical structure also makes EEEC a good candidate for understanding NLP corrections and beyond. While the resummation in the triple collinear limit is in progress, the equilateral limit could be a window to investigate symmetric trijet events in e + e − collisions. Furthermore, all these future directions can be generalized to the analysis at hadron colliders and provide new opportunities for studying Higgs phenomenology and top physics.
the European Union's Horizon 2020 research and innovation programme grant agreement 101019620 (ERC Advanced Grant TOPUP), and by the Swiss National Science Foundation (SNF) under contract 200020-204200.

A Usage of ancillary files
In the ancillary files, we provide all the results in this paper. The usage of each file is as follows.
• EEECinQCD.nb: The main Mathematica notebook. We import other files in the notebook and define a set of commands to compute EEEC. Explicitly, eeec[{x1,x2,x3}] gives the value of EEEC using the analytic formula. The options "Color" and "Parton" can be used to separate different color structures or partonic subprocesses.
-eeecEqu[x] gives the value of equilateral EEEC using the analytic formula.
• Numerical.wl: The one-fold numerical integral for EEEC in QCD, where the integrand is saved in the file EEEConefold. The main function is eeecNum.
• EEECanalyticfull: The full analytic expression of EEEC in QCD.
prefactor: The overall normalization factor.
• EEECequilateral: The analytic expression of equilateral EEEC in QCD.
bwrep: The replacement rule for Bloch-Wigner functions in equilateral EEEC.