Calculation of the transverse parton distribution functions at next-to-next-to-leading order

We describe the perturbative calculation of the transverse parton distribution functions in all partonic channels up to next-to-next-to-leading order based on a gauge invariant operator definition. We demonstrate the cancellation of light-cone divergences and show that universal process-independent transverse parton distribution functions can be obtained through a refactorization. Our results serve as the first explicit higher-order calculation of these functions starting from first principles, and can be used to perform next-to-next-to-next-to-leading logarithmic $q_T$ resummation for a large class of processes at hadron colliders.


Introduction
Parton distribution functions (PDFs) are fundamental properties of hadrons. They describe the distributions of quarks and gluons inside hadrons. Their usefulness in collider phenomenology resides in the factorization theorems, in which short range interactions are separated from long range effects. The short range interactions result in perturbatively calculable hard scattering kernels or hard functions, while the long range effects are encoded in non-perturbative or semi-perturbative objects such as soft functions, parton distribution functions and fragmentation functions. From these, it is possible to make predictions for collider observables based on perturbative calculations and experimental determination of the non-perturbative functions. This procedure has proven highly successful in the last thirty years. For most applications, the relevant factorization theorems are "collinear factorization" developed in [1][2][3]. The corresponding non-perturbative functions are so-called "collinear PDFs" or "integrated PDFs", in which only the partonic momentum component along the direction of the colliding hadron is kept, while all other components are integrated over. The collinear PDFs have gauge-invariant definitions in terms of matrix elements of nonlocal bilinear operators [4,5]. While the collinear PDFs are non-perturbative functions, their scale-dependence can be calculated perturbatively in terms of the DGLAP splitting kernels. These calculations have been performed up to 3 loops [6][7][8][9][10][11][12].
In this work, we will consider TMD factorization, where transverse parton distribution functions (TPDFs) and transverse fragmentation functions are introduced. Historically, TMD factorization frameworks were developed in three different kinds of kinematics: e + e − to hadrons, semi-inclusive deep-inelastic-scattering (SIDIS), and Drell-Yan type processes. In this work, we will mainly be concerned with hadron collider physics, and will therefore discuss TMD factorization and TPDFs for (unpolarized) Drell-Yan type processes in detail.
Consider the production of a vector boson in hadron-hadron collisions with its invariant mass Q and transverse momentum q T observed. If Q ∼ q T Λ QCD , one expects that collinear factorization is valid and the differential cross section can be factorized as (schematically) C ij (z, Q, q T , µ) are hard scattering kernels describing physics at the hard scale Q ∼ q T . The symbol ⊗ denotes convolution. Consider now another phenomenologically important region Q q T , Λ QCD . In this region, even if q T is in the perturbative domain, use of the collinear factorization formula (1.1) will lead to problems with the perturbation series due to the appearance of large logarithms of Q/q T . Therefore, one would like to factorize the two scales and resum the large logarithms to all orders in perturbation theory. Ideally, one might expect a factorization formula similar to eq. (1.1): where the hard functions H ij (z, Q, µ) describe physics at the hard scale Q, and the TPDFs B i/N (x, k T , µ) describe physics at the low scales q T and Λ QCD . However, things are not so simple. It turned out that theB functions necessarily depend on the hard scale Q through some unphysical parameter, denoted here by ξ. Moreover, depending on the regulator used and the process under consideration, an additional soft function may appear. Therefore, the correct formula looks like Note that while the individual functions depend on the unphysical parameters ξ 1 and ξ 2 , in physical cross sections these parameters are combined in a way that only the physical scale Q remains. Eq. (1.3) is not a true factorization since theB andS functions still involve the two widely separated scales Q and q T . To achieve a proper factorization of the two scales, one needs to extract the Q dependence from theB andS functions by studying their dependence on the unphysical parameters ξ 1 and ξ 2 . This procedure has various names in the literature: Collins-Soper equation in the pioneering works [15][16][17][18][19][20], rapidity renormalization group in [23,25], and refactorization in [21,26]. After this, the Q dependence can be exponentiated and one can obtain the true TPDFs. The anomalous Q dependence of the naive TPDFsB arises as follows. Similar to the collinear PDFs, the naive TPDFs can be defined as matrix elements of non-local operators. This was given in axial gauge in [4] and was rendered gauge-invariant in [5,18,19] by introducing Wilson lines. In these works, it was pointed out that by taking the gauge-fixing vectors to the light cone, or equivalently putting the Wilson lines on the light cone, one encounters singularities not regularized in dimensional regularization in the perturbative calculations. Therefore, to perform the calculations, one has to introduce an extra regulator. In [4,18,19], this was achieved by taking the gauge-fixing vectors or the Wilson lines off the light cone. Other choices of regulator are possible. For example, variations of the analytic regulator were used in [21][22][23]25]. In [24], finite imaginary parts in certain propagators were used to regulate the light-cone divergences. No matter which regulator is used, the anomalous Q dependence inevitably arises which, in the language of soft-collinear effective theory (SCET) [30][31][32], is due to the breaking of the rescaling invariance of the Lagrangian [21]. In the Collins-Soper approach [16], the explicit appearance of those singularities is circumvented by always considering the product of two TPDFs. Very recently, this approach was worked out in full generality to describe arbitrary color-singlet final states in hadronic collisions in [33].
Depending on the relative size of the two scales q T and Λ QCD , different aspects of TPDFs can be described in perturbative QCD. If q T ∼ Λ QCD , the TPDFs are fully nonperturbative and only their scale-dependence can be calculated perturbatively. In the situation q T Λ QCD , the TPDFs are semi-perturbative objects and one can further factorize the two scales. In this region, the TPDFs can be expressed as convolutions of the collinear PDFs with perturbatively calculable matching coefficient functions. These coefficient functions can be obtained via two approaches. The first is assuming the factorization formula (1.3), and extracting the coefficient functions by studying the small q T behavior of the differential cross section. This is the approach taken by [34,35], which is generalized to any color-singlet process in [33]. The second approach is starting from a gauge-invariant operator definition of the TPDFs, and straightforwardly computing the operator matrix elements. This approach is much more challenging since one directly encounters the lightcone divergences. However, the second approach, once accomplished, serves as an explicit verification of the TMD factorization framework. In [36], we derived the next-to-next-toleading order (NNLO) coefficient function for quark-to-quark transitions. Results at this order are for example relevant for a next-to-next-to-next-to-leading logarithmic (N 3 LL) transverse momentum resummation. This paper extends the calculation of the NNLO coefficient functions to all partonic channels and describes technical and methodological details. We also present several consistency checks on our results. In particular, we reproduce the process specific H (2) coefficients of [34,35] for Drell-Yan process and Higgs production, as well as the order α 2 s contributions to the DGLAP splitting kernels. This paper is organized as follows. In Section 2 we introduce our calculational framework. We provide the operator definitions of the TPDFs and the regularization of the light-cone singularities. We outline the procedure of the NNLO calculations in Section 3, with some detailed expressions collected in Appendix A and B. The main results are presented in Section 4, while several additional relations are collected in Appendix C and D. We conclude in Section 5.

