Rapidity renormalized TMD soft and beam functions at two loops

We compute the transverse momentum dependent (TMD) soft function for the production of a color-neutral final state at the LHC within the rapidity renormalization group (RRG) framework to next-to-next-to-leading order (NNLO). We use this result to extract the universal renormalized TMD beam functions (aka TMDPDFs) in the same scheme and at the same order from known results in another scheme. We derive recurrence relations for the logarithmic structure of the soft and beam functions, which we use to cross check our calculation. We also explicitly confirm the non-Abelian exponentiation of the TMD soft function in the RRG framework at two loops. Our results provide the ingredients for resummed predictions of p⊥-differential cross sections at NNLL′ in the RRG formalism. The RRG provides a systematic framework to resum large (rapidity) logarithms through (R)RG evolution and assess the associated perturbative uncertainties.


Introduction
The transverse momentum (p ⊥ ) distribution of a heavy color-neutral final state L in the process sufficiently high order provide an accurate description of the p ⊥ -spectrum. In the region where p 2 ⊥ M 2 , the distribution is peaked and the resummation of (Sudakov) logarithms ∝ ln n (p 2 ⊥ /M 2 ) is crucial for reliable theoretical predictions. The traditional approach to resum such logarithms in QCD is based on the factorization theorem devised by Collins, Soper, and Sterman (CSS) in ref. [1]. In this framework the resummation has reached the next-to-next-to-leading-logarithmic (NNLL) level [2][3][4][5][6][7][8][9]. In recent years, the factorization for p ⊥ -differential cross sections has been revisited [10][11][12] using Soft Collinear Effective Theory (SCET) [13][14][15][16][17][18]. Also a number of SCET predictions with NNLL resummation for Drell-Yan-like processes are available by now [19][20][21][22]. The SCET approach typically allows a direct control of the perturbative uncertainties by independent variations of the various unphysical factorization/renormalization scales in the problem, both in momentum and position (impact parameter) space. In principle also subleading power corrections of relative O(p 2 ⊥ /M 2 ) to the factorization theorem are systematically calculable in SCET.
Throughout this paper we will focus on the perturbatively accessible resummation region of the p ⊥ -spectrum, i.e. where Λ 2 QCD p 2 ⊥ M 2 . The leading order SCET factorization theorem generically has the schematic form dσ dp 2 (1. 2) The process-dependent (but p ⊥ -independent) hard function H contains all purely virtual QCD corrections to the partonic production of L. The hard function multiplies a convolution of one TMD soft function S(p ⊥ ) and two TMD beam functions B a,b (p ⊥ ) -one for each beam. The latter are also called TMD parton distribution functions (TMDPDFs) in the literature. 1 The beam and soft functions describe the collinear and soft radiation off the incoming hard partons, respectively. The beam functions are process-independent (universal). The soft function (for a given partonic channel) is universal within the class of processes in eq. (1.1), i.e. independent of L. Up to two loops it is the same for Drell-Yan like processes, e + e − → dijets and deep inelastic scattering [30]. It is well known that the factorization of the p ⊥ -differential cross section into separate soft and collinear sectors gives rise to so-called rapidity divergences, see e.g. ref. [10,31]. They arise from loop/phase space (PS) momentum regions, where one light-cone component, say k + , approaches infinity or zero, while the product k + k − ∼ k 2 is kept fixed. This typically happens in soft and beam functions, when the corresponding soft and collinear modes relevant to describe the process have virtualities of the same order of magnitude. The production of a final state L with measured transverse momentum p ⊥ , which recoils against the soft as well as the collinear real radiation, is a classic example for such a process. The appropriate effective field theory (EFT) framework for this kinematic situation is SCET II . By contrast, the virtuality scales of (ultra)soft and collinear degrees of freedom 1 The name TMDPDF goes back to the traditional QCD factorization formalism. We prefer to call it the TMD beam function, as in SCET it is conceptually on the same footing as the virtuality-dependent [23][24][25][26] or fully-differential (unintegrated) beam function [27][28][29] (although the latter are SCETI objects).

JHEP03(2016)168
JHEP03(2016)168 of the observable it will be beneficial to eventually have precise theoretical predictions in different resummation frameworks. This will provide important cross checks among the different approaches and also contribute to a better theoretical error estimate.
The outline of this paper is as follows. In section 2, we will briefly review the SCET definition of TMD soft and beam functions, their regularization and renormalization following refs. [10,33]. Their all order logarithmic structure is discussed in section 3. In section 4, we calculate the TMD soft function through NNLO and explicitly check its exponentiation properties and logarithmic structure. We discuss the relation of our results to results in the literature -most importantly to refs. [41,42] -in section 5. This relation together with our result for the soft function allows us to extract the NNLO TMD beam functions in the RRG scheme of refs. [10,33]. The expressions for the latter are presented in section 6. Like all our results they are given in both momentum and position space. We conclude in section 7 and collect some additional technical details in the appendices.

