Unpolarized Transverse Momentum Dependent Parton Distribution and Fragmentation Functions at next-to-next-to-leading order

The transverse momentum dependent parton distribution/fragmentation functions (TMDs) are essential in the factorization of a number of processes like Drell-Yan scattering, vector boson production, semi-inclusive deep inelastic scattering, etc. We provide a comprehensive study of unpolarized TMDs at next-to-next-to-leading order, which includes an explicit calculation of these TMDs and an extraction of their matching coefficients onto their integrated analogues, for all flavor combinations. The obtained matching coefficients are important for any kind of phenomenology involving TMDs. In the present study each individual TMD is calculated without any reference to a specific process. We recover the known results for parton distribution functions and provide new results for the fragmentation functions. The results for the gluon transverse momentum dependent fragmentation functions are presented for the first time at one and two loops. We also discuss the structure of singularities of TMD operators and TMD matrix elements, crossing relations between TMD parton distribution functions and TMD fragmentation functions, and renormalization group equations. In addition, we consider the behavior of the matching coefficients at threshold and make a conjecture on their structure to all orders in perturbation theory.


Introduction
The transverse momentum dependent parton distribution and fragmentation functions (TMDs) play a central role in our understanding of QCD dynamics in multi-differential cross sections and spin physics. Recently, factorization theorems for Drell-Yan, Vector Boson/Higgs Production, Semi-Inclusive Deep Inelastic Scattering (SIDIS) and e + e − → 2 hadrons processes, both for spin-dependent and unpolarized hadrons, has been reformulated in terms of individually welldefined TMDs [1][2][3][4], updating the pioneering works of Collins and Soper [5,6]. All these processes are fundamental for current high energy colliders, like the LHC, KEK, SLAC, JLab or RHIC, and future planned facilities, like the EIC, AFTER@LHC, the LHeC or the ILC.
So, in this work we provide a comprehensive study of TMDs at NNLO based on a direct calculation of TMD matrix elements at NNLO. In particular, our results provide an indirect confirmation of the TMD factorization theorem and the related structure of rapidity divergences. In fact we explicitly confirm that the cancellation of rapidity divergences is realized within one single TMD, and not necessarily in the product of two TMDs, which is important when studying the non-perturbative parts of these quantities.
The TMD factorization theorem at higher orders in perturbative QCD is not trivial. In fact, in the calculation one has to deal with several types of divergences (ultra-violet (UV), rapidity and infra-red (IR)), which have to be regularized and disentangled properly. The TMD factorization theorem offers a strategy to remove the rapidity divergences in order to achieve a well-defined TMDs. Recently our group has provided a direct calculation of an individual TMD at NNLO, namely the unpolarized quark transverse momentum dependent fragmentation function (TMDFF) [21], and a complete study of the structure of rapidity divergences at the same order in the soft function [22] (see also [20,23]). In this work we complete the calculation of the unpolarized quark TMDFF at NNLO, showing also the details of it and including new results for the gluon TMD fragmentation function.
On the other hand, we also calculate the unpolarized quark and gluon TMD parton distribution functions (TMDPDFs). Some properties of the TMDPDFs, like their matching onto integrated parton distribution functions (PDFs) can be found in previous works [15][16][17][18][19], where they were obtained by decomposing the product of two TMDs and did not use the fact that each TMD is per se calculable. In other words, these calculations did not fully exploit the results of the TMD factorization theorem of [1][2][3][4]. We find a complete agreement between our calculation and the results of [15][16][17][18][19], once the proper combination of collinear and soft matrix elements is considered, which represents a strong check and demonstration of the regulator-independence of the matching of each TMD onto its corresponding integrated counterpart.
In this work, we slightly go away from the standard formulation of TMD factorization, which is derived for different processes, towards a universal process-independent definition of TMDs. With this aim we introduce the process-independent TMD operators, in analogy with the parton string operators for integrated functions. So, the TMDs are the hadron matrix elements of these TMD operators. Such a reformulation suggests a new more general look on TMDs, and reveals common points between various approaches.
The formulation of a universal TMD operator is possible due to the process independence of the soft factor, that has been discussed for long time [1][2][3][4]24] and at NNLO has been explicitly demonstrated in [22]. Unlike the usual composite and light-cone operators, the TMD operator is more divergent. It contains rapidity divergences together with UV divergences. Therefore, the proper definition of TMD operator should include not only UV renormalization constant, but also a mechanism for removing the rapidity divergences. It is known that the structure of rapidity divergences within the TMDs is similar to the structure of UV divergences. The similarity is seen already for the soft factor, the logarithm of which is necessarily linear in rapidity divergences, and hence satisfies a renormalization group equation (RGE) with respect to rapidity scale. Thus, the rapidity divergences can be removed analogously to UV divergences by the "rapidity renormalization" factor, which in fact naturally appears in any formulation of TMD factorization theorem (see e.g. [13,25]).
The matrix elements of TMD operators are free from operator divergences and can have only standard IR divergences related to the external states, as checked here at NNLO for all possible unpolarized TMD operators. Using the expressions for partonic TMD matrix elements we can extract the matching coefficients of the TMDs onto their corresponding integrated functions (parton distribution functions (PDFs) or fragmentation functions (FFs)). These coefficients are free from any type of divergences and have a direct impact on phenomenological analyses. We provide the two-loop coefficients for all unpolarized TMDs, and this is the main practical result of this paper.
The regularization of rapidity divergences used in this paper is the same as in [21,22]. We use the so-called δ-regulator, which in practice is just a shift of the residue of the Wilson lines by an amount δ, to be removed at the end of the calculation. Several technical details are necessary for the proper implementation of this regulator at higher orders, which are discussed in the text. For the rest of the divergences we use the standard dimensional regularization. This particular choice of regulators simplifies significantly the calculation. For example, one can avoid a direct calculation of pure virtual contributions, which reduces the number of diagrams to be computed. The soft function presented in [22] is a key element for the NNLO calculation of all (polarized and unpolarized) TMDs. The present calculation is a confirmation of the universality of the soft function as it enters at the same footing in the calculation of both TMDPDFs and TMDFFs.
We report on the structure of the matching coefficients of the TMDs onto their corresponding integrated functions, consistently with their RGEs. Also we consider a series of technical topics which we think are interesting for the expert reader. So, we discuss the realization of the Gribov-Lipatov correspondence between TMDPDFs and TMDFFs, due to crossing symmetry. Although the computations for TMDPDFs and TMDFFs have been done independently, we have used the crossing symmetry at intermediate steps of the calculation as a check of our results. Using the obtained results we are also able to formulate a conjecture about the behavior of the coefficient at threshold at all order in perturbation theory. This result can be important to study the perturbative series for the coefficients at N 3 LO.
The paper is organized as follows. Section 2 provides the definitions of universal TMD operators, with their renormalization. We also discuss the basic structure of the small-b T operator product expansion (OPE), and its relation to the matching procedure. In Sec. 3 we discuss the regularization method, the general structure of rapidity divergences as well as the details of the calculation of TMDs. In Sec. 4, the renormalization group equations for TMDs are introduced. In Sec. 5 we present in detail the NLO computation of the TMDs and their matching coefficients. This section serves mainly as a pedagogical demonstration of all the steps of the computation and the internal structure of TMDs. The technics used for the NNLO calculation are presented in Sec. 6, while in Sec. 7 we collect all the expressions for TMD matching coefficients for both TMDPDFs and TMDFFs up to NNLO. We study the coefficient for large values of the Bjorken variables x, z in Sec. 8. Finally, in Sec. 9 we conclude. The set of appendices includes several necessary definitions, some intermediate expressions and side results which were used in the paper.