Framework
We consider the collision of two hadrons N 1 and N 2 with momenta p andp producing some color-neutral final state F of momentum q and additional unresolved remnants X Along the directions of the hadrons we specify two light-like vectors n andn with n ·n = 2.
In terms of them any 4-vector can be decomposed as where q µ ⊥ is perpendicular to both n andn. We define q 2 T = −q 2 ⊥ . We consider the differential cross section for the production of the final state F with respect to its squared invariant mass q 2 , transverse momentum q T , and rapidity y. We are especially interested in the region where the transverse momentum is much smaller than the invariant mass q 2 q 2 T . For this multi-scale problem we need to achieve the factorization of disparate scales and the resummation of the corresponding logarithms. This was done for the Drell-Yan process in the pioneering work [16]. In this paper, we will mainly follow the SCET based language in [21,26], in which the factorization formula for the Drell-Yan process can be written as where the first factor corresponds to the Born level cross section, C V is the Wilson coefficient obtained from matching the quark form factor to the effective theory. Together they form the process specific hard function. The transverse position (impact parameter) x ⊥ is the Fourier conjugate variable to q ⊥ , and x 2 T = −x 2 ⊥ . The position-space soft function S is a correlator of soft Wilson lines, and B andB are the two position-space TPDFs. In terms of x 2 T they depend on the transverse variable. Other functional dependences related to the regularization are implicit in the above expression and the −i prescription for the Wilson coefficient defining the sign of its imaginary part will be suppressed from now on. The factorization theorem holds up to power corrections in q 2 T /q 2 .

