Unpolarized Quark and Gluon TMD PDFs and FFs at N$^3$LO

In this paper we calculate analytically the perturbative matching coefficients for unpolarized quark and gluon Transverse-Momentum-Dependent (TMD) Parton Distribution Functions (PDFs) and Fragmentation Functions (FFs) through Next-to-Next-to-Next-to-Leading Order (N$^3$LO) in QCD. The N$^3$LO TMD PDFs are calculated by solving a system of differential equation of Feynman and phase space integrals. The TMD FFs are obtained by analytic continuation from space-like quantities to time-like quantities, taking into account the probability interpretation of TMD PDFs and FFs properly. The coefficient functions for TMD FFs exhibit double logarithmic enhancement at small momentum fraction $z$. We resum such logarithmic terms to the third order in the expansion of $\alpha_s$. Our results constitute important ingredients for precision determination of TMD PDFs and FFs in current and future experiments.


Introduction
Understanding the parton structures of hadron is one of the outstanding problem in Quantum ChromodyDynamics (QCD). TMD PDFs and FFs describe the distribution of parton transverse momentum inside a hadron in parton scattering or decay. Their knowledge is essential to our understanding of the confined motion of parton in nucleons [1][2][3]. Thanks to factorization and evolution [4][5][6][7][8][9][10][11][12][13][14][15], TMD PDFs and FFs also enter high precision theoretical prediction for a large variety of observables at high energy colliders. From pure theoretical point of view, TMD PDFs and FFs are also interesting since they represent light-cone correlation of quantum fields with intrinsic space-like and time-like origin, respectively. Their transparent definition in terms of light-cone correlator also allow higher-order perturbative calculation, from which one can uncover interesting analytic structure of the correlators.
Recently a first step towards N 3 LO TMD distributions have been achieved in [47] by calculating the unpolarized quark TMD PDFs, extending and significantly improving upon the methods in [26,27]. Subsequently, both unpolarized quark and gluon TMD PDFs are obtained in [48], based on an independent method from [49]. In this paper we continue our calculation in [47], and present the N 3 LO results for unpolarized quark TMD FFs, and unpolarized gluon TMD PDFs and FFs. We also present the results for unpolarized quark TMD PDFs from [47] for completeness.
Our method for calculating the unpolarized gluon TMD PDFs follows [47]. In order to obtain the results for TMD FFs, we adopt a strategy of analytic continuation proposed in [50]. The crucial observation is that although TMD PDFs or FFs are not themselves analytic function of momentum fraction, but the building blocks, space-like and time-like splitting amplitudes, are. At N 3 LO, there are four distinct contributions to TMD PDFs or FFs, namely the triple real part, the double real-virtual part, the double virtual-real part, and the virtual squared-real part. Ref. [50] shows that (a) the analytic continuation for the triple real part is trivial since it doesn't involve loop integrals; (b) the analytic continuation for the double real-virtual is also trivial at this order since the continuation of virtual loop only generate iπ terms which cancel in the sum with complex conjugate. (c) the analytic continuation of double virtual-real part and virtual squared-real part is non-trivial. But they are also simple to calculate since they only involve a one-particle phase space integral. It is suggested in [50] that one can calculate these two parts using the corresponding space-like or time-like splitting amplitudes, instead of trying to analytically continue them. Using this approach, Ref. [50] determines the complete time-like splitting functions at NNLO from the space-like counterpart, including the off-diagonal P T (2) qg , which cannot be completely fixed with previous methods [51][52][53]. With the complete NNLO time-like splitting functions, Ref. [50] also provides striking evidence for the existence of a generalized Gribov-Lipatov reciprocity relation between space-like and time-like QCD splitting functions in both non-singlet and singlet sector through NNLO. In this paper we use the idea of Ref. [50] to determine not just the splitting functions, but the full time-like TMD FFs through N 3 LO. Our results for TMD PDFs and FFs are written in terms of familiar harmonic polylogarithms [54] up to transcendental weight 5 and allow convenient analytical and numerical manipulation. We provide analytic expression for TMD PDFs and FFs in the threshold limit (x , z → 1) and high energy limit (x , z → 0). In the high energy limit, TMD FFs exhibit double logarithmic enhancement, instead of single logarithmic enhancement as the in the case of TMD PDFs. A recent discussion of the enhanced double logarithms can be found in [55]. We resum the large ln z terms in TMD FFs to the third order in the expansion of strong coupling, following the method of Vogt [56,57]. We note that a different method based on celestial BFKL equation has also been developed to resum the time-like small-z logarithms to NNLL in [55].
The structure of this paper is as follows. In Sec. 2 we give the necessary operator definition and renormalization for TMD PDFs and FFs. In Sec. 3 we present the N 3 LO results for unpolarized quark and gluon TMD PDFs and FFs. Since the expressions are lengthy, we present in the text only a numerical fit to these functions, valid to 0.1 percent in the range of 0 < x , z < 1. The full analytic expressions are provided as ancillary files in the arXiv submission of this paper. We also give the analytic expressions in the threshold limit in this section. In Sec. 4 we investigate the high energy limit of TMD PDFs and FFs, x , z → 0. We provide explicit analytic expressions, showing that TMD FFs are more singular than TMD PDFs from fixed-order point of view, namely double logs in contrast to single logs. We resum the small z logarithms for TMD FFs to the third order in the the expansion of α s . We conclude in Sec. 5.