Definitions of TMD operators
The factorization theorems for transverse-momentum-dependent cross sections are usually formulated in terms of TMDs. In this work however we follow a different strategy, namely, we focus our attention on the TMD operators. Such a consideration allows us to have a homogeneous notation and reveals the similarities between the distribution and fragmentation functions. It also allows us to formulate statements in a process-independent way.
We define the bare (unrenormalized and rapidity singular) quark, anti-quark and gluon unpolarized TMDPDF operators as follows: The repeated color indices a (a = 1, . . . , N c for quarks and a = 1, . . . , N 2 c − 1 for gluons) are summed up. The representations of the color SU(3) generators inside the Wilson lines are the same as the representation of the corresponding partons. The Wilson linesW T n (x) are rooted at the coordinate x and continue to the light-cone infinity along the vector n, where it is connected by a transverse link to the transverse infinity (that is indicated by the superscript T ). The precise definition of the Wilson lines is given in Sec. 3.
The hadronic matrix elements of the operators defined in Eq. (2.1) provide the unsubtracted TMDPDFs, as they are defined within the TMD factorization theorems [1,2,2,3]: where N is a nucleon/hadron. Here the variable x represents the momentum fraction carried by a parton from the nucleon (it also explains the TMD labelling rule f ← N ). One can see at the operator level that TMDs are like the integrated parton densities, with the only difference that parton fields are additionally separated by the space-like distance b T . The gauge connection between the parton fields follows the path uniquely dictated by the relevant factorization theorems for the given physical processes.
The definition of the operators for the fragmentation functions follows a similar pattern, with the main difference that they should be calculated on final rather than initial states. Formally, one where δ/δJ is to be understood as the state generated by the variation of the action with respect to the source J, which couples to external hadron fields. Then the unsubtracted TMDFFs are hadronic matrix elements of these operators: where again N is a nucleon/hadron. The variable z represents the momentum fraction of the parton carried into the hadron (it also explains the TMD labelling rule f → N ). The definitions of the quark TMDFFs coincide with the one coming from TMD factorization [1,2,2,3]. To our knowledge, the gluon TMDFFs were first considered in [26]. However here we find more convenient to define the normalization factor of the gluon TMDFF in analogy to the normalization of the integrated FFs.Here the normalization of TMDFFs counts the number of the physical states of a given flavor (being 2(1 − ) the number of physical gluon polarizations in d = 4 − 2 dimension). Such normalization allows the crossing relations discussed below to be fulfilled. Summarizing the expressions Eq. (2.1-2.4), the bare TMDs 1 are the hadronic matrix elements of the corresponding bare TMD operator: where the Hermitian conjugation of the states for TMDFFs indicates that these are final states to be placed inside the operator.
Unlike the usual composite or light-like operators which contain only UV divergences, the TMD operators in addition suffer from rapidity divergences. The UV divergences in the TMDs are removed by the usual renormalization factors. In order to cancel rapidity divergences one should consider both the zero-bin subtractions and the soft function. According to Soft Collinear Effective Theory (SCET) terminology the "zero-bin" represents the soft overlap contribution, that should be removed from the collinear matrix element in order to avoid double counting of soft singularities [27]. The combination of the zero-bin subtraction with the soft function has a very particular form, which is dictated by the factorization theorem, and should be included in the definition of the TMD operators as a single "rapidity renormalization factor". Therefore, we introduce the "rapidity renormalization factor" R, which completes our definition of the renormalized TMD operator. We have where Z q and Z g are the UV renormalization constants for TMD operators. The scales µ and ζ are the scales of UV and rapidity subtractions respectively. While the UV renormalization factors depend on the UV regularization method and the regularization scale µ, the "rapidity renormalization factors" also depend on the rapidity regularization method and the rapidity scale ζ. Moreover, given that the soft function is process independent (as argued with general arguments in [1][2][3][4]24] and explicitly checked at NNLO in [22]), the "rapidity renormalization factors" are also process independent. The details of the definition of the factor R are discussed in Sec. 3.
It is crucial to observe that both UV and rapidity renormalization factors, Z and R respectively, are the same for TMDPDF and TMDFF operators. That is not accidental, but the consequence of the fact that both TMDPDF and TMDFF operators have the same local structure (which makes equal the factors Z) and the same geometry of Wilson lines (which makes equal the factors R). This significantly simplifies the consideration of the operators and makes the whole approach more universal. Moreover, from the equality of renormalization factors follows that the evolution equations for TMDPDF and TMDFF are the same. The appropriate anomalous dimensions can be extracted from R and Z, see details in Sec. 4.
The renormalization of rapidity divergences needs some caution, since with some regulators the rapidity divergences can be confused with UV poles. On top of this, the particular form of zero-bin subtractions included in the factor R, is regulator dependent. Thus, in order to avoid any possible confusions we fix the exact order on how to deal with these singular factors: we first remove all rapidity divergences and perform the zero-bin subtraction, and afterwards multiply by Z's. Such an order implies that the factor R contains not only rapidity divergences, but also explicit UV divergences which are also taken into account in the factor Z. These two strategies lead to different intermediate expressions, while the final (UV and rapidity divergences-free) expressions are necessarily the same. Now, given all previous considerations, we define the individual TMDs as Such a definition implies the following relation between bare and renormalized TMDs: that follows from the TMD factorization theorem. Finally we comment that, had we chosen a different order of recombination of singularities, then we would find separate UV-renormalization factors for the soft factor and the collinear matrix element, which in turn would depend on the parameters of the rapidity regularization. Such a strategy has been recently used in [20], following the "Rapidity Renormalization Group" introduced in [13,25], which is built in order to cancel the rapidity divergences through renormalization factors from the beam and soft factors independently. In our case, however, the rapidity renormalization factor itself is constructed with the soft function. Thus, although at the end of the day the same rapidity logarithms are resummed, the definitions of TMDPDFs and thus the underlying logic is different. In Ref. [12,18,19] for TMDPDFs the soft function is hidden in the product of two TMDs.