Definition of transverse PDFs
The bare quark TPDF collinear to the n direction is represented by the gauge invariant operator matrix element [21] B q/N (z, where the sum is over all intermediate states X and the summation over the spinor indices α, β of the gauge invariant collinear quark field χ n in SCET is understood. The corresponding anti-quark TPDF is given by a similar equation with the role of the fields χ n andχ n interchanged. The TPDFs along the opposite direction, to which we refer as anti-collinear, are given by the same expressions, but with p ∼ n andp ∼n interchanged. The regularization of the rapidity divergences which we will outline in section 2.2 actually leads to a breaking of this relation. To mark this difference the anti-collinear TPDF will be denoted asB. In most aspects the discussion of these two different TPDFs is, however, completely analogous and for simplicity we therefore usually formulate it below only in terms of the collinear function. For processes initiated by gluon-gluon fusion, factorization theorems similar to eq. (2.3) hold in which gluon TPDFs are encountered. Along the n direction the latter is represented by the operator matrix element [26,37] where A µa n,⊥ is the gauge invariant collinear gluon field in SCET and the sum over the color index a is understood. Note that the gluon TPDF is a Lorentz tensor [26,37,38] in the space perpendicular to n andn. It can be decomposed into two independent components as where d is the number of space-time dimensions, g µν ⊥ is the metric tensor in the transverse space and the projection onto the two components are given by If the transverse scale is in the perturbative region, x T 1/Λ QCD , the physics of these two scales can be factorized and the TPDFs can be matched onto the collinear PDFs defined as where the sum is over all partons j. This holds up to power corrections in x 2 T Λ 2 QCD and we introduced the symbol ⊗ to denote the convolution of two functions as f (z, · · · ) ⊗ g(z, · · · ) ≡ 1 z dξ ξ f (ξ, · · · ) g(z/ξ, · · · ) . (2.10) The matching kernels I i/j and I g/j are perturbative functions. They can be extracted from eq. (2.9) and the perturbative parton-to-parton (T)PDFs B i/j , B g/j and φ i/j given by eqs. (2.4, 2.5, 2.8) with a parton j in place of the hadron N .
With the results of the matching kernels and eq. (2.9), the semi-perturbative hadronto-parton TPDFs can be obtained from the collinear PDFs as long as x T is a perturbative scale. The collinear PDFs have been extracted with high precision from experimental data by several groups. The knowledge of the matching kernels therefore provides an accurate determination of the TPDFs in the semi-perturbative domain. This not only has many phenomenological applications on its own right, but also provides necessary information for the determination of the fully non-perturbative TPDFs.
While I i/j starts at α 0 s , I g/j only starts at α 1 s . In many gluon-gluon initiated processes of interest, for example the production of a Higgs boson, the hard tensor contracting the Lorentz indices of two gluon TPDFs does not mix their two tensor structures. Then for the same level of accuracy for physical observables, the perturbative expansion of I g/j is required to one order less in α s than the expansion of I g/j . In these cases the NLO expression of I g/j which was derived previously in [26] suffices for N 3 LL precision. For this reason, the main goal of this article is to determine the NNLO corrections to the matching kernels I i/j from those of the parton-to-parton (T)PDFs, while I g/j at this order will be discussed in a forthcoming article.

Treatment of singularities
In the calculation of B i/j we have to deal with several kinds of singularities. On the one hand there are the usual ultraviolet (UV) and infrared (IR) singularities, which we regulate by dimensional regularization in d = 4 − 2 dimensions. On the other hand the functions contain extra light-cone singularities which require additional regularization. As mentioned in the introduction, there are several proposals of regulators. In our calculation, we use the analytic regulator as suggested in [22]. This amounts to introducing in the phase space integrals a factor (ν/n · l i ) α for each unresolved final-state parton with momentum l i . Here α is the analytic regulator and ν is an unphysical mass scale associated with the regulator -in a similar way as the renormalization scale µ is related to the dimensional regulator . Note that the regulating factor has to contain the same light-cone component n · l i for both the collinear and the anti-collinear region. As such, it breaks the symmetry p ∼ n ↔p ∼n between the two regions and the rescaling invariance of SCET, a fact called "collinear anomaly" in [21]. These symmetries will be restored at the end of the calculation when all divergences are removed and the limit α → 0 is taken. One good property of this scheme is that the soft function S (for Drell-Yan like processes) automatically reduces to a trivial factor of unity. For other regulators where this is not the case, one may always absorb the soft function into the two TPDFs by a redefinition (see, e.g., [5]).
The analytic regulator combined with the dimensional one suffices to regulate all singularities in the operator matrix elements. The singularities manifest themselves in the TPDFs as poles in the regulators. While the individual factors S, B andB are scheme dependent, their product is well defined and especially all poles in the regulator α cancel therein, along with the dependence on the unphysical scale ν after the limit α → 0 is taken. However, a dependence on the hard scale q 2 remains even after the regulator is dropped. The generation of the hard scale q 2 ∼ (n · p)(n ·p) through the analytic regulator can be understood from the scale ratio (ν/n · l i ) appearing along with it. For the anti-collinear region this ratio can be expressed in terms of (ν/n ·p), while for the collinear region it can be expressed in terms of (νn · p x 2 T ). In the combination of the two factors, the scale ν drops out, while the mass ratio q 2 x 2 T is left over. Extending these arguments, using the existence of the α → 0 limit of the product of two corresponding TPDFs and its independence on the scale ν, it was shown in [21] that the product can be refactorized into the form where after the cancellation of all poles in α on the left hand side the analytic regulator is set to zero and the right hand side is free of both α and ν. This defines the anomaly coefficient F and the true TPDFs B i/j which are universal process-independent functions and have the same form for the collinear and anti-collinear region. These functions still contain poles in as indicated by the label b for bare. By the operator renormalization the UV poles are absorbed into the renormalization factors Z, such that the renormalized functions B i/j (z, x 2 T , µ) and F iī (x 2 T , µ) are free of these singularities. Upon renormalization a dependence on the renormalization scale µ is introduced which is described by the renormalization group equations (RGEs) and will be discussed further below. We work in the MS scheme, which amounts to expressing the bare coupling constant as 14) and requiring that the renormalization factors contain only poles in . After renormalization, F is free of any poles, while B i/j can still contain IR poles. This signals the non-perturbative nature of the TPDFs. These IR poles are exactly the same as those in the collinear PDFs, whose renormalization takes the form Just as the functions B i/j and φ b i/j are related by eq. (2.9), the transverse and collinear PDFs are related by matching kernels via in both their renormalized and bare versions. This relation, the renormalization of B and φ as well as the result to all orders in dimensional regularization imply the renormalization of the matching kernels as In this equation the UV poles are contained in Z B i and the IR poles in φ k/j , while I i/k is free of any poles. In fact even though we do not explicitly distinguish IR and UV poles in our calculation, this equation allows us not only to extract the renormalized matching kernel, but also separately the renormalization factor Z B i and the renormalized PDFs. The separation of the last two functions can be achieved by fixing the endpoint contributions of the renormalized PDFs φ j/j from constraints on their integrals implied from momentum and quark number conservation.