SCET framework
In this paper we consider the QCD corrections to the soft and beam functions in the SCET factorization theorem for the p ⊥ -differential cross section of the process in eq. (1.1), where X now represents the inclusive hadronic final state.
We work in light-cone coordinates w.r.t. the light-like vectors n µ = (1, 0, 0, 1) and n µ = (1, 0, 0, −1) along the beam (z-) axis. An arbitrary four-momentum vector is then expressed as where p + = n · p and p − =n · p, or short p µ = (p + , p − , p ⊥ ). Euclidean (two-)vectors are denoted in boldface, p 2 ⊥ = −p 2 ⊥ . The partons initiating the hard interaction carry the (large) longitudinal momentum ω ± = Qe ±Y , where the hard scale Q = M [1 + O(p 2 ⊥ /M 2 )] and Y is the rapidity of the color singlet state L. The total center of mass energy of the colliding protons is √ s. For quarks as the incoming partons we write the factorized differential cross section as 5 The index b indicates bare quantities. Similarly for incoming gluons we have 6 The 'bare' hard function H b in SCET can be expressed in terms of (the absolute square of) a hard current matching coefficient and the combined counterterms of beam and soft functions. It is process dependent and will not be discussed any further in this paper. The factorized cross section in terms of renormalized hard, beam and soft functions takes the same form as eqs. (2.2) and (2.3). We simply replace bare by renormalized quantities. The bare soft functions are vacuum matrix elements of soft Wilson lines and are defined in momentum space as where the T andT are time and anti-time ordering operators and P is the SCET label momentum operator [15]. The soft Wilson lines are exponentials of soft gauge fields. In momentum space we compactly write them as where the sum is over all permutations of the momenta associated with the soft gluon modes [16]. We leave it implicit in the notation of the Wilson lines that the gauge fields are in the fundamental representation of SU(3) in the case of incoming quarks, or in the adjoint representation for incoming gluons. For later convenience, we have introduced an index i ∈ {q, g} to denote the soft function and all the anomalous dimensions in the fundamental (i = q) or adjoint (i = g) representation, and the generic soft function S i b (p ⊥ ). Following refs. [23,24], the bare n-collinear beam functions are defined as spin-averaged forward proton matrix elements of n-collinear quark or gluon field operators (z = ω − / √ s): 0) |p n , (2.8) and analogous for then direction with n ↔n, ω − ↔ ω + as well as the antiquark beam function. The ω ± is fixed by the delta function to the large light-cone momentum of the parton entering the hard interaction. We will often write ω ± , when the discussion is equally valid for both directions. One should keep in mind though that different collinear directions have different large components and ω + ω − = Q 2 . For more details on the relevant SCET definitions, as e.g. for the definitions of the (collinear) gauge invariant quark fields χ n and gluon fields B µ n , we refer to refs. [10,16,24]. In order to separate non-perturbative contributions we write the beam functions as a convolution of perturbatively calculable matching coefficients I ij and the usual collinear PDFs [23,44,45]. At leading order in Λ 2 QCD /p 2 ⊥ we have for the renormalized TMD beam functions

JHEP03(2016)168
Here the sum is over all partons i, i.e. quarks and anti-quarks of all light flavors as well as gluons, z denotes the large light-cone momentum fraction w.r.t. the corresponding proton momentum (P ± = √ s), i.e. z = ω ± / √ s, and µ is the UV renormalization scale. We have introduced the symbol ⊗ z for the (Mellin-) convolution (2.11) Our choice of rapidity regulator in addition introduces the rapidity renormalization scale ν as described below.

Position space
In alternative to the formulation in momentum space, the relations and results can be formulated in position (impact parameter, b) space. They are related by Fourier transformation in two dimensions:X where the ellipses denote any additional arguments. Following a common abuse of notation, we will from now on drop the tilde and use the same symbols for the b and p ⊥ space functions as the distinction is clear through their arguments.

Rapidity regulator
The beam and soft functions exhibit rapidity divergences that are not regulated by dimensional regularization. To regulate these divergences, we will use the rapidity regulator introduced in refs. [10,33] in addition to dimensional regularization. This regulator has been devised for covariant gauges and is defined by the following modification of the soft Wilson line, and a related modification of the collinear Wilson line required for the regularization of the TMD beam functions. In order to preserve gauge invariance and exponentiation, at higher loops we must not regulate the momentum of every gluon attached to the Wilson line individually. Equation (2.14) instead should be interpreted to contribute a factor of ∝ w ν η/2 |2P g3 | −η/2 for every web attached to the Wilson line as a subdiagram. The (label) momentum operator P g3 = (P − g − P + g )/2 then picks out the third component of the total (group) three-momentum flowing through the respective web. Here, a web is understood as the maximally non-Abelian piece of a Feynman diagram, which cannot be disconnected by cutting two eikonal lines. After exponentiation, the exponent of the soft function is a sum of webs only. For a precise definition of a web in the context of non-Abelian exponentiation we refer e.g. to ref. [46].

JHEP03(2016)168
For our two-loop calculation of the TMD soft function the implementation in practice works as follows: the parts of the diagrams, which are not fixed by the non-Abelian exponentiation theorem [47,48], i.e. the ones with color factors other than C 2 F , we multiply with the regularization factor (2.15) Here we identify the group momentum P g of the web with the total soft momentum crossing the final state cut. In the C 2 F pieces of the diagrams each gluon exchange individually represents a web and is regulated accordingly, see section 4.2.3. Purely virtual diagrams vanish as scaleless integrals.
Our definition of the group momentum as the total momentum passing through the web across the cut is somewhat different to what is suggested in appendix A of ref. [10]. There P g is defined to give the web's total momentum transfer from the Sn to the S n Wilson lines. The arguments given in ref. [10] proving that gauge invariance and exponentiation are preserved by the η-regulator are however equally valid for our definition of the web's group momentum. Our choice allows us to recycle previously known results for the soft one-loop current in the calculation of the real-virtual diagrams in section 4.2.1. As for the analytic regulator of ref. [34] with our definition gauge invariance is manifest, because effectively only the phase space integration is modified. Restricting P g to act on real radiation only also allows for a sensible definition of a jet function.
The η-regulator works similarly to dimensional regularization: the bookkeeping parameter w is the analogue of a renormalized coupling that will be set to one (for all values of ν) in physical (rapidity finite) results, while ν is a renormalization scale similar to µ in dimensional regularization. The rapidity divergences manifest themselves as poles at η = 0.
With this regulator, the soft and beam functions exhibit 1/η poles as well the ordinary 1/ poles from dimensional regularization and it is crucial to send η → 0 before → 0 with η/ n → 0 for all n [10]. This is because at the bare level we first have to combine soft and beam functions such that the η divergences cancel. We must then set η → 0, before we can cancel their overall 1/ poles with the ones from the 'bare' hard function, which is η independent by construction.
Requiring the bare results to be ν independent leads to RRGEs for the renormalized beam and soft functions, that enable us to resum the rapidity logarithms [10]. For this to work obviously w 2 ν η in eq. (2.15) must be ν-independent, and as a consequence The similarity to dimensional regularization also implies that with the η-regulator soft zero-bin contributions [36] in the (explicit) calculation of the beam functions amount to scaleless integrals and vanish. This is also true for the analytic regulator proposed in ref. [34] but does not hold e.g. for the delta regulator of ref. [35].
Last but not least we would like to emphasize that besides ref. [38], our computation of the NNLO TMD soft function is the only two-loop calculation and the first with an explicit demonstration of the non-Abelian exponentiation property, using the η-regulator.