The operator product expansion (OPE) at small b T
The TMDs, as a non-perturbative objects, are a highly involved functions. Any information on their behavior is important for phenomenological applications. Apart from evolution equations, QCD perturbation theory can also supply the small-b T asymptotic behavior of the TMDs, and give their matching coefficient onto their integrated collinear counterparts (see e.g. [1,2]). Such a matching is interesting because one expects it to provide a good description of x/z-dependence of TMDs in the whole region of b T , and together with a suitable ansatz for the non-perturbative contribution at large-b T , provides a reasonable phenomenological model. Also the matching represents a strong check of the theory and in this article we explicitly work it out at NNLO, both for quark/gluon distributions and fragmentation TMDs.
At the operator level the small-b T matching is a statement on the leading term of the small-b T Operator Product Expansion (OPE). The small-b T OPE is a formal operator relation, that relates operators with both light-like and space-like field separation to operators with only light-like field separation. It reads where the dimension of the transverse derivatives, ← → ∂ T = ← → ∂ /∂b T (these derivatives acts at light-like infinity, therefore the gauge field can be omitted in non-singular gauges), is compensated by some scale B T . The matching coefficients behave like where f is some function. The unknown scale B T represents some characteristic transverse size interaction inside the hadron. So for b T B T it is in practice reasonable to consider only the leading term of the OPE in Eq. (2.13), which gives the matching of the TMDs onto the integrated functions. The consideration of higher order terms is an interesting and a completely unexplored part of the TMD approach, which we do not further consider in this work. Note that the OPE onto the operators of the form in Eq. (2.14) may not be the most efficient, see discussion and alternative small-b T OPE based on Laguerre polynomials in [7,28].
For the TMDPDFs the leading order small-b T operator (i.e. the operator for the integrated PDF) is just a TMDPDF operator Eq.
while for FF kinematics one has an extra normalization factor Notice that in the equations above we have dropped a subindex 0. In this way the leading terms of the OPEs at small b T read where the symbol ⊗ is the Mellin convolution in variable x or z , and f, f enumerate the various flavors of partons. The running on the scales µ, µ b and ζ is independent of the regularization scheme and it is dictated by the renormalization group equations. Taking the hadron matrix elements of the operators we obtain the small-b T matching between the TMDs and their corresponding integrated functions, The integrated functions (PDFs and FFs) depend only on the Bjorken variables (x for PDFs and z for FFs) and renormalization scale µ, while all the dependence on the transverse coordinate b T and rapidity scale is contained in the matching coefficient and can be calculated perturbatively. The definition of the integrated PDFs are and similarly, for integrated FFs In practice, in order to calculate the matching coefficients we calculate both sides of Eq. (2.18) on some particular states and solve the system for matching coefficients. Since we are interested only in the leading term of the OPE, i.e. the term without transverse derivatives, it is enough to consider single parton matrix elements, with p 2 = 0. A study of the matching coefficients for higher-derivative operators can be found in [7,28].
3 Regularization and structure of the divergences

Explicit form of rapidity renormalization factor
In the previous Section we have defined the factor R f in a rather abstract way, as a kind of "rapidity renormalization factor". In fact, its explicit form is dictated by the TMD factorization theorem and reads where S(b T ) is the soft function and Zb denotes the zero-bin contribution, or in other words the soft overlap of the collinear and soft sectors which appear in the factorization theorem [1][2][3][4]27].
We now elaborate on this definition. The soft function is defined as a vacuum expectation value of a certain configuration of Wilson lines, which depends on the process under investigation. For example, for SIDIS it reads The Wilson lines are defined as usual The transverse gauge links T n are essential for singular gauges, like the light-cone gauge n · A = 0 (orn · A = 0), see details in Refs. [29][30][31]. In covariant gauges the transverse links are needed only to preserve the gauge invariance, but in practice do not add any contribution. Note that collinear Wilson lines W T n (x) used in TMD operators Eq. (2.1-2.4) are defined in the same way as soft Wilson lines S T n (x). However, we distinguish them since they behave differently under regularization. The zero-bin (or overlap region) subtraction is a subtle issue. In fact, the explicit definition of this subtraction significantly depends on the rapidity regularization used (see e.g. discussion in [3]). Thus, for a given regularization scheme it might be even impossible to define the zero-bin as a well-formed matrix element. Nonetheless, for any regularization scheme it has a very particular calculable expression. With a conveniently chosen rapidity regularization, the zero-bin subtractions are related to a particular combination of the soft factors. Using the modified δ-regularization, which is discussed in detail in the next Section, the zero-bin subtraction is literally equal to the SF: Zb = S(b T ). We should mention that this is not a trivial statement, and in fact, the modified δ-regularization scheme has been adapted such that this relation holds. In particular, it implies a different regularized form for collinear Wilson lines W n(n) (x) and for soft Wilson lines S n(n) (x).
So, concluding, in the modified δ-regularization that is used in this work, the expression for the rapidity renormalization factor is The relation Eq. (3.4) was first checked explicitly at NNLO in [21,22], and also confirmed for various kinematics in this work. We notice that due to the process independence of soft function [1][2][3][4]24], the factor R f is also process independent. The origin of rapidity scale ζ is explained in the next section. Let us also make a connection to the formulation of TMDs by Collins in [1]. In the JCC approach the rapidity divergences are handled by tilting the Wilson lines off-the-light-cone. Then the contribution of the overlapping regions and soft factors can be recombined into individual TMDs by the proper combination of different SFs with a partially removed regulator. This combination gives the factor R f in our notation, S(y c , yn)S(y n , yn) . (3.5) The following logical steps remain the same as with the δ-regulator.

