NNLL Resummation for Projected Three-Point Energy Correlator

The projected energy correlator measures the energy deposited in multiple detectors as a function of the largest angular distance $x_L = (1 - \cos\chi_L)/2$ between detectors. The collinear limit $x_L\to 0$ of the projected energy correlator is particularly interesting for understanding the jet-substructures, while the large logarithms of $x_L$ could potentially spoil the perturbation theory and must be resummed. As a necessary ingredient for its resummation at next-to-next-to-leading logarithmic (NNLL) accuracy, we calculate the two-loop jet functions for the projected three-point energy correlator (E3C), using direct integration method and the parameter space Integration-by-Part (IBP) method. We then present the NNLL resummation for $e^+e^-$ annihilation and an approximate NNLL resummation for $pp\rightarrow jj$ process, where the two-loop hard constant is estimated in the latter case. The convergence is improved and the hadronization effect in the collinear limit is suppressed when considering the ratio of E3C distribution to two-point energy-energy correlator (EEC). Our results show potential in precision determination of strong coupling constant using energy correlators from both $e^+e^-$ data and $pp$ data.


Introduction
Energy correlators are a class of multi-particle angle correlation functions, weighted by the particle energy.Thanks to the energy weighting, they are infrared and collinear safe observables and can be calculated in perturbation theory.The simplest energy correlator is a two-point energy correlator, or Energy-Energy Correlation function (EEC).Proposed in 1970s [1,2], EEC measures the correlation of energy deposited in two detectors as a function of the angle χ between them.In perturbation theory, the definition of EEC reads dσ [2]  d cos χ ≡ i,j where i, j run over all the final state particles, ⃗ n i and ⃗ n j are unit three-vectors that define the directions of the particles, and Q is the total energy in the center-of-mass frame.Compared with other event shape variables studied at Large Electron-Positron Collider (LEP), one advantage of EEC is its simple analytic properties.As far as we are aware of, EEC is the only event shape that can be calculated analytically beyond leading order, e.g.it's now known analytically through to next-to-next-to-leading order (NNLO) [3,4] in N = 4 super Yang-Mills (SYM) theory and through to NLO in QCD [5][6][7].
In recent years, increasing attention has been paid to generalization of EEC to N -point energy correlators, which measure the energies of the outgoing particles with N detectors at colliders and turn out to be a function of N (N − 1)/2 angles among these detectors [8][9][10][11][12][13][14].For example, the three-point energy correlator (EEEC) is defined as which gives rise to rich functional dependence on the angles and can be used to probe various properties of perturbative QCD.The LO EEEC was first computed in the triple collinear limit in Ref. [9], later genelarized to arbitrary angle dependence in both N = 4 SYM [15] and QCD [14].To reduce the dimension of the kinematic space of the measured angles without losing too much useful information, one can project the kinematic dependence into a 1D subspace, which leads to the so-called projected energy correlator [16].In momentum space, projected N -point energy correlator (ENC) is given by restricting the maximum angular distance to be x L : and for example, EEEC is then reduced to the projected three-point correlator (E3C).In this work we are mainly interested in the small angle, or collinear limit of E3C, namely x L → 0.
It is well-known in the boundary of phase space, incomplete cancellation of infrared divergences can lead to large logarithms that could possibly spoil the convergence of the perturbation theory and thus it is essential to resum these large logarithms to all orders.EEC is special as it exhibits both large logarithms in collinear limit and back-to-back limit.In this work we are interested in the large logarithms in the collinear limit, for which the most singular terms behave as α n s ln n x L at n loops.In the collinear region, EEC can be factorized into a hard function and a jet function, both of which live in the flavor space.The resummation of collinear EEC has been performed up to NNLL accuracy in both QCD [17] and N = 4 SYM [17][18][19].More interestingly, the collinear factorization can be easily generalized to three-point energy correlator [9] and even the projected N -point energy correlator [16].Previously, LL and NLL resummation has been performed in [16,20,21].To improve upon those results, it is necessary to compute the relevant jet and hard function to higher order.While the hard function is universal for them, the jet functions differ by the measurement function.One of the key new results in this paper is the calculation of two-loop jet function for projected three-point energy correlator, which is the last missing ingredient for NNLL resummation of projected three-point energy correlator in e + e − collider.
One of the main motivations for improving the theoretical accuracy of projected energy correlators comes from the possibility of determining the strong coupling constant α s by measuring the ratio of projected energy correlators [16].Measurements of strong coupling constant using classical QCD event shape observable has been actively studied for a long time, e.g.[22][23][24][25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40].In recent years, there has been increasing attention to using jet substructure observables to extract α s , such as soft-drop thrust and jet mass [41,42], see also [43] for α s determination from jet substructure by demixing quark and gluon jets.Since we are mainly concerned with the collinear limit of projected energy correlators in this paper, our results naturally provide theory input for measuring projected energy correlator within a jet, treating it as a jet substructure observable.We will show that considering the ratio of E3C and EEC can significantly reduce scale uncertainties and hadronization corrections, which makes it a good candidate for precision determination of α s using jet substructure.We also note that energy correlators have the advantage that they can be defined and calculated using charged hadrons only [16,44].Using the track function formalism [45,46], it is possible to perform precision calculation for projected energy correlators on tracks in the future.
The outline of this paper is as follows.In Sec. 2, we present the factorization theorem for ENC in the collinear limit and the RG evolution for both hard function and jet function.The desired orders required for all the ingredients to achieve NNLL resummation are briefly summarized there.In Sec.2.4, we calculate the two-loop E3C jet function.Modern multiloop techniques like IBP and differential equation (DE) are applied for both finite and contact terms.Combining all together, we are able to extract the two-loop E3C jet constants, which is the last missing piece of the NNLL resummation for collinear E3C in e + e − collision.In Sec. 3, we present the matched NNLL results for both E3C and the ratio of E3C to EEC in e + e − collision.A qualitative analysis is performed to estimate the leading hadronization correction.The resummation procedure is extended to the case of pp collision, in particular, the pp → dijet process in Sec. 4. We present the highest pertur-bative prediction given the available ingredients, the approximate NNLL, with the missing two-loop hard function constants estimated and included as an additional uncertainty.We summarize and conclude in Sec. 5.

Factorization theorem
In this subsection, we summarize the factorization theorem for the projected N -correlator in the collinear limit and describe the necessary ingredients for NNLL resummation [16].Similar to EEC, N -point energy correlator (ENC) in this limit is dominated by the logarithmic series of the largest angular distance x L dσ [N ]  dx L = where L −1 (x L ) = δ(x L ) and L j (x L ) = ln j (x L )/x L + for j ≥ 0, with standard plus distribution.We do the logarithm counting in the projected N -point energy correlator cumulant, defined as which maps At leading power, the e + e − cumulant Σ [N ] can be written in terms of a modified factorization formula in the collinear limit x L → 0 [16]: where the hard function ⃗ H [N ] ee encodes the production of a parent parton with energy fraction x with respect to the center of mass energy, and the jet function ⃗ J [N ] encodes the evolution of the parent parton into a number of collinear partons which contribute to the observable.Similar factorization formula for EEC was first obtained in [17], and checked explicitly with known NLO results in QCD [5,6] and N = 4 SYM [3,4].We note the explicit dependence on the variable x in both the jet function and the hard function.Ignoring the dependence on different quark flavor, both jet and hard functions are two-component vectors living in the flavor space, i.e. ⃗ J [N ] = {J g }, ⃗ H ee = {H ee,q , H ee,g }.We will describe their definition for both e + e − annihilation and pp collision in detail in the following subsections.We also emphasize that the factorization theorem holds for any N at leading power, though we only calculate the N = 3 case in this paper.Finally the energy weights in the distribution makes projected N -point energy correlator insensitive to the soft radiations and non-global logarithms.
In hadron colliders, the largest angular distance x L is replaced by the rapidity-azimuth distance , where X E is the set of particles that contributes to the energy weight.When the projected energy correlators are measured within a jet, as is typical for jet substructure observable, the cumulant Σ [N ] had also depends on the jet radius R 0 parameter.In the limit of R L ≪ R 0 , the modified factorization formula can be written as where p T is the jet transverse momentum.Around R L ∼ R 0 , the jet function can also depend on R 0 .However, there is no large logarithms associated with R 0 , and its dependence can be obtained from fixed-order matching.For simplicity, we will ignore the R 0 dependence in the jet function.In that case the jet function become universal between e + e − and pp collision.For pp collision, the hard function depends on the partonic scattering process, as well as parton distribution functions (PDFs).

e + e − annihilation
For e + e − , the hard function is simply the semi-inclusive hadron fragmentation function [47], which depends on the parton flavor and parton energy fraction x = 2p•q Q 2 , where q is the total momentum and p is the parton momentum.The leading order hard function follows from the born process e + e − → q q, ⃗ H (2.5) The factor 1/2 in front of the quark channel indicates for identical contribution from antiquark, since we do not dinstinguish quark and anti-quark flavor.At two-loop, the hard function can be found from the coefficient functions in [47].Similar to the hadron fragmentation function, the renormalization group evolution (RGE) for the hard function ⃗ H is simply the DGLAP equation, with P (y) being the singlet timelike splitting matrix, which is now known to three loops [48,49].While it is very difficult to derive an analytic solution for DGLAP to all orders in α s , as we will see below, our resummation only uses a α s -expanded solution (which turns out to be a very good approximation) and only requires certain moments of the hard function.
Explicitly, we will only need the regular and logarithmic moments for the hard function defined as the following [17], Here we use x N ln x = ∂ N x N and the dot on the RHS stands for the derivative.The expressions of needed hard function moments can be found in Appendix A.

pp collision
In hadronic collisions, we mainly focus on the dijet production pp → jj, which has a relatively large cross section at the LHC.Different from e + e − collider, this hard function incorporates the partonic scattering cross sections, the contribution from parton distribution functions (PDFs) and the jet algorithms for clustering the particles.Currently, to the best of our knowledge, the hard function is not know at two-loop.However, important progress are being made to compute those hard functions, e.g.[50].Similar to the e + e − case, our resummation will only need the hard function moments.
In this work we evaluate the needed moments of the hard function numerically in Madgraph5 [51,52].To investigate the sensitivity of the result to the values of α s , we used three different PDF sets: NNPDF31_nnlo_as_0112, NNPDF31_nnlo_as_0118 and NNPDF31_nnlo_as_0124 through Lhapdf [53].Each PDF set fixes also the value of α s (m Z ) and the corresponding evolution in Madgraph5.To address the fact that the hard function contains collinear divergence when resolving the energy fraction of the quarks and gluons, we use the one cut-off phase space slicing to regularize the collinear singularity, as implemented in [54].With the collinear divergent contribution singled out and calculated analytically, the remaining contributions can be evaluated numerically.The detailed discussion can be found in Appendix A.
For pp → jj, we adopt the anti-k t algorithm [55] for jet detection and use the following parameters in the calculation The two leading jets are further subject to the following cuts and cast to the corresponding p t bins for the analysis.The calculated moments need to be normalized with the cross section σ J of jet production within specific p t range.In particular, we expand H had /σ J to NLO in a s , and take the O(a 0 s ) and O(a 1 s ) as the leading and nextto-leading order results.For the purpose of phenomenological studies, we will focus on two different p t ranges: [300, 350] GeV and [500, 550] GeV.The hard function moments needed for NNLL are also summarized in Appendix A.

Jet functions
The E3C jet function, on the other hand, encodes the measurement information.From RG invariance of the modified factorization formula (2.3), the jet function satisfies a modified timelike DGLAP evolution equation (2.10) In order to write down an operator description of the E3C jet function, we first recall the collinear EEEC jet function from [9]: where W n is the collinear gluon, and P µ n⊥ form a complete set of collinear gauge invariant building blocks [56] in SCET [57][58][59][60][61].The triple collinear measurement function M EEEC is defined as with θ ij being the angle between parton i and j.Then our E3C jet function has the same form as EEEC jet function, with a replacement of the measurement function: There are two folds integration in the first line.The first one is performed in the allowed kinematic space {x 1 , x 2 , x 3 } ∈ K that will be discussed below, projecting the shapedependent EEEC jet function into a single-scale jet function.The second integration brings the differential measurement to the cumulant level.For N > 3, the measurement function takes a similar structure, with more δ functions and integrations.Perturbatively, the E3C jet function can be written as J q,g = L (α s /4π) L J (L) q,g , and we use the normalization condition 2 3 • J (0) q = 2 3 • J (0) g = 1 as in Ref. [16].The one-loop correction can be calculated from the QCD 1 → 2 timelike splitting kernel and is given by Note that the µ-dependent terms are precisely captured by the jet RGE, while the remaining constants have to come from the fixed-order calculation.One of the main result in this paper is to calculate the two-loop constants described below.

Two-loop calculation for the E3C jet function
In this subsection, we present the two-loop calculation of the E3C jet functions for both quark jets and gluon jets.Since they are universal in the small angle limit, they can be used in both e + e − collision and pp collision.
We start from recalling the definition of E3C at finite angle before taking the small angle limit.At two loops, E3C receives contributions from double-real (RR) and real-virtual (RV) as well as double-virtual (VV) corrections to q → q, from which the quark jet function can be extracted by matching to the factorization formula, (2.3).Similarly, the gluon jet function can be extracted from the NLO E3C distribution of Higgs gluonic decay H → gg.To organize the calculation, we rewrite the definition of E3C in Eq. (1.3) with the number of energy weight: where we normalize the distribution to the born cross-section in d dimension.The first line represents the contribution from nonidentical energy weights measurement and the other lines are called contact terms.If we define then in the collinear limits, they are the contact terms for δ(z z) that captures the strict squeeze limit and δ(x L ) that captures the strict triple collinear limit.The main goal of this section is to compute the collinear limit of Eq. (2.15) and extract the corresponding two-loop constants.
The lowest regular distribution of the E3C quark jet function comes from tree-level process γ * → 4 partons in electron-positron annihilation, which under the triple collinear limit, factorizes into the born process γ * → q q and the 1 → 3 splitting functions, and we will call it nonidentical energy weight term.Below we will introduce two different methods to compute this part.The traditional method is to calculate the EEEC jet function to order O(ϵ) and to integrate two angular distances x 2 , x 3 numerically by the interpolation method.The OPE singularities (sometimes called squeezed singularities) of EEEC are subtracted and integrated in d dimension separately.The second approach benefits from the parameter space IBP method [62][63][64] developed very recently.Only 7 master integrals are needed to express EEEC, allowing the precise calculation of the remaining two-fold integral.
The other two parts contribute to the contact terms and cancel the infrared divergence, which is guaranteed by the Kinoshita-Lee-Nauenberg (KLN) theorem [65,66].Similar to EEC at NLO, the measurement function in the contact terms can be treated as a nonstandard cut propagators, which allows for a generalized IBP reduction in Litered [67,68] and Fire6 [69].The master integrals then can be calculated in packages like Canonica [70] or Libra [71,72] with the differential equation method implemented.

Nonidentical energy-weight terms
We start by computing the nonidentical energy-weight contribution in the traditional approach.As discussed in Ref. [9], the inclusive jet function J ijk is related to the 1 → 3 splitting function P ijk [73][74][75] through where dΦ (3) c is the triple collinear phase space [75,76], and i, j, k run over all final-state particles.The fully differential distribution with respect to all angular distances {x 1 , x 2 , x 3 } in d = 4 − 2ϵ dimension is then written as  where G(z), F (z), H(z), • • • the shape function in ϵ expansion.The order O(1) part G(z) is computed analytically in [9] and following the same approach, we also calculate the complete result for F (z) and the z → 1 limit of H(z).We will see that these are all the needed ingredients for nonidentical part.Note that the x L dependence is defined by plus distribution, where In order to perform the integral over z, we need to figure out the integration region first.Compared with the first line in Eq. (2.15), it is straightforward to show that ) where the constant factor 6 comes from the S 3 permutation symmetry and the integration region S is given in Fig. 1.To calculate A(ϵ) numerically, we also need to subtract the OPE singularities around z → 1 at the integrand level, and evaluate its z integration analytically in d dimension.The full asymptotic expansion of z → 1 is given in the appendix C. The most singular term is proportional to The integration region S for E3C jet function.The S 3 symmetry is applied to reduce the entire region of z into 6 times the blue region.The integration range for z then becomes Here κ = ImLi 2 e i π is the Gieseking's constant living in the transcendentality-two family and η = ImLi 3 i √ is a parity-odd transcendentality-three constant.These constants are typical numbers in loop integrals, especially in trijet observable calculations.
With subtraction terms, the integral A in Eq. (2.19) up to order O(ϵ) is then written as The first term is proportional to Eq. (2.20) and it is straightforward to compute it to O(ϵ).For the second integral, we have to expand in ϵ and evaluate it numerically.To implement the interpolation method, we first change the integration variables via , such that both v 1,2 range from 0 to 1. Then we can build a 2D lattice by discretizing v 1,2 and approximate our integrand with polynomials.This allows one to perform the two-fold numerical integral directly in Mathematica.To check the stability of the integration and estimate the statistical error, we vary the lattice size and the order of polynomials and see which significant figure remains unchanged.Eventually we obtain both δ(x L ) contact term and1 x L finite term for the nonidentical energy weight contribution.The explicit expression for both quark and gluon jet function can be found in Eq. (D.1)-(D.4) in the appendix.
Alternatively, benefiting from the recent development of the IBP method in the Feynman parameter space, we can simplify the whole jet function calculation with integral reduction.First of all, recall that Eq. (2.16) takes the form (2.22)Here P is a homogeneous function of the energy fraction ω i of the final-state particles.Explicitly, it is of the form with f 1 linear in ω i , and f 2 a polynomial of ω i of degree 2. Following the idea in Ref [9], the integral d 3 J R dx 1 dx 2 dx 3 in Eq. (2.22) can be related to a Feynman parameter integral through 1 where The integral in the last line is a standard parametric Feynman integral, which can be reduced with IBP reduction [77,78] in the parametric representation [62][63][64]79] 2 .The master integrals are with the integrals I i defined by the F polynomials and α 0 = 6ϵ − 23 .The master integrals can be evaluated using the differential equation technique [84,85].For simplicity, we set µ = x 3 = 1, and introduce u and v following z = u(1 + iv).Then we construct the differential-equation system with respect to u, and derive the canonical basis [86] using Libra [71, 72] with the corresponding alphabet {u, 2u−1, x 2 , x 2 −1}.By solving the differential-equation system, we can express the master integrals via Goncharov polylogarithms (GPLs) [87][88][89].
The GPL is defined iteratively by After finishing the simplified calculation of EEEC in the collinear limit, we still need to integrate two angular distances for the projected EEEC as the previous approach.By virtue of the S 3 permutation symmetry, this amount to consider where ity corresponds to u → 0 limit, and similarly, we need to subtract the singular behavior and do the integration separately: where again we can evaluate the first integral in d dimension and expand the integrand of the second one in ϵ.
To calculate the J(u → 0), now we can directly extract the asymptotic expansion of the integral I in Eq. (2.24) from DE, in which we identify two expansion regions: (2.32) Evantually we only need to integrate the reduced master integrals in d dimension.
Regarding the second integral in Eq. (2.31), the u integral is straightforward since J(u, v) is expressed in terms of GPLs of the form G(. . ., u).However, the v integral becomes unstable in two regions v → 0 and v → ∞.To resolve this problem, we decompose the v ∈ [0, ∞] integration into three parts: and [C, ∞], with a arbitrary cut parameter C > 1.In the region ( 1 C , C), we carry out the integration numerically, with the GPLs numerically using Handyg [90].The other two regions require expanding the integrand in v (or ) and performing the integration analytically.This expansion can easily be done by asymptotically solving the differential equations satisfied by the GPLs.Eventually, we find the same result as in Eq. (D.1)-(D.4).

Contact terms
While it is convenient to calculate the nonidentical E i 1 E i 2 E i 3 part starting with the splitting functions, it is preferable to compute the full angular dependence on x L for corresponding processes (namely e + e − annihilation and gluonic Higgs decay) with energy weights i 1 , and extract the contact term from the collinear limit x L → 0. In other words, we will adopt the full matrix elements squared and compute the full phase space integral using modern multi-loop techniques, with which the collinear expansion gives dσ We start with the relevant processes in perturbation theory for two-loop jet functions, e + e − annihilation Higgs decays where V and V V denotes one-loop and two-loop correction respectively.In particular, in the x L → 0 limit, 1 → 2 processes only contribute to δ(x L )-terms (i.e., dσ shares the same structure as the original EEC, which basically follows the approach described in Ref. [5] and more detail in [6].Briefly speaking, using the Cutkosky rules [91,92], we can replace the phase-space on-shell delta functions with the cut propagators and also the EEC measurement function δ(x L − x i,j ) with where we set the center-of-mass energy Q = 1 for simplicity.After topology classification and identification as described in Ref. [6], the E 2 EC integral can be reduced to a set of master integrals I k (x L , ϵ) using IBP reduction and E 2 EC distribution can be written as a linear combination of the master integrals, Specifically, we generate the standard IBP equations using Litered [67,68], add the missing one that is associated with the EEC measurement function by hand, and do the reduction in Fire6 [69].The master integrals turn out to be the same as in NLO EEC calculation for both e + e − annihilation and gluonic Higgs decays, which can be converted into the canonical basis using the DE package Canonica [70].
In order to obtain the collinear dσ E 2 EC (x L , ϵ)/dx L , one could surely expand the differential equation asymptotically and derive the analytical expression of the master integrals in that limit.However, the fact that the most singular power of C k 's is x −8  L requires us to compute the master integrals up to O(x 7 L ) order, which turns out to be expensive and timeconsuming.This becomes worse in the higher-point energy correlator since the singular power increases as well.One antidote is to reconstruct the coefficients from DE following an ansatz on the structure of asymptotic expansion.In fact, the pattern turns out to be , where U denotes a series in x L with rational fractions of ϵ as the coefficients.
Therefore, we perform the asymptotic expansion in the following way.First of all, we solve the canonical DE at 0 < x L < 1 to transcendental-weight 5, which can be used to obtain the finite part of the contact term via Eq.(2.36).The result can be converted to Harmonic polylogarithms (HPLs) with the package Hpl [93] or even classical polylogarithms.Then we can extract the leading power x −1 L and match it to a resummed ansatz with unknown ϵ-series C 1 (ϵ) and C 2 (ϵ).The matching between fixed order calculation and the resummed structure in ϵ leads to the solution of C 1 (ϵ) and C 2 (ϵ) in ϵ expansion.Since L are defined with plus distribution similar to Eq. (2.18), now we obtain the correct O(ϵ 0 ) formula for dσ [3] E 2 EC (x L , ϵ)/dx L in the collinear limit.
The last remaining piece is dσ The computation of the self-energy correlator is much easier since its dependence on x L is factorized out by δ(x L ) and the integrals are simply standard cut integrals.The master integrals can be found in the literature, e.g.[94,95].Eventually adding dσ we obtain the complete contact terms dσ [3] C (x L , ϵ)/dx L for E3C distribution.The results are also summarized in Eq. (D.5)-(D.8).Combined with the nonidentical energy weight contributions, we find all 1  ϵ canceled and thus the infrared safety is guaranteed as expected.

Results of two-loop jet function constants
With all individual contributions at hand, the full expressions of 2-loop E3Cs in the collinear limit can be written as ) dx L (gluonic Higgs decay) . (2.39) Here a factor of 2 is added because we only consider a single jet in Sec.2.4.1.Given the treelevel hard functions, {H q , H g } = {2δ(1−x), 0} for e + e − annihilation and { H(0) q , H(0) g } = {0, 2δ(1 − x)} for the Higgs decay through the effective Hgg coupling, we can extract the two-loop jet constant directly from the δ(x L ) contribution from Eq. (2.38) and Eq.(2.39).We find that the µ dependence are in full agreement with prediction from RG evolution, providing strong check to our calculation.The µ independent part are the new results from this calculation.For the quark jet function, we get and for gluon jet functions (2.41)

Perturbative resummation
We start by defining the logarithmic order for our E3C resummation.The ingredients needed for our E3C resummation are summarized in Table 1.This includes the order of timelike splitting kernel P (y), the boundary information (hard and jet constants), the β function for running coupling as well as the fixed-order matching. 4Due to the absent of analytic method to solve the RG equation exactly, we also truncate in the number of loops of the RGE solution to the desired logarithmic order [17].
We first review the LL resummation in e + e − annihilation.Based on our resummation setting, it is safe to set x = 1 in the argument of E3C jet function in Eq. (2.10), which only affects the higher-order terms beyond LL.This leads to  1: Definition of the resummation order and their corresponding fixed-order matching.
Here, we introduce the anomalous dimension to be the moment of timelike splitting kernel Then given the boundary condition ⃗ J (0) = {2 −N , 2 −N }, we can directly write down the solution to LL jet function: Plugging both jet and hard functions into the factorization for the cumulant Σ [N ] and differentiating it with respect to x L , we obtain the LL resummed physical spectrum for E3C.
Beyond LL, the x = 1 approximation is no longer valid, and instead we have to solve the jet RGE directly.While it is difficult to obtain a close-form solution for this modified DGLAP equation, we find that a truncated solution in α s is already in good convergence.Explicitly, we assume the jet function takes the form and c i,j unknown constants, and solve both the jet RGE and β RGE order by order in α s (which is referred as expanded solution).In practice, we evaluate it numerically up to O(α 50 s ).Another advantage of using expanded solution is that we only need certain moments of the hard functions.For example, consider one term from the jet function, ⃗ J , and plug into Eq.( 2.3), we find where the three terms correspond to the standard moment, the single logarithmic moment and the double logarithmic moment of the E3C hard function.To derive the last line, we also use the following relation (2.47) In the Appendix A, we provide all the hard moments with N = 2, 3 that are required for NNLL resummation.
In this paper, we present results for the NNLL resummation of E3C for e + e − annihilation, and approximate NNLL resummation for jets from the hadronic collision process pp → jj.For e + e − annihilation, we have all ingredients needed for NNLL resummation.And since there is no accurate fixed-order data for E3C at NNLO, we will instead match the NNLL result to NLO.Regarding the dijet production, due to the absence of the twoloop hard constant, we will present the approximate NNLL resummation (which we refer as NNLL approx ), with an additional uncertainty coming from the missing two-loop hard constant.Resummation with the accurate two-loop hard function as well as the matching with fixed-order result are left as future improvements.

NNLL resummation in e + e − annihilation
With all the ingredients at hand, now we can present the NNLL resummation prediction.In this section, we first consider e + e − collision at two different energies: 250 GeV and 1 TeV.In the resummation calculation, we will use α(m Z ) = 0.118.

Resummation results
Following the discussion in Sec.2.5, our resummation is performed by perturbatively solving the jet function RG equation to order O(α 50 s ), plugging back to the cumulant factorization and finally truncating the logarithms to the desired order.In the resummation formula, we set canonical jet scale µ j = µ h √ x L in the factorization, leaving a single hard scale µ h = µ in the resummed expression.We vary the scale µ to estimate the uncertainty from higher order corrections.Regarding the observables, below we consider three cases: N = 2, N = 3 and their ratio.
The N = 2 case is precisely the EEC observable, where we directly use the result from Ref. [17], and the singular expansion has been verified against the NLO EEC fixed-order calculation.For N = 3 case, this is the main result of this paper.In Fig. 2, we first check our O(α 2 s ) expansion with the Monte Carlo program Event2.In the collinear limit, we find excellent agreement between theory and numeric result, while in the meantime, this also suggests the non-singular contribution from fixed-order calculation is negligible in this limit.
Nevertheless, the matching formula can be written as Here each term is a function of α s (µ) evaluated at the hard scale µ h = µ.In Fig. 3, we present the E3C resummation up to NNLL, matched to fixed-order.As explained above, due to the absence of NNLO data, we only match NNLL to NLO.The hard scale is chosen to be half of the center-of-mass energy µ = Q jet ≡ Q/2, the typical energy for each quark jet, and the scale uncertainty is obtained by varying the hard scale by a factor of 2. In both energies, the uncertainty band width goes down as we increase the resummation order, while at 1 TeV, we have a tighter band because the coupling α s runs slower at high energy.
At NNLL, we find a relative 4% hard uncertainty for Q = 250 GeV and 2% for Q = 1 TeV.We find large corrections as we go from LL to NNLL, as was also observed previously in [17], which emphasize the importance of higher order corrections.For higher center-of-mass energy, the convergence between different orders is improved.
To improve the convergence, we also introduce the ratio of different point energy correlators, namely [16] where µ and µ ′ are the hard scale in dσ [m]  dx L and dσ [n] dx L respectively.In particular, we focus on the ratio between fully matched E3C and EEC, i.e. ∆ 3,2 (x L ).In Fig. 4, we show the NNLL resummed ∆ 3,2 (x L ) at again Q = 250 GeV and Q = 1 TeV, and find good convergence.This implies that the ratio can be used as precision observable.For hard scale uncertainty, we use the seven-point scale variation, which amounts to varying the scales in both numerator and denominator independently by a factor of 2, to a combination of and take the envelope as the uncertainty estimation.The convergence also indicates that ENC shares similar non-perturbative behavior in the collinear limit and taking the ratio strongly suppresses the power corrections.

Hadronization corrections
In this subsection, we consider the power-suppressed hadronization corrections in the collinear limit.At present hadronization corrections cannot be computed from first principle.For simplicity, we use a phenomenological form for the leading non-perturbative power correction as suggested in [96], and fit the unknown parameters from a Monte Carlo program.This provides some insights on how to model the hadronization effect for a global fit in the future.
In general, the non-perturbative corrections in infrared-collinear safe observables are (at least) suppressed as Λ QCD /Q to some power, where Q is the hard scale of the process.Following from the LL result in Eq. (2.44), we observe that in the collinear limit, there exists a lower scale √ x L Q in the coupling, and the most important non-perturbative correction that could potentially appear is linear in Λ QCD and takes the form with an extra kinematic factor 1/x L .The sub-leading non-perturbative corrections with additional powers of Λ QCD /( √ x L Q) will become necessary down to small x L ∼ Λ 2 QCD /Q 2 , where the perturbation theory also breaks down.For the leading non-perturbative correction we are considering, such structure is in fact recovered for the EEC in the fragmentation modeling of non-perturbative radiations [1] and and analysis using renormalon or dispersive techniques [96][97][98].
As a qualitative analysis, we use the following parametrization of the leading nonperturbative correction, we verify the scaling behaviour of the non-perturbative correction in the collinear limit for both EEC and E3C distributions with Pythia8 [99], and extract the non-perturbative parameters by fitting from the difference of the hadron level and parton level predictions.Note that the issues of extracting the non-perturbative power corrections from Monte Carlo generators have been pointed out in Ref. [36].In particular, the corrections from the hadronization modeling in the Monte Carlo programs in fact unfaithfully absorb partial subleading-log contributions, as the hadronization modeling has been tuned to reproduce some collider data with limited perturbative accuracy.Therefore, in this paper we only use Monte Carlo to illustrate the impact of power correction for individual EEC and E3C distribution as well as their ratio.
For our case, we stay in the default settings of Pythia8 and obtain the following fit at the 95% confidence level.At Q = 250 GeV, we find for EEC and E3C: Λ2 = (0.956 ± 0.031) GeV , γ 2 = 0.462 ± 0.017 , Λ3 = (0.500 ± 0.040) GeV , γ 3 = 0.335 ± 0.031 .We emphasis that for too small x L value, the leading order non-perturbative approximation itself becomes invalidated.The enhancement of the non-perturbative corrections in the collinear limit must be turned off before entering the fully non-perturbative phase, where the degrees of freedom become freely interacting hadrons and a nice scaling behavior follows [20].In this qualitative analysis, we choose the lower bound of the fit range by finding the extreme point of the distributions from hadron level prediction in Pythia8.
Multiplying the extreme point by a factor of 2 gives a good estimate of the lower bound for the range where the non-perturbative correction follows the described scaling behavior.In Fig. 5, we show the relative hadronization correction from both Pythia8 and our two-parameter fit.Except the shaded region, our parameterization agrees with the Monte Carlo result and it is sufficient for understanding their structure.In Fig. 6, we include the non-perturbative correction in the matched E3C resummation, which strongly enhances the extreme collinear limit.At Q = 1 TeV, the non-perturbative correction changes our NNLL+NLO prediction by only a few percent at x L ∼ 0.1, while this modification reaches 50% at x L ∼ 10 −4 .This shows that the non-perturbative corrections for energy correlators, though being power suppressed at high energies, can become sizable even at the energy level of future e + e − colliders.However, since EEC and E3C share a close power law in the leading power correction, the enhancement is significantly canceled when considering their ratio ∆ 3,2 (x L ).As shown in Fig. 7, the leading non-perturbative correction only gives rise to roughly 4% effect at Q = 250 GeV and 2% at Q = 1 TeV for matched NNLL.This confirms that ∆ 3,2 (x L ) is insensitive to the hadronization and indeed  a good candidate for precise α s measurement.
We also investigate the impact on the final resummation results caused by the uncertainties from the two-parameter fit.The statistical error for both Λ and γ are given in Eq. (3.5) and (3.6).Fig. 8 shows the final uncertainty in the matched NNLL distribution from varying these two NP parameters.In both Q = 250 GeV and Q = 1 TeV, excluding the shaded region, the NP uncertainty is much smaller than the hard uncertainty estimated by seven-point variation.In particular, at Q = 1 TeV, the NP uncertainty is reduced to 1% in the potential fit region.Despite that, we admit that the effect of non-perturbative corrections turns to increase for such small x L region, and more accurate understanding of the non-perturbative corrections will be required to further improve the precision.The uncertainty from varying the non-perturbative parameters Λ and γ in the resummed ratio distribution ∆ 3,2 (x L ).The bottom panels stand for the relative difference of NNLL+NLO+NP with respect to its central value.

Anticipation of α s determination
In this subsection, we discuss the potential of extracting the strong coupling constant α s from measuring the resumed E3C/EEC ratio ∆ 3,2 (x L ).In literature [100,101], the back-toback limit of EEC is resummed to NNLL+NLO and has been use for α s measurement from e + e − data.Similar to other event shapes, the non-perturbative correction is significantly large in this region and require careful modeling.And how we profile the resummation and power correction has a sizable effect on the final theory uncertainty.Alternatively, we can also do the α s measurement only in the collinear limit.First of all, as we discussed in Sec.3.1, the non-singular contribution is almost zero in this limit, and thus it is safe to ignore the higher fixed-order contribution.Secondly, by considering the ratio distribution, ∆ m,n (x L ), the suppressed power corrections will lead to a smaller theory uncertainty and thus more precise α s determination.As illustration, we first investigate the sensitivity of ∆ 3,2 (x L ) when slightly changing the value of α s .In particular, we vary the value of strong coupling at Z-pole α s (m Z ) by a factor of 5%, namely α s (m Z ) = {0.112,0.118, 0.124} and compare the effect on matched resummation result.
We first consider the NNLL+NLO ∆ 3,2 (x L ) at Q = 91.2GeV with all three values of α s (m Z ).As observed in Fig. 9, the slope become sensitive to the α s in the collinear region x L = 10 −3 ∼ 10 −4 , while the relative difference with respect to α s (m Z ) = 0.118 ranges from 10% to 20%.The slope sensitivity and the cancellation of hadronization correction have made the ratio of E3C and EEC ∆ 3,2 (x L ) an advantageous observable for extracting the α s from e + e − annihilation.Similar behaviors also exist at other energies and for completeness, we present the comparison at Q = 250 GeV and Q = 1 TeV in Fig. 10.The fact that the resummed E3C/EEC ratio has larger sensitivity to α s and reduced non-perturbative corrections in the collinear limit makes it a promising candidate for the α s determination.To further improve the α s determination requires improving the resummation accuracy, matching with NNLO fixed-order correction, as well as the non-perturbative modeling.

Approximate NNLL resummation in pp collisions
In this section, we consider the dijet production pp → jj at the LHC.There are several motivations to study energy correlators in pp collisions.First of all, LHC provides unique opportunities to study energy flows correlation in QCD at extremely high energy.While the LEP or future CEPC provides a very clean environment for precise measurements, pp collisions at the LHC can produce multiple jets with very high energies (p T ≳ 500 GeV), and high angular resolution can be achieved to probe the underlying dynamics for their formation and evolution.Secondly, as we have observed in the e + e − collisions, the non-perturbative corrections for ENC have a relatively simple form compared to other event shape observables (at least in leading power), which might be easier to study nonperturbative QCD.At the same time, with multiple scales involved, pp collision can provide robust data from high energy to low energy, which is beneficial for understanding nonperturbative effects.
In this section, we still focus on improving the perturbative predictions for ENC.As in Sec. 2, the jet functions are universal across different hard processes and the new ingredients are the moments of pp hard function, both regular and logarithmic.The main complication for pp collision is that the hard function now involves convolution with PDFs and algorithmic definition of jet, allowing only numeric calculation of the hard function.
For the numerical calculation of the hard function, we adopt the anti-k t jet algorithm and choose the jet radius to be R 0 = 0.4.The complete kinematic cuts are summarized in Eqs.(2.8)-(2.9).The µ independent part of the NLO hard function are presented in Appendix.A. We observes large corrections going from LO to NLO.The µ dependent part of the NNLO hard function can be derived using the RG equation in (2.6).The µ independent part requires a genuine two-loop corrections and are beyond the scope of this work.Instead we make a simple estimate of the two-loop constant terms, and dubbed the resulting prediction approximate NNLL resummation (NNLL approx ).Specifically, we use a modified Padé approximation to estimate the two-loop hard function constants in both quark channel and gluon channel: where we vary κ in the range [0, 1/2] as a naive way to estimate our theory uncertainties on the missing two-loop constants.For the splitting function, β function, as well as the jet functions, we used the ones required by NNLL accuracy as shown in Table 1.
In Fig. 11, we show the E3C/EEC ratio ∆ 3,2 (R L ) up to NNLL approx , with the hard uncertainty estimated by seven-point variation.Due to the lack of knowledge of the genuine two-loop hard function moment, we have chosen to normalize the E3C/EEC distribution in the range of R L ∈ [0.01, 0.4] to reduce the impact from not knowing the full two-loop hard function.We find good convergence for both p t ranges: [300, 350] GeV and [500, 550] GeV.In the future, it would be interesting the compute the two-loop hard function, as well as match the resummed results to fixed order to improve the prediction around R L ∼ R 0 .

Anticipation of α s determination
Similar to e + e − annihilation, in this subsection we discuss the potential of extracting the strong coupling constant α s from the resummed ∆ 3,2 (R L ) distribution in pp → jj.In particular, we also investigate the slope sensitivity of the distribution with respect to different values of α s .For hadron colliders, we need to change the PDFs as we vary the strong coupling among α s (m Z ) = 0.118 ± 0.06.For this purpose, we use three PDF sets: NNPDF31_nnlo_as_0112, NNPDF31_nnlo_as_0118 and NNPDF31_nnlo_as_0124 when calculating the hard function using the method in [54].
As shown in Fig. 12, for each p t range, the uncertainty is significantly reduced from NLL to NNLL approx , leading to distinguishable slopes with respect to different α s .This suggests that ratios of energy correlators are good candidate for extracting α s .We note that there is larger slope variation for lower p t of the jet, in agreement with the expectation that the measurement at lower energy is more sensitive to α s due to asymptotic free nature of QCD.

Conclusion
In this paper we have performed a systematic study of resummation of projected three-point energy correlator E3C [16], and its ratio to EEC, in both e + e − collider and pp collider.We  have achieved the first NNLL accuracy for the e + e − case, and NNLL approx accuracy for the pp case.Our results show that good perturbative convergence can be achieved for the ratios of projected energy correlators.The current theoretical uncertainties are at a level of a few percent, and can be further improved in the future when the higher order ingredients become available.We have also shown that the ratio observable is sensitive to variation of α s , therefore provides a good candidate for precision α s determination using jet substructure.To achieve the above theory accuracy, one of the main new ingredients is the two-loop E3C jet function computed in this work.The calculation includes three pieces: double-real, real-virtual and double-virtual.The last two contributions only involve a single δ measure-ment function in the phase space integral and share a similar form as the analytic EEC calculation at NLO [5].Regarding the double-real emissions, which amounts to integrating the fully-differential EEEC distribution within the collinear kinematic space, we used two different approaches and find the same results.The first method is to subtract the infrared divergence in the collinear EEEC jet function, integrate it separately with d-dimension kinematic space, and expand the finite terms in ϵ.The second approach benefits from the recently developed parametric IBP, where we can also simplify the integrand with IBP reduction and calculate the integrals via differential equations.
Regarding the ENC resummation, for e + e − annihilation, we solve the E3C jet RGE (which is a modified DGLAP equation) order by order in α s with the two-loop boundary, and push the resummation up to NNLL.For pp collisions, we calculate the combined hard function moments using the method in [54] for dijet production.We present the complete NLL and the approximate NNLL resummation result, where the approximation is due to the missing of genuine two-loop hard function constant.The uncertainty is reduced compared with the previous results [16,20,21].For the fixed-order matching, we notice that the singular contribution dominates the collinear limit and the non-singular contribution from matching has only small effects in the e + e − case.Nevertheless, we perform the matching for e + e − given the fixed-order result is already available, but leave the matching with fixed-order in the pp case for the future study.
For a complete phenomenological analysis and precise α s extraction at hadron collider, there are still several ingredients needed in the future.Perturbatively, we need to compute both two-loop hard function and the NLO non-singular distribution for pp → jj, in order to achieve a full NNLL story.More over, it would be interesting to solve the RG equation exactly following [102], and compare the results with the truncation method.At the same time, for both e + e − and pp, it would be interesting to better understand the hadronization power corrections to help further reduce theoretical uncertainties.We hope that all these efforts can lead to a precision determination of α s from jet substructure in the future.
to two-loop, single logarithmic up to one-loop and the double logarithmic moments at tree level with respect to the energy fraction x: For EEC (N = 2), we have Note that the EEC hard moments are also summarized in the appendix of Ref. [17]).However, the normalization condition in [17] is different from ours, due to the scaled energy E i /(Q/2) there in contrast with E i /Q here in the definition of the jet function.For E3C (N = 3), we find For completeness, we also provide the E3C (N = 3) hard moments for the gluonic Higgs decay, which is needed for extracting the two-loop gluon jet constants.Here we use h to distinguish from the e + e − case.As one of the important checks of our calculation, we show in Fig. 13 the independence of the slicing parameter δ cut when evaluating the hard function moments using the method in [54].The values of the moments are in agreement within the numeric uncertainty for three values of δ cut across two orders of magnitude, namely δ cut ∈ {0.003, 0.03, 0.3}.Here the dot also represents the derivative with respect to N .Note that {i, j} = {q, g} and the anomalous dimension is a 2 × 2 matrix.The results for EEC (N = 2) are derived and summarized in the appendix of Ref. [17], so here we provide the expressions for E3C (N = 3).