Resummation
The differential cross section (2.3) can now be written as [21] which holds up to power corrections in q 2 T /M 2 and x 2 T Λ 2 QCD with the perturbative function The functions appearing here can be related to the quantities A, B and C ij as defined in [16]. These relations are given in eqs. (71, 72) of [21]. In eq. (2.20) each function depends only on a single physical mass scale and can be determined consistently in fixed order perturbation theory by choosing the scale µ in the vicinity of that scale such that no large logarithms are present. In a subsequent step, all functions have to be matched at the same scale µ. This is achieved by solving the RGEs for each of them, which automatically resums all large logarithms.
In terms of the logarithm the DGLAP splitting kernels P jk (z), the cusp anomalous dimension in the fundamental (adjoint) representation Γ q cusp (Γ g cusp ) and the quark (gluon) anomalous dimension γ q (γ g ), which are all listed in appendix C, the RGEs can be written as The last equation follows from the RGEs of the (T)PDFs and eq. (2.16). Eq. (2.22) takes a corresponding form for other processes; for gluon initiated processes with the anomalous dimensions Γ g cusp and γ g . This equation and the independence of the cross section on µ imply eqs. (2.23, 2.25). Also note the appearance of the hard scale q 2 in eq. (2.22) already implied a compensating dependence on this scale for the other factors. This has been found in terms of the collinear anomaly, eq. (2.11).
Since the bare functions do not depend on µ, each renormalization constant in eqs. (2.12 -2.15) obeys a RGE which exactly compensates the µ dependence of the corresponding renormalized function. Solving these equations and enforcing the MS condition, which is most conveniently done using the d dimensional coupling constant, allows us to express the renormalization constants in terms of the corresponding anomalous dimensions and the QCD β function. The results for φ i/j , Z B i and Z F i are listed in appendix D.1. Comparing these expectations with the findings in our calculation serves as a check on our results.
Provided that all coefficients in eqs. (2.20 -2.26) and the QCD β function are determined to sufficient order, any logarithmic precision goal for the differential cross section can be achieved. To obtain e.g. N 3 LL precision, the Wilson coefficient and I i/j have to be known to α 2 s , while F iī , P kj and γ i are needed to α 3 s . Moreover, Γ cusp and β are needed to α 4 s . Only some of them are known to this accuracy, which are β in [39], P kj in [11,12], γ i in [40] and for several processes also the Wilson coefficients.
The derivation of the α 2 s contributions to I i/j , as required for the N 3 LL transverse momentum resummation, is the main objective of this paper. These and F iī up to α 2 s can be obtained in the way outlined in this section from a perturbative calculation of up to n = 2. This calculation is discussed below. The main results, the NNLO matching kernels I (2) i/j , are presented in section 4. For completeness we also list further relevant perturbative results in appendix D.