JHEP03(2016)168
It can therefore be considered a valuable consistency and practicability check of the rapidity regularization method beyond one loop.

Transverse momentum in dimensional regularization
There are different schemes how to apply dimensional regularization when measuring the transverse momentum of the soft and collinear radiation, see e.g. the discussion in ref. [27]. For the calculation of the soft function as presented in this paper, we will adopt the CDR 2 scheme, where the transverse momentum, p ⊥ , of the soft and beam functions is a twodimensional vector, while all transverse loop momenta are (d−2)-dimensional vectors. 7 The delta functions to measure the transverse momentum in eqs. Due to the isotropic nature of the soft radiation, the soft function S(p ⊥ ) is independent of the orientation of the vector p ⊥ in the transverse plane. 8 We can therefore replace the delta function in eqs. (2.4) and (2.5) as where the factor 1/π on the right hand side (r.h.s. ) compensates for the contribution from the azimuthal integration, when taking the implicit sum over all soft final states, whose total transverse momentum is collected by the (label) momentum operator P ⊥ . The replacement in eq. (2.17) gives rise to a subtlety concerning the dimensional regularization of the soft function. The measurement functions on both sides of eq. (2.17) only yield identical results for the bare soft function in d dimensions if the momentum operator P ⊥ on the r.h.s. is interpreted to return only two transverse momentum components, while the residual (d − 4) components of the total soft transverse momentum crossing the final state cut remain unrestricted.
The definition of the soft function using the r.h.s. of eq. (2.17) however also admits a different scheme for the dimensional regularization. In particular we can interpret P 2 ⊥ as the square of the full (d − 2) dimensional transverse momentum of the soft radiation. In this case the scalar delta function corresponds to the vector-like measurement where also p ⊥ is promoted to (d − 2) dimensions in the definition of the bare soft functions. The latter will now however differ from the ones defined in CDR 2 , i.e. eqs. (2.4) and (2.5). Both schemes are nevertheless equally viable for the calculation of the soft function. To see this, we note that the 1/ n poles in the bare soft function originate from virtual 7 Here and in the following we refer also to the momenta of the soft partons that cross the final state cut as loop momenta. 8 An analogous argument holds for the quark beam function, but not for the gluon beam function, eq. (2.10), which has a nontrivial Lorentz (tensor) structure.

JHEP03(2016)168
ultraviolet (UV) divergences, which evade the cancellation by UV divergent real corrections because of the p ⊥ constraint. As a consequence, all NLO 1/ n divergences (n = 1, 2) are proportional to δ(p 2 ⊥ ), cf. eq. (4.35). Since the divergent one-loop virtual corrections (where P ⊥ ≡ 0) are not affected by the measurement function we find the same bare result up to O( ) terms in both schemes discussed above.
At higher orders, starting from NNLO, there are mixed real-virtual contributions. The bare expressions for the p ⊥ -dependent soft function will therefore in general depend on the scheme. We indeed find different subleading terms in 1/ . Upon renormalization these differences however cancel out and we obtain scheme independent results for the renormalized soft functions and their anomalous dimensions. We have explicitly verified this by performing the complete calculation of the two-loop soft function in both schemes, CDR 2 and using the measurement according to eq. (2.18). This also serves as a cross check of our final NNLO results.

Renormalization and logarithmic structure
In this section we briefly review the (rapidity) renormalization of the beam and soft functions according to ref. [10] and predict their logarithmic structure based on the (R)RGEs.
The renormalized soft and beam functions and their respective renormalization factors, Z i S,B , that absorb all UV and rapidity divergences are defined through where we have introduced the generic notation with the index i ∈ {q, g} also for the beam function. Lorenz indices for the gluon beam function (i = g) are understood. The symbol ⊗ ⊥ denotes the convolution in two-dimensional transverse momentum space. From the above renormalization factor Z i S we can derive the anomalous dimensions for the soft RGE and RRGE in momentum space, and analogous for the beam function anomalous dimensions. The µ anomalous dimensions actually only have a trivial p ⊥ dependence [10] as elaborated on below eq. (3.16). To make this explicit, we can write

JHEP03(2016)168
and analogous for γ i B,µ . For convenience when expressing functions as power series in α s , we define We expand the renormalized soft function and its renormalization factor as

