Three Point Energy Correlators in the Collinear Limit: Symmetries, Dualities and Analytic Results

Energy Correlators measure the energy deposited in multiple detectors as a function of the angles between the detectors. In this paper, we analytically compute the three particle correlator in the collinear limit in QCD for quark and gluon jets, and also in $\mathcal{N}=4$ super Yang-Mills theory. We find an intriguing duality between the integrals for the energy correlators and infrared finite Feynman parameter integrals, which maps the angles of the correlators to dual momentum variables. In $\mathcal{N}=4$, we use this duality to express our result as a rational sum of simple Feynman integrals (triangles and boxes). In QCD our result is expressed as a sum of the same transcendental functions, but with more complicated rational functions of cross ratio variables as coefficients. Our results represent the first analytic calculation of a three-prong jet substructure observable of phenomenological relevance for the LHC, revealing unexplored simplicity in the energy flow of QCD jets. They also provide valuable data for improving the understanding of the light-ray operator product expansion.


Introduction
Observables which probe the flow of energy in quantum field theories are interesting both from a theoretical perspective for studying the Lorentzian limits of field theories, as well as from a practical perspective, where they can be used as event shape or jet substructure observables to measure properties of QCD, and search for new physics with jets (see e.g. [1,2] for a review).
From the theoretical perspective, one of the simplest observables is the Energy-Energy Correlator (EEC), defined as [3,4] (1.1) Here E i and E j are the energies of final-state partons i and j in the center-of-mass frame, χ ij is their angular separation, and dσ is the product of the squared matrix element and the phase-space measure. The EEC can also be defined in terms of correlation functions of ANEC operators [5][6][7][8][9][10][11] E( n) = ∞ 0 dt lim r→∞ r 2 n i T 0i (t, r n) , (1.2) where the limit is taken with the retarded time u = t − r fixed, and we integrate over the advanced time v = t + r. In terms of these ANEC operators, the EEC is given by the four-point Wightman correlator for some source operator O that produces the localized excitation. This provides a connection between event shape observables and correlation functions of ANEC operators, allowing the study of event shapes to profit from recent developments in the study of ANEC operators. Conversely, the EEC provide a concrete situation for studying the behavior of ANEC operators using jets at the LHC. There has recently been interesting progress in the understanding of the EEC from a number of different directions. The EEC has been computed for generic angles at next-toleading order (NLO) in QCD [12,13], and up to NNLO in N = 4 SYM [11,14]. It has also been computed numerically in QCD at NNLO [15,16]. There has also been progress in understanding the singularities of the EEC, which occur as z → 0 (the collinear limit) and z → 1 (the back-to-back limit). In the back-to-back limit, the EEC exhibits Sudakov double logarithms, whose all orders logarithmic structure at leading power is described by a factorization formula [17,18], and whose subleading power logarithms have been resummed in N = 4 SYM [19] using recent developments in understanding the structure of subleading power corrections for this class of observable [20]. In the z → 0 limit, which will be studied I + One of the simplest observables from the theoretical perspective is the Correlator (EEC), defined as [2,3] Here E i and E j are the energies of final-state partons i and j in the center and their angular separation is ij . d is the product of the squared matrix phase-space measure. The EEC can also be defined in terms of correlation fu operators [4][5][6][7] where it is given by for some source operator O. This provides a connection between event shape correlation functions of ANEC operators allowing the study of event shape recent developments in the study of ANEC operators, and conversely, the concrete situation for studying the behavior of ANEC operators. There has recently been significant progress in the understanding of t number of di↵erent directions. For generic angles, the EEC has been comp leading order (NLO) in QCD [8,9] for both an e + e source, and Higgs dec and up to NNLO in N = 4 SYM [7,10]. It has also been computed numeri NNLO [11,12].
There has also been progress in understanding the singularities of the EEC z ! 0 (the collinear limit) and z ! 1 (the back-to-back limit). In the back-to EEC exhibits Sudakov double logarithms, whose all orders logarithmic struct by a factorization formula [13,14]. In the z ! 0 limit, which will be studie the EEC exhibits single collinear logarithms, originally studied at leading lo in [15][16][17][18][19]. Formulas describing the behavior of the EEC in the collinear lim derived in [20] for a generic field theory, and in [21][22][23][24] for the particular case limit is of theoretical interest for studying the OPE structure of non-local o phenomenological interest as a jet substructure observable.
The two-point correlator is particularly simple since it depends on a si Indeed, in a conformal field theory (CFT), its behavior in the collinear limit power law One of the simplest observables from the theoretical perspective is the Ene Correlator (EEC), defined as [2,3] Here E i and E j are the energies of final-state partons i and j in the center-of-m and their angular separation is ij . d is the product of the squared matrix eleme phase-space measure. The EEC can also be defined in terms of correlation functio operators [4][5][6][7] where it is given by for some source operator O. This provides a connection between event shape obse correlation functions of ANEC operators allowing the study of event shapes to recent developments in the study of ANEC operators, and conversely, the EEC concrete situation for studying the behavior of ANEC operators. There has recently been significant progress in the understanding of the E number of di↵erent directions. For generic angles, the EEC has been computed leading order (NLO) in QCD [8,9] for both an e + e source, and Higgs decaying and up to NNLO in N = 4 SYM [7,10]. It has also been computed numerically NNLO [11,12].
There has also been progress in understanding the singularities of the EEC, wh z ! 0 (the collinear limit) and z ! 1 (the back-to-back limit). In the back-to-bac EEC exhibits Sudakov double logarithms, whose all orders logarithmic structure i by a factorization formula [13,14]. In the z ! 0 limit, which will be studied in the EEC exhibits single collinear logarithms, originally studied at leading logarit in [15][16][17][18][19]. Formulas describing the behavior of the EEC in the collinear limit we derived in [20] for a generic field theory, and in [21][22][23][24] for the particular case of a limit is of theoretical interest for studying the OPE structure of non-local operat phenomenological interest as a jet substructure observable.
The two-point correlator is particularly simple since it depends on a single in this paper, the EEC exhibits single collinear logarithms, originally studied at leading logarithmic order in [21][22][23][24][25]. Formulas describing the behavior of the EEC in the collinear limit were recently derived in [26] for a (not necessarily conformal) gauge theory, and in [27][28][29][30] for a generic conformal field theory (CFT). The z → 0 limit is of theoretical interest for studying the operator product expansion (OPE) in Lorentzian signature, and is of phenomenological interest as a jet substructure observable for the LHC. The two-point correlator is particularly simple since it depends on a single variable, z. Indeed, in a conformal field theory (CFT), its behavior in the collinear limit is fixed to be a power law [26,29,30] where C(α s ) is a constant independent of z, γ N =4 1 (α s ) is the universal anomalous dimension of local twist-2 spin-3 operator, and we use Σ(z) to denote the cumulant [26]. In a non-conformal theory, such as QCD, one can derive a time-like factorization formula [26] Σ(z, ln Q 2 µ 2 , µ) = 1 0 dx x 2 J(ln zx 2 Q 2 µ 2 , µ) · H(x, where H is the coefficient functions for semi-inclusive fragmentation, and J is the jet function for the EEC. One finds that instead of involving a single fixed moment, the jet function involves with the timelike splitting kernel d J(ln zQ 2 µ 2 , µ) d ln µ 2 = 1 0 dy y 2 J(ln zy 2 Q 2 µ 2 , µ) · P T (y, µ) .
(1.6) Therefore, in both conformal and non-conformal theories, the dependence on the angle z is completely fixed by symmetry or renormalization group arguments, as expected for a twopoint observable.
In this paper, we study the three point energy correlation function, which we will refer to as the triple correlator, or EEEC. The EEEC is defined as (1.7) where i,j and k run over all the final-state massless partons. We will focus on the collinear limit of the EEEC, where the angles between detectors are small, x 1 , x 2 , x 3 1, but we will not assume that there is any hierarchy between the angles, so that the result will have a dependence on the ratio of angles that is not (yet) predicted by symmetries or renormalization group arguments. The complete calculation of the triple correlator for generic angles is also interesting, and will be considered in future work.
There are a number of motivations for extending the understanding of the energy correlators beyond the simplest case of the two point correlator. First, the three point correlator is the first time that a non-trivial dependence on the angles can occur, and so the result provides a more non-trivial probe of the structure of the OPE in Lorentzian signature. Secondly, the triple correlator directly probes multi-particle correlations within a jet, which is extremely interesting from the perspective of jet substructure. Observables capturing multi-particle (> 2) correlations have not been analytically computed before, and we believe that this provides a significant advance for jet substructure that will have phenomenological applications at the LHC.
In this paper we directly compute the triple correlators in QCD for both quark and gluon jets, and also in N = 4 super Yang-Mills (SYM) theory. We show that the result can be compactly written in terms of polylogarithmic functions of a complex variable that encodes the two conformal cross ratios, and we highlight a number of interesting features of the results. We also find an intriguing duality between the integrals associated in computing the EEEC in the collinear limit, and Feynman parameter integrals for loop amplitudes. This duality maps energies to Feynman parameters, and the angles of the energy correlators to dual coordinates. Due to the well established techniques for the integration of Feynman parameter integrals, this allows the result to be computed relatively easily, and we believe that it will also enable calculations of higher point correlators.
For the case of quark jets in QCD, we are able to compare our analytic result with a numerical calculation using the full e + e − → 4 parton matrix elements using Nlojet++ [31,32] and Event2 [33]. We find good agreement for all partonic channels in the singular limit, providing strong evidence both for the correctness of our results, and for the factorization properties of the EEEC in the collinear limit.
An outline of this paper is as follows. In Sec. 2, we describe our calculation of the EEEC using timelike factorization and the triple collinear splitting functions. In Sec. 3 we discuss the -4 -symmetries of the observable, the parametrization of the result, and the functions that will appear. In Sec. 4 we illustrate an intriguing duality between the integrals encountered in our calculation for the energy correlator, and Feynman parameter integrals for loop amplitudes. In Sec. 5 we present results for the EEEC in N = 4, and for both quark and gluon jets in QCD, and discuss their structure. In Sec. 7, we compare our analytic predictions with numerical calculations using the full e + e − → 4 parton matrix elements. We conclude in Sec. 8 and discuss many future directions for further understanding the energy correlators in the collinear limit, and for phenomenological applications.

Factorization in the Collinear Limit
To compute the behavior of the EEEC in the collinear limit, we will use timelike collinear factorization, and extend the factorization formula presented in [26] for the EEC. Following [26], we factorize the EEEC in the collinear limit into a hard function, H, and a jet function, J, both of which are vectors in flavor space. This factorization is shown schematically in Fig. 2. The hard function depends on the source, but is independent of the measurement, and is therefore identical for the EEEC and EEC (and more generally for an arbitrary number of correlators). The jet function, J, describes the dependence on the EEEC measurement, and will be the focus of this paper. Corrections to this formula are suppressed by powers of the angles, and are described by higher twist jet functions. In this paper we will only study the fixed order properties of the jet function, leaving a study of the all orders resummation structure of the factorization formula to a future publication.
The jet function for the EEEC for a quark (or antiquark) jet is defined as and for a gluon jet, as Here χ n and B a,µ n,⊥ are gauge invariant collinear quark and gluon fields in SCET [34][35][36][37]. The operator M EEEC implements the EEEC measurement in the collinear limit, Technically speaking, it can be written as an operator in terms of the stress energy tensor [5,6,8,10,38,39], however, here we will only need its expression in momentum space. For -5 -

Fixed by kinemat dimension ana
• Cumulant  Figure 2: A schematic of the factorization for the triple correlation function in the collinear limit. In the collinear limit, we can factorize on a single partonic state produced by the hard function, H. The measurement is then imposed on the jet function, which can be computed perturbatively for different partonic channels.
the EEEC, there must be three partons in the jet to have a non-trivial functional dependence (beyond contact terms).
To compute the jet function of Eqs. (2.1) and (2.2), we use its equivalence to a calculation of the full QCD splitting functions with the EEEC measurement function inserted. For an extensive discussion of the equivalence between jet functions and results integrated against collinear splitting functions, see [40]. We must therefore simply integrate the triple collinear splitting functions against the measurement function for the EEEC, over the triple collinear phase space.
The triple collinear splitting functions P ijk were computed in [41,42], and are summarized in App. A for convenience. The triple collinear phase space expressed in terms of the Mandelstam invariants s ij , s 123 and the momentum fractions ξ i = 2E i /Q of the partons is given by [43] dΦ (3) c = ds 123 ds 12 ds 13 ds 23 δ(s 123 − s 12 − s 13 − s 23 )dξ The angles between the partons are related to the Mandelstam variables and energy fractions used in the phase space parametrization as where in the second equality we use the collinear approximation. We can therefore write the result for the EEEC jet function at the lowest non-trivial order as This expression is also true for a generic measurement function, M, but in general it is hard to perform the integral analytically. This is due both to the constraints of the measurement function, and from the Θ(−∆) constraint appearing in the collinear phase space. However, here we will find that a remarkable simplification occurs. First, we can change variables to express the integration over the Mandelstam variables in terms of integrals over the x i , using and Now, we note that in the collinear limit, we have that and √ x 1 , √ x 2 , √ x 3 are subjected to the constraint that they must form the sides of a triangle.
Using Heron's formula, which expresses the area of a triangle, A, in terms of the lengths of its sides, we find that This implies that the measurement function automatically guarantees that the Θ(−∆) constraint is satisfied. This enables the integrals over s 123 , and the energy fractions ξ i , to be performed analytically, while the integrals over x i are fixed by the measurement function. This can be compared with other examples, for example planar flow considered in [44], where analytic integration is extremely difficult due to the Θ(−∆) constraint. In Sec. 4 we will show that the particular structure of the measurement function for the energy correlators which completely fixes the angle, but integrates over the energies allows the phase space integrals appearing in its calculation to be mapped to Feynman parameter integrals for loop amplitudes. This particular factorization of the measurement function into energies and angles is non-standard in jet substructure, but we find that it is extremely convenient and leads to many of the nice properties of the result.

Parametrization, Symmetries, Functions and Constraints
In this section we discuss our parametrization of the EEEC and its associated symmetries, and present arguments for the classes of functions that will appear. In Sec. 4 we give another way of understanding the functions that appear in the result by illustrating a relationship between the integrals for the energy correlators in the collinear limit and Feynman parameter integrals for loop amplitudes.

Symmetries on the Celestial Sphere and Parametrization
To study the structure of the energy correlators, it is convenient to switch from a parametrization in terms of energies and angles, to coordinates on the celestial sphere. We therefore parametrize the momentum of particles as Here is a rescaled energy fraction, and the complex coordinate z i are coordinates on the S 2 celestial sphere. In terms of these coordinates, the standard Mandelstam invariants are Note that in terms of standard angles, we have where the final equality holds in the small angle limit. The interesting feature of the energy correlator observables is that the measurement function fixes the coordinates on the sphere, but integrates freely (up to momentum conservation) over the energies. The result is therefore a function of the complex variables z i , and it is therefore interesting to begin by studying the symmetries acting on the z i variables. Lorentz transformations in four dimensional Minkowski space act as the global conformal group, SL(2, C), on the z i coordinates on the celestial sphere [45][46][47]. Explicitly, For a pedagogical review, see e.g. [48]. However, for energy correlators at generic angles, the Lorentz symmetry is broken by the timelike vector Q µ of the source from SO(3, 1) → SO(3). For SO(3), this leaves three generators, which act on the coordinates on the celestial sphere as translations, z → z + b , b ∈ C and rotations, z → e iθ z, θ ∈ R (see also the discussion in [29]). Due to collinear factorization, an additional symmetry is restored at small angles, which corresponds to boosts along the jet direction. Physically this arises since the jet function in Eqs. (2.1) and (2.2) depends only on the source through the momentum conserving delta function for the momentum along the jet direction δ[n · Q − Q(ω 1 + ω 2 + ω 3 )(1 + |z| 2 )]. Here we have used thatz z 1 z 2 z 3 at leading power in the collinear limit, wherez denote the jet direction on S 2 . Without loss of generality, in the rest of this paper we use the SO(3) symmetry to set the jet direction toz = 0, so that ω i = ξ i at leading power. If we apply a general SL(2, C) transform to this, then we see from Eq. (3.4) that if c = 0, then ω i → |d| 2 ω i , which can then be pulled out of the delta function. 1 This boost, which acts as ω i → |d| 2 ω i , rescales the coordinates on the sphere as z i → z i /d, and therefore acts as a dilatation. Therefore, in the collinear limit, in terms of the coordinates on the sphere, we have the following symmetry group that will act on the energy correlators These symmetries can be interpreted in two ways. First, this is the conformal group of C ∪ {∞}, with the point at infinity fixed. This has a fairly natural interpretation with the point at infinity being the source for the jet function in the collinear limit. The only missing generator for the SL(2, C) transformations on the sphere that is broken in the collinear limit is inversions. We can therefore interpret the EEEC as a three point function with these symmetries. Unlike when inversions are included, in this case the three point function is no longer completely fixed and is a non-trivial function of scale invariant cross ratios. Alternatively, and perhaps more promising from the speculative perspective of trying to develop a field theory on the celestial sphere describing the energy correlators, we can view the EEEC as a four point correlation function. For generic angles, the source momentum Q µ satisfies Q 2 > 0, and therefore lives as a point in the bulk of H 3 . 2 When one factorizes in the small angle limit, as in Fig. 2, one expands to leading power in the limit that the invariant mass of the intermediate state q 2 → 0, and therefore the source for the jet also moves to the ∂H 3 = S 2 boundary where the other massless final state particles live. This is shown schematically in Fig. 3. We can therefore interpret the EEEC as a four point correlator on the S 2 where we have already put one point at infinity (or at the very least, it possesses the symmetries of such a four point correlator).
Therefore, interpreting the EEEC as a four point correlator on S 2 , we use the standard parametrization of a conformal four point correlator, where we put the ANEC operators at {0, 1, z}, and the fourth point at ∞. This is shown in Fig. 4a. Algebraically, we then have A technically convenient property of the z parametrization is that it rationalizes the area of the triangle, which we will see appears as an argument of the transcendental functions in the result. The area of the rescaled triangle can be written as Note that for a triangle, the expression in the square root is guaranteed to be positive by Heron's formula. Finally, there is also a discrete S 3 × Z 2 symmetry that acts on the result. The Z 2 acts as complex conjugation and is equivalent to parity. The final result is invariant under the Z 2 , 3 but individual transcendental functions can be classified as being either even, or odd under this symmetry. The S 3 acts as a permutation symmetry on the three identical energy detectors. In the z variables, the S 3 acts as This is simply the anharmonic group, which is the stabilizer of {0, 1, ∞}. It is also the standard action of S 3 on the arguments of harmonic polylogarithms, which will be a fact that will be useful later when understanding the function space.
To restrict to a region in the complex plane that uniquely specifies a triangle without redundancies, we must restrict to the region shown in Fig. 4a. This region can be determined by enforcing that the sides of the triangles are ordered, x 3 > x 2 > x 1 . This is the region that is relevant for e.g. an experimental analysis. However, when expressing the result as a function of z, it is useful to view the function as defined globally in the complex plane, using the S 3 action.
With this set of coordinates in mind, we can now perform a change of variables to write the cross section as a function of x L = x 3 , and z,z. We find 4 The goal is therefore to determine the non-trivial function G(z).

Transcendental Functions
The function G(z) is obtained by integrating a rational function of z (namely the tree level splitting functions), and can only have singularities at {0, 1, ∞}, which correspond to the squeezed (or OPE) limits (see Secs. 3.3 and 6). This places strong constraints on the possible functions that can appear. In particular, it is a known mathematical result that such functions are harmonic polylogarithms [49,50] multiplied by rational functions of z. Furthermore, from the structure of the integrals the transcendentality of the functions appearing in the result is at most two.
Due to the action of the anharmonic group, the result is defined on the entire plane, and it must be single valued. There are two ways that this can occur. The first way is if the functions dependend on |z| 2 or |1 − z| 2 . The second is if the polylogarithms are single valued harmonic polylogarithms (SVHPLs) [51]. These functions are well studied, and have appeared elsewhere in the physics literature, for example in the study of the Regge limit [52], in anomalous dimensions of cusped Wilson lines [49], or in the multi-loop soft gluon emission amplitude [53]. In particular, it is well known that at weight 2, there is a unique SVHPL, the Bloch-Wigner dilogarithm, which is defined as This function has an odd parity under the Z 2 symmetry, and also has a very simple action under the S 3 action discussed before, namely [54] (3.11) Another useful identity relating to the Bloch-Wigner dilogarithm is that it can be written in terms of a unit argument In our particular case, the arguments are interpreted as the angles of the triangle. This shows that one will find polylogarithms of roots of unity when the result is evaluated for particular triangles. For example, for an equilateral triangle, one will find that the result is expressed in terms of polylogarithms of the third root of unity. It is also well known that D − 2 (z) is the area of the ideal tetrahedron in hyperbolic three space, H 3 , with points at {0, 1, z, ∞}. The celestial sphere is the boundary of H 3 , and therefore this can be interpreted as the volume of the tetrahedron formed by the three points of the energy correlators, and the point at infinity which originated the jet in the collinear limit, as shown in Fig. 3. The appearance of D − 2 (z) is therefore natural since the SL(2, C) symmetry on the S 2 acts as diffeomorphisms in the bulk H 3 . It would be interesting to understand this point better, and if a similar interpretation extends to higher points.
We will additionally need transcendental functions of |z| 2 or |1 − z| 2 . Our result is not of uniform transcendentality, and therefore we must describe functions of both weight 1 and 2. At weight 1, we have simply the logarithms log|z| 2 , log|1−z| 2 . These functions are even under Z 2 . We will also need weight 2 even functions. It is convenient to introduce the following particular variant of the dilogarithm (and its S 3 permutations) 5 where the logarithms are fixed by demanding that the result is a pure functions of either |z| 2 and |1 − z| 2 that has no double logs as z → 0, 1, ∞. This is a property of the result at this order, since there are no logarithms appearing in the OPE limit. In both QCD and N = 4 we will find that D + 2 (z) is sufficient to describe the even weight 2 part of our result, and that the functions {log 2 |1 − z| 2 , log 2 |z| 2 } do not appear independently.
We therefore expect to find the following independent functions in our result. (3.14) These functions will also appear in the permutations under the S 3 action, and will be multiplied by rational functions of z,z.
Here we have argued for these functions based purely on the form of the functions being integrated, and using known mathematical results. In Sec. 4, we will provide another way of arriving at these same functions by showing explicitly that all integrals appearing in the calculation of the N = 4 EEEC can be directly mapped to one-loop Feynman integrals for amplitudes. Here the side lengths, |z| 2 and |1 − z| 2 will be associated with Mandelstam invariants. Since the symbol [55] of amplitudes obey a first entry condition [56], namely that the first entry is a Mandelstam invariant, this guarantees that the first entry of our result is either |z| 2 or |1 − z| 2 . This is true for both D + 2 (z), D − 2 (z).

Squeezed Limits and Behavior at Infinity
A particularly strong constraint on our result comes from the squeezed (OPE) limits, namely z → 0 or z → 1. By the S 3 symmetry, these OPE limits are also related to the limit z → ∞.
For the tree level three point correlator that we are computing, we have the following OPE limits The coefficient of proportionality depends on the theory, and a more detailed discussion of the OPE limit, and the calculation of the relevant constants, will be given in Sec. 6. In particular, this result implies that all weight 1 and weight 2 functions must cancel in these limits. This power law scaling remains true beyond tree level, but it is modified by logarithms.
The behavior at infinity also strongly constrains the form of the polynomials in the result. To have higher order polynomials that describe power corrections to these limits, one must also have higher order denominators, and intricate cancellations are required to respect the required constraints. Similar cancellations were also seen in the case of the two point energy correlator in [12,13].

Collapsed Triangles
Finally, we find that another limit which plays an important role in understanding the EEC is the limit z →z, namely where the triangle collapses to a line. 6 This limit is non-singular for a generic value of (z), i.e. away from the OPE limit. Interestingly, we find that the approach to this limit largely controls the complexity of the answer, and is the primary difference between the QCD and N = 4 results.
To understand why the z →z limit plays an important role, we note that one of the two weight two functions that appears in our result, D − 2 (z), vanishes as z →z, since it is odd. Since the full result for the cross section is even, D − 2 (z) must be multiplied by an odd power of (z −z). If this power is positive, then this function completely vanishes on the line z =z. However, if the power is negative, then one must Taylor expand D − (z) in powers of (z −z) allowing this function to contribute on the line z =z.
In the N = 4 case, we find that only the following functions appear On the other hand, in QCD, we find powers up to (z −z) −11 . This leads to a very intricate set of cancellations as z →z. Furthermore, as soon as one has these large denominators in one part of the result, they immediately appear in other parts of the result. For example, assume that the function D − 2 (z)/(z −z) 3 appears in the answer. Taylor expanding the function as we see that the first term, which is now a weight 1 function must be cancelled as z →z to have a well behaved answer. This requires an equal complexity for the polynomials of the weight 1 terms. With 11 inverse powers, this generates an extremely complicated pattern of cancellations. It would be interesting to understand if one can predict before hand the powers of (z −z) that appear and to understand what generates them. In Sec. 4 we give some progress towards this for the case of N = 4, where we are able to write the result as a rational sum of Feynman parameter integrals.

Relation to Feynman Parameter Integrals
In this section we discuss an interesting relation between the integrals that are encountered when computing energy correlators in the collinear limit, and standard one loop integrals for scattering amplitudes. We will show that, at least for the three point correlator, we can interpret the energy integrals as Feynman parameter integrals, and the angles between the energy correlators as dual coordinates of Feynman graphs. This relationship is particularly interesting due to the fact that techniques for the calculation of one-loop integrals are highly developed, and their structure is much better understood.
-14 - The integrals that arise when considering the EEEC in the collinear limit take the form 7 where we have set the jet direction toz = 0 such that ω i = ξ i is the parton momentum fraction, and Here we will work in the coordinates on the celestial sphere, where To illustrate the correspondence with Feynman parameter integrals, we choose first a particularly simple term in the N = 4 splitting function (the complete splitting functions can be found in the Appendix) Writing we then find that the integral we must consider in the collinear limit is We immediately see that this resembles a Feynman parameter integral for the three mass triangle. To make this relation precise, we consider the Feynman parametrization of the integral We are then able to identify (4.8) 7 We note that in QCD it is conventional to pull out a factor of 1/(s123) 2 from the splitting functions. Due to the simpler structure of the splitting functions in N = 4, this is no longer convenient there. Here we have written it without the factor of 1/(s123) 2 pulled out.
-15 - where we have the association We therefore find an interesting equivalence between the integrals for the energy correlators in the collinear limit, which are integrals of squared amplitude over phase space, and Feynman integrals for loop amplitudes. Note that this rewrite is different from the usual application of optical theorem. This correspondence can be made even more clear if we use dual coordinates [58] In this case, we simply have a correspondence between the dual coordinates and the coordinates on the celestial sphere and between energy fractions and the Feynman parameters α i This equivalence is shown schematically in Fig. 7. We find this equivalence to be quite striking, and deserving of further investigation. From this equivalence we can already begin to easily understand the functions that will appear in the result for the triple correlator. In [59] it is advocated that the n-point integral should be considered in D n , where D n = n−2 if n is even and D n = n+1−2 if n is odd. These integrals are conjecturally of uniform transcendentality. They also have a correspondence with hyperbolic geometry. We therefore look at the three mass triangle integral in D = 4 − 2 as a base integral. The three mass triangle integral is given by (here we follow closely the notation of [60]) where D − (z) is the Bloch-Wigner function, as defined above, and Therefore, we arrive at the Bloch-Wigner function from a completely different perspective. For this particular term in the splitting function, the result of the integral for the EEEC can be expressed in terms of a three mass triangle. Naively one might have anticipated a dependence on the fourth leg that initiated the splitting. To understand how this physically occurs, we can consider the "eikonal" terms in the splitting function, taking as an example one that involves (ω 1 +ω 2 ) in the denominator. For a jet along the direction n µ = (1, 0, 0, 1), these terms can be thought of as originating from a Wilson line along the directionn µ = (1, 0, 0, −1), and therefore should introduce a dependence on the fourth direction.
To illustrate the mapping of these eikonal terms into Feynman integrals, we consider the particular contribution to the triple correlator Now, we consider the form of a box integral with an eikonal propagator, whose Feynman parametrization is given by (an implicit +i0 + for the denominators is assumed) The corresponding Feynman diagram is depicted in Fig. 6. We can now rewrite Eq. (4.13) Figure 6: One-loop box with an eikonal propagator.
-17 -in this way (4.15) We can therefore make the identification where the external momenta of the eikonal box integral are given bȳ This eikonal box integral gives rise to the the dependence on the function described previously, among other terms. Based on this result, we believe that it would be interesting to investigate the behavior of the multipoint energy correlators in the conformal fishnet theory [61]. 8 We have seen that one source of complexity for the transcendental functions involved arises from the "eikonal" terms in the splitting functions (those that involve the momentum fractions z i in the denominators). These terms are typically coming from the perturbative expansion of collinear Wilson lines. In the conformal fishnet theory, the gauge bosons are decoupled, and the splitting functions are purely functions of Mandelstam invariants. In this case there may be a simpler relation between Feynman diagrams and energy correlators in the collinear limit. Indeed, in the conformal fishnet theory, the 1 → 3 splitting function is simply given by P 1→3 ∼ 1/(s 123 ) 2 , and so one finds that the result for the EEEC can be written purely in terms of the three mass triangle (Bloch-Wigner function), and does not involve the eikonal box integral. It would be interesting to understand what happens for higher point energy correlators.
The mapping between the integrals appearing in the EEEC, and one loop Feynman integrals, in addition to providing a technical simplification for computing the result, also provides a significant number of insights into the structure of the results. In particular • It provides an understanding of the functions that will appear in the result, namely polylogarithms. Figure 7: Naive correspondence between multipoint energy correlators and dual coordinates for Feynman graphs. This correspondence is broken by eikonal terms in the splitting functions that introduce eikonal propagators into the dual Feynman graphs.
• It provides a first entry condition [56], namely that the first entry in the (1, n − 1) coproduct must be a Mandelstam. In the case of the three point correlator, the Mandelstams are mapped to zz and (1 − z)(1 −z), and we find that this condition is indeed satisfied by our result. This condition is not something that is immediately clear from the event shape perspective. It would also be interesting to understand how this condition arises from the perspective of correlation functions [9][10][11]14] or from the lightray OPE [27][28][29].
• The factor of 1/(z −z) appearing in front of the Bloch-Wigner function is identified as the leading singularity [63] of the three mass triangle diagram. The physical significance of this factor is not clear from the perspective of the EEEC.
• The result for the three mass triangle in d = 4, apart from the 1/(z −z) prefactor, is a pure function. From this perspective, the polynomials that appear in the result for the EEEC arise from integration by parts or dimension shift identities arising from factors of ω i in the splitting functions. This may provide a convenient way of understanding these polynomials. Alternatively, it may indicate that the result for the energy correlators can be most compactly written in terms of Feynman parameter integrals in shifted dimensions.
More speculatively, it is well known that one-loop Feynman diagrams have an interpretation as volumes [64][65][66][67][68][69][70][71]. It seems that at least in this simple case that we have illustrated here, the EEEC is directly realizing this dual geometry associated with the Feynman diagrams to which it can be associated. We find this is particularly interesting, and believe that it merits investigation at higher points.