Modified δ-regularization scheme
The original δ-regularization proposed in [2] consists in a simple infinitesimal shift of the i0prescriptions in eikonal propagators. However, such a rude approach appears to be not sufficient at NNLO for several reasons, e.g., the fact that at this order the zero-bin and the soft function are not equal. Therefore, in [21,22] the δ-regularization scheme was conveniently modified to overcome this issue. The modified δ-regularization is implemented at the operator level, and constructed in such a way that it explicitly preserves the non-Abelian exponentiation and the equality of zero-bin and the SF. The implementation of the regularization at the operator level grants many benefits in the analysis of the all-order structure of rapidity divergences, and allows to prove such statements as the linearity of the logarithm of the soft function in lnδ. The detailed discussion on the properties of the modified δ-regularization can be found in [22]. Here we limit ourselves to present the definitions and make the essential comments.
The modified δ-regularization scheme has to be defined at the operator level, and consists in modifying the definition of Wilson lines. So the soft Wilson lines entering the soft function in Eq. (3.2) are changed according tõ where δ ± → +0. At the level of Feynman diagrams in momentum space, the modified expressions for the eikonal propagators are written as (e.g. absorption of gluons by a Wilson line [∞ + , 0]) where the gluons are ordered from infinity to zero (i.e. k n is the gluon closest to zero). As a consequence of the rescaling invariance of the Wilson lines (that is now explicitly broken by the parameters δ ± ), the expressions for diagrams in the soft function depend on a single variable 2δ + δ − /(n ·n) = δ + δ − . The ordering of poles in the eikonal propagators, Eq. (3.7), is crucial for the perturbative exponentiation with usual properties, such as non-abelian exponentiation theorem for color-factors [32,33] or logarithmical counting [34]. As a matter of fact, within the modified δregularization, only diagrams with non-Abelian color prefactor (web diagrams) arise in the exponent. Therefore, the expression for the soft function can be written in the form where a s = g 2 /(4π) 2 is the strong coupling and C F is the Casimir of the fundamental representation of gauge group . The collinear Wilson lines appearing in the definition of the operators, Eq. (2.1)-(2.3), should be regularized in a slightly different way in order to accomplish Eq. (3.4). This is achieved by rescaling the δ-regulator with the Bjorken variables as in the case of the Wilson lines appearing in TMDPDFs, Eq. (2.1), and as in the case of the Wilson lines appearing in TMDFFs, Eq. (2.3). This rescaling is not necessary at NLO, where the contribution of the soft function is multiplied by δ(1 − x) (see details in Sec. 5), but it is necessary at NNLO and higher orders. The δ-regularized Wilson line violates the usual rules of gauge transformations. This violation is power-suppressed in δ. Therefore, throughout the calculation the δ should be considered an infinitesimal parameter, in order to avoid potential gauge-violating contributions. In most part of the calculation this is straightforward, however, the linearly divergent subgraphs should be carefully considered. A detailed discussion of this point, as well as other potential issues, can be found in [22].
The parameter ζ that appears in the factor R f is a scale that arises due to the splitting of the soft function among the two TMDs. In the calculation of the SF, one ends up with a function that depends on ln(µ 2 /(δ + δ − )). However, here the lnδ + and lnδ − represent the rapidity divergences related to different TMDs in the TMD factorization theorem. Therefore, one separates these logarithms introducing an extra scale ζ. In general one has (see e.g. [3,22] for more details) where ζ + ζ − = (p + p − ) 2 = Q 4 , with Q 2 being the relevant hard scale of the considered process. In the calculation of a single TMD (say the TMD oriented along the vector n), this operation can be effectively replaced by the substitution Here and in the following we omit the subscripts ± for the variable ζ.

Calculation of TMDs and their matching coefficients onto integrated functions
In order to calculate the leading matching coefficients of the OPE, we perform the calculation of TMD distributions on parton targets. Since at NNLO all possible flavor channels arise, we need to consider the following TMDs: where Z 2 and Z 3 are the wave function renormalization constant for quarks and gluons, respectively. During the calculation of the partonic matrix elements it is sufficient to put the momentum of the target parton at p 2 = 0. This condition is realized with p T = 0 in the momentum of target partons and restricting the light-cone momentum component p − = 0. Therefore, the momentum of the target parton is p = [p + , 0, 0 T ].
In the following, we denote by a superscript in square brackets the coefficient of the pertrubative expansions at a given order, e.g. for the partonic TMDFF f →f (z, b T , p; µ, ζ) . (3.14) The LO pertubative expansion of TMDs coincides with the unsubtracted matrix element, e.g. for quark-to-quark TMDFF, At NLO one finds The second term cancels the rapidity-divergent part from the unsubtracted expression, such that the TMD is finite when δ → 0. The last term cancels the UV divergences. After these subtractions the result remains singular for → 0 due to the collinear divergences that are part of the parton integrated FF. At NNLO the structure is richer All rapidity divergences arise and are canceled in the first line of Eq. (3.17), while in the second line we have just UV renormalization constants. In the case of TMDPDFs the perturbative expansion is the same (with the trivial substitution ∆ i → Φ i ). Finally, we calculate the matching of the TMDs onto their corresponding integrated functions. At LO the matching coefficients are trivially Comparing the matrix element at NLO we obtain C [1] f ←f = F [1] f ←f − f [1] f ←f , C [1] f ←r ⊗ f [1] r←f − f [2] f ←f , Notice the factor z 2−2 in the case of TMDFF, which comes from the operator definition in Eq. (2.17). The matching procedure in Eqs. (3.19-3.20) ensures the cancellation of the IR divergences in the matching coefficients. In our regularization scheme these divergences are regularized by dimensional regularization. That is why it is particularly important to know the dependence in Eqs. (3.19-3.20) at all orders in : one can immediately realize that the linear term in of the coefficient C [1] in combination with the single pole of f [1] contributes to the finite part of C [2] . Also, the coefficient z −2 gives a non-trivial contribution when combined with the poles of d [1,2] .

