Transverse Parton Distribution and Fragmentation Functions at NNLO: the Quark Case

We revisit the calculation of perturbative quark transverse momentum dependent parton distribution functions and fragmentation functions using the exponential regulator for rapidity divergences. We show that the exponential regulator provides a consistent framework for the calculation of various ingredients in transverse momentum dependent factorization. Compared to existing regulators in the literature, the exponential regulator has a couple of advantages which we explain in detail. As a result, the calculation is greatly simplified and we are able to obtain the next-to-next-to-leading order results up to $\mathcal{O}(\epsilon^2)$ in dimensional regularization. These terms are necessary for a higher order calculation which is made possible with the simplification brought by the new regulator. As a by-product, we have obtained the two-loop quark jet function for the Energy-Energy Correlator in the back-to-back limit, which is the last missing ingredient for its N$^3$LL resummation.

TMDPDFs and TMDFFs can be defined as hadronic matrix elements of bilinear quark or gluon field operators with a measured transverse momentum k ⊥ (or a transverse separation b ⊥ in position space). If the transverse momentum k ⊥ ∼ Λ QCD , the TMDPDFs and TMDFFs are essentially non-perturbative, and can only be extracted from experimental data or calculated using lattice methods. On the other hand, if k ⊥ Λ QCD , the TMDPDFs and TMDFFs can be related to the collinear PDFs and FFs via perturbatively calculable matching coefficients. These coefficients are known at the next-to-next-to-leading order (NNLO) for the TMDPDFs [37][38][39][40] and TMDFFs [40]. They have played an important role in a number of cutting-edge calculations, including precision predictions for the Drell-Yan process and Higgs boson production at small transverse momentum [10,29,30], and NNLO calculations for top quark pair production using the Q T subtraction method [41,42].
In this work, we revisit the calculation of the matching coefficients for TMDPDFs and TMDFFs at NNLO. We consider the quark TMDPDFs and TMDFFs in this paper, while the gluon case is left to a forthcoming article. There are several new elements in our calculation compared to those in the literature: • We employ the exponential regulator for rapidity divergences [43]. Rapidity divergences or "collinear anomalies" appear in the calculation of individual TMD functions in factorization formulas, which are cancelled in physical observables. These divergences are not regularized by dimensional regularization, and additional regulators need to be introduced [1,5,12,[44][45][46][47][48][49]. The exponential regulator has been shown to be particularly suitable in the calculation of TMD soft functions, as demonstrated in the recent next-to-next-to-next-to-leading order (N 3 LO) calculation [50]. We show in this work that the exponential regulator can also be used to calculate TMDPDFs and TMDFFs, which are more complicated objects than TMD soft functions. Our results show that the exponential regulator is a consistent rapidity regulator in both the soft and collinear sectors.
• We develop systematic calculation method based on modern techniques for loop integrals, such as integration-by-parts (IBP) identities [51,52] and differential equations [53][54][55]. Our method paves the way to calculate TMDPDFs and TMDFFs at N 3 LO.
• We obtain the bare NNLO TMDPDFs and TMDFFs up to O( 2 ), where is the dimensional regulator. They directly contribute to TMDPDFs and TMDFFs at N 3 LO upon renormalization.
• Our results for TMDPDFs agree with previous calculations [38][39][40], but we find a small discrepancy for the TMDFFs compared to those presented in Ref. [40]. We have performed several consistency checks on our results to make sure that they are correct.
• As a by-product, we obtain the NNLO quark jet function relevant for the resummation of EEC in the back-to-back limit. This is the last missing ingredient for this resummation at the next-to-next-to-next-to-leading logarithmic (N 3 LL) accuracy.
This paper is organized as follows. In Section 2 we introduce the definitions of quark TMDPDFs and TMDFFs in the context of the SIDIS process, and discuss the exponential regulator for rapidity divergences. In Section 3 and 4 we perform the calculation of the quark TMDPDFs and TMDFFs at NNLO using the exponential regulator. In Section 4 we also use the results for TMDFFs to compute the two-loop jet function for EEC in the backto-back limit. This is by itself a new result of our paper, and also serves as a cross-check of our results. We conclude in Section 5.

Kinematics and factorization
In this section, we briefly review the formalism of transverse momentum dependent factorization and introduce the definitions of TMDPDFs and TMDFFs. For our purpose, it is easiest to consider (unpolarized) SIDIS which involves hadrons in both the initial state and the final state. In SIDIS, a hadron N 1 with momentum P µ 1 is probed by a virtual photon γ * with momentum q µ and produces a jet containing a specific hadron N 2 with momentum P µ 2 . We define the kinematic invariants We introduce two light-like 4-vectors n andn satisfying n 2 =n 2 = 0 and n ·n = 2, such that we can decompose any 4-vector k µ as k µ =n · k n µ 2 + n · kn When quoting the components of a 4-vector, we use k = (k + , k − , k ⊥ ). The scalar product of two 4-vectors is given by In the hadron frame and ignoring the hadron masses, we have and we define q 2 ⊥ ≡ −q 2 T . The hadronic tensor is defined as In the region q T ∼ Q Λ QCD , the hadronic tensor can be factorized into products of hard kernels with collinear PDFs and FFs: where we have suppressed the dependence of the hard kernel on other kinematic variables.
In the language of soft-collinear effective theory (SCET) [56][57][58][59][60], the collinear PDFs and FFs can be defined as matrix elements of gauge-invariant collinear fields. For example, the bare quark collinear PDF and FF are defined by [1,59,61] where χ n and χn are the gauge-invariant collinear quark fields along the n andn directions, respectively. We have assumed dimensional regularization with d = 4 − 2 . The collinear PDF φ q/N 1 (x, µ) describes (in a sense) the probability distribution of finding the quark q with momentum fraction x inside the fast-moving hadron N 1 . The collinear FF d N 2 /q (z, µ), on the other hand, describes the probability distribution of finding the hadron N 2 with momentum fraction z inside the jet initiated by the quark q. If q T Q, however, the above picture of collinear factorization breaks down due to the appearance of large logarithms of q T /Q in the hard kernel H ij . One should instead rely on TMD factorization of the form Figure 1. Kinematics for TMDPDFs (left plot) and for TMDFFs in the parton frame (right plot).
where B i/N 1 , D N 2 /j and S ij are TMDPDFs, TMDFFs and TMD soft functions in the impact parameter space, with b ⊥ the impact parameter; whileB i/N 1 ,D N 2 /j andS ij are their counterparts in the transverse momentum space. For our purpose, we only consider i, j being quarks and anti-quarks. The quark TMDPDFB q/N 1 (x, k 1⊥ , µ) describes the probability distribution of find a quark with momentum fraction x and transverse momentum k 1⊥ inside the hadron N 1 , as depicted in the left plot of Figure 1. Naively, the bare quark TMDPDF can be defined by Similarly, the bare quark TMDFF may be defined as (2.11) Note that in the above definition, k 2⊥ represents the transverse momentum of the quark in the hadron frame (where N 2 has zero transverse momentum). In practice, it is also useful to define the TMDFFs in the parton frame where the quark has zero transverse momentum. In the parton frame, N 2 now has a non-zero transverse momentum P 2⊥ which is related to k 2⊥ by P 2⊥ = −zk 2⊥ . The parton frame quark TMDFF is then Here, P 2⊥ is defined with respect to the axis chosen such that the total transverse momentum of N 2 and X is zero. It is easy to show that (2.14) The functionF bare N 2 /q (z, P 2⊥ ) represents the probability distribution of finding a hadron N 2 with momentum fraction z and transverse momentum P 2⊥ inside the jet initiated by the quark q, as depicted in the right plot of Figure 1. From the above definitions, it is easy to see that Finally, the quark TMD soft function is given by the vacuum expectation value of a soft Wilson loop where the soft Wilson line is defined by with A s the soft gluon field in SCET.

Rapidity divergences and the exponential regulator
While the TMD factorization formula (2.9) makes some sense, the TMDPDF (2.10), TMDFF (2.11) and TMD soft function (2.16) are actually ill-defined due to the appearance of rapidity divergences which are not regularized in dimensional regularization. These divergences cancel when one combines the 3 functions in the factorization formula (2.9) to calculate physical observables. However, they also carry important information, just like the relationship between ultraviolet (UV) divergences and the renormalization group. The rapidity divergences arise due to the fact that the collinear modes and soft modes have the same typical off-shellness around q 2 T . More precisely, in the q T Q limit we have the relevant momentum regions collinear: p n ∼ Q(1, λ 2 , λ) , anti-collinear: pn ∼ Q(λ 2 , 1, λ) , soft: p s ∼ Q(λ, λ, λ) . (2.18) where λ = q T /Q 1. The effective field theory describing these modes are sometimes called SCET II . The collinear modes and the soft mode are related by a boost in the n orn direction. As a result, they cannot really be separated by a boost invariant regulator such as dimensional regularization. A brute-force separation as done in Eq. (2.9) then leads to inconsistencies manifesting themselves as rapidity divergences.
To deal with the rapidity divergences, one needs to introduce a regulator in addition to dimensional regularization. This however leads to another subtle issue. Any such regulator necessarily reintroduces a logarithmic dependence on the hard scale Q into integrals in the collinear and anti-collinear regions throughn · P 1 and n · P 2 , which was supposed to be factorized out into the hard function H in Eq. (2.9). This fact is sometimes called "collinear anomaly" or "factorization anomaly" in the literature [44,62]. Nevertheless, using the structure of the rapidity divergences, it can be shown that these Q-dependence can be extracted and exponentiated to all orders. After such a "re-factorization" The functions B q/N 1 and D N 2 /q can be regarded as the "genuine" quark TMDPDF and TMDFF which are free from rapidity divergences and are also independent of Q. The exponent function F qq is closely related to the so-called Collins-Soper kernel [1]. It has been known perturbatively to three loops [50,63]. Very recently, there are proposals to compute it non-perturbatively on the lattice [64,65]. In the literature, there are a variety of ways to regularize the rapidity divergences [1,12,[45][46][47][48][49]. In this paper, we consider the so-called exponential regulator [43] which was used to calculate the TMD soft function to the N 3 LO. We will show that it is a consistent regularization scheme also for the TMDPDFs and TMDFFs. Before discussing the exponential regulator, we briefly review the η-regulator of Ref. [46] which shares many similarities. At the next-to-leading order (NLO), the η-regulator amounts to the subsitution for the phase-space integrals over the real gluon momentum k µ , where δ + (k 2 ) = θ(k 0 )δ(k 2 ). The rapidity divergences appear as 1/η poles which can be subtracted in the same way as renormalizing the UV divergences. After the subtraction, the TMDPDFs, TMDFFs and TMD soft functions still depend on the "rapidity scale" ν. For the TMDPDFs and TMDFFs, the natural rapidity scale is ν ∼ Q, while for the TMD soft functions ν ∼ q T . The evolution equations of these functions with respect to ν can be used to exponentiate the rapidity logarithms ln(b 2 T Q 2 ) leading to the refactorization in Eq. (2.19).
While the above η-regulator is conceptually simple, it is not easy to implement in higher order calculations beyond NLO. For example, the regulator has to be carefully applied to maintain non-Abelian exponentiation in the soft sector [46,66]. In particular, when there are two real gluon emissions with momenta k 1 and k 2 , it is different to apply the regulator on k 1z + k 2z as a whole, or on k 1z and k 2z separately. Recently, a new regulator for rapidity divergences called "exponential regulator" has been proposed in Ref. [43], which leads to the same rapidity evolution equations as the η-regulator, and is easier for higher order calculations. In momentum space, the new rapidity regulator is simply multiplying each soft/collinear phase space measure by an exponential factor Note that the τ → 0 limit has to be taken after integration. Beyond NLO, when there are multiple soft/collinear partons, the regularization simple becomes Due to the exponential form, the multiple emission case naturally factorizes into products of single emissions. Therefore, non-Abelian exponentiation is manifestly preserved by this regulator. An important feature of the exponential regulator is that it leads to enormous simplification in perturbative calculations, as demonstrated by the calculation of TMD soft functions at N 3 LO in Ref. [50]. The exponential regulator also admits simple operator definitions for the TMD functions. For example, the quark TMD soft function is defined as where the rapidity regularization procedure is understood as keeping non-vanishing terms in the limit of τ → 0 (including the log τ terms which are the manifestation of rapidity divergences), and then identify the rapidity scale as ν = 1/τ . No subtraction is needed and the rapidity divergence are now renormalized. The remaining results depend on logarithms of the rapidity scale ν. Similarly, the exponentially regularized quark TMDPDF and TMDFF are defined as Note that for both the TMDPDF and TMDFF, we need to perform a zero-bin subtraction to avoid double-counting between the collinear sectors and the soft sector. The zero-bin soft function is the same as the TMD soft function Having operator definitions Eqs. (2.24), (2.25) and (2.26) for the TMD functions could be advantageous for studying non-perturbative aspects of TMD physics. In this work we focus on the perturbative part of the TMDPDF and TMDFF.

Renormalization and perturbative matching
For large impact parameter b T ∼ 1/Λ QCD , the TMDPDFs and TMDFFs are dominated by long distance contributions and are genuine non-perturbative objects. In this work, we are interested in the semi-perturbative region b T 1/Λ QCD . In this region the TMDPDF admits an operator product expansion where φ i/N is the (renormalized) collinear PDF of parton i, and I qi is a perturbatively calculable matching coefficient function describing the splitting of the parton i into the quark q. Similarly, the TMDFF can also be factorized as with perturbatively calculable coefficient functions C iq describing the fragmentation of the quark q into the parton i. The functions I qi and C iq will be the main objects we are going to study in this work. As indicated by the superscript "bare" in Eqs. (2.28) and (2.29), there are UV divergences which require renormalization. For the TMDPDF, we define the renormalization factor according to where we have used ⊗ to denote the convolution in Eq. (2.28). The matching coefficients I qi do not depend on the external state N , and can therefore be calculated with N replaced by a partonic state j = q or g. We can then extract I qi by calculating B bare q/j , performing the renormalization and subtracting the partonic collinear PDFs φ i/j . Up to the NNLO, the partonic collinear PDFs are given by ij is the LO splitting kernel and β 0 is the LO beta function. After renormalization, the TMDPDF obeys a renormalization group equation (RGE) where Γ cusp is the usual cusp anomalous dimension and γ B is the non-cusp anomalous dimension for the TMDPDF, whose perturbative expansions are collected in the Appendix. From the above equation and the famous DGLAP equation one can deduce the RGEs for the coefficient functions as Besides the normal RGE, the TMDPDF and the coefficient functions also satisfy the rapidity evolution equation [46] (2.35) The rapidity anomalous dimension γ R is known to three loops in QCD [50,63]. For our purpose, we need the first two orders which are given by 1 The renormalization equations (2.34) and (2.35) can be used to determine all the renormalization and rapidity scale dependent terms for the coefficient functions in perturbation theory. Throughout this paper, we organize perturbative expansions of various functions in powers of α s /(4π). For example Here and below we introduce two logarithms Up to O(α 2 s ), we then have Similarly for the TMDFF, the UV renormalization is given by Note that the procedure of renormalization and matching is easier to be done with the parton frame F N/q instead of the hadron frame D N/q (used in [67]). To extract the coefficient functions C iq , we calculate the bare TMDFFs with external parton states, and the partonic collinear FFs up to NNLO are given by where P T ij (z) are the time-like splitting kernels which will be presented in the Appendix. The renormalized C iq functions satisfy the evolution equations (2.43) At this point, it is worth noting that the product of the TMDPDF, TMDFF and the TMD soft function is independent on the rapidity scale ν as expected, namely where we have used It can also be shown that the µ-dependence of this product is cancelled by that of the hard function (which does not know about the rapidity divergences), such that the physical observables are independent of the renormalization scale. To see that we recall the RGEs of the hard and soft functions We note that the cancellation happens since γ B + γ H − γ S = 0 and By using the evolution equations, we can derive the scale-dependent part of C iq . Up to the NNLO we have iq (z) . (2.48) Note that here we have used the same symbol L Q as in Eq. (2.38) to denote a different meaning: which can be regarded as the crossing P 1+ → P 2− and x → 1/z.

Rapidity renormalization group and re-factorization
We now use the rapidity evolution equations of the TMDPDF, TMDFF and TMD soft function to derive the re-factorization formula (2.19). From the perturbative matching coefficients, it is evident that for the TMDPDF and TMDFF, the natural rapidity scale is ν ∼ xP 1+ ∼ P 2− /z ∼ Q, while for the TMD soft function the natural choice is ν ∼ b 0 /b T ∼ q T . In order to reconcile these different choices, we may use the rapidity RGE (2.45) for the TMD soft function to evolve it from ν where we have used γ R 0 = 0. The "genuine" quark TMDPDF and TMDFF which are free from rapidity divergences and are independent of the hard scale can then be defined as (2.52) These are essentially the functions appearing in the re-factorization formula (2.19).

Quark TMDPDF with the exponential regulator
In this section, we calculate the perturbative matching coefficients of the quark TMDPDF at NLO and NNLO using the exponential regulator. While these results are known to order 0 in the literature [38][39][40], we are able to obtain higher order terms in . The calculation with the exponential regulator is also much simpler and more systematic, which makes it possible to be extended to N 3 LO.

Quark TMDPDF at NLO
In this subsection, we briefly discuss the NLO results with the exponential regulator. While the calculation is straightforward, it illustrates the basic procedure and some interesting features of the regularization scheme.
We begin with the bare TMDPDFs before zero-bin subtraction. According to the definition in Eqs. (2.25), the TMDPDFs at NLO are given by the cut diagrams in Fig. 2. The result can be written as where we have changed to the notation that k T denotes the transverse components of k µ ⊥ , but with Euclidean signature such that The d-dimensional splitting amplitudes are given by Using the delta function for k + and the on-shell condition, we can write the exponential regulator as (3.4) At this stage we can already drop the second term proportional to τ (1−x) in the exponent, as it gives no contribution in the limit τ → 0. They might be relevant for subleading power corrections [49]. The first term involving τ /(1 − x) in the exponent provides the main service of regularizing the rapidity divergences. To see how that happens, we note that the rapidity divergence appears here as a singularity as x → 1 in the q → q splitting amplitude in Eq. (3.3). The exponential regulator provides a suppression in the x → 1 limit, and turns this singularity into a regularized distribution according to Applying the above equation to B (1),bare,unsub q/q , we find The k T integral can be easily performed with the help of the generating integral In particular, we have where ψ(x) = Γ (x)/Γ(x). Applying the above formulas to Eq. (3.6) then gives B (1),bare,unsub q/q exact in , and similarly for B (1),bare,unsub q/g . It is easy to expand the results to any order in . In particular, the O( 2 ) terms will be used for the NNLO calculations later, while the O( 4 ) terms are relevant for the calculations at N 3 LO.
We now need to subtract the collinear PDFs in Eq. (2.31) to obtain the coefficient functions I (1),bare,unsub qi , and then perform the zero-bin subtraction, where the NLO zerobin contribution is given by where L ν = ln(ν 2 /µ 2 ). After the subtraction, we find up to O( 0 ) where L Q = 2 ln(xp + /ν), and We now renormalize the UV divergences in the MS scheme with the renormalization factor and find Comparing the above form with Eq. (2.39), we can extract the (renormalization and rapidity) scale-independent part of the NLO coefficients (3.14) Remarkably, in the exponential regularization scheme, the scale independent coefficients are regular in the soft limit x → 1. As will be explicitly shown below, at NNLO there are 1/(1 − x) + distributions in the µ-independent part, but these terms are governed by the rapidity anomalous dimension and depend on the rapidity scale ν. In general this is true even at high orders in perturbation theory [40,68].

Quark TMDPDF at NNLO
We now turn to the NNLO calculations. As before, we begin with the bare TMD coefficient functions before zero-bin subtraction. At NNLO, diagrammatically there are two kinds of contributions. One is the interference of the LO amplitude with the diagrams containing one loop and one real emission, i.e., the so-called real-virtual (RV) contribution. The other is the square of the diagrams with two real emissions, i.e., the so-called double real contribution (RR). We will discuss these two contributions one-by-one in the following.

The real-virtual contribution
We adopt the light-cone gauge n·A = 0 where the relevant cut diagrams for the real-virtual contribution are depicted in Figure 3. Note that with both the exponential regulator and the analytic regulator used in [38,39], the loop integral does not need to be regularized. Therefore the treatments of the loop amplitude are rather similar. After performing the Dirac algebras and partial fractioning, there remain two classes of scalar integrals as shown in Fig. 4, which are given by where q = p − k and we make the +i prescription for all propagators implicit. The results of these integrals have already been given in [39] and we do not repeat them here. After the loop integration, the results are functions of x =n · q/n · p and k 2 T . The remaining integral over k 2 T can be carried out in the same way as the NLO calculation using Eq. (3.7).
q → q g → q Figure 3. Cut diagrams for the real-virtual contribution.

The double real contribution
At NNLO, the double real contribution is the most troublesome one to calculate. We will show that with the exponential regulator, we can apply many modern techniques for loop integrals. It is therefore possible to extend the calculation method to higher orders. We use QGRAF [69] to generate the relevant Feynman diagrams in the light-cone gauge, which are shown in Fig. 5. We then use FORM [70] to manipulate the squared amplitudes, and write them as integrals over the two cut momenta which we denote as k 1 and k 2 . We now need to apply the exponential regulator, and the integral measure then becomes It is useful to introduce an identity [71]    and rewrite the integral measure as Now the integration over k 1 and k 2 does not produce rapidity divergences and can be performed with usual techniques. Note that this fact holds also beyond NNLO where more than two cut momenta are present, due to the exponential form of the regulator. We can now use the delta function to integrate over k 2 , and rename k 1 as l. The double real contribution can then be written in the form k, l,n) , (3.19) where M qi is the squared amplitude. We will first integrate over l using the methods of reverse unitarity [72], integration-by-parts (IBP) [51] and differential equations [53][54][55]. The relevant topologies are given by (the square of) the diagrams shown in Figure 6. There are 4 topologies for the l-integrals, which are defined by where we use the subscript "cut" to label the cut propagators [72]. Integrals in each topology are further reduced to a set of Master Integrals (MIs) by IBP identities [51]. In this work, we use FIRE5 [73] and LiteRed [74] to perform the reduction. In total we have 6 MIs which can be chosen as n · l p · l , where the normalization factor We will use the method of differential equations to evaluate these MIs. For that purpose we have introduced the rescale factors w i to convert the MIs into a canonical basis [55]. They are given by where the dimensionless variables x and y are defined as The factors ω i can be easily obtained using an in-house code or the program package CANONICA [75] which implements the algorithm of [76]. Among all the MIs, F 1 , F 2 , F 4 and F 6 are easy to be evaluated in closed form where 2 F 1 is the hypergeometric function. F 3 and F 5 depend on both x and y and are more difficult to calculate. We can construct the differential equations of them with respect to y which are in the so-called canonical form [55]. Given their boundary conditions at y = 0 it is easy to solve the differential equations order-by-order in in terms of Goncharov multiple polylogarithms (GPLs). We have obtained the solutions up to weight 6, which will be sufficient for a future N 3 LO calculation.
The next step is then to perform the remaining integration over k in Eq. (3.19). The k + integral can be done using the delta function. And the k − integral can be changed to use the y variable through Note that we now have singularities at y → 0 or x → 1, which are both manifestations of rapidity divergences. These overlapping singularities often make high order perturbative calculations difficult due to the fact that the regularized integrand is often a complicated function of x and y. In our scheme, the regularization is provided by where the y → 0 and x → 1 limits are both exponentially suppressed. To perform the integration over y, we expand the above exponential regulator in terms of delta functions and plus-distributions according to Eq. (3.5) and The y integration can now be done using the package HyperInt [77] and the k T integration can again be evaluated with the help of Eq. (3.7). After the integration, the results can be expressed in terms of Harmonic PolyLogarithms (HPLs) [78] of the variable x. We use the program package HPL [79] to deal with these functions.

Final results at NNLO
Combining the real-virtual and double-real contributions, we obtain the bare un-subtracted NNLO TMDPDF. We then perform the zero-bin subtraction to remove double-counting between the collinear and soft sectors, and apply the usual α s renormalization and operator renormalization Z B q to remove the UV divergences. We have reproduced all the renormalization and rapidity scale dependent parts in Eq. (2.39), and the scale independent NNLO coefficients I (2) qi (x) are given by qq (x) , where q is a light quark flavor different from q, and we have used the shorthand notation H a 1 ,...,an ≡ H(a 1 , . . . , a n ; x) , (3.32) with H being HPLs. The p qi (x) functions are related to the DGLAP splitting kernels and are collected in the Appendix. We note that with the exponential regulator, the scaleindependent part of the TMDPDF does not involve δ(1 − x) terms, and the coefficients of 1/(1 − x) + is determined by the rapidity anomalous dimension 2γ R 1 given in Eq. (2.36). We have compared our results to those in the literature [38][39][40] and found full agreement. We have also obtained the bare NNLO TMDPDFs through to O( 2 ), which is required for a future N 3 LO calculation. Their expressions are quite lengthy and we choose to put them in an electronic file attached with the arXiv submission of this paper.

Quark TMDFF with the exponential regulator
We now turn to the quark TMDFF. Technically, it is very similar to the TMDPDF. The squared amplitudes are related via a crossing symmetry. The only subtlety is that one may perform the calculations in the hadron frame or in the parton frame. The two results should be related according to Eqs. (2.14) and (2.15). We have explicitly performed the two calculations and confirmed those relations.
In the hadron frame, one has, at a given order in perturbation theory where M iq is the squared amplitude for the q → i splitting, p is the momentum of the observed hadron in then direction, l i are loop momenta and k i are momenta of real emissions, k T,hf denotes the total transverse momentum of real emissions in the hadron frame. In the parton frame, one has instead p q k q → g Figure 7. Cut diagram for the bare q → g TMDFF at NLO. (4.5) The above integral is similar to the ones appearing in the calculation of TMDPDFs. The result is We now need to proceed with the matching procedure (2.29), where one should pay attention to the prefactor z 2−2 , which will produce logarithms of z when expanding in . We have Performing matching and renormalization as in Eq. (2.29), we then obtain Note that the scale-dependent part agrees with Eq. (2.48). We perform the calculation for the q → q fragmentation in a similar manner, where we need to use the exponential regulator for the rapidity divergences. The scale-independent coefficients at NLO are then given by where we use the shorthand notation H a 1 ,...,an ≡ H(a 1 , . . . , a n ; z) . (4.10) The NNLO calculations proceed in an analogous way, and we do not repeat the details here. The results are Again, we find that the scale independent parts do not contain δ(1 − z) terms, and the 1/(1 − z) + terms are determined by the rapidity anomalous dimension. We can convert our results to the convention of Ref. [67] and compare with the results in that work. We find that the results agree for the splitting processes q → q , q →q and q → g. However, for q → q, there is a small difference concerning a term C A C F π 4 δ(1 − z). In our framework, this term comes from the TMD soft function, which is universal for the TMDPDF and TMDFF. To address this discrepancy, we have performed several independent checks. The strongest check of our calculation comes from the calculation of the two-loop jet function of EEC in the back-to-back limit, which we shall explain in the next subsection.

Jet function for the EEC in the back-to-back limit
The EEC measures the energy correlation of two detectors in e + e − annihilation at an angle χ. The TMDFFs obtained in the last subsection can be used to calculate the jet function for the EEC in the back-to-back limit χ → π. It has been known for a long time that resummation of large logarithms for the EEC in this limit is closely related to q T resummation in the Drell-Yan process [17,80,81]. Recently, an all-order factorization formula in terms of operator matrix elements for the EEC in the back-to-back limit has been presented [35]. The factorization formula at leading power reads where z = (1 − cos χ)/2. In the back-to-back limit one has z → 1. The hard function and soft function in Eq. (4.12) is well known and can be found to two loops in Ref. [35]. The only missing ingredient for the resummation at N 3 LL accuracy is the two-loop jet function. In QCD, due to charge conjugation invariance, we have J q (b ⊥ , µ, ν) = Jq(b ⊥ , µ, ν). The quark jet function can be obtained from the second Mellin moments of the matching coefficients of quark TMDFFs [35] We expand the jet function in terms of α s as (4.14) Using the two-loop TMDFFs computed in this paper, the expansion coefficients are then given by (4.15) The µ and ν dependence of the jet function are in full agreement with the RGE and rapidity evolution equation [35]. The new result from this paper is the two-loop constant term which represents the last missing ingredient for N 3 LL resummation of EEC in the back-toback limit.
Using the two-loop jet function together with the two-loop hard and soft function, we obtain the full leading power prediction for the EEC in the back-to-back limit through two loops from the factorization formula in Eq. (4.12), including the δ(1 − z) terms, where we have set µ = Q and T F = 1/2 for simplicity, and The two-loop δ(1 − z) term is The two-loop plus distribution terms D n (z) are in full agreement with the analytical NLO calculation in Ref. [82], while the two-loop δ(1 − z) term is new. Eq. (4.19) has already been used in a previous publication to extract the δ(z) term of EEC using the energy conservation sum rule [83]. Two independent checks are made for Eq. (4.19). Firstly, EEC in the back-to-back limit obey the leading transcendental principle [84,85], which states that the maximal transcendental part of the QCD results are identical to the same quantity in N = 4 supersymmetric Yang-Mills (SYM) theory, up to trivial overall color factor [86]. In Eq. (4.19), the leading transcendental term is the ζ 4 terms. To compare with the same quantity in N = 4 SYM theory, we replace C F → C A in Eq. (4.19) and found the leading transcendental piece to be 40C 2 A ζ 4 , which is in full agreement with an independent calculation in [84]. Secondly, Besides the energy conservation sum rule, EEC in massless perturbation theory also obey a sum rule due to momentum conservation, which reads [84,87] where σ tot is the total hadronic cross section for e + e − including higher order QCD corrections. Using the analytical NLO formula of EEC for 0 < z < 1 from [82], and the end point contribution in Eq. (4.17), we explicit verify the sum rule in Eq. (4.20). The end point contributions in Eq. (4.17) are directly computed using the two-loop jet function in Eq. (4.15), which by itself reduces to moment of the TMDFFs. Therefore, the checks made for the end point contributions apply also to the TMDFFs computed in this paper, in particular to its δ(1 − z) terms.

Conclusion
In this work, we have revisited the calculation of perturbative quark TMDPDFs and TMDFFs at NNLO using a new regulator for rapidity divergences. We use the SIDIS process to set-up our calculation, while our results are universal and can be used for other processes as well. We show that the exponential regulator provides a consistent framework to carry out the calculation of the TMD soft functions, TMDPDFs and TMDFFs.
Compared to existing regulators in the literature, the exponential regulator has a couple of advantages. Firstly, the regulator can be implemented at the level of operator definitions for the TMD functions, where it manifests itself as a small shift of the spacetime coordinates. Secondly, the exponential regulator is applied to the total momentum of the extra emissions in the final state. Except for this last integration, the regulator does not change the structure of (cut)-propagators in the amplitudes. As a result, we can apply many modern techniques for loop integrals such as IBP identities and differential equations. This allows us to obtain the bare NNLO TMDPDFs and TMDFFs up to O( 2 ) in dimensional regularization, and can also be extended to a future N 3 LO calculation. Finally, the regulator can already be expanded in terms of delta-functions and plus-distributions at the integrand level, which makes the final round of integration easy to carry out.
Our results for the quark TMDPDFs up to O( 0 ) agree with the results in the literature, while our results for the quark TMDFFs have a small discrepancy with another calculation. To further check our results, we use the TMDFFs to calculate the NNLO jet function appearing in the factorization formula of EEC in the back-to-back-limit. This also serves as a new result of our paper, and is the last missing ingredient for an N 3 LL resummation. We have checked that our NNLO jet function produces the correct leading singular terms for the EEC in the back-to-back limit. This is a strong validation of our results for the TMDFFs.
Given the benefits provided by the exponential regulator, the calculation for the gluon TMDPDFs and TMDFFs (which was more difficult than the quark case using other regulators) can also be greatly simplified. We also believe that our method can be extended to the N 3 LO level. We leave these considerations to future publications.
The cusp anomalous dimension Γ cusp can be found in [88,89]. The hard and soft anomalous dimensions γ H and γ S can be extracted from the two-loop quark form factor [90,91], and can also be found in, e.g., Refs. [92,93]. Finally, the beam anomalous dimension γ B is related to γ S and γ H through γ B = γ S − γ H . And the renormalization factor for the quark TMDPDFs and TMDFFs up to O(α 2 s ) reads The QCD beta function is defined by with [94][95][96][97][98] A formula particularly useful for us is A.2 Space-like splitting functions The first order space-like splitting functions can be written as [99] P (0) qq (z) = 2C F p qq (z) + where we use the same convention as in Ref. [102].