Definition of quark and gluon TMD PDFs and FFs
In this section we give the necessary operator defintion for unpolarized quark and gluon TMD PDFs and FFs. We also give the renormalization counter terms and specify the zero-bin subtraction.
For sufficiently small b ⊥ , the TMD PDFs in Eq. (2.1) admit operator product expansion onto the usual collinear PDFs, where the summation is over all parton flavors i. The perturbative matching coefficients I bare qi (ξ, b ⊥ ) and I bare,µν 2) are independent of the actual hadron N . In practical calculations, one can replace the hadron N with a partonic state j. Furthermore, the usual bare partonic collinear PDFs are just φ bare up to power correction terms. For gluon coefficient functions, one can perform a further decomposition into two independent Lorentz structures in d = 4 − 2 dimension, where we have defined two scalar form factors, I bare gi , the unpolarized gluon coefficient functions, and I bare gi , the linearly-polarized gluon coefficient functions. They can be projected out using We focus on unpolarized TMD distributions for the current paper, and the results for the linearly-polarized gluon distribution are left for future work.

Operator definitions for TMD FFs
To specify the definition for gluon TMD FFs, it is necessary to specify a reference frame first. This is in contrast to PDFs, where the incoming hadron provides a canonical reference frame for transverse momentum. In the case of TMD FFs, two different frames can be defined, the hadron frame and the parton frame [12,26,27]. In the hadron frame, where the detected hadron has zero transverse momentum, an operator definition for gluon TMD FFs can be written down where P µ = (n · P )n µ /2 = P + n µ /2 is the momenta of the final state detected hadron. Again, we focus on unpolarized partonic TMD distributions, the projection to unpolarized gluon distributions as in Eq. (2.4) is understood from now on. In practical calculations, in particular for FF renormalization, it's also convenient to define the fragmentation functions in the parton frame, where the parton which initiates the fragmentation has zero transverse momentum. The parton frame TMDFFs are related to the hadron frame ones by [12,26] where we denote the bare partonic TMD FFs in the parton frame by F bare j/i . Our N 3 LO results for TMD FFs will be given in the hadron frame, by choosing the argument of the parton frame coefficient to be b ⊥ /z.

Renormalization counter terms and zero-bin subtraction
The TMD PDFs or TMD FFs, as well as their matching coefficients, contain both UV and rapidity divergences. We adopt dimensional regularization for the UV, and exponentional regularization [21] for the rapidity divergences. After coupling constant renormalization in MS scheme, the α s -renormalized matching coefficients still contains overlapping contributions between collinear and soft modes which is removed by a zero-bin subtraction [64]. After this, the remaining UV divergences are removed by multiplicative renormalization counter terms. After UV subtraction and zero-bin subtraction described above, the TMD PDFs or FFs still contain collinear divergence due to the tagged hadron in initial state or final state, and the remaining infrared poles are absorbed into the partonic dimensional regularized collinear PDFs 8) or the FFs ij is the (n + 1)-loop space-like splitting function [65,66], which are also known to the same accuracy in a massive environment [67,68]. P T (n) ij is the (n + 1)-loop time-like splitting function [50][51][52][53], and the symbol ⊗ is denoted as the convolution of two functions (2.10) The steps above can be summarized as the following collinear factorization formulas where B bare i/j (F bare j/i ) and S 0b (α s ) are the bare TMD PDFs (FFs) and bare zero-bin soft function, and Z B i (see in Sec. C) are the multiplicative operator renormalization constants for i = q, g. All the quantities are expressed in terms of the renormalized strong coupling α s . The zero-bin soft function is the same as TMD soft function, which is known to N 3 LO in Ref. [69]. Note also that the time-like TMD soft function is identical to the space-like one [70,71] up to this order, so we have a universal soft function up to O(α 3 s ).

N LO coefficients for unpolarized quark and gluon TMDs
In this section we give our results for coefficient functions I ij and C ij . We give only a numeric fit to these functions in the paper, but the full analytic expressions can be found in the ancillary files. For TMD FFs, we give the results for C ij with an argument b ⊥ /z, which after divided by z 2 are exactly the results in the hadron frame, see Eq. (2.7).

Renormalization group equations
The renormalized coefficient functions obey the following RG equations The rapidity evolution equations are [15,72] Expanding the perturbative coefficient functions in terms of α s /(4π), the solution to these evolution equations up to O(α 3 s ) reads, ji (x) , where in I ji we have used γ R 0 = 0 to simplify the expression and I (n) ji (z) are the scaleindependent coefficient functions. We have defined Similarly, the solution to the fragmentation coefficient functions are We stress again that due to the chosen argument, the expressions given above are for TMD FFs in the hadron frame. The anomalous dimensions appeared above are identical to those in space-like case and we have suppressed their dependence on the exact flavor. The logarithms appeared in the fragmentation coefficient functions are defined as which differ from thsoe in Eq. (3.5) only in L Q . Both space-like and time-like coefficient functions depend on the rapidity regulator being used. Rapidity-regulator-independent TMD PDFs and TMD FFs can be obtained by multiplying the coefficient functions with the squared root of the TMD soft functions S(b ⊥ , µ, ν) [26,27]

Numerical fits of the N3LO coefficients
The analytic expressions for the coefficient functions will be provided in the ancillary files along with the arXiv submission. In this section we will present their numerical fits. The coefficient functions develop end-point divergences both in the threshold and high energy limit. We first present here the results for leading threshold limit. The results for high energy limit will be discussed in next section. In the z → 1 limit, we have where γ R 1(2) are the two(three)-loop rapidity anomalous dimensions [69,73]. The relation between threshold limit and rapidity anomalous dimension has been anticipated in [25,74,75]. The explicit expressions up to three-loop read [47] We also found that threshold limit exhibits Casimir scaling up to three loops (n = 1, 2, 3),

Numerical fit for TMD PDFs
The analytic expressions of three-loop coefficient functions contain harmonic polylogarithms up to transcendental weight 5. To facilitate straightforward numerical implementation, we provide a numerical fitting to all the coefficient functions. Following Ref. [76], we use the following elementary functions to fit the results, For two loop and three loop coefficient functions, we fit the exact results in the region 10 −6 < x < 1 − 10 −6 (Numerical evaluation of HPLs are made with the HPL package [77]), and we have set the color factor to numerical values in QCD, i.e.
In more detail, we subtract the x → 0 and x → 1 limits up to next-to-next-to-leading power (x 1 and (1 − x) 1 ). Then we fit the remaining terms in the region 10 −6 < x < 1 − 10 −6 . Combining the two parts, the fitted results can achieve an accuracy better than 10 −3 for 0 < x < 1. We show below the numerical fitting with six significant digits. The full numerical fitting is attached as ancillary files with the arXiv submission. The one loop scale independent coefficient functions are given by The two-loop scale independent coefficient function are given by To present the three-loop scale independent coefficient functions, we first perform the following decompositions, qq (x), I where The numerical fitting of different color structures are given by

Small x expansion of unpolarized TMD coefficients and resummation for TMD FFs
In this section we give the results for TMD PDFs and FFs expanded in high energy limit, namely x, z → 0. A striking difference between space-like TMD PDFs and time-like TMD FFs is that there is only single logarithmic enhancement in each order of perturbative expansion in TMD PDFs, while in TMD FFs it becomes double logarithmic enhancement. We also resum the small-z logarithms in TMD FFs to Next-to-Next-to-Leading Logarithmic (NNLL) accuracy in this section.

Small-x expansion of unpolarized TMD PDFs
Using the analytic expression we obtained, it is straightforward to obtain the small-x expansion. At leading power in the expansion, the results read We note that a LL prediction for the small-x expansion has been given in [78]. After fixing a typo in that paper, we find full agreement with its LL prediction for both quark and gluon TMD PDFs at small x. 1 It would be very interesting to extend the formalism of [78] beyond LL and compare with the data presented here.

Small-z expansion of unpolarized TMD FFs
To facilitate small-z resummation for TMD FFs, we shall consider the coefficient functions in flavor singlet sector below, since non-singlet TMD FFs are at most logarithmic divergent, but not power divergent in the z → 0 limit. The flavor singlet (denoted by a superscript s ) coefficient functions can be written as a matrix, where and C ij (z) are scaleless coefficient functions as appeared in the solutions of RG equation (3.6).
In contrast to TMD PDFs, which contribute a single logarithm at each perturbative We note that while both C s qi and C s gi has double logarithmic expansion in the small-z limit, the power of leading logarithmic terms of C s qi is lower by 1 than the corresponding leading logarithmic terms of C s gi .