Anomalous dimensions of TMD operators
The renormalization group equations (RGEs) fix the scale dependence of the matching coefficients of the TMDs onto integrated functions, and follow from the very definition of the OPE, i.e. Eq. (2.19). Differentiating both sides of Eq. (2.18) with respect to the scales we obtain the RGEs for the matching coefficients in terms of the anomalous dimensions of TMD operators and integrated operators. The anomalous dimension of TMD operators is defined as Both the TMDPDF and TMDFF operators have the same anomalous dimension, as a result of the universality of the hard interactions [1][2][3]. The anomalous dimension γ f comes solely from the renormalization factor Z f . Using the standard RGE technique we obtain where AD represents the operator which extracts the anomalous dimension from the counterterm (i.e. it gives the coefficient of the first pole in with n! prefactor, being n the order of the perturbative expansion). The prefactor 2 arises from the normalization of anomalous dimension Eq. (4.1).
The flow with respect to the rapidity parameter follows from the factor R, and it is also the same for both types of operators, due to the universality of the soft interactions (see discussion in [22], also in [25]), The representation independence of non-Abelian exponentiation implies the so-called Casimir scaling of anomalous dimension D, see [22]: It is worth to mention that RGEs for TMD operators, in contrast to RGEs for integrated operators, do not mix the operators of different flavors. The rapidity anomalous dimension D f can be extracted solely from the prefactor R f [22] as where f.p. denotes the extraction of the finite part, i.e. neglecting the poles in . The singular part of the factor R is related to the renormalization factor as follows: where s.p. denotes the extraction of the singular part, i.e. the poles in . Note that these relations are independent of the regularization procedure, i.e. they hold for any rapidity regularization scheme. In the modified δ-regularization the explicit expressions for the soft function (and hence for the factor R f ) are presented in Appendices A and B. All relations in this Subsection are explicitly checked at NLO and NNLO, and the resulting anomalous dimensions, which are collected in Appendix D.2, coincide with the known values.
The consistency of the differential equations (4.1-4.3) implies that the cross-derivatives of the anomalous dimension are equal to each other ( [22,25]), (4.7) The first terms of the perturbative expansion of the cusp anomalous dimension Γ f cusp can be found in Appendix D.2. From Eq. (4.7) one finds that the anomalous dimension γ is where we introduce the notation At the level of the renormalization factors this relation allows one to unambiguously fix the logarithmic part of the factor R f , by means of the relation 2 (4.10) 2 A similar one can be found in [13,25].

RGEs for matching coefficients
The RGEs for the matching coefficients can be obtained by deriving both sides of Eq. (2.18). The only extra information which is needed is the evolution of the light-cone operator. That is given by DGLAP 3 equations where P and P are the DGLAP kernels for the PDF and FF respectively. The leading-order expressions are collected in Appendix D.2, while NLO expression can be found in [35,36].
Considering the derivative with respect to ζ we obtain the ζ-scaling for the matching coefficients The solutions of these differential equations are This defines the reduced matching coefficientsĈ andĈ, and their RGEs are where the kernels K and K are Using these equations one can find the expression for the logarithmical part of the matching coefficients at any given order, in terms of the anomalous dimensions and the finite part of the coefficient at one order lower. It is convenient to introduce the notation for the n-th perturbative order:Ĉ [n] The expressions for the anomalous dimensions, the recursive solution of the RGEs and the explicit expressions for the coefficients C and C are given in Appendix D. The known anomalous dimensions and DGLAP kernels allow to fix the logarithmic dependent pieces of the coefficients. As a result only the coefficients C

