The Energy-Energy Correlation at Next-to-Leading Order in QCD, Analytically

The energy-energy correlation (EEC) between two detectors in $e^+e^-$ annihilation was computed analytically at leading order in QCD almost 40 years ago, and numerically at next-to-leading order (NLO) starting in the 1980s. We present the first analytical result for the EEC at NLO, which is remarkably simple, and facilitates analytical study of the perturbative structure of the EEC. We provide the expansion of EEC in the collinear and back-to-back regions through to next-to-leading power, information which should aid resummation in these regions.

Introduction. The energy-energy correlation (EEC) [1] measures particles detected by two detectors at a fixed angular separation χ, weighted by the product of the particles' energies. The EEC is an infrared-safe characterization of hadronic energy flow in e + e − annihilation. It has been used for precision tests of quantum chromodynamics (QCD) and measurement of the strong coupling constant α s [2,3]. In perturbative QCD, the EEC is defined by where i and j run over all the final-state massless partons, which have four-momenta p µ i and p µ j (including the case i = j at χ = 0); Q µ is the total four-momentum of the e + e − collision and dσ is the phase-space measure. The three-vectors n i,j point along the spatial components of p i,j . The definition (1) implies the sum rule where σ is the total cross section for e + e − annihilation to hadrons. The leading order (LO) QCD prediction for the EEC has been available since the 1970s [1]: where σ 0 is the Born cross section for e + e − → qq, C F is the quadratic Casimir in the fundamental representation, and we have introduced z = (1−cos χ)/2. The cross section is strongly peaked at χ = 0 (z = 0) and χ = π (z = 1), regions that require resummation of logarithms due to emission of soft and collinear partons. At intermediate angles, higher-order corrections tend to flatten the distribution. The EEC was first computed numerically at next-toleading order (NLO) in QCD by several groups in the 1980s and 1990s, originally leading to conflicting results. Different methods were used to handle soft and collinear singularities from real radiation: phase-space slicing [4][5][6][7] subtraction methods [6,[8][9][10][11][12][13][14], or hybrid schemes [6,7,15]. Accurate numerical NLO results are available from the program Event2, based on dipole subtraction [13,14]. Quite recently, the EEC has been computed at NNLO in QCD using the CoLoRFulNNLO local subtraction method [16,17].
In perturbation theory, the EEC is singular in both the collinear (z → 0) and back-to-back regions (z → 1), as can be seen explicitly from Eq. (3). The leadinglogarithmic collinear behavior can be obtained from the "jet calculus" approach [18,19], in terms of the anomalous dimension matrix of twist-two, spin-three operators [11,19]. Resummation of the EEC in the backto-back (Sudakov) region has been performed at nextto-leading-logarithmic (NLL) and NNLL accuracy [20][21][22]. Quite recently, a factorization formula for the EEC has been derived which permits its resummation to N 3 LL [23]. Possible non-perturbative corrections to the EEC have also been investigated [24].
In N = 4 super-Yang-Mills theory (SYM), the EEC has been computed analytically at NLO in terms of classical polylogarithms [25], using an approach that bypasses the need for infrared cancellations in intermediate steps [26,27]. In the strong-coupling limit and at large N c , the EEC in N = 4 SYM can be calculated using AdS/CFT duality [28].
Despite all of this progress, the analytic computation of the EEC at NLO in QCD has remained an open problem, whose solution is desirable for several reasons. First, the analytical results can settle any remaining discrepancies between different numerical methods, and provide a benchmark for future numerical evaluations. Second, the analytical results allow extraction of the O(α 2 s ) asymptotic behavior in the collinear and back-to-back regions, not just at leading power, but any desired power. Knowledge of the subleading power corrections can be very helpful for improving the understanding of resummation at subleading power [29][30][31][32][33][34][35][36][37][38]. Third, no other event-shape variable has been computed analytically at NLO. Calculationally, the EEC appears to be the simplest such observable. Knowing it analytically at NLO marks an important step in the perturbative understanding of event-shape observables, and may pave the way for an analytic computation at NNLO. Recently, progress has been made toward computing the EEC at NLO by linearizing the measurement function [39]. In this letter, we present the first fully analytic result for the EEC in QCD at NLO. The Calculation. At LO, calculation of the EEC is straightforward, because only finite phase-space integrals need to be evaluated. At NLO, the renormalized virtual corrections contain explicit infrared (IR) poles, but no singularities from the boundary of phase space. We use the analytical one-loop amplitudes [40,41], and perform the phase-space integral directly. The real radiative corrections represent the most complicated part of this calculation, because the phase-space integrals contain unresolved soft and collinear IR divergences. We apply reverse unitarity [42,43] to write on-shell delta functions as differences of Feynman propagators with opposite signs for iε, which allows the use of integration-by-parts (IBP) equations [44,45] for multi-loop integrals. The EEC measurement function can be written in the same way, where While the application of reverse unitarity to phase-space integrals is now quite standard, Eq. (4) is special in the sense that M ij (χ) is a non-linear function of Lorentz dot products. In addition to the usual IBP equations, an extra equation, for k = 1, 2, . . ., with [δ(M ij (χ))] 0 ≡ 0, has to be added in order to fully reduce the phase-space integrals to master integrals (MIs). In our calculation, we use Qgraf [46] to generate the squared amplitudes for the LO and real NLO terms. We set all quark masses to zero, and ignore contributions from the top quark, as well as the (tiny) purely axialvector contributions in the case of e + e − annihilation via the Z boson. The color and Dirac algebra is evaluated using Form [47]. The resulting tree-level matrix elements agree fully with Ref. [40]. The squared matrix elements for the NLO real corrections, ignoring the EEC measurement function, can be divided into three integrand topologies, each consisting of nine Feynman propagators (one in the numerator). Since there are four partons in the real NLO final state, there are ( 4 2 ) = 6 different measured pairs to sum over for the EEC. Multiplying the 3 inclusive integrand topologies by the 6 pairs of measurement delta functions gives rise to 18 separate integral topologies. We use LiteRed [48,49] to generate the standard IBP equations for these integral families, and then add the additional integral relation (5) manually. We then export the resulting IBP relations to Fire [50,51] to perform the integral reduction, which leads to a total of 40 independent MIs.
We solve for the MIs by the method of differential equations (DEs) [52,53], and convert the DE systems into a canonical form [54]. Some of the DE systems can be converted to canonical form using the original variable z; for others, an algebraic change of variable to x = √ z or After identifying the appropriate variable for each integral family, the conversion to a canonical basis can be automated by the Mathematica package Fuchsia [55]. The resulting symbol alphabet, characterizing the arguments of the polylogarithms, is {1 − x , y , 1 − y , 1 + y}. Note that z, 1 − z, x and 1 + x also appear, but are not multiplicatively independent, since 1/(1 − y 2 ) = 1 − x 2 = 1 − z, etc., so we do not count them as separate symbol letters. This alphabet implies that the solution to the DEs can be written fully in terms of harmonic polylogarithms (HPLs) [56], which can be manipulated conveniently using the Mathematica package HPL [57]. Our final NLO result contains at most weight 3 HPLs, which can all be reduced to classical polylogarithms.
The most intricate part of the calculation is the determination of the constants of integration for the DEs, which requires combining several different constraints. First, we require that the leading power expansion z α of each MI in the collinear limit z → 0 has the correct power α, which can be predicted by simple power counting. We find that all the MIs in our problem have at most a z −1 pole. (Some have z 0 or z 1 as their leading behavior.) Requiring the absence of z −2 or worse poles strongly constrains the boundary constants. The second constraint is the z → −∞ limit: Before converting to the canonical basis, MIs that are pure functions of uniform transcendental weight should vanish in this limit. The third constraint comes from performing a weighted integral over z, which allows the removal of the measurement constraint, according to the integral relation Here dP S (4) is the four-particle Lorentz-invariant phasespace measure in D dimensions,ẑ ij = Q 2 p i · p j (2 p i · Q p j · Q) −1 , and I({p}) denotes a MI integrand. We choose the integers n and m to be sufficiently positive that the particular integral over z converges, and n ≤ 1, m ≤ 1 to keep the IBP reduction tractable. The integral on the left-hand side can be reduced to known inclusive fourparticle phase-space integrals [58], if we multiply the integrand on both sides by (p i · Q p j · Q) n+m . The last constraint we apply is to demand that the full NLO real corrections, after substituting in the results for the MIs, have at most a z −1 pole. This gives extra constraints, beyond the constraints applied to the individual MIs. A similar method has been applied to fix constants of integration for DEs for auxiliary EEC MIs [39]. The result. After combining the real and virtual corrections, and adding the counterterm to renormalize α s , we obtain our final result for the EEC at NLO. We write the differential distribution as where the LO coefficient A(z) has already been given in Eq. (3), and We have calculated each coefficient in the color decomposition analytically. The leading-color correction B lc reads, 3 − where the g Note that B lc contains explicit dependence on √ z through the function g 3 and its coefficient, whose product is even under √ z → − √ z. This property also holds in N = 4 SYM [25]. To describe B lc , we need just two weight 1, four weight 2, and three weight 3 transcendental functions. To express B nlc and B N f requires two more weight 3 transcendental functions. 1 Individual virtual and real contributions contain HPLs with argument However, they cancel out in the final physical result. The explicit expressions for B nlc and B N f can be found in the supplemental material for this letter. In an ancillary file, we provide computer-readable expressions for all these functions, as well as their behavior in various limits.
We have performed a number of checks on the results. First, the individual virtual and real corrections are IR divergent, but the divergent terms cancel after summing virtual and real, as required for any IR-safe observable. Second, in Fig. 1 we compare our analytical results with numerical predictions from Event2, which is based on the dipole subtraction method [13,14]. We find excellent agreement with Event2 over a large range; the apparent discrepancy in the rightmost bin is mainly due to the finite bin width used in Event2. The z → 0 and z → 1 limits of the analytical results are in perfect agreement with those predicted respectively by jet calculus [11,19] and soft-gluon resummation [22,23], as we discuss in the next section.  [13,14]. The Event2 prediction is obtained after sampling over 10 billion points, with the internal CUTOFF set to 10 −14 . Error bars represent Event2 statistical uncertainties.
Discussion. It is interesting to study the end-point asymptotic limits of the EEC, which provide useful information for resummation and for constructing more accurate parton showers. Expanding our results in the z → 0 limit gives where we have expanded through O(z 0 ). Note that individual terms in Eq. (9) are far more singular as z → 0 than is the total (11). The EEC in the z → 0 limit is dominated by collinear splitting. The leading-logarithmic term log(z)/z has been predicted [11,19] using jet calculus [18,19]. The result is expressed as a product of two 2 × 2 (quark-gluon) anomalous dimension matrices for twist 2, spin 3 operators, plus a contribution due to the running coupling. It agrees fully with the coefficient of log(z)/z in Eq. (11).
In the back-to-back limit, z → 1, we find that the expansion of B(z) to next-to-leading power reads All the terms enhanced by (1 − z) −1 were predicted previously [22], in full agreement with Eq. (12). The nextto-leading power terms are new. They will provide useful information for resumming large Sudakov logarithms beyond leading power [29][30][31][32][33][34][35][36][37][38]. We note the appearance of ζ 2 log(2) in the constant term at next-to-leading power, which originates solely from B nlc . Summary. We have presented the analytical result for the EEC in QCD at NLO. Our calculation was enabled by using the IBP equations in a novel way. The final result turns out to be rather simple; only 11 transcendental functions are required to describe the QCD results, and these functions are no more complicated than the ones in the N = 4 SYM result [25]. In contrast, the polynomial prefactors are of considerably higher degree for QCD. We have checked our results against Event2 numerically and found full agreement. We have also expanded the EEC to next-to-leading power in the collinear and back-to-back limits. The simplicity of the full NLO result provides encouragement for trying to compute the EEC at NNLO analytically. It will also be interesting to apply our method to other event-shape variables, such as the C parameter (which does appear to require elliptic functions, even at LO) [40]. We thank Marc Schreiber for extensive contributions at the beginning of this project, Alexander Smirnov for helpful instruction on the use of

