Energy correlations in heavy states

We study energy correlations in states created by a heavy operator acting on the vacuum in a conformal field theory. We argue that the energy correlations in such states exhibit two characteristic regimes as functions of the angular separations between the calorimeters: power-like growth at small angles described by the light-ray OPE and slowly varying, or “flat”, function at larger angles. The transition between the two regimes is controlled by the scaling dimension of the heavy operator and the dynamics of the theory. We analyze this phenomenon in detail in the planar N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{N} $$\end{document} = 4 SYM theory both at weak and strong coupling. An analogous transition was previously observed in QCD in the measurement of the angular energy distribution of particles belonging to the same energetic jet. In that case it corresponds to the transition from the light-ray OPE, perturbative regime described in terms of correlations between quarks and gluons to the flat, non-perturbative regime described in terms of correlations between hadrons.

1 Introduction and summary The energy correlations are among the best studied observables both experimentally and theoretically [1]. They measure the flux of energy deposited in calorimeters located at different points on the celestial sphere and carry information about the dynamics of the underlying theory. A lot of activity has recently been devoted to studying the energy correlations in QCD and in the maximally supersymmetric N = 4 SYM theory. The latter serves as a very useful toy model that one can use to develop new techniques for computing these observables in QCD. Moreover, it allows us to understand certain properties of the energy correlations that cannot be explained within the conventional perturbative QCD approach.
As an example, consider the energy-energy correlation (EEC) measuring the angular distribution of the energy of the particles that enter into two calorimeters separated by the relative angle 0 ≤ θ ≤ π [2,3]. At large total energy Q, the hadronic final states consist of collimated beams of energetic particles, or jets. In the energy correlations, jets manifest themselves as peaks located at small angles θ, as well as at finite θ corresponding to the angular separation between the jets.
For a small angle θ, the EEC describes the correlation between particles belonging to the same jet. The analysis of the experimental data shows [4] that it behaves differently at small angles, depending on how θ compares with the non-perturbative parameter θ 0 ≃ Λ QCD /Q given by the ratio of the QCD hadronization scale and the total energy, 1 Flat : EEC ∼ const, OPE : EEC ∼ 1/θ 2−γ .
(1.1) angles is determined by their relative transverse momenta Qθ. For small transverse momenta Qθ = O(Λ QCD ) the theory becomes strongly coupled. Conversely, the perturbative approach is justified for θ ≫ Λ QCD /Q. The behavior EEC ∼ θ 0 can be reproduced if one thinks about the final state as consisting of a dense cloud of hadrons weakly interacting with one another.
Describing the transition to this regime requires control over the non-perturbative QCD. A similar change of behavior of the energy-energy correlation at small angles also takes place in the N = 4 SYM theory as the 't Hooft coupling varies, but the underlying mechanism is slightly different. This theory is conformal and it looks alike at short and large distances. At weak coupling, the EEC in N = 4 SYM has the same power-like behavior at small angles as in QCD. This behavior changes as one goes to the limit of strong 't Hooft coupling constant λ = g 2 YM N c ≫ 1. Increasing the value of the coupling constant, one enhances the production of particles (gauge, gaugino and scalars) in the final state. At strong coupling, the final state in N = 4 SYM consists of an infinite number of soft particles whose energy is distributed homogeneously on the celestial sphere. As a consequence, for λ → ∞ the energy-energy correlation does not depend on the angle, EEC(θ) = 1 + O(1/λ). As in the case of QCD, the change of behavior of the EEC as a function of the angle is associated with the presence of a large number of particles in the final state. To describe the transition in detail, one needs to know the energy correlation in N = 4 SYM for an arbitrary 't Hooft coupling, which is out of reach at present.
In this paper, we point out that there exists another physical mechanism of transitioning between the two characteristic regimes (1.1). For simplicity we restrict our attention to conformal field theory (CFT), but the basic mechanism should be applicable to theories with dimensionful scales, such as generic four-dimensional gauge theories including QCD. Because the transition is driven by the large number of particles produced in the final state, we can create such a final state by exciting the vacuum with a "heavy" operator O H (x) carrying a large scaling dimension ∆ H ≫ 1. Physical intuition suggests that for ∆ H → ∞ the state |H(q)⟩ = d 4 x e iqx O H (x)|0⟩ should contain arbitrarily many soft particles which are largely uncorrelated. As a consequence, for ∆ H → ∞ the energy correlations are expected to be angle-independent, up to corrections suppressed by a power of 1/∆ H .
In order to formulate this property, it is convenient to introduce the energy flow operator E(n) which measures the energy flux in the direction specified by a null vector n µ = (1, ⃗ n) with ⃗ n 2 = 1 [12][13][14]. In the rest frame of the source, for q µ = (Q, ⃗ 0), its expectation value in the state |H(q)⟩ is determined by the total momentum, ⟨E(n)⟩ H = Q/(4π). It does not depend on the choice of the source operator and yields the total energy after integration over the celestial sphere, d 2 ⃗ n ⟨E(n)⟩ H = Q. where E(n) = E(n)/⟨E(n)⟩ H is the normalized energy flow operator with ⟨ E(n)⟩ H = 1. To study their behavior in the limit ∆ H → ∞, we find it convenient to apply the cumulant expansion ⟨ E(n 1 ) . . . E(n k )⟩ H = 1 + i<j ⟨ E(n i ) E(n j )⟩ c + i<j<m ⟨ E(n i ) E(n j ) E(n m )⟩ c + . . . , (1.3) where ⟨. . .⟩ c denotes the connected correlation. Here the terms ⟨ E(n i ) E(n j )⟩ c depend on the single relative angle between the ith and jth detectors; ⟨ E(n i ) E(n j ) E(n m )⟩ c depend on the three relative angles between the ith, jth, and kth detectors, etc. Starting from four detectors, for k ≥ 4, the expansion (1.3) involves non-linear terms, the simplest being ⟨ E(n i ) E(n j )⟩ c ⟨ E(n m ) E(n q )⟩ c , see e.g. [15] and Appendix A. We conjecture that when the source becomes heavy, ∆ H ≫ 1, each subsequent term in (1.3) is further suppressed compared to the previous one, namely The values h k and the precise form of ⟨ E(n 1 ) . . . E(n k )⟩ c depend on the theory and the physical state in question. Notice that in (1.4), when taking ∆ H → ∞, we keep the angles between the detectors, as well as the number of detectors, fixed. The same property was observed before in strongly coupled gauge theories with gravity duals [5]. As mentioned above, in this case the multi-particle final state is not generated by considering a heavy source, but by the strongly coupled dynamics of the theory. The role of ∆ H is played by the 't Hooft coupling λ, and h k = k 2 . In contrast, we would like to emphasize that the relation (1.4) holds in the limit of large ∆ H for an arbitrary coupling (including the free theory!).
In this paper we study the energy correlations (1.2) in planar N = 4 SYM. To define the heavy state |H(q)⟩, we choose the operator O H (x) to be a half-BPS scalar operator of the form O H (x) = tr[ϕ K (x)]. Its scaling dimension is protected from quantum corrections, ∆ H = K, and the heavy state limit corresponds to K → ∞. Another advantage of this choice is that the two-point correlation ⟨ E(n 1 ) E(n 2 )⟩ c , defining the leading angular-dependent contribution in (1.3), can be computed explicitly both at weak and strong coupling for arbitrary K ≥ 2. This computation is done most efficiently by the method based on correlation functions and Mellin transforms, developed in [5,[16][17][18][19][20].
Using the explicit expression for ⟨ E(n 1 ) E(n 2 )⟩ c , we can verify that its dependence on the relative angle cos θ = (⃗ n 1 ⃗ n 2 ) is a function of the scaling dimension of the source K. We find that for K = 2, the energy-energy correlation at weak coupling in N = 4 SYM has a shape very similar to that in perturbative QCD. Namely, it is peaked around the end points, θ = 0 and θ = π, and is flat in between. This shape corresponds to a final state containing two jets, one per each scalar field in the definition of the source operator tr[ϕ 2 (x)]. As K increases, we observe that the peak at θ = π disappears and the function flattens out for 0 < θ < π. Moreover, ⟨ E(n 1 ) E(n 2 )⟩ c vanishes as 1/K in the limit K → ∞. To the lowest order in the coupling we find, up to corrections suppressed by powers of λ and 1/K, ⟨ E(n 1 ) E(n 2 )⟩ c = λ 8π 2 K 3 z + + 2Li 2 (1 − z) − 6 log z − 2ζ 2 − 13 2 + 5 2 δ(z) + . . . , (1.5) where z = (1 − (⃗ n 1 ⃗ n 2 ))/2 = sin 2 (θ/2) in the rest frame of the source and 1/z + is the plusdistribution. The result (1.5) is in agreement with our expectation (1.4) for ∆ H = K and h k = k − 1.
The relation (1.5) holds at weak coupling for K ≫ 1. The situation changes at strong coupling for λ ≫ K ≫ 1. In this limit, we find that the energy-energy correlation does not depend on K and it vanishes as 1/λ, (1.6) The relations (1.5) and (1.6) illustrate two different mechanisms of producing slowly-varying energy distributions. At weak coupling, the leading correction is controlled by the scaling dimension of the source. At strong coupling, it comes from the evolution of the state and is controlled by the coupling constant.
Below we show that at large K and finite coupling constant the energy-energy correlation in N = 4 SYM has the following behaviour at small angles z = sin 2 (θ/2) ≪ 1 where C and γ depend on the 't Hooft coupling. In the case of QCD we expect that the EEC should have the same behavior in the OPE region (1.1), for θ > θ 0 . An important difference compared to N = 4 SYM is that the "kinematical" large parameter K in (1.7) is replaced in QCD by a dynamical one defined by the multiplicity of the produced hadrons within a jet. Above we focused our discussion on energy correlations in QCD and its conformal cousin N = 4 SYM. It is interesting to ask what happens for event shapes, or matrix elements of light-ray operators, in a general interacting CFT when the source becomes heavy. 3 We present some evidence that the clustering structure (1.4) is valid for energy correlations in any CFT. Let us briefly explain the reason for that. According to the state-operator correspondence, a local operator O H (x) with scaling dimension ∆ H can be associated with an energy eigenstate E cyl = ∆ H /R in a CFT defined on a cylinder R × S d−1 , where R is the radius of the sphere. A heavy operator with large scaling dimension ∆ H ≫ 1 corresponds to a highly excited energy eigenstate. In an interacting theory it is expected to look thermal when simple enough observables are considered [22] (the statement known as the Eigenstate Thermalization Hypothesis). The clustering structure of energy correlations (1.4) in a heavy state is then related to the clustering of correlation functions of stress-energy tensors in the thermal state as the separation between the operators becomes large [23].
The paper is organized as follows. In Section 2 we compute the multi-point energy correlations ⟨E(n 1 ) . . . E(n k )⟩ H in a free theory and show that they satisfy (1.4) in the limit of large scaling dimension of the source. In Section 3 we consider the energy-energy correlation ⟨E(n 1 )E(n 2 )⟩ H in planar N = 4 SYM at weak coupling. We choose the source operator to be the simplest half-BPS operator built out of K scalar fields and examine the dependence of the energy correlation on K at the leading order in the coupling constant (Born level). In Section 4 we compute the one-and two-loop perturbative corrections to the energy-energy correlation for arbitrary K and show that they verify the relation (1.4) at large K. In Section 5 we present an approach that allows us to compute the O(1/K) correction to (1.4) directly, without going through the details at finite K. In Section 6 we obtain the contact terms necessary to describe the singular behavior of the energy correlations at the end points, for finite K and for K → ∞. In Section 7 we compute the same observable at strong coupling including the first stringy correction. In Section 8 we present arguments in favor of (1.4) in a generic interacting CFT. The paper contains several appendices and ancillary files.

Energy correlations in the free scalar theory
The energy correlations for a heavy source are expected to have the general form (1.4). In this section, we show that the relations (1.4) are indeed satisfied in the free scalar theory.
We consider a free massless scalar field ϕ(x) and choose the source operator to be O K (x) = ϕ K (x). It creates the state |H(q)⟩ = d 4 x e iqx O K (x)|0⟩ containing K massless scalar particles with the total momentum q µ (with q 0 > 0 and q 2 > 0). The differential cross-section describing the probability to find K particles in a final state with on-shell momenta p i (with p 2 i = 0) is given by where dLIPS(p) is the Lorentz invariant phase space measure of a particle with momentum p, with δ + (p 2 ) = δ(p 2 )θ(p 0 ). The total cross-section is given by In this representation, the phase space integral factorizes into a product of two-point Wightman functions where the '+i0x 0 ' prescription ensures that the Fourier transform is different from zero for p 0 > 0. The calculation shows that .
Notice that σ K (q) grows as a power of q 2 with the exponent being a linear function of the weight (or scaling dimension) K of the operator O K . Let us now examine the k−point energy correlation in the K−particle final state described by the differential distribution (2.1). It is given by where the weight factor w n (p) selects the particle in the final state moving along the null direction n µ = (1, ⃗ n) (with ⃗ n 2 = 1), The combinatorial factor K! /(K − k)! in the numerator of (2.6) is due to the Bose symmetry of the scalar particles. It counts the total number of events in which k out of K particles enter the calorimeters. The diagrammatic representation of the relation (2.6) is shown in Figure 1. Substituting (2.1) into (2.6) we can express ⟨E(n 1 ) . . . E(n k )⟩ as an integral over the phase space of the K particles. The integral over the undetected (K − k) particles gives rise to the total cross-section σ K−k (q) whereq µ = q µ − k i=1 p µ i is the total momentum of these particles. Then, ⟨E(n 1 ) . . . E(n k )⟩ is given by the integral over the phase space of the detected k particles, Using (2.7), this expression can be written as a k−fold integral over the energies w i = p 0 i of the detected particles, The integral (2.9) involves the ratio of the invariant mass of the undetected particles and the total energy,q Notice thatq 2 /q 2 < 1 and, therefore, for large K the dominant contribution to the integral (2.9) comes from the integration over the regionq 2 /q 2 = 1 + O(1/K), or equivalently ω i = O(1/K). Changing variables as we find in the limit K → ∞ with k held fixed where z ij are dimensionless variables . (2.13) In the rest frame of the source, for q µ = (Q, ⃗ 0) and n i = (1, ⃗ n i ), these variables are related to the relative angles between the detectors on the sphere, z ij = (1 − cos θ ij )/2. Combining together (2.9) and (2.12) we find that the energy correlations in the large K limit take a remarkably simple form, (2.14) The integral in this relation describes an ensemble of k noninteracting particles whose energies ω = εQ/(2K) are distributed according to the law dP (ε) = dε ε 2 e −ε . The first term in the brackets in (2.14) is independent of the angular variables z ij defined in (2.13). It describes a homogeneous distribution of the energy on the celestial sphere. The angular dependence comes from the second term, which is suppressed by a factor of 1/K and involves a quadratic polynomial in the energy variables ε i . To higher orders in 1/K, the coefficient of 1/K p is given by a polynomial in ε i of degree 2p. As we show below, it contributes to the connected part of the p−point correlations in (1.3). The calculation of the integral in (2.14) yields Switching to the normalized operators E(n i ) on the left-hand side of (2.15), the term of order 1/K on the right-hand side is identified with the sum of the two-point connected correlations It is straightforward to compute the subleading terms in (2.15). Bringing them to the form (1.3), we can determine the higher-point connected correlations, e.g.
The relations (2.15), (2.16) and (2.17) are valid for n i ̸ = n j and they do not take into account the contribution of contact terms localized at n i = n j . The presence of such terms can be detected as follows. By definition, ⟨E(n 1 ) . . . E(n k )⟩ measures the energy flux on the celestial sphere in the directions specified by the unit vectors ⃗ n 1 , . . . , ⃗ n k . The integral of E(n)n µ over the unit sphere ⃗ n 2 = 1 yields the total momentum of the final state. As a consequence, the energy correlation has to satisfy the sum rule One can check that the leading O(K 0 ) term in (2.15) verifies this relation but this is not the case for the O(1/K) term. The reason for this is that the integral in (2.18) receives contributions from contact terms localized at coincident directions n i . A contact term proportional to the product of L delta functions on the sphere, δ(n 1 , n 2 ) . . . δ(n 1 , n L+1 ), comes with the combinatorial factor K k / K k−L ∼ 1/K L . To the leading order in 1/K, only the contact terms with L = 1 contribute. They are proportional to 1≤i<j≤k δ (2) (⃗ n i , ⃗ n j In distinction to (2.18), it is not sensitive to contact terms. We have checked that ⟨ E(n i ) E(n j ) E(n m )⟩ c in (2.17) satisfies (2.19). We conclude this section by reiterating our main conjecture (1.4): The connected correlations in the free theory have the expected heavy weight behavior with ∆ H = K and h k = k − 1 (with k = 2, 3, . . . ).

Energy correlations in N = 4 SYM
In the previous section, we demonstrated that the energy correlation in a final state consisting of a large number K of free scalar particles is given by simple expressions like (2.16) and (2.17). The question arises, how does the interaction between the particles in the final state affect this result?

Energy-energy correlation
In order to address this question, we consider the energy correlations in a particular fourdimensional gauge theory -the maximally supersymmetric N = 4 Yang-Mills theory with gauge group SU (N c ). For the sake of simplicity we shall concentrate on the two-point correlation ⟨E(n 1 )E(n 2 )⟩. To create a final state with a large number of particles produced, we excite the vacuum by acting on it with a half-BPS operator of the form It is built of six real scalar fields X I with Y I being an auxiliary null complex vector of SO(6) defining the orientation of O K (x) in the isotopic R−space. This operator has the R−charge of the [0, K, 0] representation of SO(6) ∼ SU (4), as well as conformal weight (or scaling dimension) K, the latter being protected from quantum corrections by N = 4 superconformal symmetry. It creates K scalar particles out of the vacuum. Going to the limit K → ∞, we encounter the final state discussed in the previous section. The two-point correlation function of the operators (3.1) is protected from quantum corrections and is given by the product of K free scalar propagators (2.4) multiplied by an SU (N c ) color factor. In the planar limit we have Notice that the operators are not time ordered. To simplify the formulae, we put (YȲ ) = 1 in what follows. Then, the total cross-section for producing an arbitrary number of particles in the final state is given by (see (A.4) in [16]) where σ K (q) is defined in (2.5). Like (3.2), it is protected from quantum corrections. The relation (3.3) expresses the total cross-section as the Fourier integral of the two-point Wightman correlation function (3.2). It is equivalent to (2.3) after inserting the sum over the final states between the operators. In a similar manner, the energy correlations (2.6) also admit a representation in terms of correlation functions involving two scalar operators (3.2), as well as the energy flow operators E(n i ) (for the definition see Appendix A). The latter play the role of calorimeters detecting the particles in the final state.
For instance, the two-point energy correlation is given by Lorentz symmetry allows us to fix the form of the correlation, up to an arbitrary function F K (z) of the Lorentz covariant angular separation z ≡ z 12 between the detectors (recall (2.13)), Here θ 12 is the angle between the unit vectors ⃗ n 1 and ⃗ n 2 in the rest frame of the source, q µ = (Q, ⃗ 0). The kinematical factor on the right-hand side of (3.5) has Lorentz weight (−3) under the independent rescaling of the detector vectors n µ i , matching the corresponding weight of the energy operators E(n i ) on the left-hand side (see (A.1)-(A.3)). The factor (q 2 ) 4 in the numerator corresponds to the scaling dimension (+1) of E(n i ).
For K = 2 the energy correlation (3.5) has been studied both at weak and strong coupling [5,[16][17][18][19][20][24][25][26][27]. Our goal in this paper is to extend these results to K ≥ 3 and to understand the properties of the energy correlation in the limit of large K. We show below that at large K the function F K (z) satisfies (3.7) Substituting this relation in (3.5) we reproduce (1.4) for ∆ H = K and h 2 = 1.
In distinction with the total cross-section, the function F K (z) is not protected in N = 4 SYM and depends on the coupling constant g 2 YM , as well as on the rank of the gauge group N c . In the planar limit, for N c → ∞ with the 't Hooft coupling constant λ = g 2 YM N c held fixed, it admits a weak coupling expansion, where F In the Born approximation, the energy correlation (3.5) is given by (2.9) evaluated for k = 2. Matching (3.5) with (2.15) we find that in the limit of large K the scaling function F (0) K (z) is given by As was explained in the previous section, this relation describes an ensemble of K noninteracting particles whose energies are distributed according to the law (2.14). Turning on the interaction between these particles, one generates additional correlations of their energies. They are described by the functions F (ℓ) K (z) in (3.8) which we compute below for ℓ = 1, 2. We will show that at large K the loop corrections scale as F (ℓ)

Sum rules
The function F(z) satisfies nontrivial conditions that follow from the sum rules (2.18).
We recall that the integral of the energy flow operator over the celestial sphere yields the total energy-momentum of the final state, d 2 n n µ ⟨E(n)⟩ K = q µ . Combined with (2.18), this yields the following relations in the rest frame of the source, (3.10) Substituting the general expression (3.5) of ⟨E(n 1 )E(n 2 )⟩ K and using (3.6), we find that F K (z) has to satisfy the sum rules [6-8] These relations hold for an arbitrary coupling constant λ.
Replacing F K (z) with its perturbative expansion (3.8) and matching the coefficients of the powers of λ on both sides of (3.11) we find that the Born approximation F (3.12) The sum rules for the loop corrections F (ℓ) It is straightforward to verify that (3.9) satisfies the sum rules (3.12). Notice that the function accompanying 1/K in (3.9) gives zero contribution to the sum rules (3.12).

Born contribution
In this subsection, we compute the Born contribution to (3.8) at finite K and study the transition from K = 2 to the large K behavior (3.9).
In the previous section, we used the conventional amplitude approach to obtain the integral representation (2.9) for this contribution. Let us show how the same result follows from the representation (3.4) of the energy correlations in terms of correlation functions.
The energy flow operator E(n) is given by an integral involving the stress-energy tensor of N = 4 SYM, see (A.1) in Appendix A. The latter is bilinear in the scalar fields X I constituting the source operators (3.1). As a consequence, the calculation of the four-point correlation function on the right-hand side of (3.4) in the Born approximation reduces to performing Wick contractions between the K scalar fields in the operators O K (x) (source) andŌ K (0) (sink) and the two pairs of scalar fields in the energy flow operators E(n 1 ) and E(n 2 ). 4 In this way, in the planar limit the correlation function ⟨O K (x)E(n 1 )E(n 2 )Ō K (0)⟩ (0) can be factorized into a product of (K − 2) propagators connecting the source and sink and the same correlation function evaluated for K = 2, where the scalar propagator D(x) is given by (2.4) and the superscript '(0)' indicates the Born approximation. The relation (3.14) is represented diagrammatically in Figure 2. The combinatorial K−dependent factor in (3.14) counts the number of contributing diagrams. Substituting (3.14) into (3.4) we find that ⟨E(n 1 )E(n 2 )⟩ K is the convolution of the Fourier transform of D K−2 (x) and ⟨E(n 1 )E(n 2 )⟩ (0) K=2 . The former is given by the function σ K−2 defined in (2.5) and the latter is given in (3.22) below. This leads to  where n 1 τ 1 and n 2 τ 2 are the momenta of the particles entering the calorimeters. Replacing σ K−2 with its expression (2.5) and changing the integration variables as τ i = ω i q 2 /(2(qn i )), we find that (3.15) takes the expected form (3.5) with where the integration is restricted to the region 1 − ω 1 − ω 2 + zω 1 ω 2 ≥ 0. The result can be expressed in terms of a hypergeometric function, This relation is valid for 0 < z < 1. The values of z = 0 and z = 1 correspond to the situation where the two detectors are located, respectively, atop of each other or antipodal to each other. In general, F K (z) contains contact terms proportional to δ(z) and δ(1 − z). We discuss them below in Section 3.4.
The function F Figure 3 for several values of K. For z → 0 it approaches a finite value, For z → 1 it grows logarithmically for K = 3 and stays finite for K ≥ 4, (3.20) In the large K limit we get instead in agreement with (3.9). We repeat that the z−dependence only appears at the level of the O(1/K) corrections. We observe that the function F (0) K (z) flattens out as K increases and becomes constant for K → ∞. This property is in agreement with our expectation (1.4) that the correlations between the particles in the final state are suppressed at large K.

Contact terms for K = 2
As was mentioned above, the relation (3.17) is valid for 0 < z < 1. One of the reasons for this is that F (0) K (z) is expected to contain contact terms localized at z = 0 and z = 1. The easiest way to reveal their presence is to apply the sum rules (3.12).
For K = 2 the naive answer F (0) K=2 (z) = 0 in (3.18) does not satisfy the sum rules (3.12). In reality, in this case the final state consists of two particles that move back-to-back in the rest frame of the source. As a consequence, the energy-energy correlation is different from zero only if the two calorimeters are atop of each other on the celestial sphere (z = 0) or antipodal (z = 1). This means that it is given by the sum of two contact terms, F with coefficients that can be determined from the sum rules (3.12), Substituting this relation into (3.5), one reproduces the known expression for the energy correlation ⟨E(n 1 )E(n 2 )⟩ K=2 in the Born approximation.
where a small parameter ϵ > 0 was introduced to exclude the contribution of the end-points z = 0 and z = 1. In close analogy with the case K = 2, the latter is expected to be given by a sum of contact terms δ(z) and δ(1 − z). Comparing (3.23) with (3.12), we observe that the contact terms should contribute zero to the first sum rule in (3.12) and 3/(K + 1) to the second sum rule. Put together, these conditions fix the coefficients of the contact terms. The resulting expression for the function F As compared with (3.17), it is valid for 0 ≤ z ≤ 1. At large K the contact term in (3.24) agrees with (3.9). Combined with (3.5), the relation (3.24) yields the expression for the energy-energy correlation in planar N = 4 SYM, to the lowest order in the coupling constant. For finite K, the function (3.24) is peaked around the end-point z = 1 indicating that, like in the two-jet final state, the most of the energy in the final state is deposited at the calorimeters located at two antipodal points on the celestial sphere. Going to the limit K → ∞ we find that, in agreement with (1.4), the energy-energy correlation ceases to depend on the angular separation of the calorimeters. This means that the jets disappear and the energy is homogeneously distributed on the celestial sphere.

Energy correlations at weak coupling
In this section, we study the corrections to the energy correlation (3.5) in the planar N = 4 SYM theory due to the interaction between the particles in the final state. At weak coupling, the Feynman diagrams contributing to (3.5) can be obtained from those shown in Figure 2 by adding interaction vertices. In the planar limit, the interaction can only occur between adjacent scalar particles. Depending on the choice of these particles we can distinguish three cases of interaction between: (i) undetected particles, (ii) one detected particle and undetected ones, and (iii) two detected and undetected particles. The first and the second class of diagrams constitute the total cross-section (2.5) and the single energy correlation ⟨E(n)⟩. Because both quantities are protected from quantum corrections in N = 4 SYM, the above mentioned diagrams do not contribute to the unprotected energy correlations (3.5). We are therefore left with the Feynman diagrams in which two detected particles interact among themselves and with other undetected particles. To the first few orders in the 't Hooft coupling constant, examples of such diagrams are shown in Figure 4.
As follows from the form of these diagrams, at order O(λ ℓ ) the interaction can affect (ℓ + 1) particles at most. At order O(λ) the two detected particles can only interact with each other. At order O(λ 2 ) there is an additional possibility for them to interact with one undetected particle (spectator). At order O(λ ℓ ) the number of spectators cannot exceed (ℓ − 1).
For ℓ < K the number of interacting particles ℓ + 1 is smaller or equal to the number K of particles produced (this is equivalent to the absence of wrapping corrections to the energy correlation (3.4)). The remaining (K − ℓ − 1) particles propagate freely from the source to the sink. In close analogy with (3.14), their contribution to the correlation function ⟨O K (x)E(n 1 )E(n 2 )Ō K (0)⟩ is just a power of the free scalar propagator, The combinatorial factor on the right-hand side gives the ratio of the number of Wick contractions of scalars in the two correlators in the planar limit. We would like to emphasize that the relation (4.1) is only valid for ℓ ≤ K − 1. In particular, in the limit K → ∞ it holds to any given order of the weak coupling expansion.
At one and two loops the relation (4.1) reads Notice that the first relation holds for K ≥ 2 and the second one for K ≥ 3.
Comparing the last two relations with (3.14) we observe an important difference. At large K the expression on the right-hand side of (4.2) and (4.3) is suppressed by a factor of 1/K as compared with (3.14). The reason for this is that the leading O(K 3 ) contribution to the Born level result (3.14) comes from the diagrams shown in Figure 2 on the right panel. These diagrams do not contribute to (4.1) in the planar limit as explained above. In application to the energy correlation (3.4), (3.5) this implies that the loop corrections to the scaling function F K (z) are suppressed by a factor of 1/K as compared to the Born contribution (3.9).

Recurrence relation for the perturbative corrections
The relations (4.2) and (4.3) can be used to compute the correction to the energy correlation (3.8) at one and two loops. Let us introduce the auxiliary functions where ℓ = 1, 2 and the scalar propagator D(x) is given by (2.4). They satisfy the differential equation Below we show that it yields a recurrence relation for the functions F (ℓ) K (z) (for ℓ fixed) that allows us to efficiently determine the loop corrections to the energy correlation (3.8).
We combine equations (3.14) and (4.1) together with (3.4) and (3.5) to get , where σ K (q 2 ) is given by (2.5). Next, the relation (4.5) leads to differential equations for the functions F (0) where ℓ ≥ 1 and D is the second-order differential operator The relations (4.7) and (4.8) are valid for K ≥ 3 and K ≥ 2 + ℓ, respectively. One can check that the Born approximation result (3.17) verifies (4.7). In fact, we could have obtained (3.17) (up to an overall normalization) by solving the recurrence relation (4.7) starting from F

One-loop correction
In this subsection we show that the recurrence relation (4.8) can be effectively used to compute the energy correlation function F that is valid for 0 < z < 1.

Boundary conditions
The solutions to the second-order differential equations (4.8) are defined up to the zero modes of the operator (D − K(K − 1)), This freedom can be fixed by imposing boundary conditions. Notice that the last term on the right-hand side of (4.11) grows as z −K−3 for z → 0, thus producing a divergent contribution to the sum rules (3.13). To avoid it, we put c 2 = 0. According to (4.10), the function F K=2 (z) behaves as O(1/z) at small z. Examining the differential equation (4.8) for K = 3, 4, . . . , one finds that its solution has to have the following asymptotic behavior at small z, For K = 2 we find from (4.10) that a 2 = 1 and b 2 = 0.
To fix the coefficient c 1 in (4.11) we examine the behavior of the function F (1) K (z) around z = 1. In the Born approximation, we deduce from (3.20) that F (0) K (z) grows logarithmically for K = 3 and it approaches a finite value for K ≥ 4. Requiring F (1) K (z) to stay finite for K ≥ 4 as z → 1 implies c 1 = 0. This condition can be derived from the sum rules (3.13) as follows. Multiplying (4.8) by z and integrating by parts, we obtain (4.13) The boundary term at z = 0 vanishes due to (4.12). This derivation relies on the sum rules (3.13), i.e. we assume that for any K ≥ 3 the relation As we show below, it is satisfied for K ≥ 4 only. We therefore deduce from (4.13) that F (1) K (z) should stay finite at z → 1 and K ≥ 4. This leads to c 1 = 0 because the term proportional to c 1 in (4.11) grows logarithmically as z → 1.
For K = 3 one can check using (4.10) that In this case, in order to satisfy the sum rule (3.13), one has to add to (4.10) a contact term localized at z = 1. As a result, for K = 3 the right-hand side of (4.13) scales as 3 2 log 2 ϵ leading to F (1) The term proportional to c 1 in (4.11) yields a subleading contribution in this limit. To fix the coefficient c 1 we can use the sum rule (3.13). 7 We can now turn to solving the differential equation (4.8) supplemented with the boundary conditions specified above. Our strategy is as follows: we first solve explicitly the differential equation for F   Solution for K = 3 Solving the differential equation (4.8) for K = 3, we replace F (1) K=2 (z) with its expression (4.10) and pick the solution that satisfies (4.12). The resulting expression for F (1) K=3 (z) is given by a linear combination of special functions of different transcendental weights accompanied by seven polynomials of z, (4.14) Here the function L(z) contains only terms of the maximal logarithmic weight 3: 8 For K > 3 the sum rule (3.13) leads to c 1 = 0. 8 The terms in the second line are not included in the polynomials c [3] 5 and c [3] 7 in order to keep their coefficients rational.
The coefficients c [3] m (z) (with m = 1, . . . , 7) are polynomials of degree 2: The coefficients c [3] 5 and c [3] 7 depend linearly on the zero-mode coefficient c 1 in (4.11). As explained above, its value is fixed by the sum rule Figure 5. For z → 0 it has the asymptotic behavior in agreement with (4.12). For z → 1 we find As expected, this function grows as a power of log(1 − z).
It is straightforward to verify that the function (4.14) does not satisfy the second sum rule in (3.13), It is needed to regularize the contribution of the pole in (4.17) to . More details about the contact terms can be found in Section 6.

General solution for K ≥ 4
The solution of the recurrence relation (4.8) for any K ≥ 4 takes the same form as (4.14), where c m (z) (with m = 1, . . . , 7) are polynomials of degree K − 1. They are determined recursively by solving (4.8) and requiring F (1) K (z) to be finite for z → 1 and to satisfy (4.12) for z → 0. In particular, the polynomial c (4.20) Instead of presenting explicit formulas for the remaining polynomials at each weight K, in the next subsection we provide a generating function for all F Figure 5 for several values of the weight K. For small z they behave as for several values of the weight K. The limiting curve for K → ∞ is labeled 'asympt', see (4.27). Owing to the normalization prefactor, all plots originate from 3 at z = 0, see (4.21). For z → 1, the one-loop correction has a pole at K = 2, a logarithmic singularity at K = 3 and it is finite for K ≥ 4, see (4.10), (4.18) and (4.22). Its value at z = 1 decreases with K and approaches − 7 2 − π 2 3 ≈ −6.7899 in the limit K → ∞, see (4.29).

Generating function
It is convenient to introduce an auxiliary variable t and combine the one-loop functions F (4.23) Replacing F (1) K (z) with its expression (4.19), we find that each of the seven polynomials c Then the generating function (4.23) takes the form We have found closed-form expressions for the functions G m (z, t). They are given by linear combinations of classical polylogarithms of weights up to three, decorated with rational terms. More precisely, G 1 , which resulted from the summation of the polynomials (4.20), is a rational function; G 2 , G 3 , G 4 contain logarithms and rational terms; G 5 , G 6 contain dilogarithms, logarithms and rational terms; and G 7 is the most complicated function containing tri-logarithms, dilogarithms, logarithms and rational terms. The arguments of the polylogarithms depend both on z and t in such a way that the generating function G(z, t) is given by 2dHPL functions [28].
We provide an explicit expression for G(z, t) in the ancillary file. We can apply the generating function (4.25) to show that F Combining this relation with (4.23), we expect that G(z, t) should scale as − log(1 − t)φ (1) (z) for t → 1. Indeed, examining the generating function G(z, t) for t → 1 we reproduce the expected behavior and identify the function We therefore conclude that, as announced earlier, the one-loop correction to the energy correlation (3.8) is suppressed at large K by a factor of 1/K, as compared with the Born contribution (3.9). For z → 0 the relation (4.26) takes the expected form (4.12) with It also agrees with (4.21). For z → 1 we have, in agreement with (4.22), In the previous relations we keep only the leading term at large K.

Two-loop correction
The two-loop corrections F (2) K (z) satisfy the recurrence differential equations (4.8). These equations are valid for K ≥ 4 and they allow us to express F K=3 (z). The latter can be determined by extending the relation (4.8) to the case K = 3 with a properly defined inhomogeneous term F (2) aux (z).

Two-loop solution of the recurrence relation
Below we show that F Here a (w) i,K (z) and b K (z) are polynomials in z of degree (K − 1) and (K − 2), respectively, with rational K-dependent coefficients. The second line in (4.30) is relevant only for K = 3.
The functions L (w) i ( √ z) are pure polylogarithms of transcendental weight w. They are given by multi-linear combinations (with rational coefficients) of harmonic polylogarithms (HPL) [29] with argument √ z and zeta-values ζ 2 , ζ 3 , ζ 5 , which are graded by the transcendental weight.
The functions L  • Weight-one functions Here we omit the arguments of the HPLs for the sake of brevity, H a,b,... ≡ H a,b,... ( √ z) [29].
The remaining higher-weight combinations L (w) i have a similar form. They are defined in the ancillary file.
The polylogarithmic combinations of HPLs of transcendental weights up to four, L (w) i with 0 ≤ w ≤ 4, can be expressed as classical polylogarithms. However, this does not apply to the weight-five combinations L (5) i , so we prefer to present all weights using the same HPL notations.

Solution procedure
In this subsection we present some details of the derivation of (4.30). Namely, we first compute F 1 + 2zL To find the function F K=3 , we have to solve the differential equation (4.8) for K = 3. This may seem contradictory, since earlier we declared that at two loops equation (4.8) is valid only for K ≥ 4. In reality, we can extend its validity to K = 3 by making the following observation. As explained in Appendix B, the correlator (B.5) takes a factorized form for all K ≥ 3 with the common function F (2) K≥3 given in (B.9). As compared with the function F (2) K=2 that defines the correlator for K = 2, it is given by a different linear combination of the same two-loop conformal integrals. We recall that the factorized form (B.5) is at the origin of the differential recursion (4.8) (see (4.4) and (4.5)). However, in the case at hand the starting point of the recursion is not the 'physical' F (2) K=2 (z) from (4.34) but a different 'auxiliary' function F (2) aux (z) that is defined below in (5.19).
In this way we arrive at the differential equation This is a particular case of the general recursion (4.8), specified to ℓ = 2 and K = 3, but with a different inhomogeneous term F aux (z), not to be confused with F K=2 (z). In other words, we have extended the relation (4.8), initially derived for K ≥ 4, to the case K = 3.
To find the explicit expression for F 1 + L 2 + zL 1 + zL . (4.36) At two loops, the differential equations (4.35) and (4.8) can be solved recursively for K = 3, 4, . . . in the same manner as in the one-loop case. The general solution takes the form (4.30). The polynomials a (w) i,K (z) and b K (z) in (4.30) can be found by substituting the ansatz (4.30) into (4.35) and (4.8) and by taking into account that the polylogarithmic functions L where c i,K (z) and b K (z).
The two boundary conditions needed to fix the solution of the second-order differential equations are the absence of z −2−K singularities at z → 0, and the finiteness of F  K=3 (z) has a logarithmic singularity at z = 1 (see (4.42) below). As before, we use the first sum rule in (3.13) to fix the remaining freedom in the solution for F (2) K=3 (z). The sum rule requires the evaluation of integrals involving HPLs. In Section 6.3, we describe a semi-numerical implementation of the sum rules for the two-loop energy correlation (4.30) that yields exact values of the remaining unknown.
We have solved explicitly the recurrence relations for the polynomials a (w) i,K (z) and b K (z), up to K = 100. In the anciliary file we provide a Mathematica code that constructs the solutions recursively. The procedure is completely algebraic and it only requires solving a system of linear equations with rational coefficients at the K-th step. Its solution is then fed in the linear system of the (K + 1)-th step, etc.
Unlike the one-loop case, we have not attempted to find a closed-form expression for F K (z) for arbitrary K. Instead, we used the obtained solutions for several small values of K to deduce the asymptotic behavior of F (2) K (z) as z → 0 for any K, where The expression on the right-hand side of (4.38) contains a pole at z = 0. To make F K (z) integrable at the origin, one has to add to (4.38) the contact terms localized at z = 0. These terms are derived in Section 6.3 for any K. The two-loop EEC F (2) K=2 , calculated in [18], has a similar asymptotic behavior at z → 0, see (6.12). 9 Let us note that in order to impose the boundary conditions it is sufficient to work out only the leading terms in the expansions z ∼ 0 and z ∼ 1 of the HPL combinations L (w) i , and they do not depend on K. Then fixing the boundary conditions is completely algebraic.
As compared with the one-loop result (4.21), the two-loop correction (4.38) is enhanced at small z by a factor of log z. Adding together (4.21) and (4.38), we keep the most singular terms at each loop order to get aF (1) where a = λ/(4π 2 ). This relation is in agreement with the expected small z behavior of the energy correlation which follows from the operator product expansion on the celestial sphere (see (8.7) below).
Here γ(λ) = λ/(2π 2 )+O(λ 2 ) does not depend on K and coincides with the anomalous dimension of the twist-two operators of spin 3. For z → 1, the two-loop correction to the energy correlation F For K = 2, the function F K=2 (z) has a stronger singularity at z → 1, namely it grows as O(log 3 (1 − z)/(1 − z)) at the end point, see (6.12).
The functions (K + 1) z F K (z) are plotted in Figure 6 for several values of K. The additional factor of (K + 1)z is inserted to soften the singularity at z = 0, see (4.38), and to ensure that the function stays finite as K → ∞. As follows from (3.13), the integral of (K + 1)zF (2) K (z) over the interval 0 < z < 1 has to vanish. This explains why the functions shown in Figure 6 change sign at some value of z for K ≥ 3. For K = 2 the function takes negative values for 0 < z < 1 and its integral vanishes after one takes into account the contact term proportional to δ(1 − z). We recall that such contact terms are absent for K ≥ 3.
We observe that, in agreement with (4.41), the curves in Figure 6 look alike for small z. The situation is different for finite z, the functions (K +1)F (2) K (z) become more and more flat as K increases. We have observed the same pattern at one loop. We recall that the z−dependence describes the angular distribution of the energy on the celestial sphere. The flattening of the z−dependence of the energy correlation at large K implies that the energy distribution becomes more homogeneous.

Heavy sources at weak coupling
In the previous section, we computed the energy correlations at weak coupling for an arbitrary weight K of the source operators and then examined their asymptotic behavior at large K. In  Figure 6. The two-loop correction to the energy correlation (K + 1) z F K (z) for several values of the weight K of the source. The limiting curve for K → ∞ is labeled 'asympt'. Owing to the normalization prefactor, all plots scale as 3 log z + O(z 0 ) at z → 0, see (4.38). For z → 1 the two-loop correction develops a pole at K = 2, scales as a power of log(1 − z) at K = 3 (see (4.42)) and approaches a finite value for K ≥ 4. The latter decreases with K and is given by 49 144 + 181 180 ζ 2 + 85 6 ζ 3 + 7 4 ζ 4 ≈ 20.9176 in the limit K → ∞, see (5.24) below.
this section, we present another approach that allows us to obtain this asymptotic behavior directly, without going through the details at finite K. It is based on the method developed in [16][17][18]20] for the calculation of event shapes with sources of weight K = 2. The main advantage of this method is that it is applicable to the energy correlations (3.4) for an arbitrary coupling constant.

Mellin approach
The energy flow operator E(n) in (3.4) is given by the stress-energy tensor integrated along the light-ray defined by the null vector n µ (see (A.1)). As a consequence, the energy correlation (3.4) is related to the four-point function ⟨O K (1)T (2)T (3)O K (4)⟩ involving two heavy source operators and two stress-energy tensors. In the maximally supersymmetric N = 4 SYM theory, the stress-energy tensor belongs to the same supermultiplet as the simplest half-BPS operator It proves convenient to use the Mellin representation for this function, where the integration contours run parallel to the imaginary axis to the left of the origin, Re j i = −δ with 0 < δ < 1. The Mellin amplitude M K (j 1 , j 2 ) depends on the coupling constant. The expansion coefficients M (ℓ) K of M K (j 1 , j 2 ) at weak coupling in the planar limit, are not all independent. As we have already discussed in the beginning of Section 4, and provide more details in Appendix B (see (B.15)), the ℓ-loop perturbative corrections with K ≥ ℓ + 1 are all the same, Applying the N = 4 superconformal Ward identities, one can express the correlation function in (3.4) in terms of the function (5.2) or equivalently the Mellin amplitude M K (j 1 , j 2 ). The result is (see [19,25] for details) where the function G K (γ) is given by a Mellin integral involving the same amplitude M K (j 1 , j 2 ) as in (5.2), and γ is a dimensionless Lorentz scalar variable invariant under the independent rescaling of the null vectors n i , γ = 2(xn 1 )(xn 2 ) x 2 (n 1 n 2 ) . (5.7) The differential operator acting on G K (γ) in (5.5) stems from the relation between the correlators with different spins mentioned in the beginning of this subsection.
Substituting (5.5) into (3.4) and going through the calculation we can obtain the expression for the energy correlation ⟨E(n 1 )E(n 2 )⟩ K . It takes the expected form (3.5) with the function F K (z) given by the double Mellin integral Here the Mellin amplitude M K (j 1 , j 2 ) defines the four-point correlation function (5.2) and depends on the coupling constant. The function K K (j 1 + j 2 , z) is called the detector kernel. It is independent of the coupling constant and is given by The derivation of this relation can be found in Appendix C, see (C.8) for s 1 = s 2 = 2. We see that the dependence on the angular variable z in (5.8) comes only through the detector kernel. We remark that the relation (5.8) holds for any value of the coupling constant. One can check that the function F K (z) defined in (5.8) satisfies the differential equation (4.8), irrespectively of the choice of the Mellin amplitude M K (j 1 , j 2 ), provided that the latter does not depend on K. For example, this is the case at weak coupling in the planar limit for sufficiently large K, namely K ≥ ℓ + 1, see (5.4).
In the special case K = 2 the kernel simplifies drastically and it is given by (B.17) (see [16]- [19]). In the next subsection, we discuss the asymptotic behavior of the kernel (5.9) at large K. As we show below, the main advantage of the Mellin representation (5.8) is that it allows us to demonstrate directly that the corrections to the energy correlation scale as 1/K for a finite coupling constant. In particular, at weak coupling This relation generalizes (4.26) to any loop order.

Leading term of the large K expansion
For arbitrary K the detector kernel (5.9) is given by a sum of three terms, each containing a hypergeometric function. At large K these functions can be replaced by 1. In this way, we obtain that, to the leading order in 1/K, the detector kernel (5.9) takes the form where the leading coefficient function κ(j, z) is given by This relation is conveniently rewritten in terms of a fourth-order differential operator, It coincides with the analogous operator in (5.5) after replacing γ with 1/z. Inserting (5.13) into (5.8) and changing the integration variable j 2 to j = j 1 + j 2 , we arrive at the relation (5.10) with the function φ (ℓ) (z) given by To illustrate the relation (5.14), we repeat the calculation of the function φ (ℓ) (z) at one loop, i.e. for ℓ = 1. The one-loop Mellin amplitude is given by (B.14). After its substitution in (5.14) the integration over j 1 can be done immediately leading to Here in the second relation we closed the integration contour in the left half-plane and picked the residues at the poles j = −1, −2, . . . . Applying the differential operator (5.13) we recover the expression (4.27), obtained in the previous section from the generating function.

Integral relation between the kernels at large K
The differential operator in (5.13) allows us to establish a remarkably simple relation between the detector kernels (5.9) for K = 2 and for K → ∞. We recall that the kernel K K=2 (j, z) is given by (B.17). Applying the identity we get from (5.13) Notice that the Mellin parameter j enters the right-hand side of (5.17) through the argument of the K = 2 detector kernel. Then, substituting (5.11) and (5.17) into (5.8) we can swap the Mellin integration with that in (5.17) to get where F (ℓ) aux (τ ) is the convolution of the Mellin amplitude of the correlation function (5.2) and the K = 2 detector kernel, At one loop, the Mellin amplitude M K is independent of K and is shown in (B.14). As a consequence, the function F In summary, we have calculated the leading term in (5.10) in three different ways: (i) extracting it from the generating function (4.23); (ii) doing the Mellin integrations in (5.14); (iii) doing the one-fold integration in (5.18).
At two loops, the Mellin amplitude M K takes different forms for K = 2 and K ≥ 3, see (B.15). In the former case, the function F

Two-loop corrections to the large K asymptotics
In this subsection, we use the relation (5.18) at ℓ = 2 and compute the function φ (2) (z) defining the leading term of the large K asymptotics of the energy correlation at two loops. The evaluation of the integral in (5.18) is much simpler than the double Mellin integral (5.14). As explained above, the function F According to (4.36), the function F (2) aux (z) is a combination of classical polylogarithms of transcendental weight up to three. It can also be expressed in terms of harmonic polylogarithms [29], whose arguments define an alphabet of three letters, Note that the two-loop correction (4.30) involves the same HPL alphabet. The integration of F  Table 1. Maximal transcendental weight of the polylogarithms in the expression for the energy correlation F (ℓ) K (z) at ℓ loops for finite K and in the limit K → ∞.
we can write where g Note that (5.21) is invariant upon the reflection √ z → − √ z that is equivalent to the complex conjugation x → x * . Substituting (5.21) into (5.18), we apply the differential operator (5.13) and obtain a closedform expression for the function φ (2) (z), to be found in the ancillary file. We recall that this function defines the leading large K behavior (5.10) of the energy correlation at two loops. It is interesting to note that, due to the appearance of the new letter in the HPL alphabet, its expression involves a larger set of functions than the two-loop energy correlation F (2) K (z) at finite K. Also, the maximal transcendental weight of these functions varies with K, as summarized in Table 1.
The dependence of zφ (2) (z) on the angular variable 0 < z < 1 is shown in Figure 6 (see the curve labelled 'asympt'). For z → 1 it approaches a finite value, For z → 0 it grows as log z/z, We verify that this relation is in agreement with the large K limit of (4.38).

Subleading corrections
So far we have discussed only the leading term in the large K expansion of the energy correlation (5.10), The differential equation (4.8) allows us to find the subleading coefficient functions φ where the differential operator D is defined in (4.9). Substituting the expansion (5.26) into this equation and comparing the coefficients on both sides, we derive a recursion relation for φ n (z) with n ≥ 2. The solution is For z → 1, the functions (5.26) approach a finite value. At small z, replacing φ (ℓ) (z) ∼ (log z) ℓ−1 /z (see (4.41)) in (5.28), it is straightforward to verify that the ratio φ (ℓ) n (z)/φ (ℓ) (z) goes to (−1) n−1 for z → 0. Therefore, the subleading terms in (5.26) effectively modify the coefficient in front of the leading term to 1/(K + 1) where the dots denote the subleading corrections suppressed by powers of λ and z. At two loops, the relation (5.29) is in agreement with (4.40).

Contact terms
The expressions for the energy correlations derived in the previous sections are valid only for 0 < z < 1. To extend them to the end points z = 0 and z = 1, we have to add contact terms. As was explained above, such terms are needed to ensure that the sum rules (3.11) are satisfied. Furthemore, the obtained expressions for the energy correlations have non-integrable singularities at the end points z = 0 (for K ≥ 2 and K → ∞) and z = 1 (for K = 2) and require a careful treatment.
In this section we present two complimentary approaches to computing the contact terms. We start with the simplest example of the one-loop correction to the energy correlation F (1) K=2 (z) given by (4.10) and proceed to the case of arbitrary K at one loop. In Section 6.3 we compute the contact terms at two loops.

Sum rule approach
At K = 2, the one-loop correction to the EEC (4.10) is singular at the end points z = 0 and z = 1. Following Section 4 in Ref. [7] (see also [6,8]), we define the regular (i.e. integrable) part by subtracting the non-integrable singularities from (4.10), The resulting function F (1) reg K=2 (z) is integrable in the interval 0 ≤ z ≤ 1. Next, we compensate the subtracted singular terms by adding them up but now interpreted as (integrable) plusdistributions (for the definitions see (D.1) and (D.2)). In addition, we allow for contact terms at the end points with arbitrary coefficients, The above procedure is the most general regularization of the energy correlation at the end points consistent with the sum rules (3.11). Indeed, subtracting the poles makes the function integrable. The subsequent addition of the subtracted term in the form of the plus-distributions does not modify the energy correlation (4.10) for 0 < z < 1, at the same time maintaining the integrability. At this stage the regularized expression (6.2) still contains the arbitrary coefficients C 1 and C 2 , which is the usual ambiguity in singular distributions. They can be determined from the sum rules (3.13). Substituting (6.2) into the sum rules (3.13) we get for K = 2 leading to C 1 = −1 and C 2 = −ζ 2 . Finally, the complete one-loop expression for the energy correlation at weight K = 2 is where the regular part F (1) reg K=2 (z) is given by (6.1). For higher weights K ≥ 3, the contact terms for the energy correlation F (1) K (z) can be found in the same way. We recall that the function F (1) K≥3 (z) is singular for z → 0 (see (4.21)) but it is finite at z → 1 for K ≥ 4 (see (4.22)). Then we define the regular function F (1)reg K (z) by subtracting the non-integrable singularity at z = 0, where the first term on the right-hand side involves the generating function (4.23). To get the complete one-loop result, we add to F (1)reg K (z) the plus-distribution term [1/z] + and calculate the coefficient of the contact term δ(z) with the help of the sum rules (3.13). In this way, we find for K ≥ 3 Let us emphasize that only the contact term at the end point z = 0 is required, and its coefficient is fixed by one of the sum rules in (3.13). The other sum rule is not sensitive to the contact term in (6.6), and it should be automatically satisfied. This is a useful cross-check of our EEC calculation, At large K, we apply (4.27) to find the regular part of φ (1) (z), The sum rule (3.13) enables us to calculate the contact terms, We see that the distribution terms in (6.10) agree with those in (6.6) at large K.

Mellin approach
In this subsection, we follow Refs. [26,27] to compute the contact terms for the function φ (1) (z) which is given by the Mellin integral (5.15).
We computed the integral in (5.15) by closing the contour from the left and picking the poles at j = −1, −2, −3, . . .. The integrand in (5.15) involves D(z −j+1 ). It is given by the sum of three terms of the form z α with α = −j − 2, −j − 1, −j which we now treat as distributions z α + , see (D.3). The key observation is that this distribution has a contact term in its Laurent expansion, see (D.4). This creates an extra pole under the Mellin integral at α = −1, whose contribution accounts for the contact term in the energy correlation. Of the three values of α listed above, only α = −j−2 is inside the integration contour. Replacing z −j−2 (6.11) in agreement with (6.10). The energy correlation F (1) K=2 with its two contact terms is treated similarly.
We see that the two approaches -the sum rule method of fixing the contact terms and the approach based on the careful calculation of the Mellin integral -give identical results for the one-loop energy correlation. However, technically the Mellin approach may be more difficult to implement beyond one loop. Indeed, there the function M (ℓ) (j 1 , j − j 1 ) in (5.14) is given by a 2(ℓ − 1)−fold Mellin integral. The sum rule approach is more efficient, provided that the function F (ℓ) (z) is known explicitly for 0 < z < 1. We apply it at the two-loop level in the next subsection.

Two-loop contact terms from the sum rules
Let us start with K = 2. The two-loop function F (2) K=2 (z) (4.34) was calculated in [18] for 0 < z < 1. Like the one-loop case, this function has non-integrable singularities at z = 0 and z = 1. Promoting them to plus-distributions, adding the contact terms δ(z) and δ(1 − z) and using the sum rules (3.13) to fix their coefficients, we obtain Here the regular function F K=2 (z) by subtracting the poles at z = 0 and z = 1.
For K ≥ 3, the two-loop function F (2) K (z) contains a non-integrable singularity at z = 0 calculated in (4.38). For z → 1, it is regular for K ≥ 4 and has an integrable logarithmic singularity for K = 3. As a result, only the contact term δ(z) is required. As before, we promote the poles at z = 0 to plus-distributions to get where c 1/z,K is given in (4.39). To find the coefficient of the contact term c (2) δ,K , we substitute (6.13) into the sum rules (3.13). We implemented the integration in (3.13) numerically and carried out the calculation for K ≤ 15. We interpolated the HPL expressions (4.30) by their generalized series expansions in the vicinity of z = 0 and z = 1, and integrated numerically the interpolation formulae achieving a precision of at least 40 digits. This level of precision was sufficient to reconstruct c (2) δ,K with K ≤ 15 as rational linear combinations of the numbers {1, ζ 2 , ζ 3 , ζ 4 } using the PSLQ algorithm. Exploiting these results, we arrived at an expression for c (2) δ,K valid for any weight K, where H (2) n is the n-th generalized harmonic number of degree two. We have also checked numerically that the second sum rule in (3.13) is satisfied for K ≤ 15, Similarly, we use (5.25) supplemented with the sum rule (3.13) to calculate the contact terms for the two-loop function φ (2) (z) defining the leading large K asymptotics (5.10) of the energy correlation, We would like to emphasize that, unlike (6.13), this relation was derived analytically.
In the large K limit the function (6.13) takes the expected form (5.10). Comparing the coefficients of the contact terms in (6.13) and (6.16) we obtain the consistency relation We verified using (6.14) that it is indeed satisfied. Analogous consistency relations hold for c 1/z,K in (4.39) at large K.

Event shapes at strong coupling
After exploring the dependence of the energy correlations on K at weak coupling λ ≪ 1, in this section we perform the same analysis at strong coupling λ ≫ 1. We also consider event shapes in N = 4 SYM with detectors other than energy calorimeters.
The energy correlations at strong coupling were analyzed in [5], where their computation was mapped to the propagation of a probe particle (dual to the source operator) on a shock wave background (dual to the energy calorimeters). In particular, it was found that in theories which are dual to matter minimally coupled to gravity in AdS, the energy correlations do not depend on the angle, This universality is the expression of the equivalence principle in the bulk since the energy correlations in this approximation effectively measure the Shapiro time delay, which in general relativity does not depend on the type of particle in question. The leading stringy correction to the EEC was computed in [5], It was also found there that the connected energy correlations obey This relation is analogous to (1.4) where the role of ∆ H is played by the coupling constant λ.
In this section we rederive and generalize these results starting from the Mellin space representation of the four-point functions ⟨O K O 2 O 2ŌK ⟩ at strong coupling. We consider different event shapes measuring the correlations of various conserved charges (not just energy, see also footnote 22). In close analogy with (3.5) we introduce ⟨J s 1 (n 1 )J s 2 (n 2 )⟩ K = (q 2 ) s 1 +s 2 2(4π) 2 (qn 1 ) s 1 +1 (qn 2 ) s 2 +1 F s 1 ,s 2 K (z) , (7.4) where the flow operator J s (n) is built out of the conserved current of spin s. The energy-energy correlation studied in the previous section corresponds to F 2,2 K . For s = 1 we get the charge detector. The corresponding flow operator J s=1 (n) is given by the light-ray transform of the spin one conserved current (see (A.4)). For s = 0 we get the scalar detector J s=0 (n), it is given by the light-ray transform of the half-BPS operator of dimension ∆ = 2 (see (A.4)). More details about the correlations (7.4) can be found in appendix C.
In this section we compute F s 1 ,s 2 K (z) in the supergravity approximation and derive the leading stringy correction to it, where the dots denote subleading corrections. We find that in the K → ∞ limit the event shapes (7.4) do not depend on the angle and are given by the product of the one-point functions ⟨J s 1 (n 1 )⟩ K ⟨J s 2 (n 2 )⟩ K . We also analyze the event shapes away from the large K limit, and we find that for event shapes other than the energy-energy correlation the dependence on K is nontrivial. For the EEC, in contrast to the weak coupling results where the suppression factor is ∼ 1/K, now the leading correction is suppressed by 1/λ. In our analysis we consider K's that do not scale parameterically with λ, see e.g. [32]. Taking K ∼ λ 1/4 , which corresponds to a source dual to short massive string modes, we do not expect to observe any change to (7.1) and (7.2) above. For K ≳ √ λ, in which case the source is described by a big classical string [33], we expect to get back the 1/K suppression of the leading correction to the energy correlation [34]. It would be interesting to check this explicitly.

Supergravity approximation
The four-point function of the half-BPS operators ⟨O K O 2 O 2ŌK ⟩ has been actively studied at strong coupling recently, starting from [35,36]. We will be only interested in stringy corrections here, leaving quantum gravity corrections aside.
The stringy corrections to the ⟨O K O 2 O 2ŌK ⟩ Mellin amplitude were considered in [37]. To utilize their results, let us notice that the Mellin variables (s, t) used in that paper are related to (j 1 , j 2 ) used in the present paper as follows 12 The Mellin amplitude M K (s, t) in [37] is related to our M K (j 1 , j 2 ) as follows, Let us notice that while the normalization of the scalar correlator ⟨O K O 2 O 2ŌK ⟩ does not have an intrinsic meaning, it is important for us because we use superconformal Ward identities to relate it to the correlators with conserved currents which have canonical normalization, see appendix C for details. This explains the factor 1 K in the formula above. The supergravity result for the Mellin amplitude M K (s, t) takes the form .
It is then easy to compute various event shapes by plugging this Mellin amplitude into the master formula (C.7). We get the following results It is interesting to notice that while the energy correlations do not depend on K, in agreement with the arguments of [5], the results for other event shapes require simply polynomial corrections. We can also consider the large K → ∞ limit. According to the general clustering arguments at the beginning of this paper all the event shapes effectively cluster, in other words they become angle-independent, with the correction to clustering being of order 1/K. Let us comment on the power of K that appears in (7.10). As the two-point function clusters, the power is dictated by the square of the corresponding one-point event shape. This in turn is given by the light transform of the corresponding three-point function which in our case scales as ∼ K. The light transform of a detector operator with quantum numbers (∆, J) produces an extra factor ∆ 1−J H . This is why starting from the three-point function which scales as K, we get K 4 for scalar-scalar correlations (detectors of spin J = 0), K 2 for charge-charge correlations (spin J = 1), and K 0 for energy-energy correlations (spin J = 2).

Stringy correction
The computation of the stringy corrections to event shapes is subtle because the event shapes are sensitive to the high-energy limit of the scattering amplitudes in AdS. In particular, computing the higher derivative stringy corrections to the Mellin amplitude and plugging them into the formula (C.7) produces divergent results.
To solve this problem we adopt the approach used in [24] based on the Borel re-summation of the leading high-energy corrections to the Mellin amplitude. The basic observation is that these are controlled by the scattering of strings in flat space, see also [5]. One can then use the formulas that relate the Mellin amplitude to its flat space limit to predict the form of the relevant series [38].
The leading stringy correction to the Mellin amplitude takes the form [37] M stringy If we try to plug this formula into the generating function for the event shapes, we find that the t integral takes the form dt 2πi M stringy K (s, t) and is divergent. Therefore we need to re-sum the stringy corrections. More precisely, rescaling t → √ λt we notice that infinitely many terms in the stringy expansion of the Mellin amplitude become of the same order dt t 2n λ 3/2+n ∼ 1 √ λ , which is indeed the expected leading stringy correction to the event shape [5].
In [24], using the flat space limit of the amplitude, it was found that for K = 2 the relevant series takes the form n=0,n−even c n x n , c n = 1 4 Γ(6 + n)ζ(3 + n) 2 n+1 . (7.12) For general K, using [37], we get instead n=0,n−even c n x n , c n = 1 4 Γ(4 + n + K)ζ(3 + n) 2 n Γ(1 + K) . (7.13) To perform the Borel sum, we change the summand c n x n → cn(σx) n Γ(n+1) . We then substitute ζ(3 + n) = ∞ k=0 1 (k+1) 3+n and perform the sum over n. As a result, we get the following integral when evaluating dt where σ is the integration of the Borel transform. To evaluate this integral we notice that if we rescale x → x σ , which assumes σ ̸ = 0, we get zero. Therefore we can limit the integral over σ to an infinitesimal interval around the origin and drop e −σ from the integrand.
In this way we get for the integral of the re-summed Mellin amplitude Using this result we can compute the stringy corrections to the event shapes: As an extra consistency check, we verify that the charge-charge correlation F 1,1 stringy (z) satisfies the relation 1 0 dz F 1,1 stringy = 0, which follows from the conservation of the total charge. Similarly, 1 0 dzF 2,2 stringy = 1 0 dzzF 2,2 stringy = 0 due to the energy-momentum sum rules (3.11). Putting K = 2 in (7.16), we recover the results of [17] that, due to the superconformal Ward identities, all F s 1 ,s 2 stringy (z) are proportional to each other up to a power of z. This simple relationship between event shapes does not hold for K > 2.

Clustering in CFT
In this section we present some arguments in favor of (1.4) in a general CFT. We start by discussing the clustering of correlation functions involving the stress-energy tensor in heavy states. It is analogous to familiar clustering of correlation functions in QFT as the spatial separation between operators becomes large. We then discuss the leading non-universal correction to the disconnected result which depends on the details of the theory and the heavy operator in question. We then proceed to event shapes which is the main topic of this paper, and we conjecture that the clustered structure of the stress-energy tensor correlators is not affected by the light-ray transform. Finally, we discuss the special case of heavy states in planar theories c T ≫ ∆ H ≫ 1, where c T is defined via the two-point function of stress tensors ⟨T µν T ρσ ⟩ ∼ c T . 13 The statements in this section are based on some basic physics intuition and we do not attempt to prove them rigorously starting from the CFT axioms [40].

Clustering of local operators in a heavy state
In QFT, the correlation functions of local operators factorize when the mutual separation between the local operators goes to infinity [41]. In CFT a closely related property is expected to hold for the correlation functions of local operators in heavy states. We define Let us also introduce the notation T (x) ≡ T µν (x)z µ z ν for the stress-energy tensor operator contracted with a null polarization vector z µ . We keep the dependence on the polarization implicit to avoid cluttering. We can now formulate the clustering of the stress-energy tensor operators in a heavy state as follows Let us motivate why (8.2) might be true. We consider a CFT on R × S d−1 with the heavy operator defining a state of energy E = ∆ H R and the physical distances between the operators being L ij = R|x i − x j |. The key point is that we can naturally associate to the limit ∆ H → ∞ an infinite volume or thermodynamic limit by sending at the same time R → ∞, see e.g. [42]. Such a limit is not uniquely defined because it requires specification of which quantities are kept fixed as we take the limit.
It is natural to keep certain local densities fixed as we take the limit. For our purposes we can keep the energy density ϵ = E R d−1 = ∆ H R d fixed. It is also natural to keep the physical distances between the operators L ij fixed. Notice that it corresponds to the multi-OPE limit |x ij |→ 0 in the space of cross-ratios. On general grounds, we then expect to get nontrivial correlation functions in the infinite volume R d−1 in an excited state characterized by a characteristic scale L ϵ . In CFT, due to the absence of dimensionful parameters, such a correlator depends on the dimensionless ratios L ij Lϵ . We expect that for L ij Lϵ ≪ 1 we get nontrivial correlations in the excited state, whereas for L ij Lϵ ≫ 1 the correlator clusters as the relative spatial separation between the operators goes to infinity. 14 Because the stress-energy tensor measures the energy density, the resulting one-point function ⟨T ⟩ ϵ ̸ = 0, and the disconnected piece indeed provides the leading answer to the correlator at large distances. 13 See, for example, [39] for the precise definition. 14 For the closely related finite temperature CFT correlators the clustering property was discussed in [23].
There is an obvious generalization of the argument above for CFTs with global symmetries. Let us assume for simplicity a CFT with U (1) global symmetry and denote the corresponding conserved current J(x) ≡ J µ (x)z µ . We can then consider a source of large charge and apply the same argument to argue that where this time we keep the charge density q = Q R d−1 fixed as we take the limit and we use the fact that ⟨J⟩ q ̸ = 0 because the conserved current measures charge density.
We expect similar arguments to hold for any local operators and not just conserved currents. However, in this case the non-vanishing property of the one-point function is not guaranteed.
Let us next discuss the leading correction to the disconnected result discussed above. For this purpose it is convenient to consider the four-point function. At the level of the original correlation function on the plane the thermodynamic limit discussed above corresponds to the following limit where we have put all operators on a 2d plane with complex coordinates (z,z) and we set them The separation between the operators on the right-hand side after taking the limit is x 2 = ww. The extra factor ∆ −2 H has a kinematic origin and comes from mapping the correlator on the plane to the one on the cylinder, see [42] for details. The two-point function ⟨T µν (w,w)T ρσ (0, 0)⟩ ϵ stands for the two-point function in the microcanonical ensemble with energy density ϵ for the theory on R 1,d−1 . Locally, it is the same as the finite temperature correlator where the temperature is fixed to correctly reproduce the energy density ϵ.
As we take the operator insertions apart, in an interacting theory the leading correction to the disconnected result is expected to decay exponentially, see e.g. [43], where m th is the so-called thermal mass. It could also happen that the effective theory that emerges in the limit ∆ H → ∞ is gapless, in which case the leading correction to the disconnected piece is a power. For the thermal case it happens for example in the case of the free scalar theory, and for the large charge case in the case of the three-dimensional O(2) model [44,45]. In addition to the connected corrections in the macroscopic state we can imagine having corrections related to the finite curvature of the sphere, e.g. terms that go to zero as Lϵ R . To summarize, while the leading result (8.2) is universal, the corrections to it depend on the details of the theory and the nature of the heavy state. As such, they cannot be computed on general grounds and require a detailed case-by-case analysis.

Event shapes
Going from the formulas for correlation functions of local operators to the calculation of event shapes includes the extra step of integrating over the detector times, or, equivalently, performing the light transform. When we do that, the separation between the detector operators |x 12 | becomes arbitrarily small [46]. The small separation regime includes the light-cone limit and the Regge limit. It could then happen that while the clustering holds for local operators at fixed separation, it fails for the light-ray operators. As a simple example, consider the chargecharge correlation related to (8.3). As we perform the light-transform, the one-point charge correlation produces a finite result, whereas the two-point (or higher-point) function in general will diverge in the Regge limit. 15 Therefore, the analog of (8.3), in general will not hold for charge correlations.
The situation is better for the energy correlations which are finite, or IR-safe, observables. In this case, the contribution of the Regge limit to the light transform produces a finite result. The question that remains then is whether this contribution is suppressed, compared to the factorized result in the limit of the heavy source ∆ H → ∞. We do not know how to show it in general, but based on the explicit example analyzed in the paper and the intuition that the Regge limit is naturally suppressed for energy correlations, we conjecture that in a general CFT Celestial clustering: lim To understand this better in a non-perturbative setting it is interesting to consider the light-ray OPE, which captures the behavior of the energy correlations at small angles. For simplicity, let us start with the two-point energy correlation ⟨E(n 1 )E(n 2 )⟩. The leading OPE contribution at small angles z = sin 2 (θ/2) ≪ 1 takes the form (including the 1 from the disconnected contribution) where O + 3 is the spin three light-ray operator of positive signature belonging to the stress-energy tensor Regge trajectory and γ + 3 is its anomalous dimension. As follows from (8.7), the contribution of O + 3 to the energy-energy correlation diverges at small angles for γ + 3 < 2. Assuming that this relation is satisfied, a necessary condition for (1.4) to hold is We then see from (8.7) that there is an emerging characteristic angle defined as 1 It plays the same role as θ 0 in QCD. Namely, for θ ≲ θ * the second term on the right-hand side of (8.7) dominates and the energy-energy correlation is in the OPE regime. For θ ≳ θ * this term is subdominant and the energy-energy correlation is in the flat regime. Notice that as the angle between the detectors decreases, the transition in QCD is OPE-to-flat, whereas for CFT energy correlations in heavy states it is flat-to-OPE.
Similarly, considering the multi-collinear limit of energy correlations, see e.g. [8,49], we get a spin k + 1 light-ray operator Let us now see how this comes about in the simple multi-particle model of the final state. Imagine that the source carries total energy E which is shared between K particles. In this case the operator O + J schematically measures the (J − 1)-th moment of the energy distribution In the heavy limit we expect that K(∆ H ) → ∞, which creates the hierarchy between the different connected contributions. For example, in the free scalar theory the source ϕ K creates a K-particle state, so that h k = k−1.

Planar theories
In the discussion above we kept c T fixed as we took ∆ H → ∞. We expect that the suppression of connected energy correlations takes place in planar theories as well, where c T ≫ ∆ H ≫ 1. In this case, however, the argument about clustering in the thermodynamic limit does not apply and a different argument is needed. As we have explicitly shown in this paper using the example of half-BPS operators in N = 4 SYM at large charge K, the suppression of connected correlations is different at weak (1.5) and at strong (1.6) coupling.
At weak coupling, perturbatively in λ, the connected energy correlations are suppressed by 1 K . The mechanism is precisely the one described above, where the heavy operator creates a multi-particle state of weakly interacting particles with the number of particles ∼ K. We expect this picture to be qualitatively correct for λ ≪ 1 ≪ K.
At strong coupling, perturbatively in 1 λ , the connected energy correlations are suppressed by 1 √ λ . In this case, the strongly coupled dynamics leads to a copious production of particles [50], with the characteristic energy controlled by λ and not by K. We expect this picture to be correct for 1 ≪ K ≪ λ.
More generally, it would be interesting to exlpore different scalings between K and λ [32] (in particular, see Figure 1 in that paper): K ∼ 1 (SUGRA); K ∼ λ 1/4 (short massive strings); K ∼ λ 1/2 (big classical strings). For example, the latter case is related to the so-called Frolov-Tseytlin limit [33,34], where the heavy state describes a classical "fat" string in AdS and to leading order the correlators are again simple, given by the one-point functions on this classical solution. It would be also interesting to study energy correlations in the recently introduced large charge 't Hooft limit [51].
We recall that in the maximally supersymmetric N = 4 Yang-Mills theory the energymomentum tensor T µν together with the R-symmetry current J µ and the half-BPS operator of weight two O 2 are members of the same supermultiplet. By analogy with the energy flow (A.1) one can define the R-charge and scalar flow operators (the R-symmetry indices are not displayed) where J + ≡n µ J µ (x)/(nn). We denote collectively these three flow operators as J s (n) where the spin s takes the values s = 0, 1, 2.
The definition (A.1) of the flow operator involves two steps: (i) we send the detector, i.e. the projection T ++ of the energy-momentum tensor, weighted with the factor x 2 + , to spatial infinity; (ii) we integrate over the entire working time of the detector −∞ < x − < ∞. These manipulations are done with each insertion of the energy-momentum tensor into the Wightman correlation function of the half-BPS scalar source and sink of weight K, In close analogy with the k−point correlation functions, the correlations (A.5) can be decomposed into connected pieces (cumulants) ⟨E(n 1 ) . . . E(n k )⟩ c . It is convenient to do it by introducing the energy correlations via a generating functional, (A.6) Here J i are sources coupled to the energy flow operators E(n i ) and the generating functional Z(J) is defined by In particular, Note that the expressions for the connected correlations of four or more operators involve nonlinear terms. Introducing the normalized energy flow operators E(n i ) = E(n i )/⟨E(n i )⟩, one finds that (A.6) is equivalent to (1.3).
In this paper we are mainly interested in the case of two energy flow operators, the socalled energy-energy correlation (EEC). 16 For illustration purposes, in this appendix we show the details of the calculation of a single energy insertion at Born level [5], 17 The starting point is the three-point function of two half-BPS operators of weight K and one energy-momentum tensor. In Euclidean space it is given by where we have set x 3 = 0, as indicated in (A.9). The last term on the right-hand side of (A.10) ensures that δ µν (G E ) µν (1, 2, 3) = 0. We perform a Wick rotation to the Wightman function (G W ) µν in Minkowski space-time by inserting the prescription x 2 ij → x 2 ij − iϵx 0 ij for i < j. Next we project the free vector indices withn µ , according to the definition (A.2) (this removes the trace part). Using the decomposition (A.3) in the limit x 2+ → ∞ we obtain . With this we take the limit We observe two poles in the variable x 2− located on the opposite sides of the real axis. Closing the integration contour in, say, the upper half-plane, we get Notice that the auxiliary light-like vectorn has dropped out of the right-hand side. The last step in the procedure is the Fourier transform in (A.9) which gives In Section 3.3 we deal with the two-point correlation ⟨E(n 1 )E(n 2 )⟩ at Born level. The main contribution comes from the first diagram in Figure 2. It factorizes into two one-point expressions of the type (A.12), with two different direction vectors n 1,2 : (A.14) Introducing Schwinger parameters for the Fourier integral of the function in the second line, we arrive at Eq. (3.15) in the main text. There exist other diagrams where the two detectors exchange one propagator. They give rise to a contact term, as explained below.
x 1 x 4 x 2 x 3 x 1 x 4 x 2 x 3 Figure 7. Cross-talk diagrams for Born-level QSC. The blob with a cross and the grey blob denote the charge and scalar detectors, respectively.
The origin of the contact term δ(z) at Born level Let us consider a simple example of an event shape, the Born level QSC (or charge-scalar correlation), defined in (7.4) with s 1 = 1, s 2 = 0. It can be obtained from the Born level correlator ⟨O (1) where O =φϕ is a scalar operator made of a charged scalar field ϕ(x), and J µ = iφ ← → ∂ µ ϕ is the U (1) current. We want to turn the operators at points 2 and 3 into detectors according to (A.4). The procedure becomes singular if the two detectors can exchange a particle ('cross talk'). At Born level this means to connect the detector points by a propagator. The result, as we show below, is the contact term Fourier transform x 14 → q turns this into the expected contact contribution q 2 (n 2 q) 2 (n 3 q) δ(z) (with z as defined in (2.13)) to the event shape QSC.

(A.22)
This simple example explains the origin of the contact term δ(z) in the event shapes at Born level. For a more general discussion see [27].

B Four-point correlation functions of half-BPS scalar operators
In this appendix, we summarize the properties of the four-point correlation functions in N = 4 SYM that we use in computing the energy correlations.
The four-point correlation functions of interest are where O K (x, Y ) are half-BPS operators defined in (3.1). They are annihilated by half of the N = 4 supersymmetry charges and their scaling dimension ∆ = K is protected. The half-BPS scalar operator (3.1) is the lowest component of a short N = 4 supermultiplet. Among the operators (3.1) with different weights K, the one with K = 2 plays a special role. The corresponding N = 4 supermultiplet becomes ultrashort [52]. It contains all the conserved currents of the theory, including the SU (4) R-symmetry current of spin one and the energy-momentum tensor of spin two, hence the name "stress-energy supermultiplet". In virtue of N = 4 superconformal symmetry, the four-point correlation functions of the operators belonging to the same supermultiplet are related to each other by Ward identities.
In application to the energy correlations, we encounter the four-point correlation function ⟨O K (1)T µ 1 ν 1 (2)T µ 2 ν 2 (3)O K (4)⟩ where the half-BPS operators define the source/sink and the stress-energy tensors describe the calorimeters. Applying the N = 4 superconformal Ward identities, this correlator can be expressed in terms of the four-point function of scalar operators ⟨O K (1)O 2 (2)O 2 (3)O K (4)⟩. Below we summarize the properties of the latter.
The correlator G K22K splits into the sum of two terms 2) The first, Born-level part is a rational function of the space-time coordinates and a polynomial in the auxiliary variables Y , independent of the coupling. The second, coupling dependent correction part involves non-trivial functions originating from Feynman integrals. The expression for the Born term G 0 K22K can be obtained by Wick contracting the scalar fields belonging to the four half-BPS operators. It is a polynomial of degree (K + 2) in the free scalar propagators The coefficients of this polynomial depend on the rank of the gauge group N c . In the planar limit, for N c → ∞, the connected part of G 0 K22K is given by where the K−dependent combinatorial factor counts the number of Wick contractions.
The interacting (coupling dependent) part of (B.2) takes the following form in the planar limit (see [53] and references therein) Here R is a universal rational prefactor carrying SU (4) weight 2 and conformal weight 1 at each point, As explained in [17], putting Y 2 = Y 3 is equivalent to projecting the product of the two operators where the functions F (ℓ) K can be expanded over a basis of conformal ℓ−loop four-point integrals. At K = 2 the expansion (B.8) was derived in [54][55][56][57] up to ten loops. For K ≥ 3 the analogous expansion is known up to five loops [53,58]. 21 The conformal integrals are independent of K but they appear in the expression for F K (x, λ) with K−dependent coefficients.
The explicit expressions of the function in (B.8) for ℓ = 1, 2 are . (B.10) Combining together (B.8) and (B.9), one can verify that the function F K=2 is invariant under the exchange of any pair of points in (1,2,3,4), whereas F K>2 is symmetric in the two pairs of points (1,4) and (2,3) separately, in agreement with the Bose properties of the correlation function (B.2). As follows from (B.9), the one-loop correction in (B.8) is independent of K. At two loops, the functions F (2) K=2 and F (2) K>2 are given by different linear combinations of the same conformal integrals. The three-loop results of [53] show that the functions F K are different for K = 2 and K = 3, and again become the same for all K ≥ 4. The proposed generalization of [59] indicates that a similar pattern continues at higher loops. Namely, the function F (ℓ) K ceases to depend on K for K ≥ ℓ + 1. This feature allows us to study the event shapes recursively, see Section 4.3.

C Detector kernel of the Mellin representation
In this appendix, we present some details of the derivation of the Mellin representation (5.8) and (5.9) of the energy correlation. We recall that the operator E(n) in (3.4) measures the flux of the energy in the direction specified by the null vector n µ and it is built out of the energy-momentum tensor. If the underlying theory contains a conserved current of spin s, 22 we can define the flow operator J s (n) in an analogous manner. For spin s = 2 it coincides with E(n) whereas for arbitrary spin s it measures the flow of the corresponding conserved charge. 23 We generalize (3.4) and define the correlation of charges with different spins ⟨J s 1 (n 1 )J s 2 (n 2 )⟩ K = σ −1 tot (q) d 4 x e iqx ⟨O K (x)J s 1 (n 1 )J s 2 (n 2 )Ō K (0)⟩ , (C.1) where the total cross-section σ tot (q) is given by (3.3).
In N = 4 SYM the charges J s are members of the stress-energy multiplet. As a consequence, the correlation functions ⟨O K (x)J s 1 (n 1 )J s 2 (n 2 )Ō K (0)⟩ with different spins s 1 and s 2 are related to each other by N = 4 superconformal Ward identities. The general solution to these identities was derived in [19], ⟨O K (x)J s 1 (n 1 )J s 2 (n 2 )Ō K (0)⟩ = 2 s 1 +1 i s 1 +s 2 (n 1 n 2 ) s 1 +1 where γ is defined in (5.7) and the function G K (γ) is independent of the spins s 1 ≥ s 2 .
The relation (C.2) was derived in Ref. [19] for sources of weight K = 2. As discussed around (B.5), the generalization to arbitrary weight K can be achieved by simply attaching (K − 2) additional scalar propagators stretched between the source/sink operators, thus introducing the dependence of G K (γ) on the total weight. Going through the analysis in [19] we see that the origin of the relation (C.2) has to do with the supermultiplet buildup at the detector points. At the source/sink points we stick to the lowest weight state of the supermultiplet, so the result where ϕ (k) (0) denotes the k-th derivative, and similarly for the distribution (1 − z) α + . Note that z −1 + is not the value of z α + at α = −1. The distribution z α + admits the following Laurent series expansion in the vicinity of α = −1: