Anti-$k_T$ jet function at next-to-next-to-leading order

Jets constructed via clustering algorithms (e.g., anti-$k_T$, soft-drop) have been proposed for many precision measurements, such as the strong coupling $\alpha_s$ and the nucleon intrinsic dynamics. However, the theoretical accuracy is affected by missing QCD corrections at higher orders for the jet functions in the associated factorization theorems. Their calculation is complicated by the jet clustering procedure. In this work, we propose a method to evaluate jet functions at higher orders in QCD. The calculation involves the phase space sector decomposition with suitable soft subtractions. As a concrete example, we present the quark-jet function using the anti-$k_T$ algorithm with E-scheme recombination at next-to-next-to-leading order.

Jets constructed via clustering algorithms (e.g., anti-kT , soft-drop) have been proposed for many precision measurements, such as the strong coupling αs and the nucleon intrinsic dynamics. However, the theoretical accuracy is affected by missing QCD corrections at higher orders for the jet functions in the associated factorization theorems. Their calculation is complicated by the jet clustering procedure. In this work, we propose a method to evaluate jet functions at higher orders in QCD. The calculation involves the phase space sector decomposition with suitable soft subtractions. As a concrete example, we present the quark-jet function using the anti-kT algorithm with E-scheme recombination at next-to-next-to-leading order.

I. INTRODUCTION
Jets have been the center of attention at many frontiers of high-energy particle and nuclear physics. At the Large Hadron Collider (LHC), for many years, jets have served as indispensable tools for precision test of the Standard Model (SM) and new physics searches at the TeV scale. Measurements at the LHC and the Relativistic Heavy Ion Collider (RHIC) demonstrate that jets can also be unique probes of non-perturbative dynamics, such as collinear parton distribution functions (PDFs) [1] as well as transverse momentum dependent ones (TMDPDFs) [2], the intrinsic spin of the nucleon [3][4][5] or the hot medium effects of the quark-gluon plasma [6][7][8]. These efforts are supported by theory developments of novel jet substructures [9,10], and enable accurate extractions of the SM parameters [10,11] and the three-dimensional tomographic images of nucleons [12][13][14][15][16][17][18] out of the jets. All these studies will receive a further boost at the future Electron-Ion Collider (EIC) [19,20], including the possibility to consider polarization degrees of freedom.
The theoretical foundation of such precision measurements roots in the factorization of the cross section dσ with N exclusive jets with small jet radius R, whose schematic form is [21][22][23] where F encodes the parton distributions of the colliding nucleons, such as the collinear PDFs or TMDPDFs, in case transverse momentum dependence is kept. The traces are over the color matrices. H is the hard function that describes the hard interaction which initiates the process while S G includes the wide angle radiations 2 n µ J , with n J a light-like vector along the jet axis, to split into m particles with momenta {p 1 , . . . , p m } and Θ jet alg. enforces that all the m collinear particles are clustered within one jet.
The measurements of jets and jet substructures are formulated in the angular convolution ⊗ Ω of the energetic collinear jet function J m and the collinear-soft function S cs,m , where m is the energetic collinear parton multiplicity. The angular convolution guarantees the proper inclusion of the non-global logarithms (NGLs) [24] if present. We remark that the multi-emission feature of the jet cross section in Eq. (1) is also shared by the forward scattering in the small-x framework [25]. The necessity for increasing theoretical precision for jet cross sections, both at fixed order perturbation theory and with resummation, requires the evaluation of each component in the factorization Eq. (1) at least to next-to-next-to-leading order (NNLO).
On the other hand, many theoretical predictions of the jet measurements are currently based on the simplified formula [13][14][15][26][27][28][29][30] where one decouples the angular correlation between J m and S cs,m and averages each function over its solid angle in Eq. (2) to have J = m J m Ω and S cs = S cs,1 Ω . As indicated by the factor e L ngl , the leading NGLs are resummed through the exponential. In order to reliably justify the logarithmic accuracy of the simplified form, a direct NNLO calculation of the factorization is also desired. For many processes, the hard and the soft functions are already known to two loops or even beyond, see, e.g. [31]. However the two-loop information for the jet function is still missing and therefore limits the theoretical accuracy of the resummation to the next-to-leading logarithm (NLL). Compared with the other components in the factorization, the computation of the jet function beyond next-to-leading order (NLO) is dramatically complicated by the presence of the recursive clustering procedure, which is the essence of the widely used algorithms for jets and jet substructure, such as anti-k T [32] with either the E-scheme or the winner-take-all (WTA) recombination scheme [33,34] and soft-drop grooming [9].
In this work, we propose a method to evaluate J m or J at NNLO efficiently using sector decomposition for the phase space together with suitable designed soft subtractions. We demonstrate the feasibility of the method by presenting the first anti-k T quark-jet function at O(α 2 s ) using E-scheme recombination as a concrete example.

II. NLO
We take the well-known NLO calculation [15,26] of the quark-jet function as a warm up to illustrate our approach. The NLO calculation involves evaluating the integration of the matrix element for q a → q i g j over the collinear two-body phase space, which reads in d = 4−2 dimensions where Z α is the α s renormalization factor, which at NLO is s ij = 2p i · p j and z is the momentum fraction carried by the gluon g j andz = 1 − z. The first line of Eq. (4) is the phase space measure in the collinear limit and we note that the matrix element in second line of Eq. (4) is nothing but the leading order (LO) splitting function 8παsµ 2 sij P q→qg . The θ-function originates from the requirements of the jet algorithm to cluster two partons into one single jet with radius R. Here where ∆η ij and ∆φ ij are the rapidity and azimuthal angle differences between parton i and j, respectively. p T is the jet transverse momentum and p i,T and p j,T are the transverse momenta for partons i and j. The narrow jet approximation R 1 has been used.
If we introduce the variables we can write where L = log µ p T R . All the singularities are given by x i → 0 and have been isolated in the factor so that the remaining terms in the second line of the above equation are finite as x i → 0.
With this manipulation, to evaluate the integral, we can first expand x −1−ai i using the Laurent expansion and then perform the integration either numerically or analytically to find the coefficients of the -poles and the finite terms. This reproduces the well-known NLO bare quark-jet function We note that at NLO, all k T -type jet algorithms give the same results.
Although we present here only the x i -integrated jet function, our approach is actually fully exclusive in x i . There are no difficulties to generate distributions differential in x i , as can be seen from Eq. (9), especially in the angular variable x 1 =s ij . The implications of the possibility to provide exclusive results, especially their role in resumming non-global logarithms, will be discussed later.
A similar strategy will be used to compute the NNLO jet functions but the calculation is more involved, as we will explained in detail in the rest of the work.

A. Real-virtual contribution
Now we proceed to calculate the jet function at NNLO. We start with the real-virtual corrections. The phase space of the real-virtual corrections is identical to the NLO calculation. The matrix element for the jet function is given by the one-loop correction to the splitting kernel, which is well documented in [35][36][37][38][39]. The calculation for the real-virtual correction is straightforward and gives with and B. Real-real contribution

Matrix elements
We turn to the calculation of the real-real corrections. The matrix element involved is nothing but the tree level a → ijk splitting kernel, where s ijk = s ij + s ik + s jk and z i is the momentum fraction carried by parton i, with z i + z j + z k = 1. The explicit form of the splitting function P a→ijk (z i , z j , z k ) can be found in [40] (see also [41]). For the quark-jet we need q1q2q3 , P (ab) g1g2q3 , and P (nab) g1g2q3 , (15) and explicit expressions are listed in the appendix A. The matrix element |M| 2 contains structures which will develop singularities in if one or more of the s ij 's vanish upon integration over the phase space. To extract the real-real corrections, our strategy is to write the phase space integration as where F is regular as x i → 0. Hence we can first perform the Laurent expansion in x −1−ai i to isolate thepoles and compute their coefficients and the remaining constant terms either numerically or analytically. We emphasize again that Eq. (17) allows not only for the integrated jet function but also the differential distributions in x i 's. To arrive at the expression in Eq. (17), we need a proper parameterization of the phase space.

Phase Space Parameterization
In the collinear limit, the three-body phase space can be shown to take the form [39,42] where the Gram determinant is given by with z k = 1 − z i − z j . Now we introduce the angular variable between two arbitrary partons a and b and further make a variable transformation using one of thes abs with t ∈ [0, 1] to find Ifs ik happens to be in the denominator of the matrix element, the variable change in Eq. (21) will develop a line singularity, which can be mapped tos ij =s jk by a non-linear transformation [43,44] and thus which gives Here we have introduced the variable x 5 through Without loss of generality, we can assume z i ≤ z j and s ij ≤s jk , and in this case we can then introduce 1 1] to reach the phase space parameterization Other orders can be obtained in the same manner. In the third line of the phase space parameterization, the Laurent expansion in Eq. (9) can be performed to extract the -poles and the finite contributions. On the other hand, the last line of Eq. (28), when multiplied with the matrix element, will be free of any singularities as x i → 0 or x 2 → 1. It is suitable for a numerical evaluation, as long as appropriate subtraction terms are constructed as will be explained in the following section.

Jet algorithm and subtraction terms
In this work, we are always interested in k T -type jet algorithms, which involves the comparison of the metrics for some parameter α and α = 1 defines the anti-k T jet algorithm. In the small jet radius limit, or in the collinear limit, we can approximate p T,i = z i p T up to power corrections. Furthermore, if we demand all three partons to be grouped into the same jet, we will require in the last stage of the clustering, where we have assumed that the partons i and j are clustered first into the pseudo-particleĩj, which then clusters with parton k to form the final jet.
Within the E-scheme recombination, the fourmomentum of the pseudo particleĩj is given by p μ ij = p µ i + p µ j and thus From this, it is easy to see that the requirement in Eq. (30) of all three partons being clustered into one jet can be written as valid in the small-R limit. Here again,s ij ≡ sij zizj In the earlier history of the clustering, in order to have i and j to be clustered first, we need as required by the jet algorithm. This maps to the conditions for the following cases: • case I: • case II: • case III: Other cases can be obtained by switching the order between i and j.
The jet clustering conditions can be expressed as constraints in the phase space integrals. For case I, we have where dΦ 3 is the three-body phase space given in the previous section and c ijk is given in Eq. (32). We note that the θ-functions in the first line due to the clustering condition are dangerous, in the sense that they will change the pre-determined fractional power of in the exponent in Eq. (28), associated with the z j → 0 poles and therefore invalidate the Laurent expansion in Eq. (9).
To understand the issue, we consider a simple example by comparing the integrals and Apparently we can safely Laurent-expand the integrand in I 1 to perform the integration. This naive expansion is not allowed in I 2 , since the correct power of for x 3 is not −a 3 but −(a 3 + a 2 α) , which can be easily seen by integrating out x 2 . To avoid this complication, we introduce a subtraction term, to remove all possible the z j → 0 singularities by taking the z j → 0 limit in the clustering condition to find where we have used the fact that as z j → 0,s ij → 0 as forced by the θ-functions, and therefores jk →s ik . With the subtraction, the integral (dΦ I 3 − dΦ I 3,sub. )|M| 2 will be free of the z j → 0 poles. The term dΦ I 3,sub. |M| 2 alone still contains the z j → 0 singularity and therefore the potential fractional power modification due to θ s ij < z 2α js jk . But we will see later that when combined with other subtraction terms introduced soon, all such dangerous θ-functions go away.
We have, for case II and case III, the phase space integrals and respectively. Here for later use, we have relabeled subscripts i, j and k to stick to the case that z i ≤ z j ≤ z k , but with i, k clustered first in case II and j and k first in case III. This is always allowed since eventually we sum over all possible orders of z i 's and all possible clustering histories. Following the same logic as in case I, we construct the subtraction terms for case II and case III respectively, to remove the z j → 0 singularities. The phase space integral with the jet algorithm implemented reads and we note that by construction, it is easy to check that where all θ-function constraints involving powers of z α j cancel in the sum. Therefore, the phase space integration does not induce any modifications of the fractional powers in as z j → 0. However at this stage, the phase space parameterization in Eq. (28) is still not sufficient to isolate all singularities within |M| 2 . Particularly in III i=I dΦ i 3,sub. , a further ordering between the variables x 2 and x 4 is required to fully disentangle the poles generated by s ijk = z i z ksik + z j z ksjk + z i z jsij → 0. To avoid such an ordering, we introduce instead a subtraction term in the matrix element to regulate the single soft limit, as x 4 → 0, of |M| 2 . Thus, a natural choice for the subtraction term |M| 2 sub. is which is nothing but the product of the LO eikonal factor and the LO splitting kernel. The final phase space integration suitable for a numerical evaluation is thus given by Similar subtraction counter terms can be constructed for the WTA scheme recombination as well as other jet algorithms such as the Cambridge-Aachen algorithm and thus used to compute the soft-drop groomed jet function. Further applications to obtain the semi-inclusive jet function [45,46] at O(α 2 s ) are also straightforward.

Results and validations
With the setups developed in the previous sections, we find the double-real correction to the anti-k T jet function to be where we find and K rr N F T F = − 1 6 3 − 7 9 2 + 0.1067(3) + 17.230 (2) . (48) Here and below we always include in brackets the numerical error on the last digit.
To validate our results, we have performed independent calculations and found full agreement. Also we note that parts of our calculations associated with dΦ 3,sub. can be done (semi-)analytically and for these parts, we have checked our numerical computation against the (semi-)analytic results to find full agreement. Some of these checks are presented in the appendix C.

Renormalization group equation
The renormalization group equation (RGE) for the quark-jet function serves as strong independent check on the calculation. The leading poles up to −2 of the double-real correction in Eq. (46) combined with the realvirtual one in Eq. (11) can be predicted by solving the RGE up to the logarithmic term α 2 s L 2 and taking the NLO quark-jet function from Eq. (10) as input. Details are given in the appendix B.
As a result, the RGE predicts the leading three poles normalized to α 2 s /(2π) 2 for the C 2 F term as which all agree exactly with Eqs. (12) and (47), in the case of the −2 pole within numerical errors as shown in the second line in Eq. (51). For the C F N F T F contribution, solving the RGE we find the leading poles which fully agrees with what we obtained in Eq. (48).
As for the C F C A term, the naive RGE alone can predict the −4 and −3 poles correctly, however the −2 term is affected by the non-global contribution, for which we find that an additional − π 2 12 2 correction is needed to produce the correct pole and the logarithm at α 2 s L 2 . The same situation has been encountered in the Sterman-Weinberg cone jet in [21] and for the hemisphere mass distribution in [47].
Once every piece is taken into account, the poles predicted for the C F C A contribution are where the agreement is within numerical errors. For the −2 term it is at O(1‰) as shown in the second line in Eq. (56).

IV. TWO-LOOP JET FUNCTION AND DISCUSSIONS
Now we present the anti-k T bare quark-jet function at two loops where J (1) bare has been given in Eq. (10). The new two-loop result reads and with additional O(α 2 s ) contributions coming from the α s renormalization of J (1) bare in Eq. (10). The renormalized jet function J is then found to be where we have set the scale µ = p T R. The associated numerical errors originate entirely from the NNLO realreal corrections. Before concluding we make two comments on our calculation and the result. First, although we present m=3 are already encoded in our calculation, which can be easily seen from the fully differential nature of the phase space sector decomposition and the relation that J = m J m Ω . To demonstrate such feasibility, we plot in fig. 1, the distribution ins ij for the 0 part of the jet function J (0) m=3 with the P (id) q1q2q3 contribution only. A detailed study of the angular-dependent jet function will be left for future work. The results for the individual J m will be more suitable for resumming NGLs and for improving the logarithmic accuracy of the predictions.
Second, our two-loop result confirms that the factorization theorem in Eq. (3), as used in many phenomenological studies of jets, TMDs and in spin physics, is only valid up to the single logarithmic level (O e α n s L n ) due to the existence of the NGLs, which is generic for many exclusive jet processes. Beyond this order, the factorization in Eq. (1) or alternative approaches [48] should be used for predictions. Our calculational framework can be useful for realizing the resummation in Eq. (1). We leave this for future studies.