SUPPLEMENTAL MATERIAL
In this supplemental material, we provide detailed analytic formulae for the remaining coefficients in the color decomposition of the EEC at NLO, as well as for an identical-quark piece that cannot be isolated by color alone. We describe the z → −∞ limit of the various color components, and provide the z → 0, 1 and −∞ limits of the identical-quark piece. The z → 0 limit of this piece can be interpreted using jet calculus at next-to-leading logarithm. Finally, we validate our analytic results for each color component numerically against Event2.

Analytical results for remaining color coefficients
The most complicated part of this calculation involves the real corrections. We show representative cut diagrams for real corrections with final states qqgg (Fig. 2a), qqq q (Fig. 2b), and qqqq ( Fig. 2c and 2d). Note that the class of diagrams in Fig. 2d vanishes for a virtual photon thanks to Furry's theorem, because there are three vector couplings attached to each fermion "loop" and the measurement function is invariant under charge conjugation [40]. (At the Z pole, the purely axial-vector contributions represented by Fig. 2d vanish for the first two generations in the massless quark limit, while the third generation contribution is tiny [59].) The subleading color coefficient B nlc appearing in Eq. (8) is given by where we have introduced two more weight 3 transcendental functions, The contribution to Eq. (8) that is proportional to the number of light quark flavors, N f , is given by 3 + 8z 3 − 66z 2 + 71z + 7 60(1 − z)z 5 g 4 + 2 50z 4 − 100z 3 + 66z 2 − 16z + 1 g 1 .
It is straightforward to take the limits of B lc , B nlc and B N f as z → 0 and z → 1, in order to obtain Eqs. (11) and (12), as well as further subleading powers if desired.
Another limit that we can study with the exact NLO result is z → ∞, where the limit can be taken in any radial direction on the complex plane. Although this limit requires an analytic continuation out of the physical region, the analytic properties may still prove useful for understanding the EEC as a limit of a four-point correlator [25][26][27][28], and perhaps for constructing it at higher perturbative orders. In fact, we find that the LO result and the NLO coefficients are quite suppressed in this limit. In the NLO case, the suppressed behavior requires cancellations between many different terms in the exact result. The leading-power terms as z → ∞ have the form, B nlc (z) = 1 z 3 Here the branch cut is chosen along the negative real axis for both log(−z) and √ −z. In the the remainder of this subsection, we discuss the contribution from identical-quark exchange terms of the type shown in Fig. 2c. In the real corrections, both the qqgg final state and the qqqq final state with identical quarks contribute to the subleading color coefficient. We can write the real corrections to B nlc as There are no virtual corrections to the identical-quark interference terms at this order. Hence B qqint is by itself IR finite and gauge invariant, although it is just one piece of the color component B nlc . It has the following form, . (21a) The one additional weight 2 and four additional weight 3 transcendental functions introduced to describe B qqint are defined as is also known as the Bloch-Wigner function [60]. It is single-valued (real analytic) in z = ir,z = −ir; here we only need it on the imaginary z axis. Although g (2) 5 and g 9 contain logarithms and classical polylogarithms with complex arguments, they are real for z ∈ (0, 1). While g (3) 9 is even under r → −r, g 5 is odd under this transformation. On the other hand, the square roots in the prefactor for g (2) 5 in Eq. (21a) ensure that only integer powers of z and 1 − z appear in the expansions around these limits.
Interestingly, the identical-quark interference terms are more complicated than the final results for B nlc in Eq. (13a), in the sense that five more transcendental functions are needed to fully describe them, including g with their complex polylogarithmic arguments. In the full NLO result, the coefficients of these functions cancel between the qqgg cuts and the qqqq cuts.
It also worth pointing out that while the identical-quark interference contribution is IR finite by itself, it does have a leading-power singularity (but not leading-log) in both the collinear and back-to-back regions. Expanding the result (21a), we find the asymptotic z → 0 limit to next-to-leading power (NLP) is given by Remarkably, the coefficient of 1/z in Eq. (23a) can be reproduced using the jet calculus approach at next-to-leadinglogarithm [61]. While the leading-log log(z)/z contribution in the full NLO QCD result comes from iterating two 1 → 2 splittings, q → qg and g → gg, and accounting for running-coupling effects [11], the next-to-leading logarithms also involve the 1 → 3 splitting process. The identical-quark exchange terms in the triple-collinear splitting process q → qqq are somewhat unique in not requiring the subtraction of any iterated 1 → 2 splittings. Assigning longitudinal momentum fractions x 1 , x 2 and x 3 = 1 − x 1 − x 2 to the three quarks, these exchange terms can be described by the function given in Eq. (B.7) of Ref. [61]. The relative transverse momenta have already been integrated over to obtain E(x 1 , x 2 ). Assuming that all angles between the three final-state partons are comparable, we multiply E(x 1 , x 2 ) by x 1 x 2 +x 2 x 3 + x 3 x 1 in order to account for the energy weighting in the EEC, and integrate over the remaining longitudinal momentum fractions, Accounting for the constant prefactors in Ref. [61], we reproduce the 1/z term in Eq. (23a). For other components of B(z), in which there is also a log(z)/z term, the 1/z term is subleading and the integration of the analogous functions in Ref. [61] does not reproduce the correct result, presumably because the assumption of comparable angles between the three partons fails. It may still be possible to extract the 1/z term in this case by using the full dependence of the triple-collinear splitting process on the two-particle invariants p i · p j [62,63].
In the opposite, back-to-back region, the expansion of B qqint to NLP is given by In the limit as z → ∞, the identical-quark exchange terms have the same falloff as the other NLO color components:

Validation of individual color components
For validating the color components of the NLO result, we use a slightly different color basis. We decompose the NLO coefficient B in terms of the following color structures: Each of these coefficients can be computed separately using Event2. In Figs. 3a, 3b, 3c we compare our analytic results with numerical results from Event2. We sample Event2 at over 10 billion points, and set its internal CUTOFF to 10 −14 . Overall we find excellent agreement, except perhaps for B C F in the range −0.5 < cos χ < 0, where large statistical fluctuations are seen in the bottom panel of Fig 3a. Numerical evaluation of B C F in this region is challenging, as there are large cancellations between the virtual and real corrections. A close-up of this region is shown in Fig. 4. Finally, in Fig. 5 we plot the identical-quark interference terms, B qqint .