The Energy-Energy Correlation in the back-to-back limit at N$^3$LO and N$^3$LL$^\prime$

We present the analytic formula for the Energy-Energy Correlation (EEC) in electron-positron annihilation computed in perturbative QCD to next-to-next-to-next-to-leading order (N$^3$LO) in the back-to-back limit. In particular, we consider the EEC arising from the annihilation of an electron-positron pair into a virtual photon as well as a Higgs boson and their subsequent inclusive decay into hadrons. Our computation is based on a factorization theorem of the EEC formulated within Soft-Collinear Effective Theory (SCET) for the back-to-back limit. We obtain the last missing ingredient for our computation - the jet function - from a recent calculation of the transverse-momentum dependent fragmentation function (TMDFF) at N$^3$LO. We combine the newly obtained N$^3$LO jet function with the well known hard and soft function to predict the EEC in the back-to-back limit. The leading transcendental contribution of our analytic formula agrees with previously obtained results in $\mathcal{N} = 4$ supersymmetric Yang-Mills theory. We obtain the $N=2$ Mellin moment of the bulk region of the EEC using momentum sum rules. Finally, we obtain the first resummation of the EEC in the back-to-back limit at N$^3$LL$^\prime$ accuracy, resulting in a factor of $\sim 4$ reduction of uncertainties in the peak region compared to N$^3$LL predictions.


Introduction
In the quest of understanding the nature of the strong interaction, the study of QCD radiation produced in high-energy electron-positron collisions provides a powerful lens into the behavior of Quantum Chromodynamics (QCD). Among the observables that one can consider to perform precision studies of QCD radiation, the Energy-Energy Correlation (EEC) [1] stands out for its simplicity, and has been an important benchmark observable to test QCD and to extract the strong coupling constant α s , which has been explored both at LEP and at SLAC [2][3][4][5][6][7][8][9]. The EEC is an e + e − event shape to measure the energy-weighted angular distance between any pair of particles in an event [1], and as such is one of the earliest examples of an infrared and collinear (IRC) safe observable. It is defined as The differential cross section is distribution-valued in the small-angle limit (sometimes also referred to as collinear or forward limit) χ → 0 and in the back-to-back limit χ → π. By expressing the EEC in terms of z the singular points of the distributions are mapped onto z → 0 and z → 1, respectively. The differential cross section at Born level is given by where α s is the strong coupling constant. In this article we will direct attention towards the singular limits of the EEC. Consequently, it is useful to split the cross section into three different contributions.
where dσ 0 dz and dσ 1 dz contain all terms that behave as 1/z and 1/(1 − z), respectively. More precisely, they are expressed in terms of (plus) distributions to regulate these divergences, while dσreg.
dz is a regular function of z that is holomorphic in the entire unit interval, z ∈ [0, 1]. The EEC can also be expressed in the variable (1.5) The gluon and photon-induced EEC was computed through order α 2 s in refs. [14,15] in terms of classical polylogarithms as a function of z. Here, we note that expressing the differential cross section in terms of the variable x allows us to represent it in terms of harmonic polylogarithms [33] with argument x and indices {−1, 0, 1}. We provide analytic formulae for the EEC expressed in x including all distribution valued terms through α s as ancillary files together with the arXiv submission of this article. The recent computation of the EEC at O(α 3 s ) in N = 4 sYM in ref. [19] finds evidence that at this order the analytic formula contains elliptic functions and is no longer expressible in terms of HPLs.
At Born level the EEC vanishes for z = 0, 1. Consequently, in the context of the fixed-order calculations reported above it is customary not to include the distributional behavior as z → 0, 1 and count O(α s ) as the leading order (LO). Accordingly O(α n+1 s ) contributions are counted as N n LO. Such fixed-order calculations become unreliable in the singular limits z → 0 and z → 1, where large logarithms ln(z) and ln(1 − z) can spoil perturbative convergence. In these limits the EEC needs to be resummed to all orders in perturbation theory to retain predictive power.
In the forward limit, this resummation was performed at leading logarithmic (LL) accuracy a long time ago [34]. Recently, a factorization theorem was derived in this limit and the resummation was improved in ref. [24]. In the back-to-back limit, the resummation was carried out at next-to-next-to-leading logarithmic (NNLL) accuracy [7,35] based on transverse-momentum dependent (TMD) factorization in e + e − [36][37][38][39]. Recently, ref. [40] proofed that the all-order factorization for the EEC in the back-to-back limit indeed follows from TMD factorization, and presented first results at N 3 LL acurracy. In N = 4 sYM, the factorization of the EEC in both its forward and back-to-back limit has also been explored up to four loops by using the operator product expansion for light-ray operators [41], and by relating the EEC to four-point correlation functions of conserved currents [23]. Note that these factorization theorems contain the full distributional structures, and thus their fixed-order expansions start at O(α 0 s ). Hence, in the context of factorization theorems one counts O(α n s ) as N n LO accuracy. In this paper we calculate the full singular structure of the EEC in QCD in the back-toback limit at N 3 LO, i.e. O(α 3 s ). This is achieved by calculating the EEC jet function at the same order, which is the only unknown ingredient of the factorization theorem derived in ref. [40]. As an application, we resum the EEC in the back-to-back limit at N 3 LL accuracy, i.e. N 3 LL resummation combined with the full N 3 LO fixed-order boundary condition. This constitutes the highest level of accuracy obtained for an event shape sensitive to QCD radiation to date. We show that thanks to the resummation of large logarithms up to N 3 LL and the inclusion of fixed-order boundary terms at N 3 LO, we achieve a ∼ 4-fold reduction of uncertainty compared to previous results obtained at lower accuracy. We thoroughly discuss different schemes to estimate the uncertainties due to missing higher order corrections both in the boundary terms as well as in the anomalous dimensions and compare them with the ones used in the literature for this observable. Adopting a scheme in line with the ones previously used in the literature we obtain a 0.5% uncertainty at the peak. Using a more conservative scheme, which also estimates uncertainties from soft physics, we obtain a 4% uncertainty for our result at N 3 LL . Furthermore, using a momentum sum rule we also analytically calculate the N = 2 Mellin moment of the regular part at N 3 LO, which will be an important check of the full result once it becomes available. The jet functions calculated in this work are necessary ingredients to describe the singular behavior of the TEEC in the large angle limit at O(α 3 s ) as well as to the resummation of large logs at N 3 LL and N 4 LL accuracy for this observable.
This paper is organized as follows. In section 2 we briefly review the factorization theorem for the EEC in the back-to-back and show how to extend it to gluon-induced Higgs decays, clarifying its nontrivial helicity structure. We also obtain the three-loop jet functions from our results for the transverse-momentum dependent fragmentation functions (TMDFFs) calculated in the companion paper [42], 2 and present the full singular structure of the EEC at N 3 LO in QCD in the back-to-back limit. In addition, we validate the conjecture that in the back-to-back limit, the leading transcendental terms in QCD match those in N = 4 supersymmetric Yang-Mills ref. [23]. In section 3, we exploit the fact that the EEC obeys a set of sum rules to analytically obtain the N = 2 Mellin moment of the bulk of the EEC distribution at NNLO in QCD. In section 4, we carry out the resummation of the EEC at N 3 LL accuracy to illustrate the improved perturbative accuracy compared to previous results. We conclude in section 5.
2 The EEC in the back-to-back limit 2.1 Factorization of the EEC in the back-to-back limit In this section, we briefly review the factorization of the EEC in the back-to-back limit z → 1. In the case of electron-positron annihilation, i.e. quark-initiated EEC, this was first derived by Collins and Soper [36,37] and Kodaira and Trentadue [38,39]. Recently, it was also formulated in ref. [40] using Soft-Collinear Effective Theory (SCET) [44][45][46][47], which clarified the role of nontrivial fixed-order boundary terms in the factorization formula. So far, no details have been given in the literature on the Higgs EEC. To fill this gap, we briefly review the derivation of the quark EEC in section 2.1.1, and then show how the same strategy can be applied to the gluon EEC in section 2.1.2, where additional complications arise from a nontrivial Lorentz structure.