Perturbative calculation
Once the perturbative results for the parton-to-parton (T)PDFs to sufficient order in the strong coupling and the two regulators are determined, the extraction of the final results according to eqs. (2.11 -2.18) is straightforward. We therefore only discuss the former in more detail. We begin with the collinear case.
Since the relevant matrix elements (2.4, 2.5, 2.8) contain solely collinear fields and the purely collinear SCET Lagrangian has the same form as the full QCD Lagrangian, we can use QCD Feynman rules to evaluate them. In a general gauge, any number of gluons can couple to the Wilson lines contained in the gauge invariant fields (χ, A) and the associated vertices lead to denominators with momentum components projected to then direction. A special gauge is the light cone gauge withn chosen as the light cone vector. In this gauge the Wilson lines reduce to factors of unity, but one still finds then dependent denominators -this time introduced through the gluon propagators. We will focus our discussion to this gauge, although we also performed the calculation in Feynman gauge as a cross check.
In our regularization scheme, the perturbative corrections to the bare collinear PDFs lead to scaleless integrals vanishing in dimensional regularization, such that their all order results are given by eq. (2.17).
For the transverse PDFs the additional scale x ⊥ is present. The corresponding expressions essentially correspond to the square of matrix elements as in figures 1 and 2 where the momentum p − k of the parton coupling to the gauge invariant field can be off-shell. Calling the momenta of the emitted partons l i , i = 1, . . . n r , and their sum k = i l i , the phase space factor takes the form wherek z =n · [k − (1 − z)p]. The last factor arose from the t integral in eqs. (2.4, 2.5), the exponential from the x ⊥ dependence of the gauge invariant fields and the α dependent factors arise from the analytic regularization. It is essentially these factors which lead to difficult integrals, where many standard calculational methods become inapplicable. Another complication is the presence of light-cone propagators due to the use of light-cone gauge (or alternatively the presence of Wilson lines).
For the anti-collinear case, the arguments are completely analogous. Relabeling p ∼ n ↔p ∼n, one finds the same expressions as for the collinear cases, the only change is the appearance of the analytic regulator which now enters in eq. (3.1) as (ν/n · l i ) α . Using this relabeling, in the following we can discuss the collinear and anti-collinear cases in parallel.
Discussing the individual contributions up to NNLO, we first observe that for n r = 0 emitted partons, k = 0 and the x ⊥ dependence is lost. Setting α = 0, no scale dependence remains in these cases and dimensionless integrals are found. Hence, the bare TPDFs receive no contributions from purely virtual corrections. Then to obtain the corrections up to α 2 s to the trivial LO results B , the only cases we have to consider are the real NLO corrections as well as the double real and the virtual real NNLO corrections. Their amplitudes correspond to the diagrams in figures 1 and 2 with appropriate placement of partons as well as diagrams obtained from shrinking individual lines to points. The full contributions are obtained from appropriate combinations of these amplitudes with their Hermitian conjugates. Sums and averages over color and spin of external partons are understood. Therein the factors in eqs. (2.4, 2.5, 2.7) contracting the two gauge invariant fields lead to the factorn αβ /2 if they are (anti)-quarks, and −zn·p g ⊥µν if they are gluons.
In this sense, the NLO contributions correspond to the square of diagram 1(a). The two different 1-loop amplitude topologies with unspecified partons are depicted in figure 1(b,c). For the virtual-real contribution, these diagrams or their versions with a shrinked propagator are combined with the NLO diagram of figure 1(a). The double real diagrams without specified partons are given in figure 2. For all three amplitude topologies one propagator carries momentum p − k. The other momentum is either p − l, p − k + l or k depending on the amplitude topology. By shrinking the propagator with the second momentum, we receive the same additional amplitude subtopology from all of them. For the double real NNLO contribution, these diagrams are combined with each other. We use QGRAF [41] to generate the amplitudes and FORM [42] to manipulate them.
The NLO contributions can be solved in closed form. Having used the δ distributions, the only integral required is (3. 2) The corresponding results are given in eq. (D.12). Using appropriate parametrizations, we will identify this integral as subset of the integrals of the two NNLO cases -the virtual-real and double real corrections, which we discuss in the following two subsections. Expressing the bare coupling constant via eq. (2.14) by the renormalized one, introduces powers of the MS factor and an additional NNLO contribution stemming from the NLO contribution multiplied by the α 1 s term in the renormalization factor

Virtual-real contribution
The calculation of the virtual-real diagrams is straightforward. We first perform the integrals over the loop momenta. Using partial fraction decomposition and shift of momentum, we can reduce the scalar loops integrals to two generic types: where q = p − k and in all the propagators an imaginary part of −iδ is implicit. These integrals can be calculated using standard techniques. Taking I VR   1 as an example, we first use a Feynman parameterization to combine the propagators and perform the integration  over the loop momentum. We are then left with a multi-dimensional integral over the Feynman parameters: . The remaining integrals are not difficult to carry out and the results can be written in closed form in terms of hypergeometric functions. The final forms of the two integrals are Having performed the loop integrals, the remaining integrals over k are similar to those at NLO and can be readily evaluated using eq. (3.2).

Double real contribution
For the combined double real emission diagrams we find integrals of the form where then·k integral has already been performed using δ(k z ) such thatn·k = (1 − z)n·p .
In the argument of the squared amplitude |M | 2 we have listed all possible scalar products that can appear. We introduce a variable change y = k 2 T /(n·kn·k), and the n·k integral becomes To evaluate the l integral, we boost into the rest frame of k, such that the vectors can be parameterized as and the scalar products are given bȳ We also define D 5 = y and D 6 = 1 − y. From the above equations, we see that whenever D 8 and D 9 both appear, we can use a partial fraction decomposition to get rid of one of them. We can also use partial fraction decompositions for the pairs {D 1 , D 3 } and {D 2 , D 4 }. However, it is not always possible to get rid of these due to the analytic regulator.
Inserting the above parameterizations into eq. (3.8), we see that the k T dependence is power-like, and the k T integral can be easily performed using eq. (3.2). Performing also the l 0 and | l| integrals using the delta functions, we finally arrive at integrals of the form where {a i } = {a 1 , a 2 , a 3 , a 4 , a 5 , a 6 , a 7 , a 8 , a 9 } is the collection of the powers of denominators. Either {a 1 , a 3 } or {a 2 , a 4 , a 5 } will contain the analytic regulator α, and in general we will calculate the integrals as a power series in α and .
There are two situations where we need to keep the α regulator in the integrals. One is if the integral itself is divergent for α → 0. Another is if the integral is multiplied by the first term in the expansion of In the latter case the expansion of the integral at z = 1 is needed to α 1 . It proves useful to distinguish the two cases a 8,9 = 0 and a 8,9 = 0. For all integrals with a 8,9 = 0 that we encountered, neither of the two conditions above apply, and therefore we can always drop the α regulator in them. For integrals with a 8,9 = 0, we can use the freedom of parameterizing n andn to exchange a 1 ↔ a 2 and a 3 ↔ a 4 , and always bring the α dependence to a 1 and a 3 . For both cases, we can then use a partial fraction decomposition and the symmetry between l and k − l to reduce a 4 to 0. In the end, what we need to calculate are: I RR (a 1 , a 2 , a 3 , 0, a 5 , a 6 , a 7 , 0, 0) with α in a 1 , a 3 and possibly also in a 5 as well as I RR (a 1 , a 2 , a 3 , 0, a 5 , a 6 , a 7 , a 8 , 0) and I RR (a 1 , a 2 , a 3 , 0, a 5 , a 6 , a 7 , 0, a 9 ) without α regulator.
The corresponding calculations are further outlined in appendix A. The solutions to most of the relevant integrals can then be obtained straightforwardly, while the remaining solutions are listed in appendix B.