Recurrence relations
The (R)RG structure of the soft and beam functions is conveniently discussed in impact parameter (b) space. Upon the Fourier transformation in eq. (2.12), convolutions of the type in eq. (3.3) turn into ordinary products. For notational convenience we define All b space expressions given below can be straightforwardly transformed back to p ⊥ space. 9 The two-dimensional Fourier transform of powers of the logarithm in eq. (3.10) can be expressed in terms of plus distributions L T n (p ⊥ , µ) defined in appendix A, where we also give corresponding translation tables.
The RGE and RRGE of the soft and beam functions in b space read 14) The FO cross section must be independent of µ and ν order by order in α s . The anomalous dimensions must therefore fulfill the consistency relations where γ i H,µ denotes the anomalous dimension of the hard function. Besides µ the latter only depends on the hard scale Q. In eqs. (3.11) and (3.13) we therefore anticipated that, due to eq. (3.15) and the fact that in SCET II the soft and collinear sectors are located JHEP03(2016)168 on the same invariant mass hyperbola, the µ anomalous dimensions γ i S,µ and γ i B,µ are also independent of b [10]. Given the well-known Sudakov form of γ i H,µ we can thus infer the structure of γ i S,µ and γ i B,µ from eq. (3.15), The universal cusp and the soft and beam non-cusp anomalous dimensions only depend on µ through a s . The coefficients Γ i n are known up to three loops [49,50]. We From eqs. (3.16), (3.17) and because derivatives commute, we have Given the simple relation in eq. (3.16), we define a common ν anomalous dimension Equation (3.21) allows us to determine the structure of γ i ν . Expanding it we find for the n-th coefficient where the β i denote the usual QCD β-function coefficients with β 0 = (11C A − 4T F n f )/3. Upon integration we thus obtain the recurrence relation Analogous to this derivation, the RGE and RRGE for the soft function in eqs. (3.11) and (3.12) imply differential equations for its expansion coefficients. Integrating them, we arrive at 10 We can choose the scales µ S and ν S such that these integration constants are free of large logarithms and even become b-independent numbers e.g. by setting them to their canonical values, given below in eq. (3.38). Also note that for the canonical value of µ S we have ln µ 2 /µ 2 S = L b . In this way all RG and RRG logarithms (∝ L n b ) are generated by the recursion. We can therefore use eqs. (3.24) and (3.25) to predict the logarithmic structure of the ν anomalous dimension and the soft function, respectively. Through NNLO we find perfect agreement with the results of our explicit two-loop calculation of the soft function in section 4.4.
The structure of the beam function matching coefficients I ij and J gj in eqs. (2.9) and (2.10) can analogously be derived from the RGE and RRGE for the beam functions, eqs. (3.13) and (3.14). To do so, we first note that the PDFs obey the DGLAP equations where the splitting functions can be expressed as power series in a s , ij (z) functions in the notation used here as well as convolutions among them, can e.g. be found in the appendices of refs. [25,26]. From eqs. (3.13), (3.14) and eq. (3.27) we obtain the following RGE and RRGE for the matching coefficients: We expand Just as for S i , we obtain a recurrence relation for the coefficients I ij by integrating the differential equations for these coefficients following from eqs. (3.29) and (3.30):

JHEP03(2016)168
with the integration constant gj (z) become b independent and are pure functions of z. We will use the recursion relations to derive the logarithmic structure of all beam function matching kernels through NNLO. Note that to this order all ingredients of eq. (3.32) but the I (n) ij (z) are known, once we have determined the soft (non-cusp) and ν anomalous dimensions from the direct calculation in section 4.
By matching our expressions onto the beam functions calculated in another (rapidity regularization) scheme, we are able to fix the missing coefficients I

Resummation
In addition to using the RGE and RRGE of S i and B i to obtain the logarithmic structure at each FO in α s , we can also use them to resum the logarithms to all orders in α s . The RGE and RRGE for X = S i , B i both have the form with the generic solution is independent of ν, the solutions of the RRGEs simplify to where only a single logarithm of ν appears in the exponent of the resummation factors. Accordingly, the maximum power of the rapidity logarithms at O(α n s ) in the cross section is n. This is in contrast to the series of double (Sudakov) logarithms generated by the µ anomalous dimensions which themselves depend on µ and originate from the interplay of soft and collinear singularities. Since the ordinary PDF f i is independent of ν, the analogue of eq. (3.37) holds for the matching kernels I ij in eqs.

Calculation of the soft function
In this section we compute the soft function in eq. (2.4) to NNLO. We perform the calculation in momentum space and with the Wilson lines in the fundamental color representation.
Our results are straightforwardly generalized to an arbitrary color representation due to Casimir scaling and translated to impact parameter space. At two loops, we separately treat the contributions with non-C 2 F (i.e. C A C F and C F T F n f ) and C 2 F color factor. For the latter we explicitly verify non-Abelian exponentiation in the presence of the η-regulator, see section 4.2.3. Throughout this paper, we use the MS (UV) renormalization scheme and Feynman gauge for the calculations. The MS relation between bare and renormalized couplings is The bare results presented in this section are obtained for the measurement eq. (2.17) in the CDR 2 scheme. The renormalized results are independent of this scheme, see discussion in section 2.
3. An extended discussion of several aspects is given in ref. [51].

NLO
The NLO soft function in the RRG formalism has been calculated in the adjoint representation in ref. [10]. The relevant one-loop diagrams are shown in figure 1. In agreement

JHEP03(2016)168
Color Factor Diagrams Table 1. Table of diagrams (2nd column) and the color factors they involve (1st column). The corresponding diagrams are displayed in figure 1, 2 and 3. After renormalization, the one-loop diagram contributes to both C F C A and C F T F n f through the β 0 term in eq. (4.1).
with ref. [10] we obtain for the bare result in the fundamental representation. In the second step we performed the integration and expressed g b through the renormalized coupling a s using eq. (4.1). This directly reveals the NNLO contribution from the one-loop calculation, which is proportional to β 0 . In section 4.2.3 we will show that the square of eq. (4.2) (in a convolutional sense) determines the C 2 F part of the NNLO result.