Quark EEC
In the back-to-back limit, the EEC can be related to identified hadron production, e + e − → h 1 h 2 + X, at small transverse momentum q T of the dihadron system [36,37,40], Here, z 1,2 = (P 1,2 · q)/q 2 describe the longitudinal momenta carried away by the hadrons, with the total momentum transfer given by q µ = p µ e + + p µ e − . 3 The factorization of dihadron production at small q T was derived in seminal works by Collins and Soper [36,37] (see also ref. [48]), and can be written as Here, summation over all quark flavors q is kept implicit, the hard function H qq encodes virtual corrections to Born process e + e − → qq,D h/q is the fragmentation function encoding the probability to obtain the hadron h from the fragmentation of a quark q, andS q is the TMD soft function. 4 Note, that in the literature, the soft function is often combined with the fragmentation function asD h/q S q , whereas we keep it explicit. The scale ν is the so-called rapidity renormalization scale, which is closely related to the Collins-Soper scale ζ. Note that the fragmentation functionsD h/q and the soft functionS q both depend on the chosen rapidity renormalization scheme. This scheme dependence cancels in the combinationD h/q S q , and consequently also in the combinationD h 1 /qDh 2 /qSq in eq. (2.2).
For more details on the fragmentation functionD h/q , we refer to ref. [42]. By combining eqs. (2.1) and (2.2), we obtain the EEC factorization theorem in the back-to-back limit as stated in ref. [40], In a frame where the hadrons are aligned along back-to-back lightlike directions, i.e. p µ 1 = p − 1 n µ /2 and p µ 2 = p + 2n µ /2 with n 2 =n 2 = 0 and n ·n = 2, these evaluate to z1 = p − 1 /q − and z2 = p + 2 /q + , up to corrections suppressed by qT . In this frame, qT is the transverse component of the momentum transfer q µ . 4 The TMD soft function is universal between e + e − and pp processes [49], i.e. it is independent of the direction of the contained Wilson lines, and thus we do not further distinguish these. For a detailed discussion of the equivalence of the TMD soft function in the context of the EEC, see ref. [40,50].
where H qq andS q are the same hard and soft functions as in eq. (2.2), J 0 (x) is the 0-th Bessel function of the first kind, and the jet functions J q are defined as the first moments of the fragmentation functions, The jet function J q has the same dependence on the rapidity renormalization scheme as the fragmentation functionD h/q . Similar to eq. (2.2), this scheme dependence cancels in the combination J q JqS q that arises in eq. (2.3). In analogy to TMD factorization, where it is often customary to absorb the soft function in the TMDPDF or TMDFF, one can construct a manifestly rapidity-regulator independent jet function as where the Collins-Soper scale ζ 2 = Q as in TMD factorization. Using eq. (2.5), eq. (2.3) can be equivalently written as For more details on this construction in the context of TMDPDFs and TMDFFs, see refs. [42,51]. In the following, we will restrict our discussion to J q , as the corresponding results for J q can be trivially obtained using eq. (2.5).
To simplify the jet functions, we use that for perturbative b T Λ −1 QCD they can be perturbatively matched onto the collinear fragmentation functions d h/q as whereK qq is a perturbative matching kernel. Applying eq. (2.7) to eq. (2.4), we obtain where we used the momentum sum rule of the fragmentation function, Thus, the EEC jet function is free from nonperturbative hadronic matrix elements, which makes the EEC much less susceptible to nonperturbative effects than the q T distribution itself. We note that perturbative power corrections to eq. (2.3) can be systematically studied using the operator formalism of SCET [22,52,53] and involve the treatment of rapidity divergences beyond leading power [54,55]. Nonperturbative power corrections to the EEC have been explored in refs. [56,57].
The factorization for the EEC in the back-to-back limit has been explored in the literature since a long time. In refs. [36,37], Collins and Soper derived the q T factorization theorem in eq. (2.2), which albeit using a different notation already contained hard and transverse-momentum dependent fragmentation functions, with the soft function absorbed into the TMDFFs. They also showed how to resum large logarithms by solving evolution equations in the unphysical scales µ and ν, which we will discuss in detail in section 4. Using the relation between z and small q T , they also obtained a resummed formula for the EEC. However, it was only provided at LO in the matching, where the hard, jet and soft functions all evaluate to unity. Similarly, Kodaira and Trentadue presented q T factorization with TMDFFs that are matched onto collinear FFs, but only provided formulas for the resummed EEC spectrum valid at LL and NLL, where the hard and jet functions again evaluate to unity [38,39]. Nontrivial fixed-order terms in the factorization formula were first pointed out in ref. [35], which however did not separate between hard, jet and soft functions, such that their hard function is equal to the product of H qq J q JqS q when evaluated at fixed order. 5 However, for the purpose of resummation, it is important to distinguish the hard function, which describes physics at the high scale µ ∼ Q, from the jet and soft functions which describe physics at the low scale µ ∼ 1/b T . (In the context of TMD factorization, the TMDFFs and soft function are often combined.) This will be addressed in section 4.4. This separation was first achieved in ref. [40], and we closely follow their conventions.