Resummation of small-x logarithms for unpolarized TMD FFs
In this subsection, we shall derive the all-order resummation at NNLL accuracy (resummation of the highest three logarithms) in Eq. (4.7), following an idea proposed in [56].
To this end, we begin with the unrenormalized version of collinear factorization formula in Eq. (2.11) for singlet TMD FFs (see (4.6) for the definition of singlet combination), where the convolution is in z. Note that F s i/j (z, ) is a quantity to which the usual strong coupling renormalization, zero-bin subtraction and operator renormalization have been performed. However renormalization of collinear FFs has not been performed, which is why we still keep the dependence in (4.12). From now on we concentrate on the scaleindependent part of the coefficient functions by setting all the scale logarithms to zero. We can do this because the scale logarithms depends either on anomalous dimension, which has no z dependence, or on time-like splitting functions, whose small-z behavior is known to NNLL accuracy [57]. It proves convenient to work in Mellin-N space also, which is defined as 13) where N = N − 1. Small-z logarithms becomes poles inN under Mellin transformation, (4.14) In Mellin space the unrenormalized collinear factorization formula in Eq. (4.12) becomes , where γ T (N ) is the time-like singlet splitting function in Mellin space. Its complete NNLO results can be found in [50], see also [51][52][53].
The crucial observation of [56] is that unrenormalized collinear functions have specific singular behavior in the limit of small-z in dimensional regularization. In the case of TMD FFs, we can write down an general ansatz at small z, where c (1,l,n) gi is the leading term in the expansion and small-z expansion, whose knowledge correspond to LL resummation as labeled in (4.18), and similarly for other terms. Precisely, for C s gi (z) the LL series correspond to α n s ln 2n−1 z terms, while NLL correspond to α n s ln 2n−2 z, and NNLL to α n s ln 2n−3 z. For C s qi (z) the corresponding power of ln z is lower by 1. We have verified this general ansatz through explicit N 3 LO calculation from its operator definition in (2.6) and (2.7). In Mellin space the corresponding ansatz reads Equations (4.18) or (4.19) provides a way to resum all the large logarithms of z. Specifically, if one knows all the c 1,l,n gi for all l and n, then one can do LL resummation for F s g/i , and similarly for NLL and NNLL resummation. In this paper instead of working out the constants for all n, we provide results for n up to 15, which is sufficient for phenomenological purpose.
We note that the neglected terms in (4.19) are higher orders in and in N . To facilitate easy extraction of the constants, it is convenient to define a "small-z" weight: 20) such that c gives us n conditions and is enough to determine the n unknowns at LL accuracy. To acheive NNLL accuracy, one only needs two more power of expansion, that is we need to know F s(n) g/i (N , ) to order −n+2 . Similar analysis shows that we also need to know F s(n) q/i (N , ) to order −n+2 to achive NNLL accuracy.
Having this important information in hand, let us contrentrate on the right hand side of Eq. (4.15), which generates the necessary poles. The dimensionally regularized partonic FFs d s (N , ) in MS scheme is purely divergent in and can be determined easily by solving Eq. (4.17) order by order in α s , .