NNLO
The two-loop diagrams contributing at NNLO can be divided into two categories: realvirtual (single real emission) and real-real (double real emission) graphs. They are displayed in figures 2 and 3, respectively, and their color structure can be read off table 1. As previously stated, we ignore purely virtual diagrams since they vanish for our choice of regulators.

Single real emission
We use the known result for the one-loop soft current [52] to calculate the C F C A part of the sum of real-virtual diagrams (R) in figure 2. In the end, the calculation differs from the one-loop calculation in eq. (4.2) only by an overall factor in the integrand:

Double real emissions
The two-loop diagrams for the TMD soft function are actually the same as for the dijet (e + e − hemisphere) soft functions in refs. [53][54][55]. Only the measurement (delta) function differs. In ref. [54] the explicit expressions for the unintegrated amplitudes of the double real emission diagrams shown in figure 3 are given in Feynman gauge. We can readily use these amplitudes for our calculation. We just have to implement our measurement and rapidity regulator according to sections 2.3 and 2.2 and finally integrate as described below. Similar to ref. [54], we divide the double real emission diagrams into five groups as displayed in figure 3. With the right choice of integration variables, all the real-real diagrams can be expressed through the same family of two-loop integrals. They can therefore be calculated with the same method, which we now describe in detail.
The two integrations are over the d-dimensional real particle momenta k 1 and k 2 . With our regulator and measurement, it turns out to be convenient to use the integration variables = k 1 and k = k 1 + k 2 . The phase space integration then takes the generic form where the two δ + (k 2 i ) ≡ θ(k 0 i )δ(k 2 i ) originate from cutting the propagators of the two real particles. The expressions for the unintegrated amplitudes A i for i = {I, T , G, H, Q} can directly be taken from appendix B of ref. [54].
Considering the k integral first, we find a set of convenient integration variables with the d − 2 dimensional k ⊥ and A major advantage of this variable choice is that only the regulator factor |v| −η depends JHEP03(2016)168 on v. Hence the integration over k's light-cone components reduces to To perform the integral in eq. (4.4), we boost to the rest frame of k, where the involved vectors can be parametrized as angles θ i . All angle dependent scalar products that can appear in the amplitudes are then Note that none of these factors depends on the other d−3 spherical coordinates. Hence, the corresponding angular integral can be performed trivially. It is also convenient to define D 5 = y , (4.16) Implementing the above parameterizations in eq. (4.4), the k ⊥ dependence is power-like and we can easily perform its integral: The 0 and | | integrals can be done with the delta functions in eq. (4.4) and we finally arrive at integrals of the form I RR (a 1 , a 2 , a 3 , a 4 , a 5 , a 6 ) ≡ 1 0 dy 1 2π This integral family is studied in refs. [42,56] and we can use their results for the relevant integrals in our case:

JHEP03(2016)168
Note that eq. (4.22) is not valid for 1 ≤ a 1 + a 2 ≤ 2. 12 In this case we must instead use Now, matching the amplitudes onto the I RR functions in eqs. (4.20)-(4.23) gives the final results for the double real diagram groups defined in figure 3: (1, 1, 0 The common prefactor is Of the eight distinct I RR functions in eqs. (4.24)-(4.27), five reduce to combinations of ordinary gamma functions and are thus straightforwardly expanded in η and . The remaining three I RR functions contain nontrivial hypergeometric functions which are expanded using the HypExp package [57,58]. We list the corresponding results in appendix B.

Non-Abelian exponentiation
In the calculation of the soft function we have so far neglected any C 2 F pieces in the diagrams. Instead of computing it explicitly, we can obtain the C 2 F part from the non-Abelian exponentiation theorem [47,48]. According to this theorem, the soft function takes the form of an exponential of webs in position space. We therefore write the bare soft function as where the factor of 1/(2π) 2 is a normalization factor. We can identify the O(a s ) term in the exponent with the Fourier transform of the NLO bare TMD soft function computed in eq. (4.2). 12 The 3F2 hypergeometric function arises from an integral over a 2F1. For the values 1 ≤ a1 + a2 ≤ 2 one integrates over a singularity and the 3F2 is therefore not well defined. Using an Euler transformation on the 2F1, the singular part can be extracted to arrive at a well defined integral representation of a 3F2.

JHEP03(2016)168
The non-Abelian exponentiation theorem states that s q(1) (b) is purely C F , while s q(2) (b) only contains the color factors C F C A and C F T F n f . Hence, the total C 2 F part of the soft function is equal to the [s q(1) (b)] 2 term in eq. (4.29). This allows us to obtain the C 2 F piece by simply squaring the NLO result. Fourier transformation to momentum space turns the square into a convolution, and the full C 2 F part of the soft function is It is straightforward to see that the total C 2 F part comes exclusively from the I diagrams in figure 3. With a rapidity regulator for each real gluon momentum, the C 2 F piece of I takes the form Simply performing the convolution integral in eq. (4.30) we find Consequently the (non-web) C 2 F piece of the real-virtual diagrams (b) and (c) in figure 2 must be zero. In fact, one can see this already at the integrand level. In the sum of the diagrams the virtual part factorizes due to the eikonal identity and gives a vanishing scaleless integral as an overall factor. With real gluon momentum and virtual gluon momentum k, the two-loop integral explicitly evaluates to The calculations in eqs. (4.31)-(4.33) explicitly show that the η-regulator preserves non-Abelian exponentiation at two loops. We stress that the separate regulation of each web (in this case gluon) momentum is crucial for this. Likewise the non-Abelian exponentiation holds for the RRG renormalized TMD soft function.

Results in momentum space
The final result for the bare NNLO soft function in the fundamental representation is All pieces have been calculated in the previous sections. The bare one-loop result can be found in eq. (4.2) and in ref. [10]. The C F C A part of the NNLO real-virtual diagrams,

JHEP03(2016)168
R C F C A , is written in eq. (4.3). As shown in section 4.2.3, we can use the non-Abelian exponentiation theorem to express the C 2 F part (= I C 2 F ) in terms of the NLO expression in eq. (4.2). The double real terms, I C F C A , T , G, H and Q, are given in eqs. (4.24)- (4.27).
For the expansion of the bare soft function in and η we make use of eq. (A.1). The expanded momentum space results will therefore involve the plus distributions L T n (p ⊥ , µ) defined in eqs. (A.2) and (A.3). We give tables for Fourier transformations and the relevant convolutions among the plus distributions in appendix A. In the following we will usually suppress the arguments of the L T n . After expanding, we renormalize the soft functions according to eq. (3.1) and find the NLO counterterm and soft function coefficients in eq. (3.8): They agree with the ones of ref. [10] (for i = g). 13 Similarly, we obtain the NNLO coefficients 13 Up to misprints in eq. (5.62) and eq. (5.63) in the arXiv and journal version of their article, respectively.

JHEP03(2016)168
Here we have generalized our results to an arbitrary color representation according to Casimir scaling and parametrized them in terms of the anomalous dimension and integration constants introduced in section 3.1 for conciseness. The explicit values for the constants γ i S n , γ i ν n and S i n extracted from the above results are given below in eqs. (4.48)-(4.53). In S i we have set w = 1. In the renormalization factor Z i S we must however keep the η-dependent w as it is crucial for the derivation of the RRGE, see eq. (2.16). With eqs. (3.4) and (3.5) we thus find for the anomalous dimensions will in turn be relevant at higher orders. Such terms are of course also present in the position space results and the results for the beam functions given below. For the sake of brevity we will suppress them in the following.
We note, that our momentum space results depend on the renormalization scales µ and ν only through L T n , ln µ ν and a s . Correspondingly, the position space results depend on µ and ν only through L b , ln µ ν and a s . The scale dependence is dictated by the RGEs and RRGEs. In section 3.1, we have made this explicit in terms of the recursion relations eqs. (3.24) and (3.25). They are in perfect agreement with our results for the soft function and its ν anomalous dimension up to two loops.

Results in position space
Fourier transforming the p ⊥ -dependent soft function in the previous section to impact parameter space or using the recurrence relation in eq. (3.25) gives

JHEP03(2016)168
Similarly, one can Fourier transform the ν anomalous dimension coefficients in eqs. (4.41) and (4.42) or use the recurrence relation in eq. (3.24) to arrive at The µ anomalous dimension coefficients according to eqs. (3.6) and (3.17) have already been written in eqs. (4.39), (4.40). The constants that are not predicted by the RGE and RRGE are fixed by our explicit calculation of the soft function to where we have introduced the color factor C i , which equals C F (C A ) for i = q (i = g). We note the striking similarity of the (non-cusp) constants γ i S 1 and γ i ν 1 . Given the scheme dependence of γ i S 1 , we refrain however from investigating this further at the time.

Comparison to the literature
Several different formulations of the SCET TMD factorization theorem can be found in the literature. They differ by the treatment of the rapidity divergences, i.e. the way they are regulated, and how the individual soft and collinear functions are defined. Of course formulations in momentum as well as in position space are possible and directly related by Fourier transformation. While the combination of the TMD soft and beam functions is unambiguous, the three parts individually are scheme-dependent. In our case, this is reflected by the dependence on the unphysical scale ν. Varying ν allows a systematic quantification of effects from higher order terms related to the rapidity logarithms and therefore contributes to a reliable estimate of the overall perturbative uncertainty of the cross section. In section 5.1 we study the relation of our scheme to the one of refs. [11,21]. The NNLO factorization ingredients for the latter have been calculated in refs. [41,42]. Recently, the authors of ref. [43] have computed the soft function in the scheme of ref. [12] to NNLO. The result still contains explicit poles in and their rapidity regulator. The corresponding beam functions are only known to NLO so far. 14

JHEP03(2016)168
In refs. [11,21] the bare (b space) soft function is, by construction and to all orders in perturbation theory, identical to one. The analog to the renormalized/resummed soft function is absorbed into the collinear functions (TMDPDFs) and the 'collinear anomaly' exponential as we will see below. In refs. [41,42] the corresponding TMDPDFs and the collinear anomaly exponent are determined to NNLO. Identifying the exact relation to our scheme together with our NNLO results for the soft function allows us to extract the NNLO expressions for our beam functions. We will work this out in detail below. The results for the rapidity renormalized NNLO beam functions we present in section 6.
TMD distributions have actually first been considered in traditional QCD. Again, the final result for the differential cross section should agree with the SCET predictions (up to higher order terms). The relation of the SCET scheme we are comparing to the direct QCD approach of refs. [1][2][3][4][5][6][7][8][9], has been discussed in refs. [11,42].

Extraction of the TMD beam functions in the RRG scheme
In section 3.1, we have obtained recurrence relations that can be used to derive the logarithmic structure of the soft and beam functions to all orders. We have also computed the complete NNLO soft function and thus know by virtue of eq. (3.15) all anomalous dimensions at this order. Hence, the pieces we are missing are the boundary terms I (2) ij (z) and J (2) gj (z) in the beam function matching coefficients, see eq. (3.32). We now show how to determine them from the results for the NNLO TMD beam functions in refs. [41,42], which are based on the framework of refs. [11,21]. The comparison of the two formalisms is most conveniently performed in position space, where the TMD convolutions turn into ordinary products and the L T n distributions into logarithms L b . The translation between position and momentum space can easily be carried out with the help of tables 2 and 3.
To perform the comparison on the level of perturbative functions, we deconvolve the universal PDFs f i from the collinear functions and divide out the common hard function. The equality of the (perturbative) physical cross section then implies (2π) 6 where k, i, j ∈ {g, q,q, q },ḡ = g. 15 There are additional equations for k = g, where one or both of the pairs (I gl , I g/l ) are exchanged for (J gl , I g/l ) associated with the other tensor structure in eq. (2.10). These cases are understood in the discussion below. The expression in the first line of eq. (5.1) are the combined soft and collinear parts from the renormalized version of our factorization formulas in eqs. (2.2) and (2.3). The second line corresponds to the equivalent expressions in the scheme of refs. [11,21], where the first term is a resummed expression generated by the 'collinear anomaly' with the respective coefficient F k and the I i/j are the TMDPDFs. The details can be found in 15 The names of the variables in refs. [41,42] have been translated to our notation as x 2 T = b 2 and L ⊥ = L b . They also use a different definition of P (n) ij (z) which is related to ours by adding a factor 2 −n . JHEP03(2016)168 section 2.2 of ref. [42]. Note especially eq. (2.1) therein. The factor (2π) 6 arises because the normalizations of the three functions on both sides of eq. (5.1) differ by a factor of (2π) 2 each.
Strictly speaking eq. (5.1) only holds once both sides are expanded to a fixed order in α s (or for the hypothetical exact all-order solutions of beam and soft functions). If we include the resummation factors, according to eq. (3.35), and the analog for the I i/j at a finite logarithmic order, in general the two sides of eq. (5.1) will disagree by terms beyond that order. The reason is related to the resummation of the rapidity logarithms. The lack of the ν argument on the r.h.s. of eq. (5.1) indicates that for resummed beam and soft functions exact agreement of both sides is only achieved for a specific choice of the renormalization scales ν ± B and ν S , which turn out to be the canonical scales in eq. (3.38). Note that the ν dependence exactly cancels between the resummation factors of soft and beam functions, cf. eqs. (3.36) and (3.37).
Switching off the (MS) µ evolution by setting µ B = µ S = µ we can easily identify the exponential on the r.h.s. of eq. (5.1) with the product of the soft and beam function resummation factors. Concretely we have .
The last equality requires that we set ν S and ν ± B to the canonical values in eq. (3.38) and Exploiting the symmetry between the n-collinear andn-collinear beam functions, we can now extract the unknown coefficients iteratively from eq. (5.1) order by order in α s for each combination of partons k, i, j. Let us sketch two possible ways to do so: (i) We set all scales, i.e. µ B = µ S = µ as well as ν S and ν ± B , equal to the respective canonical values in eq. (3.38). We then divide both sides of eq. (5.1) by the exponential factor in eq. (5.2). In this way we have removed all logarithms from eq. (5.1) and the soft and beam function coefficients on the r.h.s. reduce to S i n and I gj (z), respectively. We can thus directly determine the missing coefficients I (ii) We expand eq. (5.1) to O(α 2 s ). Then the dependence on ν (and of course on µ S , µ B , ν S , ν ± B ) on the l.h.s. drops out. The only µ-dependent logarithm on both sides of the equation is therefore L b . We now write function and the symmetry between the matching coefficients of the n-andn-collinear beam functions. In fact it reveals the physical origin of the hard scale Q 2 = ω + ω − in the combined soft-collinear sector as an interplay of the rapidity related scales. We can thus uniquely identify the complete matching kernels I ki (z n , b, µ, ω ± /ν).
Method (i) is somewhat easier to carry out, but relies on the recursion relation in eq. (3.32). In contrast, method (ii) does not require any knowledge about the logarithmic structure of the beam functions. Instead the comparison of the result to eq. (3.32) serves as a strong consistency check of both frameworks and our result for the soft function. Another strong check is that our result for γ i ν confirms (5.3). We have independently verified that both methods work perfectly. We present the results for the I ki (z n , b, µ, ω ± /ν) in the next section.
Finally we emphasize again that the complete factorization of logarithms according to eq. (5.4) and the freedom to choose the ν scales is an advantage of using the RRG formalism. Eventually this provides the complete set of relevant scale variations as a tool to estimate the uncertainty in the resummed cross section.

TMD beam functions
As described in section 5.1 we use our soft NNLO results and eq. (5.1) to extract the NNLO TMD beam functions in the RRG scheme from refs. [41,42]. Here we present the results.

Results in position space
In b space we parametrize the matching kernel of the gluon beam function as ( We first give explicit expressions for the beam function matching coefficients for n = 0, 1, 2 according to the recurrence relation in eq. (3.32). We start at LO, i.e. O(α 0 s ): We then proceed iteratively and find gi (z) , (6.5)

JHEP03(2016)168
and The splitting functions P (n) ij (z) in our notation can be found in refs. [25,26]. The (Mellin) convolutions can be performed easily e.g. with the Mathematica package MT [59].
The non-cusp anomalous dimension γ i B can be obtained by analyzing the µ dependence of eq. (5.1). We find with γ i as defined in ref. [42]. This yields the following constants in the expansion according to eq. (3.19) γ q B 0 = 6C F , (6.10) (3)) . (6.12) From the procedure described in section 5.1 we finally obtain the constants in eqs. (6.4)-(6.7): 14) q q (z) = I (2) q /q (z, 0) , (6.16) q/q (z, 0) , (6.17) Explicit expressions for the I i/j (z, 0) are given in section 4 of ref. [42]. The J (2) gi matching constants have not been calculated in any scheme so far. For the case of scalar and vector boson production (including their decays) the J (2) gi are however not required for predictions with NNLL accuracy [4,42]. The reason is that the Lorentz structure associated with the J gj in eq. (2.10) is orthogonal to g µν and does not occur at LO, cf. eq. (6.3).
Note that due to charge conjugation and flavor symmetry, the above expressions determine the beam functions for all partons. Respecting flavor symmetry we only wrote q and q = q to denote quarks of unspecified flavor. For two quarks of flavor a and b we can write e.g. I qaq b = δ ab I qq + (1 − δ ab )I q q . From charge conjugation symmetry we have Irs = I rs and I rs = Ir s , with r, s ∈ {g, q, q }. Up to two loops also Iq q = I q q .