Possible Generalization to Higher Points
Although we will not consider higher point correlators in this paper, it is interesting to speculate whether this correspondence can provide a guide for higher point correlators. First, we must note that it may not be true that at higher points one can always map the energy correlator integrals to one-loop Feynman diagrams. However, it is generically true that one will have integrals of the form perhaps with some Gram determinant constraints. Even if the integrals do not take precisely the form of Feynman diagrams due to the structure of the splitting function P 1→N , we believe that thinking about the integrals in relation to Feynman parameter integrals may be useful in two ways. First, and most practically, regardless of the exact form of P 1→N , these integrals will have similar form to Feynman parameter integrals, and therefore many of the techniques can be carried over. Second, and more conceptually, in the three point case, the geometry, and its simple parameterization were suggested by the dual Feynman graph, and could be understand largely just from the 1/s 2 12···N factor. It is therefore possible that this term will continue provide a guide at higher points for what variables are convenient for parametrizing the result.
If we assume P 1→N = 1, as would occur if one had a scalar theory with an n + 1 point interaction, then, as mentioned above, we have a correspondence between the distances on the sphere between the points of the energy correlator and the dual coordinates (4.20) In particular, this identification gives the sides of the energy correlator shapes as the masses of the external legs, and the chords as Mandelstam invariants. We note that there are of course 3n − 10 invariants plus the masses, and there are 3n − 10 chords. Unlike the case of a Wilson loop that has been more thoroughly studied, these are not null polygons. However, higher point integrals with external masses have been studied, see e.g. [72,73]. It is also interesting to note that the dilatation symmetry for the energy correlators would correspond to a dual conformal symmetry [58,74,75] of the corresponding Feynman diagram. An intriguing test case will be for the four point EEEC. The one loop four mass box integral can be written in terms of the Bloch-Wigner function [76] (see also [77] for a detailed discussion) if one uses the variables defined as where is the Kallen function. It will be interesting to see to what extent this holds true for the four point energy correlator, and whether this duality determines the geometry. A first step in this direction will be interesting to understand the structure of the 1 → 4 splitting functions in QCD or N = 4 to see to what extent different terms can be manipulated into the form of Feynman diagrams. The 1 → 4 quark splitting function appeared recently in [78].