V. CONCLUSIONS
In this work, we have developed a method to calculate the jet functions, which arise in the factorization of cross sections with exclusive jets of small radius R, with a jet clustering procedure up to NNLO. Explicit results at NNLO have been presented, for the first time, for the quark-jet function with the anti-k T jet algorithm. The new two-loop contribution to the anti-k T quark-jet function provides the missing input for the evaluation of the factorization theorem in Eq. (1) to NNLO and for resummation beyond NLL accuracy. This will substantially improve the precision of predictions for observables with jets as hard probes. The computational framework is not limited to anti-k T E-scheme jets but readily applicable to many other jet observables up to NNLO, such as the semi-inclusive jet function, the WTA jet and the soft-drop groomed jet radius R g and the jet momentum fraction z g . In addition, the obtained two-loop quark-jet function will serve as an important building block when using the factorization in Eq. (1) for the development of a colorful q T -subtraction scheme, which is applicable to both, calculations at fixed order in perturbation theory up to NNLO and to parton shower simulations.
In summary, given the broad interest in jets from both the high-energy particle and the nuclear communities, we expect the results to have numerous phenomenological applications at colliders, such as the LHC and a future EIC for the study of jet cross sections, TMDs and in addressing the spin structure of the nucleon.
with the abbreviation For the final state with identical quark flavorsq 1 q 2 q 3 , we need where P (id) q1q2q3 accounts for the interference term The splitting function for the quark decay into g 1 g 2 q 3 is given by where the Abelian contribution is Using the fact that d log µ = β −1 [α s ] dα s , with β[α s ] = −2α s n=0 β n αs 4π n+1 , the RGE can be solved to find Expanding the solution to O(α 2 s ), we find the logarithmic terms in the jet function at α 2 s contain where J (l) 0,ren are the renormalized jet functions at NLO (l = 1) and at NNLO (l = 2) with µ = p T R and → 0. On the other hand, if we assume the bare jet function have the form and hence the NLO renormalization factor of the jet function is It is straightforward to check that J 0,ren . From these, the logarithmic contributions can be obtained via The coefficients of the -poles J −n with n = 4, 3, 2 serve as a check, although due to the existence of the non-global logarithms, the −2 -pole of C A C F contribution will be affected and will deviate from the RGE prediction. The term J (2) −1 allows for the extraction of the two-loop expression for the anomalous dimension of the quark-jet γ 1 , which can be determined as With quark-jet function computed to NNLO in Eq. (61) and all anomalous dimensions known to the relevant order we can also predict the logarithmic structure up to the next-to-next-to-leading logarithm at O(α 3 s ) based on the naive RGE, which reads 0,ren + where sub-leading logarithmic terms require additional inputs from higher order calculations. Finally, we note that in dimensional regularization the -poles of the quark-jet function display the a pattern of exponentiation analogous to the quark form-factor, cf. e.g. [50], with the cusp anomalous dimensions determining the leading poles and the anomalous dimension of the quark-jet γ Jq governing the sub-leading ones.
The analytic computation of some of the phase space integrals associated with Eq. (43) is straightforward, specifically for integrals like For (a) we directly obtain while for (b), we find By splitting up the integral into two regions,s ik ≤s jk ands ik >s jk , we have fors ik ≤s jk , we lets ik =s jk u 2 , where we have used 2 F 1 1, Fors jk <s ik , we lets jk = u 2s ik and obtain 1 0 dz i dz j (4π) 4−2 (z i z j z k ) 1−2 θ(z i < z j < z k ) The results of (a) and (b) in Eqs. (C3) and (C8) can be used to validate the numerical calculations. For instance, for the term proportional to 1 s12s13 in the splitting function P we have the following result for z 1 < z 2 < z 3 , by using (b) in Eq. (C8) and after integrating over z 1 and z 2 , For z 2 < z 3 < z 1 , using (a) in Eq. (C3) and integrating over z 2 and z 3 , we find The same check can be done for the splitting function P g1g2q3 . For instance for the contribution proportional to we find for z 1 < z 2 < z 3 , Similar other checks have also been performed but are not listed here.