Results
Combining the contributions to the NNLO result (D.14) expanded up to the finite terms in the regulators α and , carrying out the refactorization (2.11), identifying the matching kernels (2.16) and renormalizing them and the anomaly coefficients (2.18, 2.13), we obtain the final results. In the FORM module [43] associated with this article, we provide the full set of results in digital form. Here, we present only the parts which are free of scale logarithms and obtained for µ = µ x ≡ 2e −γ E x T . These are F (2) iī (L ⊥ = 0) and I (2) i/j (z, L ⊥ = 0). The corresponding expressions at µ = µ x , containing powers of L ⊥ , can straightforwardly be obtained from these expressions as explained in section D.2.
The NNLO anomaly coefficients result in accordance to [21] into The NNLO matching kernels are expressed in terms of harmonic polylogarithms H an ≡ H an (z) introduced in [44], ζ values and functionsp ij related to the lowest order DGLAP splitting kernels P . The gluon-to-gluon kernel is given by The quark-to-gluon kernel reads while the gluon-to-quark kernel is obtained as The matching kernel of a quark evolving to a quark of the same flavor is given by +p qq (z) 8H 0,1,0 + 4H 0,1,1 − 4H 1,0,0 + 8H 1,0,1 + 8H 1,1,0 + 3H 0,0 + 8H 0 + 24ζ 3 + 2(1 + z)H 0,0,0 + (3 + 7z)H 0,0 + 4(1 − z)H 0,1 For a quark evolving to a quark (or anti-quark) of different flavor, it reads instead while for a quark evolving to an anti-quark of the same flavor it is obtained as q /q (z, 0) . In a slightly different notation, we reported these results already in [36,45]. All other splitting kernels I (2) i/j are related by charge conjugation or flavor symmetry to these results. The charge conjugation symmetry implies the equality Iī / = I i/j and to respect the flavor symmetry we introduced above only a quark q of unspecified flavor and a quark q of different flavor. Moreover, the relation Iq /q = I q /q holds up to NNLO. As a check of our results, we also considered other combinations of partons and found agreement.

Relation to q T -resummation in the Collins-Soper framework
In [34,35] the hard-collinear coefficient functions for Drell-Yan and Higgs production were calculated within the framework established in [17,20] up to NNLO+NNLL. A processindependent formulation of this framework for q T -resummation is derived in detail in [33]. The same framework is also used as construction principle for a subtraction scheme [46] for fixed-order NNLO calculations.
Our results are obtained in a completely different approach to q T -resummation, based on a different factorization into individual contributions. Consequently, the building blocks of the resummed cross section can not be compared one-by-one between the approaches, since they are scheme-dependent. Both approaches must agree on the scheme-independent expression for the resummed cross section, as we will verify explicitly below.
In eq. (6) of [33], the differential cross section is expressed in a factorized and resummed form, which contains the hard factor [H F C 1 C 2 ]. For qq initiated processes the latter is given by the product for gg initiated processes it is given by the following contraction of tensors [38] where the dependence on Laplace-space variables and coupling constants has been omitted for clarity. In our language, H F corresponds to the square of the Wilson coefficient C F , which arises on matching QCD on the effective field theory. The process-independent factors C 1 and C 2 correspond to the collinear and anti-collinear matching kernels I, respectively. However, there is no one-to-one correspondence, since these expressions are scheme-dependent.
Nevertheless, their product related to the physical cross section by eq. (6) of [33] and our eqs. (2.19, 2.20) is well defined after carrying out the convolution in the momentum fractions z 1 and z 2 . In [34,35] this is given by for Drell-Yan and Higgs production respectively. From the process-dependent Wilson coefficients and our results on the process-independent matching kernels, we can determine these H functions as