(4.22)
where we have traded d ln µ 2 for da s with the help of the (4 − 2 )-dimension beta function, Assuming d s (N , ) = 1 + ∞ k=1 a k s d s k and γ T (N ) = ∞ k=0 a s k+1 γ T k , the solutions can be worked out order by order. For example, at first four orders we have As explained above, to NNLL accuracy, we need to determine the coefficients of the pole terms in F s(n) g(q)/i to order −n+2 . Since the matching coefficients C s g(q)i (N , ) in (4.15) must be finite in the limit of → 0, we need to determine the coefficients of the pole terms in d s n also to order −n+2 . Also recall that d s n are function of N and we are only interested in the small N limit. For NNLL resummation we only need to keep the terms in d s n with "small-z" weight up to 2 − 2n, see Eq. (4.20). We also note that the lowest weight term in γ T n is 1 − 2n, corresponding to the leading 1/N 2n−1 pole. Using these information it can be shown that the required inputs for d s n to achieve NNLL accuracy are β 0 and γ T 0,1,2 . One can check that this is indeed the case frm Eq. (4.24), where explicit examples for n ≤ 4 are shown.
Having known the general form of counter terms d(N , ) to any order in α s , and the fact that the coefficient function C s g(q)i is finite in the → 0 limit, we can solve Eq. (4.15) recursively to get the coefficient for the required pole terms in F s g(q)/i , order by order in α s . In summary the input data we need to achive NNLL accuracy are (4.25) Note that although we have the explicit results for F s (3) g(q)/i from direct calculation, they are not needed for predicting the small-z logarithms at NNLL accuracy. Rather they can be used as a check for our resummation.
Following the approach outlined above, we have determined the coefficients of poles in F s(n) g(q)/i up to −n+2 for n ≤ 15. From these coefficients we solve for the LL to NNLL constants defined in (4.18). Expanding the results in gives us the resummed NNLL series truncated at order α 15 s , which should be sufficient for phenomenology. We have checked that a truncation of perturbative series at α 15 s leads to a less that 1% relative uncertainty. The analytic expressions for the truncated resummed perturbative series can be found in the ancillary files of this paper.
The series can be resummed analytically, leading to an closed form expression for LL results, It's interesting to note that Eq. (4.27) coincides with that of the transverse coefficient functions for semi-inclusive e + e − annihilation [56,79]. In Fig. 1 and 2 we plot the fixed-order coefficient functions and with different orders of small-z resummation. We use N f = 5 throughout the calculations. We note that even at N 3 LO, the effects of resummation is important for z < 10 −2 .

Conclusion
In summary we presented calculations for the unpolarized quark and gluon TMD PDFs and FFs at N 3 LO in QCD. The unpolarized quark TMD PDFs at N 3 LO have already been reported in Ref. [47]. The rest of the calculations are new. Unpolarized quark and gluon TMD PDFs have also been calculated in [48] recently using an independent method, whose results are in full agreement with ours 2 . Our calculations for TMD PDFs are based on the method proposed in [47], which in turns is based on decomposition of light-cone correlators into phase space integration of collinear splitting amplitudes of different multiplicities. The advantage of this decomposition is that it allows better understanding of the analytic continuation property of TMD PDFs and FFs [50]. Using the analytic continuation prescription of [50], we successfully obtained the N 3 LO TMD FFs from corresponding TMD PDFs, without the need to compute everything from scratch. Our results open the avenue for precision phenomenology of TMD physics at N 3 LO in perturbative QCD.
We also provide threshold and high energy asymptotics of TMD PDFs and FFs through N 3 LO by expanding the corresponding analytic expressions. The high energy (small-z) limit of TMD FFs features double logarithmic enhancement, in contrast to the single logarithmic enhancement of TMD PDFs. We resum the small-z logarithms through NNLL accuracy using a method proposed in [56]. The resummation leads to better behaved perturbative convergence for z < 10 −2 .
Our method of calculation is general and is not limited to unpolarized distribution. Works towards to polarized TMD distributions at N 3 LO are in progress. Note added: A few days after this paper appeared, an independent calculation for unpolarized quark and gluon TMD FFs at N 3 LO was submitted to arXiv [80]. During private communication, the Authors of [80] uncovered a minor error in one of their routine for analytic continuation. After fixing it, they found full agreement with our results.

A QCD Beta Function
The QCD beta function is defined as with [81]  where the coefficients for quark are given by Since cusp and soft and rapidity anomalous dimensions exhibit Casimir scaling, the corresponding anomalous dimensions for gluon could be obtained by multiplying in above with C A /C F . The beam anomalous dimensions do not exhibits Casimir scaling, thus should be list separately. The beam anomalous dimensions for quark are 3) The beam anomalous dimensions for gluon are The cusp anomalous dimension Γ cusp can be found in [65]. The beam anomalous dimension γ B is related to the soft anomalous dimension γ S [82] and the hard anomalous dimensions γ H [83][84][85] by renormalization group invariance condition γ B = γ S − γ H . The rapidity anomalous dimension γ R can be found in [69,73]. Note that the normalization here differ from those in [69] by a factor of 1/2.

C Renormalization Constants
The following constants are needed for the renormalization of zero-bin subtracted [64] TMD PDFs through N 3 LO, see e.g. Ref. [26,27]. The first three-order corrections to Z B and Z S are Keep in mind that the anomalous dimensions appeared above depends on the flavor, they should be replaced by the corresponding values in Sec B. We also remind the reader that the renormalization constants are formally identical for TMD PDFs and TMD FFs, the logarithms appeared above should be replaced by their corresponding values in each case, and we have with b 0 = 2 e −γ E for both TMD PDFs and TMD FFs. For TMD PDFs, while for TMD FFs, L Q = 2 ln P + z ν . (C.4)