B β-function RGE and running coupling
The well-known QCD β-function is written as where the coefficient up to three loops are given by [103][104][105][106] At one-loop, the β-RGE can be solved exactly.At two-loop and beyond, there are different solutions.In terms of L ≡ ln µ 2 Λ 2 QCD , a expanded solution can be written as: Here we can obtain the two-loop running coupling for NLL resumation by setting β 2 = β 3 = 0 and three-loop running coupling for NNLL by only β 3 = 0.Alternatively, one can iteratively solve the RGE order by order in a formal expansion parameter ϵ ∼ βn β 0 , with n ≥ 1.For NLL, the two-loop running coupling is written as and at three loops for NNLL (B.5)For the resummation in this paper, we use the iterative solution (the latter one) and set the coupling at Q = 91.2GeV to be the world average value α s (m Z ) = 0.118.

C Squeeze limit of EEEC jet functions
In this section, we provide the perturbative data for the squeeze limit of the EEEC jet function in Eq. (2.17), which is needed for E3C jet function calculation.Given the conformal parameterization, the squeeze limits correspond to z → 0, 1, ∞, related by a S 3 symmetry.Without loss of generality, we provide the z → 1 limit for the shapes function up to O(ϵ 2 ).In the quark jet, we find for G(z) and for F (z): as well as the H(z): Here the red stands for the most singular term, which contributes to 1 ϵ divergence in the E3C jet function calculation.For the gluon jet, we also find

D Result of two-loop E3C jet function calculation
We list the individual results for the two-loop jet function calculation in Sec.2.4.As we discussed above, the calculation is reorganized as nonidentical energy weight contribution and contact terms.For the nonidentical energy weight in Sec.2.4.1, we find for the quark jet where λ is the effective Hgg coupling 5 [107].These results are then used to extract the two-loop jet constants.

E Fixed-order expansion
In this section, we provide the singular expansion of projected energy correlator up to NNLO O(α 3 s ) in e + e − annihilation.This can be achieved by expanding our resummed distribution with canonical scale µ = Q.For EEC, we find 1 σ 0 dσ [2]  dx

Figure 2 :
Figure 2: The comparison of fixed-order result and the singular expansion from resummation.The left panel shows good agreement between LO expression and the O(α s ) singular expansion of E3C resummed result.The difference between them is the non-singular result.The right panel gives the corresponding contributions at NLO, with the numerical full NLO prediction from Event2.

Figure 3 :Figure 4 :
Figure 3: The resummed E3C distribution up to NNLL+NLO, multiplied by a factor x L (1 − x L ), for e + e − collision at 250 GeV (left top panel) and 1 TeV (right top panel).Uncertainty bands are obtained by varying the hard scale µ around the nominal value Q jet = Q/2 by a factor of 2. The lower panels show the relative scale uncertainty of the NNLL+NLO distribution around the central value.

Figure 5 :
Figure 5: Comparison of Pythia8 result and the fitting using Eq.(3.4).The blue curves are the difference of hadron-level and parton-level distribution for EEC and E3C, at both 250 GeV and 1 TeV.The red curves are our fitting result with parameters from Eq. (3.5) and(3.6).The shaded region stands for the range where the parametrization of the leading non-perturbative correction is no longer valid and should be excluded from the fit range.

Figure 6 :
Figure 6: The E3C distribution to NNLL+NLO including the non-perturbative (NP) hadronization corrections estimated with Pythia8 data, with different plot ranges for collision energies at 250 GeV (left panel) and 1 TeV (right panel).

Figure 7 : 2 L
Figure 7: The E3C/EEC ratio to NNLL+NLO with non-perturbative (NP) hadronization corrections.The non-perturbative hadronization corrections are estimated as above by Λ/(x 3/2 L Q) for both the EEC and E3C distributions with the coefficients fitted from Pythia8.

Figure 8 :
Figure8: The uncertainty from varying the non-perturbative parameters Λ and γ in the resummed ratio distribution ∆ 3,2 (x L ).The bottom panels stand for the relative difference of NNLL+NLO+NP with respect to its central value.

Figure 9 :
Figure 9: Left panel is the matched NNLL ratio distribution ∆ 3,2 (x L ) with different strong coupling constants: α s (m Z ) = {0.112,0.118, 0.124} at Q = 91.2GeV.The uncertainty is the hard scale variation.Right panel shows the relative deviation of all three bands with respect to the central prediction of α s (m Z ) = 0.118.

Figure 10 :
Figure 10: Left panels are the matched NNLL ratio distribution ∆ 3,2 (x L ) with different strong coupling constants: α s (m Z ) = {0.112,0.118, 0.124} at Q = 250 GeV and Q = 1 TeV.Right panels are the relative deviation of all three bands with respect to the central prediction of α s (m Z ) = 0.118.

Figure 12 :
Figure 12: Normalized NLL (upper) and NNLL approx (lower) resummation result for E3C/EEC ratio at pp → jj with α s (m Z ) = 0.118 ± 0.06, i.e., varied by about 5%, for two different jet p t ranges [300,350] GeV and [500,550] GeV.Lower panels show the relative difference from the result at the central scale with α s (m Z ) = 0.118.

Figure 13 :
Figure 13: NLO hard function moments for N = 2 (left), and N = 3 (right), with p t ∈ [300, 350] GeV and δ cut ∈ {0.003, 0.03, 0.3}.The lower panels show the relative variation of the moments compared with the average value for three δ cut .The error bars represent the Monte-Carlo numeric uncertainty given by Madgraph5.

Table 2 :
The following table gives the hard function moments for pp → jj calculated in Madgraph5 in two different p t ranges: [300, 350] GeV and [500, 550] GeV, needed for the resummation of both EEC (N = 2) and E3C (N = 3).Values for hard function moments in pp collision for different p t ranges.The NLO corrections turn out to be significant.
Regarding the contact term in Sec.2.4.2, for e + e − annihilation, we have the sum of E 2 EC and E 3 C with the gluonic singular term r g (µ, Q, ϵ)r g (µ, Q, ϵ) = C A T F n f − 2 δ(x L )f g (µ, Q, ϵ) + 1 x L C F T F n f − For the case of gluonic Higgs decays, we normalize the E3C into the form where the LO E3C is from Eq. (2.40)-(2.41).