5)
H H gg←jk z, α s , log where each function is evaluated at a value of the renormalization scale for which no large logarithms arise, which is the invariant mass of the produced final state and µ x = 2e −γ E /x T , respectively. In the second line, I µν g/j is the gluon matching tensor which is related to B µν g/N in eq. (2.5) in a completely analogous way as I g/j is related to B g/N [26,45]. It can be decomposed into the two independent components I g/j and I g/j analogously to eq. (2.6).
is the hard tensor. For Higgs production it has the explicit form with the Wilson coefficients arising on first integrating out the top quark and then matching to SCET. To determine eq. (4.6) to NNLO, the NLO results of I g/j are required which we calculated finding results in accordance with [26]. The resulting expressions for the H coefficients are found in full agreement with the results in [34,35], and constitute a fully independent validation of them in a completely different calculational approach.

Further checks
Below we describe further observations and checks confirming our results for the matching kernels I (n) and anomaly coefficients F (n) with n ≤ 2 . We first observe that these functions depend only on the scale logarithm L ⊥ and the momentum fraction z. As required by consistency, no dependence on the analytic regulator α or the associated scale ν remained, but they canceled in eq. (2.11), where moreover all dependence on the hard scale q 2 had been refactorized from the resulting functions. This not only confirms our results but also the consistency of the whole framework and explicitly demonstrates the applicability of the analytic regulator of [22] in high order calculations.
Moreover, in our results no poles in the dimensional regulator remained, but they could consistently be removed by renormalization (2.13, 2.18), where the exact renormalization factors had been implied already by their RGEs in terms of known functions and are listed in section D.1. We also explicitly confirmed that F iī (L ⊥ , α s ) and I i/j (z, L ⊥ , α s ) themselves obey the RGEs (2.23, 2.24) and that their L ⊥ dependent terms can be reconstructed through the relations in appendix D.2 from the results listed here and the expressions in appendix C and D.3. These points are yet another strong confirmation of our results.
Furthermore, we did not only perform the calculation in light cone gauge as described in this article, but also in Feynman gauge finding identical results. This not only serves as test to our calculation, but also explicitly demonstrates that the individual factors in our framework are gauge invariant.
In addition to that, we compared our results to literature: we could explicitly confirm the expressions for the anomaly coefficients and the NLO matching kernels as given in [21,26].

Conclusions
In this paper, we have derived perturbative QCD corrections to all parton-to-parton TPDFs at NNLO. Our calculation is based on a gauge invariant operator definition [21,26] with an analytic regulator [22]. We demonstrate for the first time that such a definition works beyond the first non-trivial order, and that it provides a fully complementary approach to q T -resummation in the CSS framework [16,17,20,33]. From our calculation, we extract the coefficient functions relevant for q T -resummation at N 3 LL accuracy. Our results can be applied to any process yielding a colorless final state, provided the NNLO virtual corrections are known. They confirm the recent structural findings in [33], while working with a completely different methodology [21,22,26] based on SCET. Combined with the work of [47], our results could also be applied to the transverse momentum resummation in tt production. For gluon-gluon initiated processes with a general spin structure, in addition to the results presented here, N 3 LL transverse momentum resummation may require the NNLO corrections to the second tensor structure of the gluon TPDFs, which we will present in a separate article. We documented our calculation in detail, and validated our results with numerous non-trivial checks, including an independent re-derivation of the secondorder contributions to the hard factors H for Drell-Yan and Higgs production that were obtained previously in [34,35]. A digital form of our results is provided in [43].

Acknowledgments
We would like to thank Matthias Neubert, Guido Bell and Massimiliano Grazzini for useful discussions. T.L. would like to thank Peking University for the hospitality during part of this work has been completed. This work was supported in part by the Schweizer Nation