NLO computation
The calculation of TMDs at NNLO is a complex task. Technically it is convenient and safe to split it in several steps, and perform intermediate checks. In order to illustrate the procedure and also for Figure 1. Diagrams contributing to soft factor at NLO S [1] . The complex conjugated diagrams should be added pedagogical reasons, in this Section we present the NLO calculation of TMDs and their matching coefficients, with attention to some important details. Here and below Feynman gauge is used for the calculations. The Feynman diagrams for the bare TMDFFs at NLO are drawn in Fig. 2. The bare TMDPDFs are given by the same diagrams, but interpreting the external lines as the initial states and the momentum p as incoming. For the final/initial gluons we choose the polarization plane perpendicular to p µ and n µ . Thus, the possible diagrams with final/initial gluons radiated by the Wilson lines are zero, and are not shown in Fig. 2.
The only physical Lorentz-invariant scale present in the calculation is b 2 T , because the target parton is massless, p 2 = 0, and has no transverse components. The scale b 2 T appears only in the diagrams with left and right parts connected by gluon/quark exchange. Therefore, the pure virtual diagrams (i.e. diagrams without any cut propagator) are zero 4 . The only piece of the virtual diagrams which is relevant for our purposes, is the UV-divergent part, that enters the operator renormalization constants Z q and Z g . The pure virtual diagrams are independent of the kinematics and the operator, which implies that the renormalization constants Z q and Z g are the same for PDF and FF operators and independent of z and x. At NLO, the pure virtual diagrams are diagrams A in Fig. 2 (for quark-to-quark and gluon-to-gluon sectors), as well as diagram A for the soft factor in Fig. 1. Calculating the ultaviolet limit of the virtual diagrams we obtain Here it is important to preserve the previously defined order of subtraction of divergences (see Footnote ??). So, according to our definition we first recombine the rapidity divergences and then the UV-divergences.
Diagrams B and C in Fig. 2 provide the quark-to-quark matrix elements, where B B B = b 2 T /4. We have similar expressions for the other flavor channels. One can see that the expressions in Eq. (5.2) are connected by the relation The validity of this relation to all orders in perturbation theory can be proven in a diagram-bydiagram basis, comparing the expressions in both kinematics. In fact, if we do not remove any regulator, the TMDPDFs and the TMDFFs are related to each other by the crossing symmetry x ↔ z −1 . This is the generalization of the well-known Gribov-Lipatov relation between PDF and FF for the TMD operators. We have where the factor N arises from the difference of the operator normalization. Comparing the operator definitions we find N q,q = N g,g = 1, Before combining the collinear and soft matrix element, we develop our results in the limit δ → 0. This step allows to pass from the analytical functions Eq. (5.2) to the distributions, where the singularity at z, x → 1 is regularized. Within δ-regularization this step can be done using when the functions are regular at x, z → 0. In the case that the functions are singular at x, z → 0 (i.e. TMDFF and gluon distributions) we extract an extra factor of z as ∆(z, δ) = 1 z (z∆(z, 0)) + + δ(z) perturbation theory by analogous contributions from diagrams with real quark/gluon exchanges, see proof in [22].
The powers of δ are irrelevant for our calculation and are dropped. In the limit δ → 0 the expressions in Eq. (5.2) are Let us make a comment on the small-δ expansion in Eqs. (5.6)-(5.7). This operation breaks the analytical properties of the calculated functions in the complex plane of x, z. Therefore, at this stage of the calculation, one brakes the crossing relation between PDF and FF kinematics in Eq. (5.4). Indeed, the distributions in Eq. (5.8) are not analytical functions of x and z and can not be analytically continued to each other straightforwardly. That could be done using some regularization method, e.g. by restoring the δ-regularization parameter. This is a simple exercise at NLO but becomes involved at higher orders, see e.g. corresponding analysis for DGLAP kernels in [37]. In practice it results simpler to calculate the TMDPDFs and the TMDFFs independently, without using this analytical continuation property. In order to complete the calculation of the TMDs we have to include the contribution of the soft factor, which is computed at NLO and NNLO in [22]. At NLO the soft function is given by the diagrams shown in Fig. 1. The expression for these diagrams, by means of the substitution in Eq. (3.12) is where the color prefactor depends on the representation of Wilson line: C K = C F (C A ) when a quark (gluon) is the initiating parton. Combining Eq. (5.8) and Eq. (5.9) one can immediately check the exact cancellation of the rapidity singularities in the limit δ → 0 (represented by λ λ λ δ ) between unsubtracted TMDs and the soft factor.
Expanding in and combining together all the pieces of the TMD matrix element according to Eq. (3.16), we obtain where p qq (x) = (1 + x 2 )/(1 − x). This is the final expression for the TMD partonic matrix elements. They are free from the rapidity and UV divergences, as predicted by the TMD factorization theorem [1,3,4]. The final expressions for unsubtracted TMD at NLO for all other flavor channels are similar and are collected in the Appendix A. In Eq. (5.10) one recognizes the -pole, which is part the corresponding integrated functions. In order to complete the matching between the TMDs and integrated functions we need to calculate the matrix elements of the integrated operators. The diagrams contributing to these matrix elements are all zero, due to the absence of a Lorentz-invariant scale. Therefore, the only non-zero term is the UV renormalization factor, which can be deduced from the DGLAP kernel. So, for quark-to-quark channel we have The matching prescription of Eq. (3.19) allows to derive the coefficients where -poles are exactly cancelled. The final matching coefficients for TMDPDFs at LO are and the rest are zero. At NLO we find where the definitions of functions p(x) are given in Appendix A, Eq. (A.1). Hereafter we follow the notation/convention of [35,36], so that the piece of the coefficients divergent at x, z → 1 should be understood as "plus"-distribution. The expression for C [1] q←q has been already obtained in many articles (see e.g. [2,7,8,[12][13][14][15][16][17][18][19]38]), the expression for C [1] q←g is also well-known [7,8,17,19] (note that there is a misprint in [19]), and the expression for C [1] g←q and C [1] g←g have been obtained in [17,19].
The matching coefficients for TMDFFs at LO are and the rest are zero. At NLO we find The functions p(z) are related to the one-loop DGLAP kernels and are defined in Eq. (A.1). The coefficient C q→q has been calculated in [4,8], and C q→g agrees with the one calculated in [8]. The coefficients C [1] g→q and C [1] g→g are presented here for the first time. One can see that the expression for the matching coefficients for TMDFFs have an extra ln(z) in comparison to TMDPDFs. This contribution comes from the difference in the normalization factor z −2 , see Eq. (2.17). This logarithm is the main source of difference between the TMDFF and TMDPDF matching coefficients. At higher orders the effects of the z −2 normalization factor are more involved.

NNLO computation
The diagrams that contribute to the unsubtracted TMD matrix elements can be generically classified in pure-virtual diagrams (i.e. diagrams with no cut propagator), virtual-real diagrams (i.e. diagrams with one single cut propagator) and double-real diagrams (i.e. diagrams with two single cut propagators). Alike NLO case, pure-virtual diagrams are zero due to absence of a Lorentzinvariant scale. Virtual-real and double-real diagrams are proportional to B B B 2 . In total, there are about 50 virtual-real diagrams and about 90 double-real diagrams.
The generic expression for virtual-real diagrams is (for TMDFF kinematics) where for brevity we drop the i0-prescription of propagators. The function F contains all "plus"components of the momenta and the parameter δ, while the function f contains only scalar products of momenta. The discontinuity of the propagator is The generic form of a double-real diagram in the same notation takes the form 3) The functions f can be re-expressed via the propagators, and so the diagrams can be split into several integrals with i a i = 3 for virtual-real diagrams and i a i = 2 for double-real diagrams.
In order to decouple the functions F from the scalar loop integrals we introduce the auxiliary unity factor With the help of this trick the dependance of functions F on k + and l + can be re-written as a function of z and ω, and all numerators simplify. The integration over the loop-momenta is straightforward and all non-zero integrals appearing in the calculation are presented in Appendix C.
In this way, we are left with a set of one-dimensional integrals over ω. The evaluation of these integrals is technically the most difficult part of the calculation. Most part of these integrals are evaluated in terms of Γ-functions and their derivatives, while several are expressed through hypergeometric functions (and one integral in g → g and g → q channels that has been expressed via Appell function F 1 ). All diagrams are calculated in d = 4 − 2 dimensions.
During the evaluation of the integrals we have used that we need only their asymptotic behavior at δ → 0. In order to find the small-δ limit we expand the eikonal propagators in Mellin-Barnes contour integral around δ = 0. Then we calculate the integrals over ω, and close the contour over the closest to zero poles. If an integral has a singularity at z → 1 it should be regularized by means of a "plus"-distribution (see Eqs. (5.6-5.7). The final expression for a diagram takes the generic form It is important to mention that the functions f 2 and f 3 exactly cancel in the sum of all diagrams, which we have checked explicitly. The unsubtracted TMDPDFs can be calculated in the same manner.
Due to the symmetry of the operators, the expressions for TMDPDFs and TMDFFs satisfy the crossing relation Eq. (5.4) in a diagram-by-diagram basis. However, since we consider only the leading contribution at δ → 0, the diagram-by-diagram crossing is violated, due to the fact that IR singularities and rapidity singularities get different phases during the procedure of analytical continuation. In the sum of diagrams all terms f 2,3 cancel, and one can check the crossing relation in Eq. (5.4) without any special tricks. Since we calculated TMDPDFs and TMDFFs independently, such a relation grants a very strong check for our results (however we have not compared the δcontribution for q → q and g → g channels, for the reasons explained earlier).
Having the expressions for the unsubtracted TMDs, we multiply them by the rapidity and UV renormalization factors. At NNLO this procedure is given by Eq. (3.17). The expressions for the factor Z [1] and soft factor S [1] are given in previous Section. The NNLO expression for the soft factor has been obtained in [22] and is given in Eq. (B.1). At this stage we also perform the expansion in of the expressions. To perform the renormalization procedure in Eq. (3.17) we have to calculate the operator renormalization constants at NNLO Z [2] q and Z [2] q . They are given by the UV-part of pure-virtual diagrams. In our calculation we have not calculated these constants explicitly, but found them by demanding the cancellation of UV poles. We obtain the following expressions: As was discussed in Sec. 4.1, most part of the UV counterterm should be related to the factor R and to the known anomalous dimensions. Finally, we perform the matching procedure as in Eq. (3.20). The integrated matrix elements are zero due to the absence of a Lorentz-invariant scale and are given solely by their UV renormalization counterterm. They can be deduced from the DGLAP kernels, and given by f [2] f ←f = The obtained matching coefficients are free from any kind of divergences. The results of the calculation are presented in next Section.

Expressions for matching coefficients
In this section we present the expressions for the finite part of the small-b T matching coefficients. The logarithmic part can be restored by using the RGEs and is explicitly given in Appendix D.1. For completeness we present LO, NLO and NNLO finite parts together.

TMD parton distribution functions
The LO matching coefficients are and all other flavor configurations are zero at leading order.
The NLO matching coefficients are Here, the coefficient C q←q (x) is the coefficient with all possible mixing flavor channels. Thus q can be any quark or anti-quark, even of the same flavor as q. In other words, the matching coefficient for, say, u ← u is given by the sum C q←q + C q←q , as well as the matching coefficient for u ←ū is given by C q←q + C q←q .
The NNLO matching coefficients are 3) These matching coefficients were first calculated in [15][16][17] by a direct calculation of a crosssection, and in an SCET framework in [19,20]. Our results agree with these previous calculations once the proper combination of collinear and soft matrix elements is considered.