Higgs EEC
We now discuss the extension of the EEC factorization to gluon-induced processes, e.g. the Higgs EEC in e + e − → H → gg + X, which has not yet been given explicitly in the literature. As in the quark case, one relates the EEC in the back-to-back limit to TMD factorization for identified hadron production, e + e − → h 1 h 2 + X, which for gluon-induced processes reads To the best of our knowledge, this formula has not yet been explicitly given in the literature, but it follows immediately from the similar structure of TMD factorization at hadron colliders [58][59][60]. Similar to eq. (2.2), H,D h/g andS g denote the hard, fragmentation and soft function, respectively. For more details on the definition of the gluon fragmentation function, we refer to ref. [42].
The key difference between eqs. (2.2) and (2.10) is the Lorentz structure of the hard function and the fragmentation functions, reflecting the transverse polarization of the fragmenting gluons in the factorization limit. Since the gluon fragmentation functionsD ρσ h/g only depend on one Lorentz vector, namely b µ ⊥ = (0, b T , 0), their most general decomposi-tion is given byD where for brevity we suppressed the scales.
In this work, we will only be interested in the Higgs-initiated gluon EEC. In this case, the scalar nature of the Higgs boson implies a trivial Lorentz structure of the hard function, Inserting this into eq. (2.10) and using eq. (2.11), we obtain i.e. we simply encounter the sum of two factorized expressions, one for the polarizationindependent contribution and one for the polarization-dependent contribution. (Note the factor 1/2 from contracting eq. (2.11) with itself.) We can now apply the same steps as in section 2.1.1 to obtain the factorization formula for the Higgs EEC as (2.14) As before, the gluon jet functions are related to the matching coefficients of the gluon fragmentation functions, whereK gi andK gi are the matching kernels ofD h/g andD h/g onto the gluon fragmentation function d h/g , with the matching taking the same form as eq. (2.7). Since the polarized TMDFF K gi starts at O(α s ), the same holds for the polarized jet function J g , implying that it contributes to eq. (2.14) starting at O(α 2 s ). As in the quark case, the gluon jet functions J g and J g depend on the chosen rapidity regulator, and manifestly regulator-independent jet functions can be constructed as in eq. (2.5), see section 2.1.1 for more details.
Note that our factorization theorem for the Higgs EEC disagrees with the statement in ref. [61], where the gluon jet function is claimed to be the linear combination J g + J g . In this case, the J g would already contribute to the cross section at O(α s ). However, by comparing our results for the Higgs EEC in the back-to-back limit with the fixed-order calculation of ref. [15], we confirm that this is not the case, and find perfect agreement with the prediction from eq. (2.14).

Jet function for the back-to-back limit at N 3 LO
Using eqs. (2.8) and (2.15), the EEC jet function can be easily obtained from the N 3 LO result for the TMDFF calculated in our companion paper [42]. This result has been obtained using a recently developed method for the expansion of cross sections around the collinear limit [62].
Throughout this paper we will present results computed with the exponential regulator [63]. As a renormalization scheme we employ the one of ref. [58] by relating the regulator τ to the inverse of the rapidity renormalization scale ν. This is standard procedure for calculations in the exponential regulator and further details on it can be found in ref. [63]. The choice of this regulator is due to the fact that it allows to regulate rapidity divergences while only retaining information on the total momentum of the real radiation, which is at the basis of the framework of ref. [62] that we have employed for this calculation. The rapidity-regulator independent combination J i defined in eq. (2.5) can be obtained by combining our results with the N 3 LO soft function computed in refs. [62,64]. While we only discuss results for J q,g in the following, for completeness we also provide J q,g for quarks and gluons in the ancillary files. Using J i , it is trivial to obtain the N 3 LO jet function in any rapidity regularization scheme for which the soft function is known at the same order, which so far is only the case for the exponential regulator used here.
To present our results, we expand the jet functions as where the n-th order coefficient J only depends on the logarithms is fully encoded by the jet function RGEs, which due to its definition in eq. (2.4) are identical to the RGEs of the TMDFF, where the anomalous dimensions have the all-order expressions Here, Γ i cusp (α s ) is the cusp anomalous dimension,γ i J (α s ) is the jet function noncusp anomalous dimension, which is identical to the noncusp anomalous dimension γ i B (α s ) of the TMD beam function, andγ ν (α s ) is the rapidity anomalous dimension. Explicit results for these in our notation are collected in ref. [65]. For more details on these RGEs, see section 4.1.
We have explicitly checked that our result for the jet function obeys eqs. (2.18) and (2.19). Solving these equations order-by-order in α s , one easily obtains the logarithmic structure of the J (n) i (L b , L Q ) in terms of the anomalous dimension and the constant piece of the jet function, which we define as (2.20) In appendix A, we provide the full fixed-order structure of the J (n) i through N 3 LO, and for brevity in the following we only present the constant terms j (n) i through N 3 LO. The complete N 3 LO jet functions are also provided as ancillary files with this submission.

Quark jet function
The quark jet function has already been calculated using the exponential regulator in ref. [61] at NLO and NNLO, with which we find full agreement with their results. For completeness, we repeat these results, The new result of this paper is the three-loop coefficient, which is given by

Gluon jet function
In the gluon case, we have to consider both the polarization-independent jet function J g and the polarization-dependent jet function J g .
At NLO and NNLO, the finite terms for J g can be obtained from the calculation of the unpolarized gluon TMDFFs of ref. [66], The new result at three loops can be obtained from our calculation of the TMDFF in ref. [42]. We obtain The polarization-dependent jet function J g starts at O(α s ), but for Higgs production only interferes with itself, see eq. (2.14). Hence, it is sufficient to know J g at O(α 2 s ) to obtain the Higgs EEC at N 3 LO. For completeness, we report their results from ref. [66],

The EEC in the back-to-back limit at N 3 LO
Combining our result with the hard function from ref. [67] and the soft function from ref. [64], we have all ingredients to obtain the EEC in the z → 1 limit at N 3 LO. The Bessel integral in eqs. (2.3) and (2.14) can be easily evaluated analytically, see appendix B.
To present our results, we expand the EEC in the back-to-back limit as where we have divided out the Born partonic cross sectionσ 0 . In eq. (2.26), L h = ln(Q/µ) is the only remaining logarithm of µ, and its structure is entirely governed by the β function. For brevity, we only report the nonlogarithmic terms (which is the same as the entire result for the choice µ = Q), but provide the full result in an ancillary file with this submission.