A Calculation of double real integrals
In this appendix we outline the determination of the integrals I RR ({a i }) defined in eq. (3.12) which appear for the double real emission. As explained in section 3.2, we distinguish three relevant subsets of integrals.
We first consider the integrals with a 8,9 = 0. It is convenient to define integrals of the form The full integrals are then given by If one of its arguments is 0, the I RR 1 integrals can be readily calculated to be If furthermore a 7 = 0, the remaining integral over y can be carried out, and the result is We also have I RR (a 1 , 0, a 3 , 0, a 5 , a 6 , a 7 , 0, 0) = I RR 1 (a 1 , 0, a 3 ) I RR 2 (a 5 , a 6 , a 7 ) .
For more generic cases, we change variables to which allows us to rewrite the integral as ¿From this representation, it is obvious that if a 2 ≤ 0, the integrand can be expanded and written in terms of powers of u, 1 − u, v, 1 − v, y and 1 − y. The integrals over u and v then lead to some Γ functions, while the powers of y and 1 − y can be absorbed into a 5 and a 6 . The remaining y integral can then be performed with the help of eq. (A.1). For a 2 > 0, we first perform the v integral to get For each of the two integrals above, we change variable from u to the last argument of the hypergeometric function, which we call t. We then apply 2 F 1 (a, b; c; t) = (1−t) c−a−b 2 F 1 (c− a, c − b; c; t) and insert the resulting expression into eq. (A.2) to arrive at a 2 , a 3 , 0, a 5 , a 6 , a 7 , 0, 0 From here, the remaining integrals in general cannot be performed in closed form, and a series expansion in α and is required. These expansions are documented in the next appendix.
We now turn to the cases where a 8 > 0 or a 9 > 0. As mentioned above, we can always drop the analytic regulator α for these integrals. Therefore we can always reduce a 4 to 0. Following the same procedure as before, cases with a 2 ≤ 0 can be performed straightforwardly. For a 2 , a 8 > 0 we obtain a 1 , a 2 , a 3 , 0, a 5 , a 6 , a 7 , a 8 . (A.10) The main complication here is that may both appear. This prevents us from changing variable to the last argument of the hypergeometric function, since regardless of whether we substitute u or y, either D 7 or D 8 will become very complicated in terms of the new variable t. We therefore now consider specific cases. For a 7 = 0 one obtains The representation of I RR (a 1 , a 2 , a 3 , 0, a 5 , a 6 , 0, 0, a 9 ) is essentially the same as above, with D −a 8 8 replaced by D −a 9 9 . The relevant cases, where in addition a 7 > 0, are a 7 = 1, a 1 , a 3 = 0 and either a 8 or a 9 = 1. We then partial fraction decompose D 7 with D 8 or D 9 , respectively. After changing variables from either u or y to the last argument of the hypergeometric function, which we call t, and if relevant renaming y to u, one obtains I RR (0, a 2 , 0, 0, a 5 , a 6 , 1, 1, 0 and an even more involved version of this for a 9 = 1. At intermediate steps, an additional regulator is introduced in a 5 which however does not lead to poles in the final result for the integral.

B List of double real integrals
In the previous appendix we described the methods of calculating the double real integrals. Some integrals can be represented in an exact form in terms of hypergeometric functions The QCD β function is Higher order coefficients of Γ i , γ i and β can be found in [11,39,40], respectively. The DGLAP splitting functions are where the first order coefficients are with the functionsp The second order coefficients can be obtained from the results in [9,10] and we do not repeat these expressions here. The coefficients up to third order are given in [11,12]. By (. . .) + we denote the plus prescription with support on [0, 1] regulating the pole at z = 1. To express our results, we also usep ij (−z). In those cases the plus prescription is dropped.

D Further results
In this section, we collect a number of results which either appear in intermediate steps or have been given in the literature already. Due to the flavor and charge conjugation symmetry of QCD, the number of independent functions reduces and below we use a notation, where q refers to a quark of unspecified (but same) flavor and q to a quark of different flavor. Up to NNLO we moreover have Bq /q = Bq /q and corresponding relations for the other functions.

D.1 Renormalization factors
As explained in section 2, the renormalization factors are related by their RGEs to the anomalous dimensions and QCD β function listed in appendix C. Here we list the resulting expressions for the perturbative coefficients according to eq.
while the coefficients for Z F i in eq. (2.13) are given by

D.2 Dependence on scale logarithms
The L ⊥ dependence of F iī and I i/j can be recovered from their values at L ⊥ = 0 by solving the RGEs (2.23, 2.24). More explicitly, we can expand these functions in both α s and L ⊥ according to for the coefficients defined above with n, l ≥ 0. Coefficients on the right hand side with l outside the range specified in eqs. (D.5, D.6) are understood to vanish. The coefficients with values (n , l ) are thus expressed in terms of coefficients with lower values of n and l and the QCD parameters listed in appendix C. Applying these equations recursively, one can remove all terms with l > 0 on the right hand sides of these equations as can be shown easily by induction. Phrased differently, the functional dependence on L ⊥ of F iī and I i/j can be recovered from their values at L ⊥ = 0 , F The NNLO results have been given in eq. (4.1) and the terms containing L ⊥ are implied by eq. (D.7). The full results agree with [21]. The renormalized matching kernels up to NLO are g/g (z, 0) = − C A ζ 2 δ(1 − z) , g/q (z, 0) = 2C F z , q/g (z, 0) = 2T F z(2 − z) , q /q (z, 0) , I The NNLO results have been presented in section 4. The full results are available in the FORM module [43] associated with this paper.

D.4 Bare TPDFs
For the logarithms associated with the analytic regulator we identify L a = log ν n·p and L c = log νn·p x 2 T 4e −2γ E (D.11) for the anti-collinear and collinear region, respectively. Then the exact NLO results for the bare TPDFs of the collinear and anti-collinear region are given by B (1) i/j (z, x 2 T , µ, ν) = e αLc+ L ⊥ e −( +2α)γ E Γ(− − α) Γ(1 + α) (1 − z) α f f (2,2) q/q (z, s) = C F C A δ(1 − z) − q /q (z, s) , f (2,2) q/q (z, s) = O(α 0 ) . (D.17) Using the results listed above, it is a straightforward exercise to confirm the cancellation of all poles in the analytic regulator up to NNLO in α s on the left hand side of eq. (2.11) as well as the cancellation of the related scale ν and the associated generation of the hard scale q 2 ∼n·p n·p by the difference of the two logarithms in eq. (D.11).