TMD fragmentation functions
The LO matching coefficients are the rest flavor configurations are zero at LO.
The results for the quark sector were first presented by us in [21] 5 . The mixed flavor and gluon contributions are presented here for the first time. Moreover, to the best of our knowledge, the NLO expressions for gluon TMDFF are also presented for the first time.

Matching coefficients at threshold
Using our NNLO results for TMDs we observe that it is possible to find a recurrence in the behavior of the matching coefficients for x, z → 1.
In order to establish the idea, we recall that the matching coefficients for processes where collinear factorization applies, behave as a k s (lnx) 2k−1 /x [42,43]. These corrections are dominant for x, z → 1 and have to be resummed for phenomenological applications (threshold resummation). In the case of TMDs, by analyzing the structure of divergences one may expect the leading behavior to be at most like a k s (lnx) k−1 /x. In fact, in the case of TMDs the singular behavior at x, z → 1 should also be universal, due to the universality of the soft function. This statement can be seen in the following way: in the regime x, z → 1 the real soft gluon exchanges in Feynman diagrams are dominant, and are the source of the rapidity divergences in TMD operators. The rapidity divergences are removed by the R f factors which are universal for both PDF and FF kinematics. Thus, the leading x, z → 1 behavior of TMDs should be the same. At the same time, the leading asymptotic term of the integrated functions is independent of the kinematics and goes like ∼ Γ cusp /(1 − x) + [44], so that we expect the behavior of the matching coefficients to be also universal in the threshold limit.
At two loops, the leading singular behavior at x, z → 1 should be the same for gluons and quarks (up to a trivial change in the color factor), since it is produced solely by the convolutions of one-loop soft subgraphs, which are the same for quarks and gluons. Indeed for TMDs F f ←f and D f →f we observe from our results that where dots denote the less dominant contributions and collinear poles. The sub-leading contribution, proportional to 1/(1 − x) + or 1/(1 − z) + , is different for gluons and for quarks and depends on l ζ , as expected. We observe that the difference between the gluon and quark channels, as well as the dependence on ζ, disappear after the matching procedure. In fact, we obtain a simple expression for the leading term at x, z → 1: where the dots denote the contributions with δ-functions and the non-singular terms at x, z → 1.
The values of d (2,i) can be found in Appendix D.2. In Eq. (8.2), the µ-scale dependent terms follow from the RGE, while the coefficient for the finite part is peculiar, because it is directly connected to the perturbative expansion of the D f function, which governs the evolution of the TMDs. If one then extrapolates a similar behavior to an arbitrary loop order, we can make a conjecture for the leading term at x, z, → 1 for the finite part of the TMD matching coefficients, based on one-and two-loop calculations: The terms proportional to L k µ can be deduced from the general formulas of Appendix D.1 in the threshold limit. Notice that according to this conjecture, and using the recent result for d (3,0) obtained in [41], one can give an estimate of these coefficients at threshold at N 3 LO.
As a final remark, we notice that, since the soft function enters the polarized TMDs on the same footing as unpolarized TMDs, a similar result can be obtained for all of them.