Quark EEC
The results through NNLO were already given in ref. [61], with which we fully agree, and which are repeated here for completeness: Here, we introduced the standard plus distributions The new result at N 3 LO reads Here, following the notation of ref. [67] the factor N F,V arises from virtual diagrams where the exchanged vector boson couples to a closed quark loop, and hence this contribution is not proportional to the charges of the Born process. For the simplest case of photon exchange, it is given by N F,V = ( f e f )/e q , where e q is the flavor of the external quark. Note that d abc d abc = (N 2 c − 4)(N 2 c − 1)/N c and N r is the dimension of the fundamental representation, hence N r = 3 in QCD.

Gluon EEC
The singular structure of the Higgs EEC in the back-to-back limit is given by The logarithmic structure matches the fixed order calculation of ref. [15], while to the best of our knowledge, the δ(z) term has not appeared in the literature before. Note that the polarized jet function J g appears for the first time in the δ(z) coefficient at O(α 2 s ), as it vanishes at tree level and interferes only with itself.
For the singular structure at N 3 LO we find (2.31)

The maximally transcendental limit from N = 4 sYM
It is conjectured that the leading transcendental term of the EEC in QCD (with C F , C A → N c ) in the back-to-back limit is given by the EEC in the back-to-back limit in maximally supersymmetric Yang-Mills theory (N = 4 SYM) [18]. This relation has been observed to hold up to O(α 2 s ) [14,15,24,61] and it is known not to hold beyond leading power [19,22]. Here, we validate this conjecture at O(α 3 s ), which also provides us with an independent check of the leading transcendental limit of the three-loop jet functions.
In N = 4 SYM, the back-to-back limit of the EEC is given by a (simplified) conformal version of the factorization formula in eq. (2.3) [18,23,41], Here,z ≡ 1−z, and Γ cusp (a) and Γ(a) are the cusp and the collinear anomalous dimensions [68][69][70][71][72]. In ref. [23] it has been shown that the boundary function H(a) can be obtained from the OPE coefficients of twist-two operators at large spin up to 3 loops [73,74] and it reads .
If the principle of maximal transcendentality holds at this order, eq. (2.34) should predict the leading transcendental terms of the δ(1 − z) coefficient of the EEC in QCD through O(α 3 s ). 6 We find perfect agreement between eq. (2.34) and the leading transcendental terms of the quark and gluon EEC in eqs. (2.29) and (2.31), confirming that the conjectured principle of maximal transcendentality holds for the EEC in the back-to-back limit through O(α 3 s ). Conversely, if we instead assume that eq. (2.34) predicts the maximal transcendental limit of the EEC in QCD, we can use it to predict the leading transcendental term of the EEC jet function at O(α 3 s ), as the required hard and soft functions are already known at this order [64,67]. We obtain which precisely agrees with the leading transcendental limit of the quark and gluon jet functions given in eqs. (2.22) and (2.24). Since the jet functions are obtained from a weighted integral of the corresponding TMDFF, summed over all contributing partonic channels, this also provides a remarkable cross check on the N 3 LO TMDFF matching kernels calculated in the companion paper [42].
3 Sum rules and the N = 2 Mellin moment of the EEC at O(α 3 s ) An interesting property of the EEC is that it obeys the sum rules [23,24,41] 1 where σ is the inclusive cross section for e + e − → hadrons (or Higgs decay to hadrons), which is known up to O(α 4 s ) [75][76][77]. These sum rules are due to energy and momentum conservation, respectively. The first sum rule can be derived by noting that The second follows from where i p µ i = q µ = (Q, 0, 0, 0) in the rest frame of the source. 7 The presence of these sum rules puts interesting constraints on the EEC distribution and allows one to use information about one region to extract nontrivial information on the distribution away from that region [23]. In particular it is useful to divide the EEC distribution in different terms [23,24] This isolates the endpoint singular terms with respect to the regular part of the distribution, such that dσ dz reg. is finite as z → 0 and z → 1 and the plus prescription for φ 0 is taken to act at z = 0, while that for φ 1 at z = 1. The functional form of the leading term of the 7 Note that the choice of frame is already implicitly taken in defining the EEC observable in terms of the angle χ a,b which is clearly frame dependent. An alternative, and equivalent, way of defining the observable in a frame independent way is to let z = pa·p b 2pa·qp b ·q . It is easy to check that, in the rest frame of the source, this gives back the definition in eq. (1.2). EEC in the z → 0 and z → 1 limit in QCD is known at all orders due to the factorization theorems presented in ref. [24] and ref. [40]. In particular, order by order in α s one finds where all the coefficients c (k) n,m with n ≤ 3 are known [24,40]. In particular, the coefficient c is the coefficient of L m (z) in eqs. (2.29) and (2.31) for e + e − and Higgs, respectively. Expanding the total cross section and the boundary constants as we can make use of the sum rules and obtain relations order by order in α s . We can write the relations at O(α n s ) in compact form as . (3.11) In appendix C we collect the expressions for the terms entering eq. (3.9), eq. (3.10) and eq. (3.11). The different ingredients can be obtained with wildly different techniques and, as it is often the case, some are easier to obtain than others. Therefore, the power of these sum rules lies in the ability of extracting information about the endpoints from the bulk of the distribution and vice versa. This idea has been applied in ref. [23] to obtain the endpoint behaviors of the EEC in N = 4 up to three loops. 8 In QCD, it has been applied in ref. [24] to extract V 0 , i.e. the two loop δ(z) coefficient, using the two loop δ(1−z) coefficient V coefficients, is known [24,40], but neither the regular part of the distribution nor the δ(z), δ(1 − z) terms are known. In this work we have calculated the EEC jet function at O(α 3 s ), which is the last missing ingredient to obtain V 1 , i.e. the three-loop coefficient of δ(1−z), which we have presented in eqs. (2.29) and (2.31) for e + e − and Higgs, respectively. 8 The contact terms of the EEC in N = 4 were also obtained in ref. [41].
Combining our new results with eq. (3.10), we obtain the N = 2 Mellin moment of the EEC in the bulk at O(α 3 s ). For e + e − we obtain where we adopt the same convention for the d abc d abc color structure as in eq. (2.29).
We stress that to calculate the N = 2 Mellin moment using eq. (3.10) we need the O(α 3 s ) logarithmic coefficients also in the forward (z → 0) limit, whose factorization theorem was obtained in ref. [24]. However, at the time of the publication of ref. [24], the NNLO timelike splitting function P (2) qg available in the literature was not correct. A new result for P (2) qg was more recently obtained in ref. [27] and we have confirmed this result in our calculation in ref. [42]. Since P (2) qg is part of the singlet time-like splitting kernel matrix governing the logarithmic structure of the EEC in the z → 0 limit [24], the new result for P (2) qg modifies the small angle limit of the Higgs EEC at O(α 3 s ) and beyond. The result in e + e − is not modified since in that case P (2) qg contributes to the z → 0 logarithmic coefficients only starting at O(α 4 s ). In the Higgs case, using the correct splitting function we obtain The N = 2 Mellin moment of the EEC for e + e − and Higgs at O(α 3 s ) in eqs. (3.12) and (3.13) constitute two new results. In particular it is the first piece of information on the EEC in QCD analytically at NNLO in the bulk of the distribution. As a matter of fact, at NNLO the EEC is only known numerically in QCD for e + e − [7,10] while no results at this order at all are available for the EEC in gluon-initiated Higgs decay. Noting that the full result for the EEC at this order is known in N = 4 [19], we can check if the N = 4 result indeed constitutes the leading transcendental terms of the QCD result. While the result of ref. [19] is expressed in terms of a two-fold integral and therefore cannot be directly checked analytically, in ref. [23] the N = 2 Mellin moment of the distribution away from the endpoints has been obtained and reads (3.14) We see that this result has no uniform transcendentality, but the leading transcendental piece exactly matches our result in QCD after taking C A , C F → N c . Therefore we confirm that the principle of maximal transcendentality holds also at NNLO for the N = 2 Mellin moment of the distribution in the bulk.