Results in momentum space
For completeness we here give the NNLO beam function results in momentum space. We obtain them by Fourier transforming the position space beam functions in section 4.4 using table 3. Care has to be taken in the transformation of the coefficients J ij as the second Lorentz tensor in (6.1) itself depends on b. We again express the matching coefficients in eqs. (2.9) and (2.10) as a power series in a s : 16 ij (z, p ⊥ , µ, ω ± /ν) a n s , (6.21) gi (z, p ⊥ , µ, ω ± /ν) a n s , with gi (z, p ⊥ , µ, ω ± /ν) = 0 , (6.24) gi (z)L T 0 , (6.26) 16 Note the factor of (2π) 2 difference to eq. (3.31).

Conclusions
We have calculated the TMD soft function in the RRG scheme of refs. [10,33] to NNLO. To regularize rapidity divergences we have employed the η-regulator [10,33]. We have explicitly demonstrated that this regulator preserves non-Abelian exponentiation of the TMD soft function at two loops. This represents a valuable consistency and practicability check of the regularization method. We present the new result in momentum (p ⊥ ) as well as in position (b) space in eqs. (4.38) and (4.45), respectively. We also obtain the soft two-loop anomalous dimension for the RGE and RRGE. The corresponding expressions in eqs. (4.49) and (4.50) exhibit an interesting similarity.
Based on the equality of cross section predictions and with the two-loop soft function at hand we have furthermore extracted the NNLO TMD beam functions in the RRG scheme from the results of refs. [41,42]. We have checked that our expressions for the soft and beam functions match the logarithmic structure predicted by recursion relations we have derived from the corresponding RGEs and RRGEs. The results for the beam function matching kernels and anomalous dimension are collected in section 6. All expressions for the soft and beam function coefficients are also available in electronic form upon request to the authors.
Our results represent the universal ingredients in the SCET factorization theorem for the peak region of the transverse momentum distribution of a color-neutral final state at the LHC with RRG resummation at NNLL . The perturbative uncertainties of such resummed cross sections can be systematically studied by independent variations of the different involved renormalization/factorization scales associated with the RG as well as the RRG. The phenomenological analysis at NNLL (+NNLO) of the transverse momentum spectra for processes like Drell-Yan or Higgs production is left to future work.