Conclusions
In this paper we present a comprehensive study of the unpolarized TMDs at NNLO. To make it as general as possible, we have introduced the TMD operators, such that TMDs are matrix elements of these operators. We find that the understanding of the TMDs benefits from such a language, as provided by the present work. In fact, in these terms, it is possible to introduce a common formalism to describe the universality of soft interactions, the parallelism between the renormalization of UV divergences and rapidity divergences in the TMDs, and their matching onto integrated functions. In addition, the consideration of any TMD can be performed without an explicit reference to any given process.
The TMD operators are involved objects, which contain both rapidity divergences as well as UV divergences, and thus are different from usual light-cone operators. The rapidity divergences can be absorbed by "rapidity renormalization factors", alike the usual UV divergences. The explicit form of "rapidity renormalization factors" is obtained from the factorization theorems for semi-inclusive DIS, Drell-Yan and e + e − → 2 hadrons [1][2][3][4]. It is important to note that the "rapidity renormalization factors" are the same for all kind of TMD processes, for distribution and fragmentation kinematics and that, together with the UV renomalization, they give direct access to the RGE for TMD operators. That completes the analogy with UV renomalization and it allows to construct a universal TMD operator.
One of the main outcomes of the paper are the matching coefficients of all the unpolarized TMDs onto their integrated analogues. According to the operator language they are the Wilson coefficients for the leading term in the small-b T operator product expansion, as explained in the text.
The calculation of TMD matrix elements needs a rapidity regulator in addition to a UV regulator. For that we have used the (modified) δ-regularization [21], in which the form for the "rapidity renormalization factor" is especially simple. The presented matching coefficients for the TMDPDFs agree with the results of [18][19][20], once the proper combination of collinear and soft matrix elements is considered. Here, instead, we have provided a method that realizes the cancellation of rapidity divergences within a single TMD, and we have checked this fact explicitly at NNLO. The results for quark TMDFF were partially presented in [21], while here we provide the complete results, which include also the gluon TMDFFs, that were unknown. All these matching coefficients are necessary for accurate phenomenological studies, and allow to consider exclusive and inclusive processes on the same level of theoretical accuracy.
The performed calculation has a complex structure which involves the calculation of TMD matrix elements, integrated matrix elements, TMD soft factor and TMD renormalization constants at NNLO. Some of these ingredients are already known at NNLO. So, the TMD soft factor has been presented by our group in [22] and the integrated matrix elements can be related to DGLAP kernels in our regularization scheme. The TMD matrix elements and TMD renormalization constants have been computed at NNLO in this work for the first time with the δ regulator. During the calculation of TMD matrix elements we have used many checks, which include: a check of logarithmic parts by RGEs, an independent extraction of anomalous dimensions and a check of the crossing relations between TMDPDF and TMDFF. The regularization method described and implemented here, together with the results for the master integrals, can be useful also for the study of polarized TMDs.
In addition, we have studied the limit of large x,z and found that the behavior of the matching coefficients is universal for all unpolarized TMDs. It is naturally controlled by the anomalous dimension D f , which allows us to make an all order conjecture on the leading contribution at large x, z.
The obtained matching coefficients are necessary in order to pursue phenomenological studies at N 3 LL accuracy. Some recent developments towards this goal can be found in [23]. There, although using a different regulator for rapidity divergences [41], the authors have assumed the structure of rapidity divergences which has been explicitly checked in the present work. Together with the large-x conjecture presented in this work, it opens the door to a very precise estimate of these perturbatively calculable contributions. The phenomenological applications of these results will be exploited in future works. We expect all these efforts to be necessary in order to have a unified picture of Drell-Yan, semi-inclusive DIS and e + e − → 2 hadrons.
Note added: while this article was under submission G. Lustermans, W. J. Waalewijn and L. Zeune [46] confirmed the threshold behavior of the coefficient obtained in Section 8.
where B B B = b 2 T /4. The singularities at x → 1 are understood as "plus"-distribution. The unsubtracted TMDFFs are The matrix elements of the integrated functions are given solely by their UV counterterms. For the PDF kinematics they are For the FF kinematics we have The expression for the NLO soft factor is For the completeness of exposition we also present the renormalization constants for fields

B NNLO expressions
In this appendix we present all side expression used for NNLO calculation.
The soft factor at NNLO has been calculated in [22]. We present the NNLO contribution to the exponent Eq. (3.8). The -expansion of NNLO soft factor reads for quark (gluon) soft-factor. Here, the logarithm l δ is ln µ 2 /|δ + δ − | , while after substitution Eq. (3.12) it reads The constants d (n,k) are given in Sec. D.2. The NNLO TMD operator constants are calculated in Sec. 6 and reads The NNLO field renormalization constants are [39] Z

C Results for integrals
In this appendix we present the loop integrals that are used to calculate the TMD PDF and TMD FF at NNLO. The parameter ω is introduced to in order to resolve the k + and l + dependance as explained in Sec.6.

C.1 Integrals for virtual-real diagrams
The scalar integrals for the virtual-real diagrams has generally the form (C.1) The corresponding integral in PDF kinematics reads For our observable only integral with sum of indices equal to 3 contribute. The momentum p has no transverse component and p 2 = 0. Due to it, the integral with decoupled virtual loop are zero. For example: Integrals with negative index can be rewritten using identity The only non-zero integrals with positive indices are where [dx] = δ(1 − x 1 − x 2 − x 3 )dx 1 dx 2 dx 3 . We leave the integral over Feynman parameters in F 10101 , since it is convenient first over ω with the help of δ-function. There are two another integrals that appear in calculation and can be reduced to the previous cases F 00111 (ω) = F 10101 p + + k + p + − ω ,

C.2 Integrals for double-real diagrams
The scalar integrals for the double-real diagrams have generally the form (C.6) The components k + and l + must be integrated with the help of δ-functions as explained in Sec. 6 and do not participate in the loop-integration (that is indicated by d − 1-dimensional integral).

D Recursive relations from RGE and anomalous dimensions
In this Appendix we collect all the expressions necessary for the application of the RGEs, as well as the explicit expressions for the logarithmical part of the matching coefficients.

D.1 Recursive form of RGEs
The derivation of RGEs is given in Section 4. The ζ-dependence of the matching coefficients can be explicitly solved by Eq. (4.13). The µ-dependence is then given by the RGEs in Eq. (4.14). For practical purposes it is convenient to rewrite the RGE application in a recursive form. We use the notation of Eq. (4.16). Then the equation of the logarithmic dependent part of the TMDPDF matching coefficient reads h←f (x) . (D.1) The same logarithmic part of the TMDFF matching coefficient reads f ←f = δ f f δ(x) The expressions for TMDFF matching coefficients can be obtained from these ones by changing the directions of the arrows and replacing DGLAP kernels as P → P/z 2 . Explicit expression for these equations can be found in a supplementary file [40].

E Alternative form of matching coefficients
For practical purposes, it is convenient to write the matching coefficients as overall "plus"-distributions.
In this Appendix we rewrite the expressions for the matching coefficients in such a form. Only the flavor-diagonal coefficients need to be rewritten in this way, since the non-diagonal channels are integrable at z, x → 1.
The NLO expressions read