Resummation of the EEC at N 3 LL accuracy
In this section, we use our results to obtain for the first time the EEC spectrum in the back-to-back limit resummed at N 3 LL accuracy. In section 4.1, we review how the large logarithms ln(1 − z) can be resummed to all-orders by solving the renormalization group equations (RGEs) of the hard, jet and soft functions. Details of its implementation are presented in section 4.2, before we show our numeric results in section 4.3.

Renormalization group evolution
The hard, jet and soft functions entering eqs. (2.3) and (2.14) obey the same RGEs as the hard, beam and soft functions in TMD factorization, The anomalous dimensions with respect to µ have the all-order form Here, Γ i cusp (α s ) is the cusp anomalous dimension, γ q (α s ) and γ g (α s ) are the quark and gluon anomalous dimensions, andγ i J (α s ) andγ i S (α s ) are the jet and soft noncusp anomalous dimensions, respectively. In eq. (4.2), the notation γ[α s (µ)] indicates that their scale dependence only arises through the strong coupling constant α s (µ), and we distinguish the noncusp anomalous dimensions from the full anomalous dimension by the number of arguments. Note that the jet anomalous dimension is identical to that of the TMD beam and fragmentation functions, and in particularγ i J (α s ) ≡γ i B (α s ). µ independence of the cross section implies that The rapidity anomalous dimension itself obeys an RGE, which follows by commutativity of applying both d/d ln µ and d/d ln ν to either J i orS i . The solution to eq. (4.4) can be written as where the integral resums logarithms ln(µ/µ 0 ), leaving only logarithms ln(b T µ 0 /b 0 ) in the boundary term that is evaluated at fixed order, as indicated. This generalizes eq. (2.19) to arbitrary boundary scales µ 0 . Eq. (2.19) is obtained by choosing the canonical scale µ 0 = b 0 /b T , which eliminates all large logarithms in the boundary term such that it can be reliably calculated in fixed order. The boundary term at this particular choice is commonly abbreviated asγ As for the µ anomalous dimensions, the boundary term is distinguished from the full γ i ν (b T , µ) by the number of arguments. We note that the rapidity anomalous dimension becomes nonperturbative for b T Λ −1 QCD , irrespective of whether µ is perturbative or not. In our numeric analysis, we will simply freeze out µ 0 to avoid the Landau pole, but note that it would be very interesting to extract nonperturbative contributions toγ ν using EEC data from LEP. Another promising approach is to calculate it using lattice QCD as suggested recently [78][79][80], and first exploratory results have recently been obtained in refs. [81,82], see also ref. [83] for first estimates of its large-b T asymptotics.
By solving eq. (4.1), one can evolve the hard, jet and soft functions from their natural scales to the common scales (µ, ν). This two-dimensional evolution is independent of the chosen path by virtue of eq. (4.5), and we choose to first evolve in virtuality µ before evolving in rapidity ν. The resummed cross section in eq. (2.3) is then given by and similar for the gluon case shown in eq. (2.14). Eq. (4.7) is manifestly independent of the overall rapidity scale ν, while the dependence on µ cancels due to eq. (4.3). 9 By choosing the canonical resummation scales as the fixed-order boundary terms in the first line of eq. (4.7) are free of large logarithms and can be reliably evaluated in fixed-order perturbation theory, while all large logarithms are explicitly exponentiated. The logarithmic accuracy of the resummed cross section in eq. (4.7) is classified by the exponentiated logarithms. For example, LL refers to exponentiating all terms α s L 2 , requiring only the one-loop cusp anomalous dimension and beta function. The fixed-order boundary terms H, J and S must be chosen at an order such that they cancel all logarithms obtained from expanding the exponential in fixed order. In practice, one often chooses the boundary terms at one order higher as this significantly reduces residual scale dependencies, which is referred to as prime counting [86]. For completeness, table 1 lists the ingredients required up to N 3 LL .
At N 3 LL , we need to implement the hard, jet and soft function at three loops. The Drell-Yan hard function is known at O(α 3 s ) from the quark form factor [67,[87][88][89][90][91][92][93][94], and is explicitly given in ref. [67]. 10 The soft function has been calculated at O(α 3 s ) in ref. [64] and confirmed in ref. [51], while the EEC jet function at three loops is the remaining ingredient provided in this paper. The resummation also requires knowledge of the fourloop cusp anomalous dimension [68,70,[95][96][97][98][99][100][101][102] (see ref. [101] for a complete list of partial four-loop results). The three-loop quark and gluon anomalous dimensions γ q,g are known from the corresponding form factors [87][88][89][90][91][92]103] (see ref. [102] for the result at four loops). The beam and soft noncusp anomalous dimensionsγ i B,S were first obtained by consistency with the invariant-mass dependent jet function and threshold soft function, all of which are known at three loops [95,96,[104][105][106][107], and were confirmed by explicit calculations 9 The integrals over µ are often implemented using an approximate analytic solution which leads to a small residual dependence on µ [84,85]. To avoid this effect, we have implemented all integrals and the running coupling constant exactly, but we have checked that the difference to the analytic approximation is negligible. 10 Note that starting at three loops, the exchanged vector boson can couple to closed quark loops, whose couplings differ from those of the Born process. This is the origin of the NF,V piece in ref. [67]. For our numerical illustration later, we simply set NF,V = 0. in refs. [51,64,108]. The rapidity anomalous dimensionγ ν (α s ) is known at N 3 LO [64,109,110]. Finally, resummation at N 3 LL accuracy also requires knowledge of the QCD β function at four loops [111][112][113][114]. Explicit expressions for all these anomalous dimension through N 3 LO in the notation used in this paper are collected in refs. [65,115].

Resummation scales and perturbative uncertainties
We choose the canonical resummation scales as Here, we employ a local b * prescription to freeze out the virtuality scales to avoid the Landau pole at large b T , with The functional form in eq. (4.9) is identical to the b * prescription of refs. [36,37], but following ref. [116] we only modify the resummation scales rather than globally replacing b T by b * T (b T ), as this would induce a global power correction O(b 2 T /b 2 max ). Note that we always choose (variations around) the canonical resummation scales in eq. (4.9). In a detailed phenomenological study, one would smoothly turn off the resummation when the power corrections to eq. (4.7) become comparable to the terms predicted by the factorization theorem. Here, we refrain from doing so, as we only intend to illustrate the impact of the new three loop results on the resummation in the regime where canonical resummation is justified.
To estimate perturbative uncertainties, we follow the procedure developed in ref. [117] and separately consider a fixed-order and a resummation uncertainty, and in addition consider an uncertainty from our nonperturbative prescription. Since these sources are considered uncorrelated, the individual uncertainties are then added in quadrature, which we apply symmetrically around the central prediction.
The fixed-order uncertainty is estimated by varying all scales except µ 0 in eq. (4.9) by a common factor of 1/2 or 2 and take ∆ fo to be the maximum deviation. This probes all fixedorder boundary terms in the first line of eq. (4.7), but does not affect the exponentiated logarithms in the second line, and thus is akin to a standard fixed-order scale variation.
The resummation uncertainty is probed by individually varying all scales by a factor of 1/2 or 2 around their central value, constrained such that the arguments of all exponentiated logarithms in the second line of eq. (4.7) are varied up or down by a factor of at most 2. We do not vary µ H , whose variation is already covered by the fixed-order uncertainty, resulting in 35 variations in total [117]. These variations probe the cancellation of the exponentiated logarithms with those in the fixed-order boundary conditions, and thus are interpreted as a resummation uncertainty. Since these variations can be highly correlated, we define ∆ res as the envelope of all variations.
Finally, we vary b 0 /b max = 0.5, 2 GeV to probe the uncertainty ∆ np of our nonperturbative prescription.

Numerical results
We illustrate the impact of our new results by numerically studying the EEC spectrum differential in the angle χ, for photon-induced hadron production up to N 3 LL . We always work on the Z-pole with Q = m Z = 91.1876 GeV, and choose α s (m Z ) = 0.118 and evolve it according to table 1. The resummed spectra are evaluated using an implementation in SCETlib [118]. Figure 1 shows the breakdown of the total uncertainty into the different sources, namely fixed-order (∆ fo , blue), resummation (∆ res , green) and nonperturbative (∆ np , orange) uncertainties, for all resummation orders except LL. In all cases, the smallest uncertainty is ∆ np from varying b max , and only becomes relevant for the very precise N 3 LL prediction. The resummation uncertainty ∆ res is always larger than ∆ fo , but both are of comparable size at NNLL and higher. As is clear from figure 1, the overall uncertainties reduce drastically with increasing the resummation accuracy, i.e. as NLL → NNLL → N 3 LL.
A striking feature is the huge reduction of uncertainties when going from N n LL (left panel) to N n LL (right panel), which has a much bigger impact than simply increasing the resummation order from N n LL to N n+1 LL. This illustrates the importance of prime counting, i.e. including the fixed-order boundary terms in the resummation at the same perturbative order as the resummation (confer table 1).
We remark that we encounter much larger uncertainties than observed in the NNLL results in refs. [7,8,35]. Our uncertainties are dominated by variations of the low scales, as they probe α s (1/b * T ) and thus become large at large b T . Due to applying these uncertainties symmetryically, the resulting uncertainty bands even become negative up to NNLL, but greatly improved beyond NNLL. We will address this in more detail in section 4.4.  Figure 2 shows a comparison of the resummed EEC spectrum at various orders. In the left panel, we compare NLL through NNLL , while in the right panel we compare NNLL through N 3 LL . In both cases we also plot the nonsingular distribution defined as where dσ 1 /dz is the leading-power cross section in the back-to-back limit as predicted by the factorization formula in eq.  Thus, dσ nons /dz shows the impact of the fixed-order corrections beyond leading power. In the left plot, the dot-dashed black line shows the neglected correction from the NLO fixedorder results calculated in ref. [14] (with respect to the O(α 2 s ) leading power result), while in the right plot we show the nonsingular distribution at O(α 3 s ) using the NNLO fixed-order results calculated numerically in ref. [7]. 11 As expected, these fixed-order corrections are neglibile for χ 160 • , thus justifying employing canonical leading power resummation in the shown range. One notable feature in figure 2 is that going from NLL to NNLL actually increases the size of the scale variations, except in the region χ ∼ 170 − 177 • . Nevertheless, the central value at NNNL (blue dashed) is much closer to the central value at NNLL (red solid) than at NLL (green dotted). We observe a similar pattern when comparing NNLL and N 3 LL in the right panel, even though the effects are less pronounced at this order. Overall, the central values show very good convergence beyond NNLL, with greatly reduced uncertainties at N 3 LL of about ±4% at the peak, compared to about ±15% at N 3 LL.
In figure 3, we overlay the N 3 LL and N 3 LL resummed spectra of figure 2 with the fixed order NNLO numerical calculation of ref. [7] as well as the experimental measurement of the EEC [4] from the OPAL collaboration at LEP. We observe that resummation effects improve substantially the behavior of the spectrum in the peak region compared to the fixed order calculation. We also note that pushing the resummation to N 3 LL accuracy is crucial to achieve a perturbative control on the theory prediction that is competitive with the experimental uncertainty on this observable. However, it is important to point out that for a realistic comparison with experimental data and for the extraction of the strong coupling constant from this event shape, it is essential to include hadronization effects [2][3][4][5][6][7][8][9]56], as resummation effects alone, even at N 3 LL , cannot fully bridge the gap between partonic fixed order results and data. For completeness, we also show in figure 3 11 We thank Gábor Somogyi for private communication of these results -25 - . Comparison of the N 3 LL and N 3 LL resummed EEC partonic spectrum with the fixed order NNLO numerical calculation of ref. [7], the LEP data from the OPAL collaboration [4], and the analytic leading power spectrum at O(α 3 s ) 12 from eq. (2.29).
the analytic leading power spectrum in the back-to-back limit at O(α 3 s ) 12 (in red) from eqs. (2.27) and (2.29). In the peak region, this result perfectly agrees with the fixed order NNLO numerical calculation of ref. [7], thereby again justifying resummation in this region. Note that the difference between those two results is exactly the non-singular distribution plotted in the right panel of figure 2.

Comparison to literature
Resummed predictions for the EEC in the back-to-back limit have previously been reported in refs. [7,8,35] at NNLL, based on the approach developed in refs. [36][37][38][39]. To compare our formalism to theirs, we start from eq. (4.7) and choose the resummation scales as i.e. we distinguish only an overall high scale µ h and low scale µ l , but always evaluate the rapidity anomalous dimension at its canonical scale µ 0 = b 0 /b T . With these choices, 12 Note that, as explained in the Introduction, in the context of factorization theorems it is customary to count O(α n s ) contributions as N n LO. Therefore, the leading power spectrum for z → 1 at O(α 3 s ) is referred as the EEC in the back-to-back limit at N 3 LO, hence the title of section 2.3. However, since in figure 3 we included the leading power fixed order result mainly for comparison with the numerical calculation of ref. [7], which counts O(α n+1 s ) corrections as N n LO as it is usually done for fixed order calculations of event shapes, we refrained from mixing the notation and therefore labeled it as NNLOLP in the legend.
eq. (4.7) can be rewritten as The coefficients A q and B q in eq. (4.15) are given by where we remind the reader thatγ i is the boundary term of the rapidity anomalous dimension. Notably, the A q coefficient differs from the cusp anomalous dimension starting at O(α 3 s ), which already contributes at NNLL [119]. In eq. (4.15), the first line contains the fixed-order boundary terms, which at canonical scales are free of any logarithms and only depend on α s (µ h ) = α s (Q) and α s (µ l ) = α s (b 0 /b T ). The second line in eq. (4.15) contains the Sudakov form factor that exponentiates the large logarithms. The third line only contributes when µ h = Q or µ l = µ 0 , i.e. when scales are not chosen exactly canonically, and thus can be used to assess resummation uncertainties by separately varying µ h and µ l . We note that this procedure is not quite as refined as the one introduced in section 4.2, where we separately vary all resummation scales. Also note that sinceγ i ν = O(α 2 s ), the second term in this exponential first contributes at O(α 3 s ). To compare eq. (4.15) to the results used refs. [7,8,35], we now explicitly choose the canonical scales, using which eq. (4.15) reads dσ dz =σ For comparison, the resummation formula given in ref. [7] reads 13 where σ tot is the total hadronic cross section. Comparing eqs. (4.15) and (4.18), we first notice that both formulas contain the same Sudakov form factor. 14 However, as already 13 Their expansion of the Sudakov form factor contains an implicit µR dependence, which formally cancels with the µR dependence of the overall hard function. 14 Note that ref. [35] did not use the correct N 3 LO result for A q , which was first obtained in ref. [119].
discussed at the end of section 2.1.1, their result only contains a combined functionH instead of separating physics at the high and low scales into a hard function and jet and soft functions, respectively. It obeys Here, each term is evaluated at canonical scales, and thus only depends on the scale through the running coupling. Eq. (4.19) is to be understood as a reexpansion in α s (b 0 /b T ) = α s (Q) + O(α 2 s ). Due to eq. (4.19), both eqs. (4.17) and (4.18) recover the correct fixedorder expansion of the EEC in the back-to-back limit. However, only eq. (4.17) yields the correct NNLL result, as eq. (4.18) does not contain the correct boundary terms.
As remarked earlier, while the original works in ref. [36][37][38][39] did not yet contain separate hard and jet functions, as they do not yet contribute a the NLL accuracy they work at, the existence of these functions can already be seen in the q T factorization used in those works to derive the EEC factorization in the back-to-back limit. For instance, ref. [39] also explicitly mentions corrections to the TMDFF in α s (1/b T ).
Another crucial difference lies in the estimation of perturbative uncertainties. In our approach, we vary all resummation scales, which in particular separately varies the high scale µ h ∼ Q and the low scale µ l ∼ 1/b T . This reflects that in the back-to-back limit, the EEC contains two parametrically different scales, and varying both resummation scales probes the physics at both scales. In contrast, in refs. [7,8,35] effectively only the hard scale µ R is varied, while the low scale is always kept canonical at µ l = µ 0 = b 0 /b T . This completely neglects the variation from the last line in eq. (4.15), and thus largely underestimates the perturbative uncertainties. Note that due to the lack of a jet and soft function evaluated at the low scale, variations of the low scale can not even cancel formally in eq. (4.18), in contrast to variations of the hard scale µ R . Also note that refs. [7,8,35] avoid the Landau pole by deforming the integration contour into the complex plane, rather than freezing out the scale as done in our analysis.
The above observation already explains why the perturbative uncertainties observed in our analysis are larger than those seen in refs. [7,8,35], as we cover a larger set of scale variations. Moreover, since variations of the low scale µ l ∼ 1/b T probe the strong coupling at much larger values than variations of the high scale µ h ∼ Q, it is not surprising that the former are in fact the dominant uncertainties. To validate this, we compare three methods of estimating uncertainties: while all other scales are kept as in eq. (4.9). This roughly mimics the procedure in refs. [7,8,35].
3. ∆(µ l ): We only consider three variations, while all other scales are kept as in eq. (4.9). These are the only scale variations where the hard scales µ H,J are unchanged, and the low scales are only varied down.
We also keep the rapidity scale ν S = b 0 /b * T canonic, as it does not enter the running coupling.
The results are shown in figure 4, in a similar pattern as in figure 1 such that one can easily compare the two figures. At NLL, we see that the variations of both the high scale (orange) and low scale (blue) are significant, giving rise to a very large overall uncertainty (green). In contrast, at all higher orders we clearly see that it is indeed the low-scale variation ∆(µ l ) (blue) that dominates the total uncertainty, while the high-scale variation ∆(µ h ) and the remaining scale variations are almost negligible.
In ref. [35], the uncertainty of the peak at NNLL was given by roughly ±8%, compared to ∆(µ h ) ∼ 12%. By only considering the variation ∆(µ h ), our uncertainty at N 3 LL (N 3 LL) reduces to only ±0.5% (±2%), compared to our more conservative uncertainty of about ±4% (±15%) when using all scale variations. In both approaches, the uncertainty reduces by a factor of about 4 when going from N 3 LL to N 3 LL , illustrating the importance of including the N 3 LO boundary terms computed in this work. However, we stress again that only varying the high scales neglects important uncertainties from soft physics, and thus the ∆(µ h ) variation alone is not sufficient to obtain a robust estimate of theory uncertainties.
We close by remarking that it was already remarked in ref. [39] that the EEC is quite sensitive to the high b T region, and nonperturbative model functions were introduced in order to achieve agreement with CELLO data. This is consistent with our observation that variations of the 1/b T scales yield the dominant uncertainties.

Conclusions
In this work we have calculated the full singular structure of the Energy-Energy Correlation (EEC) in the back-to-back limit at O(α 3 s ) in QCD, including contact terms. Our work applies both in the case of e + e − annihilation as well as in gluon induced Higgs decays. To obtain these results we have calculated the quark and gluon jet functions for the EEC in the back-to-back limit at N 3 LO, which were the last missing ingredients for the factorization theorem in this limit at N 3 LO.
The computation of the jet functions relies on the calculation of the kernels of the transverse-momentum dependent fragmentation functions at N 3 LO in our companion paper [42] which have been obtained using a recently developed method for the expansion of cross sections around the collinear limit [62]. We checked that the logarithmically enhanced terms of our calculation match those predicted by the rapidity renormalization group (RRG) evolution.
By comparing the leading transcendental part of our results we show that the principle of maximal transcendentality, which states that a quantity obtained in N = 4 SYM constitutes the leading transcendental term of the same quantity in QCD, holds for the EEC in the back-to-back asymptotic up to O(α 3 s ). In particular, we show that the leading transcendental part of the EEC in e + e − is identical to the one for the EEC in Higgs decay to gluons and that both match the result in N = 4 [19,23,41], not only for the logarithmic part but also for the contact terms. This also provides a non-trivial cross check on both the quark and gluon jet function calculations of the δ(1 − z) constants in addition to the one on the logarithmic parts coming from the RRG evolution.
Leveraging on the fact that the EEC obeys a set of non trivial sum rules [23,41], we used our newly calculated result for the contact terms at O(α 3 s ) as well as the logarithmic enhanced contributions in the small angle limit [24] to obtain the N = 2 Mellin moment of the EEC distribution in the bulk at NNLO in QCD analytically. This constitutes the first piece of analytic information on the EEC distribution at this order in QCD away from the endpoints, both in the case of e + e − annihilation as well as in gluon induced Higgs decays.
Finally, we have carried out the resummation of the EEC in the back-to-back region at N 3 LL accuracy. This is the first time an event shape observable is resummed at this level of accuracy and, more generally, this constitutes the highest level of resummation for any infrared and collinear safe observable in QCD to date. We show that by performing the resummation at N 3 LL we obtain a reduction of uncertainties by a factor of ∼ 4 in the peak region compared to previous results obtained at lower accuracy. We thoroughly discuss different schemes to estimate the uncertainties due to missing higher order corrections both in the boundary terms as well as in the anomalous dimensions. We compare these different schemes with the ones used in the literature for this observable. Adopting a scheme in line with the ones previously used in the literature we obtain a 5 per-mille uncertainty at the peak. Using a more conservative scheme, which includes a significant contribution from non-perturbative regions, we obtain a 4% uncertainty for our result at N 3 LL . We point out that the accuracy for lower order results is severely affected by the choice of scheme and that varying the low-energy scale gives a dramatically larger estimate of uncertainties.
The recent progress in understanding Energy-Energy correlators is very promising and we believe it shows that they will play a crucial role in improving our understanding of the strong interaction in the years to come. For example, as the EEC has been often used to determine the strong coupling constant [2,3,[5][6][7][8][9], it would be interesting to leverage the high level of perturbative control we gained on this observable thanks to the N 3 LL resummation, to improve the extraction of α s , complementing the extractions based on other event shape observables in e + e − such as thrust and C-parameter [120][121][122][123]. In addition, it could also be used to extract nonperturbative corrections to the rapidity anomalous dimensions. We expect that further studies of perturbative [22] and non-perturbative [56,57] power corrections will be important to improve the theoretical control on this observable.
It will also be interesting to explore the application of the techniques used in this work and its companion paper [42] to higher point energy correlators in QCD where recent progress has been obtained [28,124,125]. eq. (2.16) are given by Here, the Γ i n and γ i J n are the O[(α s /4π) n ] coefficients of the cusp and jet noncusp anomalous dimensions, respectively, where we remind the reader that the jet noncusp anomalous dimension is identical to that of the TMD beam function, γ i J n ≡ γ i B n . Explicit expressions for these anomalous dimensions in our conventions are collected in ref. [65]. The corresponding fixed-order expansion of the polarized gluon jet function J g can be obtained from eq. (A.1) by dropping all terms that do not contain an explicit factor j where the plus prescription on the right hand side now acts with respect toz. For example, the first few Bessel transforms read I 0 = 1 2 δ(z) , where L h = ln(Q 2 /µ 2 ). Starting from n = 3, the R (n) 2 terms start to induce ζ values.

C Sum Rules Ingredients
Here we collect all required results for the quantities in eq. (3.8).
First, we note that V The expressions for V (n) 0 , the coefficients of δ(z), can be in principle extracted via the factorization formula for the EEC in the collinear limit of ref. [24]. Here, we have extracted them by using the sum rule in eq. (3.9) and analytically integrating the regular part of the distribution obtained from the fixed order calculation of refs. [14,15]. For e + e − we obtain Finally, for the expressions for the total cross section we take the results of ref. [77]. For e + e − they read 15 15 Note that for the color structure d abc d abc of the singlet part we adopted the same convention as in eq. (2.29), which is different from the one adopted in ref. [77].
For Higgs, they are given by