A Plus distributions
When expanding the bare expressions in η and , we make use of the distributional identity for µ 2 , p 2 ⊥ > 0, where we define the L T n (p ⊥ , µ) for n ≥ 0 in terms of the usual plus distributions 17 L T n (p ⊥ , µ) ≡ For an equivalent definition of the L T n (p ⊥ , µ) and more details about their properties and generalizations we refer to ref. [10]. By definition and with our choice of boundary condition 18  where D κ = {p ⊥ : |p ⊥ | ≤ κ} denotes a disc in p ⊥ -space around the origin with radius κ, which is set to µ in the first and an arbitrary positive value in the second equation. When dealing with the RGE structure of the theory, it is often convenient to work in position space where the plus distributions become ordinary logarithms and convolutions turn into ordinary products. In the following we give explicit expressions for the relevant Fourier transformations as well as for the convolutions of the L T n .
Fourier transformations. In table 2 we give the results of Fourier transforming the L T n with n ≤ 4 to b space in terms of L b = ln b 2 µ 2 e 2γ E /4 . The inverse Fourier transformations are collected in table 3. 17 Here we introduce the superscript T to distinguish the p ⊥ -dependent plus distributions from the notation Ln(x) = [θ(x)/x ln n x]+ as introduced in ref. [62]. It holds L T n (p ⊥ , µ) = (−1) n 2πµ 2 Ln(p 2 ⊥ /µ 2 ) for n ≥ 0. 18 This corresponds to λ = 1 in ref. [10]. Table 2. Fourier transforms of plus distributions L T n in terms of L b = ln b 2 µ 2 e 2γ E /4 . b-space p-space Convolutions. Here we list the convolutions of plus distributions relevant for this work. They can be derived from expanding eq. (F.23) in ref. [10]: . In most cases we encounter, one of the first three indices of the hypergeometric function is a non-positive integer, while the last two indices are no integers. These hypergeometric functions then trivially reduce to rational functions in the regulators. For the three non-trivial cases, we used the Mathematica package HypExp [57,58] to expand the hypergeometric function and obtained

JHEP03(2016)168
The 3 F 2 appearing in the last integral could not be expanded directly with HypExp. Using the integral representation of 3 F 2 , expanding the 2 F 1 in the integrand, and then integrating all the terms separately leads to the stated result.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.