Analytic Results for the Three-Point Correlator
In this section we give the result for the non-trivial shape dependence of the EEEC for N = 4 SYM, as well as for quark and gluon jets in QCD. We then briefly discuss the differences between the QCD and N = 4 results.

Result in Terms of Loop Integrals
Using the relation with Feynman parameter integrals, we can write down a compact representation of the result for the EEEC in N = 4 in terms of three-mass triangle integral, the three-mass box with an eikonal propagator integral, and a rational term that incorporates the OPE limit. Explicitly, we have where , and the permutations refer to the S 3 action in Eq. (3.8), including the corresponding Jacobian (see Eq. (5.14) for an explicit example). The triangle and box integrals appearing here are UV and IR finite, and can be defined in integer dimension, d = 8 and d = 10. One can also use the powerful toolkit of Feynman integral calculation, such as Integration-By-Parts (IBP) reduction [79]. To do that, it is necessary to continue the dimension to d = 8 − 2 and d = 10 − 2 , respectively. The reason is the the master integrals obtained from IBP will generically contain UV or IR divergences, but the physical sum do not. We give the explicit results for these integrals in the appendix A.3.
We find the simplicity of this result to be quite remarkable. The unusual dimensions of the integrals in this result mean that when they are expressed in terms of four dimensional integrals, which are pure functions, complex polynomials are generated. This also leads support to the idea that the result can be easily bootstrapped if the correct basis of functions is chosen. In particular, in Eq. (5.1), with the exception of the OPE term, R(z), which is fixed by the behavior in the OPE limit, the coefficients of the other functions are rational.
Due to the remarkable simplicity of this result (compare with the expanded result in Eq. (5.2)), we believe that it is worth better understanding this relation between the integrals for the energy correlator and Feynman parameter integrals in shifted dimensions. Large polynomials also appear in the expression for the two-point correlator in QCD [12,13] (and to some extent in N = 4 [11,14]), and it would be interesting to understand if they can also be expressed in terms of Feynman parameter integrals in shifted dimensions. It will also be extremely interesting to understand to what extent this relationship extends to higher point correlators, or beyond the collinear limit.
Finally, we wish to emphasize another feature of the integrals appearing in the result, namely that due to the fact that the observable is infrared and collinear safe, the loop integrals appearing in its result must be infrared finite. Such integrals are much simpler than their infrared divergent counterparts, and many powerful techniques exist for their study. In particular, finite integrals involving similar functions to those appearing here obey a number of interesting differential equations (see e.g. [58,[80][81][82]), and it would be interesting to understand if these could be applied in the current case. Either way, any relationship between event shape observables and infrared finite loop integrals is extremely intriguing and worthy of further study.

Result in Terms of Polylogarithms
The N = 4 result can also be written directly in terms of the z variable, and the transcendental functions appearing in Sec. 3. We have The transcendental functions D ± 2 (z) are as defined in Sec. 3.2. The polynomial appearing in front of the Bloch-Wigner function is given by One can easily verify that this result obeys all the constraints discussed in Sec. 3. We will discuss the OPE limit in more detail in Sec. 6.
We see that the complexity of the result resides in the rational prefactors of the different transcendental functions, which are automatically accounted for by the representation in terms of Feynman integrals in shifted dimensions. This illustrates that these standard transcendental functions do not seem to be the optimal basis for expressing the result. Nevertheless, we still find it interesting to understand the structure of the polynomials Eq. (5.2) without reference to the Feynman integral perspective, so we explore their properties a bit more in this section.
We have chosen to write them in an explicit form in terms of the z andz variables at the expense of the S 3 symmetry. Here we show that taking into account this S 3 symmetry, these polynomials are in fact quite minimal. To understand how to write in terms of S 3 invariants, we recall that we defined the function G(z,z) as We can re-write this in a more symmetric way as where we recall from Eq. (3.7) that we defined and This shows that the function that has trivial transformation properties under the S 3 group is However, due to the non-trivial action of the S 3 symmetry, given in Eq. (3.8), S 3 invariant polynomials can take a non-trivial form. First consider the OPE limit. In this limit, we have The completion of this into an S 3 symmetric polynomial is Interestingly, this is exactly the weight 0 term in the full result in Eq. (5.2), once we account for the conversion between F (z,z) and G(z,z), as given in Eq. (5.8).
We also find that the weight 1 terms can be expressed as a sum of two S 3 invariant terms where we have written r 1 = |z| 2 and r 2 = |1 − z| 2 to shorten the expression. In particular, we have which is the weight 1 term in our result for the EEEC, showing that this polynomial is in fact fixed by only two coefficients.

(5.13)
Here we see that we need one S 3 invariant at each order in the A expansion.
We also note that while the result for the EEEC is not of uniform transcendental weight, the weight 0 and weigh 1 terms are minimal, in that they are fixed in terms of the weight 2 terms by the requirement of satisfying the correct behavior in the OPE limit, as well as in the collapsed triangle limit. We believe that this makes it promising that this result could be entirely bootstrapped. However, to completely bootstrap the result would require a better understanding of the polynomials that can appear in the result, in particular, it would be nice to have an a priori prediction for their order.

QCD Jets
Having presented, and discussed in some detail the N = 4 result, which takes a particularly simple form, we now present the results for both quark and gluon jets in QCD. These take a more complicated form, however, this complexity is a feature of the rational prefactors, not the transcendental functions. Despite some effort, we were unable to simplify the structure of the polynomials, which we believe in QCD are inherently complex and encode the structure of power corrections in the QCD splitting functions. Calculations in this section were performed using the program HyperInt [83].
-24 -To simplify the structure of our results, we write the rational functions in terms of the variables r = zz, s = z +z, t = z −z , and we use the following shorthand for the transcendental functions 2 = π 2 .
Note that the functions g 2 , g 2 and g 2 are related by the S 3 permutation symmetry. The results of this section can equally well be written in terms of loop integrals, however, unlike the case of N = 4, this does not provide significant insight, and so we will write the results in terms of the transcendental functions given above.

Discussion
Having presented the results for both QCD and N = 4, we see that there is an enormous difference in the complexity of the rational prefactors. This can also observed by comparing the analytic calculations of the two-point EEC at generic angles for Higgs decays [13] and e + e − annihilation [12] with the N = 4 result [11]. However, the situation is even worse for the EEEC. In particular, in the case of quark and gluon jets in QCD, we find factors of 1/(z −z) 11 , as compared with 1/(z −z) 3 in N = 4. It is interesting to speculate on the reason for this remarkable simplicity in N = 4. There are several features of our results that we find of particular interest in this regard. First, even the pure gluon result in QCD is very complicated, despite the fact that we are working at tree level, so that the underlying amplitude is equivalent to the pure gluon amplitude in N = 4. Second, we can look at the result in N = 1 with an adjoint gluino by setting C A → C F , and n f → C A in the QCD result. In this case, we are also unable to observe any particular form of simplicity in the result. We are therefore led to believe that this simplicity is associated with a particular feature of the splitting functions in N = 4, when occurs when the amplitudes are squared and summed over particle species (i.e. at the cross section as opposed to amplitude level). Indeed, it is well known [84,85] that for the 1 → 2 splitting functions in N = 4, the particular splitting functions for a given particle species are not any simpler than in QCD, but when summed over all final states, all the quantum corrections cancel, and one gets precisely the classical result. To our knowledge, there has been limited study of the summed splitting functions for higher point splittings in N = 4. However, for the 1 → 3 splitting function used for the the calculation of the EEEC, we find a particularly simple result We believe that it would be interesting to understand this in more detail, in particular, if there is some emergent classicality if one is measuring energy type observables, as opposed to considering quantum scattering amplitudes. We also believe that it would be particularly interesting to understand the consequences of the integrability of N = 4 at the level of the cross section, which has received relatively little attention compared to its study at the level of amplitudes. As a final comment on the structure of the results, we note that there does not seem to be a naive form of the principle of maximum transcendentality for the multi-point energy correlators. While the principle of maximal transcendentality holds for the twist-two anomalous dimensions at least through to three loops [86][87][88][89], and therefore controls the scaling of the energy correlators, it does not seem to persist for the non-trivial cross ratio dependence. This is perhaps not surprising since these contributions are in no sense eikonal, but it could be interesting to understand better the extent to which it is broken.

The Squeezed Limit of the Three-Point Correlator
In this section we briefly study the squeezed (or OPE) limit of the collinear three-point correlator, which is illustrated in Fig. 8. For a generic collinear three-point correlator 3 1, the squeezed limit is defined as further taking the limit x 1 x 2 = x 3 = x L + O(x 1 ). One can also take the limit in any other ordering, but the result is the same due to bosonic symmetry of the energy measurement operator. Using the analytical results computed in this paper, one can directly verify that the triple correlator in the collinear squeezed limit scales as where the constant C depends on the theory under consideration. In the language of ANECs, the collinear squeezed limit corresponds to sequential light-ray OPE limits [29]. It can also be easily understood in terms of the momentum space factorization picture for EEC proposed in [26], which is the approach we take here. We start from a generic triple correlator, without imposing the collinear limit. Then the limit of x 1 → 0, x 2 → x 3 is captured by the factorization formula, 2) The jet function j f i (u 2 |z| 2 ) is the two-point energy correlator at an angle |z| 2 , in a jet initiated by a parton with flavor f i , and momentum fraction u = 2p f i · q/q 2 , where q is, e.g., the centerof-mass momentum of e + e − collision, and p f i is the momentum of the jet. It is related to the cumulative jet function for the EEC by To leading logarithmic (LL) accuracy, the dependence on the momentum fraction u in the jet function is irrelevant, and we can simply set it to 1. Then the integral over the hard function can be interpreted as a two-point energy correlator, but with unequal energy weighting on the two points, This function is not collinear safe due to the E 2 i weighting. However, the collinear divergence will be canceled by the boundary contact term in the triple correlator. For our purposes, we only need the result in the limit x 3 → 0. In this limit, the hard function further factorizes into where H f k (v) is the semi-inclusive fragmentation coefficient with momentum fraction v, which also appears in the context of EEC, and j f i /f k (x 3 v 2 ) is the two-point correlator with unequal energy weight for a jet initiated by f k with momentum fraction v. Again, when working at LL accuracy we can set v = 1, and H q (v) = δ(1 − v). 9 At the first order in α s , j f i /f k can be computed schematically by where dPS (2) is the collinear two-body phase space measure, P f i f j ←f k is the splitting function for f k → f i + f j , y is the momentum fraction of f i with respect to f k , and s = ( The difference as compared with the jet function for the EEC is the y 2 weighting, and the tagging of a parton flavor f i . We therefore find that, to LL accuracy in the collinear squeezed limit, the triple energy correlator factorizes as We can now use this formula to compute the squeezed limit for a quark jet, a gluon jet, and a jet in N = 4 SYM. We first recall the EEC jet function. At O(α s ) and LL accuracy, it is simply given by (j q (|z| 2 ), j g (|z| 2 )) = α s 4π (1, 1) The evolution matrix on the right hand side of Eq. (6.8) is the third Mellin moment of puresinglet time-like splitting matrix. For the LL jet function of the two-point correlator with 9 We do not distinguish q andq in e + e − due to the charge conjugation invariance of QCD.
-35 - For a quark jet, we then find find the collinear squeezed limit to be (6.10) For a gluon jet, we find (6.11) These results are in perfect agreement with expanding the the results in Sec. 5.2.1 and 5.2.2. One can also convert the results above to N = 1 SYM, by setting C F = C A = n f = N c . We find which agrees with the expectation of N = 1 supersymmetry. For N = 4 SYM, all the flavor dependence disappears, and one only needs a universal jet function for both EE and E 2 E correlation, We then obtain the collinear squeezed limit in N = 4 SYM, (6.14) Again, this results is in complete agreement with the expansion of Eq. (5.2). It will be interesting to extend this analysis to higher point correlators where there are multiple iterated squeezed limits, which should placed strong constraints on the structure of the result.

Numerical Checks
Particularly in QCD, our result is fairly non-trivial, and also relies on collinear factorization for a multi-particle correlation. Since this is beyond what has traditionally been studied, we -36 -  can use numerical codes for the full e + e − → 4 parton matrix elements to check that our analytic results do indeed reproduce the correct singular behavior for the EEEC observable. This check can be performed separately for each partonic channel, which provides a strong cross check on our analytic calculation. In this section we use both Nlojet++ [31,32] and Event2 [33] to check our results.
We start from numerically verifying the case where the three point span an equilateral triangle, which can be obtained from our analytic calculation as (see e.g. Eq. (5.4) and (5.14)) We begin by verifying the behavior for the color summed result. In Fig. 9 we show the analytic result for the EEEC for an equilateral triangle configuration in the collinear limit, as a function of the side length, x L . This is compared with numerical results from both Nlojet++ and Event2. In both cases, we see agreement between the numerical result and our analytic calculation as x L → 0. This shows that our factorization formula correctly reproduces the singular behavior in the collinear limit.
To test this agreement in more detail, we can decompose the cross section by color structure  Figure 10: A comparison of the analytic results for the EEEC in an equilateral triangle configuration with numerical results from EVENT2, separated by color channel. Results are shown for the C 2 F , C F C A and C F N f T f channels, following the decomposition given in the text. Good agreement is observed for each channel.  Figure 11: The identical quark interference term compared with Event2 for an equilateral triangle configuration. and separately compare our analytic predictions with EVENT2 for each of the color channels C 2 F , C F C A and C F N f T f . This comparison is shown in Fig. 10, where we again see good agreement for all the partonic channels. This provides a strong check on our result.
Another gauge invariant subset of contributions that we can test are the identical quark -38 - (qqqq) interference terms. These terms were studied in detail for the EEC in [12]. They are a subleading color contribution, with color factor C F (C A − 2C F ). Although these terms are numerically small, they can be isolated in Event2 and provide an additional check on our calculation. In Fig. 11, we show a comparison of our analytic result with Event2 for the identical quark interference term, again finding good agreement as x L → 0.
To further examine the shape dependence of the EEEC, we show in Fig. 12 three separate cases, where the three point span acute, obtuse, and right triangle, respectively. In all cases, we find good agreement between our analytical results and Monte Carlo.
The combination of all these different tests provides strong support for the correctness of our calculation, and for the factorization of the EEEC in the triple collinear limit presented in this paper. It would be interesting to extend our calculation to the next order in α s to study the radiative corrections to the EEEC, and compare with the predictions of factorization.

Conclusions and Future Directions
In this paper we have presented a calculation of the three point energy correlator in the collinear limit. This extends recent calculations of the two-point correlators, and is particularly interesting, since the three point correlator is the first time where non-trivial shape dependence appears. Our calculation uncovered a number of interesting properties of the energy correlators, which suggest that they can be computed to higher points and at higher loops in both QCD and N = 4. Perhaps most intriguingly, we found a duality between the integrals involved in the calculation of the energy correlator, and Feynman parameter integrals. This allowed the N = 4 result to be expressed in a single line, in terms of triangle and box integrals. We presented analytic results for N = 4, and for both quark and gluon jets in QCD, and contrasted the results in the two theories. For the case of QCD quark jets, we compared our results with a numerical calculation using the full e + e − → 4 parton matrix elements, finding agreement for all partonic channels in the collinear limit.
Since this is the first calculation of a multipoint energy correlator, and due to the interesting nature of our results, we believe that it opens up several directions for future study, some of which we describe here: All Orders Resummation of Multi-Point Correlators: In this paper we have only considered the fixed order calculation of the non-trivial shape dependence for the triple correlator. At higher perturbative orders, one will encounter logarithms, log(x L ). As for the case of the two particle correlator, these are single collinear logarithms, and it would be interesting to understand their all orders structure for a generic multi-point correlator. This is interesting theoretically, and is also essential for phenomenological applications at the LHC.
Jet Substructure: One of the primary goals of jet substructure is the understanding of the correlations of energy flow with a jet [1]. Much of the interesting information is encoded in multi-point correlators, however, these have until now evaded analytic calculation due to their apparent complexity. Here we have shown that these multiparticle correlations can indeed be analytically calculated, and that the results are much simpler than naively expected. It would be interesting to perform a more phenomenological study of these observables as relevant for the LHC. It would also be interesting to understand better their relation to more standard jet substructure observables, such as the energy correlation functions [90,91] or energy flow polynomials [92]. Indeed, the moments of some of the most well used jet substructure observables, namely D 2 [93][94][95][96] and N 2 [91] are directly related to three point energy correlators. We also believe that the multipoint correlators are interesting as precision probes of QCD at the LHC, and it would be interesting to measure them directly.
Higher Point Correlators and Bootstrap: The simplicity of our results suggest that, particularly for N = 4, the extension to higher point correlators is feasible. An interesting feature of the higher point energy correlator observables, is that due to their dependence on multiple variables, one can study their behavior in a number of different kinematic limits, much like amplitudes. It would be interesting to understand if one can set up a bootstrap program, similar to that which has been successful for scattering amplitudes [97][98][99][100][101][102][103]. This would require a better understanding of the space of functions that describe the energy correlators perturbatively, as well as better constraints on the rational functions appearing in the results.
Relation to Feynman Parameter Integrals: An interesting observation made in this paper was the relation between the calculation of the energy correlators in the collinear limit, and Feynman parameter integrals for loop amplitudes. It would be interesting to understand if this extends to higher point energy correlators, as well as to higher orders in the perturbative expansion. This relation is interesting in its own right, but is also technically extremely useful due to the wealth of knowledge regarding the structure of loop amplitudes, as well as the well developed technology for performing Feynman parameter integrals. We also believe that this relation is promising for giving a geometric interpretation to the calculation of the energy correlators, given that such an interpretation is known for amplitudes (at least to one loop) [64][65][66][67][68][69][70][71].
Full Angular Dependence: It would be interesting to understand whether the calculation performed here in the collinear limit could be extended to derive an analytic result for the full angular dependence of the EEEC, and to understand the class of functions that appear in the result. For the EEC, there exist powerful sum rules which constrain the endpoint behavior of the result, and it would be interesting to understand the sum rules for the EEEC, and also for higher point correlators.
Effective Theory on the Celestial Sphere: More speculatively, it would be interesting to understand if there exists a field theory on the celestial sphere that can be used to compute the energy correlators directly. There has recently been much work in this direction at the amplitude level (see e.g. [104,105]). Much like the Mellin transform used in these studies, the integral over the energy weighting appearing in the definition of the energy correlator observable turns the result into a conformal primary on the sphere. One advantage of the energy correlator observables as compared with the case of amplitudes is that they are necessarily infrared finite at higher orders. By computing higher point functions, one can build up this speculative theory on the celestial sphere, and study in detail the structure of its OPE. The three point function provides the first non-trivial data in this direction.
Consequences of Integrability for Cross Section Level Observables: There has been significant progress in understanding the consequences of integrability in N = 4 for amplitudes and anomalous dimensions [106]. However, there has been relatively little work on understanding the consequences of integrability on infrared safe event (or jet) shape observables, which are the physical observables that are measured in collider experiments. We believe that the energy correlators provide an interesting playground to study this question.
We hope that the simplicity of the results presented here motivate more work towards understanding the structure of the energy correlators, both on the more formal side, as well as in phenomenological applications to jet substructure, and we hope for a fruitful interplay between the two. We believe that there is considerable unexplored simplicity in the substructure of jets that could have a significant impact on analyses at the LHC, and in turn, we hope that the high quality jet data at the LHC can be used to better understand Lorentzian quantum field theories. (A.2)

A.2 N = 4 SYM
We extract the universal 1 → 3 splitting function in N = 4 SYM by taking the three-particle collinear limit of the squared amplitudes for 1 → 4 decay for a source operator trF 2 . We find a particularly simple answer after summing over all flavors, It would be interesting to better understand the simplicity of this result.