Transverse Momentum Dependent Gluon Distribution within High Energy Factorization at Next-to-Leading Order

We discuss Transverse Momentum Dependent (TMD) gluon distributions within high energy factorization at next-to-leading order in the strong coupling within the framework of Lipatov's high energy effective action. We provide a detailed discussion of both rapidity divergences related to the TMD definition and its soft factor on the one hand, and rapidity divergences due to high energy factorization on the other hand, and discuss common features and differences between Collins-Soper (CS) and Balitsky-Fadin-Kuarev-Lipatov (BFKL) evolution. While we confirm earlier results which state that the unpolarized and linearly polarized gluon TMD agree in the BFKL limit at leading order, we find that both distributions differ, once next-to-leading order corrections are being included. Unlike previous results, our framework allows to recover the complete anomalous dimension associated with Collins-Soper-Sterman (CSS) evolution of the TMD distribution, including also single-logarithmic terms in the CSS evolution. As an additional result we provide a definition of $k_T$-factorization, i.e. matching of off-shell coefficients to collinear factorization at next-to-leading order within high energy factorization and the effective action framework. We furthermore establish a link between the QCD operator definition of the TMD gluon distribution and a previously derived off-shell TMD gluon-to-gluon splitting function, which is within the present framework obtained as the real 1-loop correction.


I. INTRODUCTION
Transverse momentum dependent (TMD) parton distribution functions (PDFs) [1][2][3] are objects of increased interest, since they allow to provide a more precise kinematic description of partonic scattering processes already at the leading order (LO) of perturbation theory. This is in particular true, if the observable of interest is not entirely inclusive. In that case, TMD PDFs provide an important advantage over a description based on collinear parton distributions. TMD PDFs arise naturally in processes which are characterized by a hierarchy of scales. With Q the scale of the hard reactions, TMD PDFs are at first defined for the hierarchy Q q T Λ QCD with q T the transverse momentum of the parton and Λ QCD the QCD characteristic scale of a few hundred MeV. The QCD description of such events gives then rise to the so-called Collins-Soper-Sterman (CSS) [4][5][6] resummation formalism. A different kinematic hierarchy in which TMD PDFs arise is provided by the perturbative Regge or low x limit √ s M Λ QCD , where √ s denotes the center of mass energy of the reaction and x = M 2 /s. While the resulting high energy factorization [7][8][9] does not primarily address the description of transverse momenta of final states, the ensuing formalism naturally factorizes cross-sections into transverse momentum dependent coefficients, so called impact factors, and transverse momentum dependent Green's function, which summerize logarithms in the center of mass energy. In particular Balitsky-Fadin-Kuraev-Lipatov (BFKL) evolution [10][11][12][13][14][15], as well as its non-linear extensions, keep track of * martin.hentschinski@udlap.mx transverse momenta along the evolution chain.
Both kinematic limits have a region of overlap, characterized through the hierarchy √ s M q T Λ QCD , which is of particular interest due its sensitivity to the emergence of a semi-hard dynamical scale in the low x limit, the so-called saturation scale [16]. The relation of both frameworks has been explored in a series of publications [17][18][19][20][21][22][23][24][25] and is currently used for a wide set of phenomenological studies see e.g. [26][27][28][29][30][31]. A somehow orthogonal approach has been put forward in [32][33][34]: instead of studying the region of overlap of both kinematic regimes, the goal has been to derive TMD evolution kernels which are meant to achieve a simultaneous resummation of both Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) and BFKL logarithms [35,36]. Such an approach seems to be of particular interest for Monte-Carlo applications such as [37][38][39][40][41]. While currently only real splitting kernels have been derived, which reduce in the regarding limits to the respective real DGLAP and BFKL kernels, the relation to CSS resummation is at the moment less clear, however the gluon-to-gluon splitting kernel could be shown to reduce to the Ciafaloni-Catani-Fiore-Marchesini kernel [42,43] in the soft limit.
In the following we aim at solving various open questions related to the previously mentioned studies. In particular we will give a next-to-leading order (NLO) study with respect to an expansion in the strong coupling constant α s of the gluon TMD in the high energy limit. While a related study has been already presented in [22] within the Color-Glass-Condensate approach, we will investigate this problem within the context of Lipatov's high energy effective action [44,45]. Starting with [46][47][48][49][50] the systematic determination of perturbative higher order corrections has been worked out for this framework, while in [51] equivalence with the Color-Glass-Condensate formalism has been demonstrated, including an re-derivation of the Balitsky-JIMWLK evolution, see also [52,53] for reviews; for further recent studies based on this framework see also [54][55][56][57][58][59][60][61]. In [62] a systematic framework for the determination of next-to-leading order corrections at cross-section level has been worked out. For the present study we will further extend this framework to include asymmetric factorization scale settings as required for a matching to collinear factorization, i.e. k T -factorization [8]. While this framework is currently limited to the dilute regime, i.e. 2 (reggeized) gluon exchange at the level of the cross-section, it has the great advantage that it allows to systematically study different choices of factorization parameter and schemes, which will be of particular use for the further exploration of the relation between BFKL and Collins-Soper (CS) [4,5] evolution, initiated in [20]. In particular we will determine systematically the NLO coefficients which relate the QCD operator definition of the unpolarized and linearly polarized gluon TMD PDF with the unintegrated gluon density of high energy factorization, which in turn will allow us to recover the complete CSS resummation scheme in combination with BFKL evolution, following closely related calculations based on collinear factorization [63,64]. We expect our result to be useful for a precise description of final states with small transverse momenta within high energy factorized cross-sections, at a similar level of accuracy as descriptions based on collinear factorization.
Another aspect of our result relates to the derivation of TMD splitting kernels in [32,34]. While the original derivation was based on a combined implementation of the Curci-Furmanski-Petronzio formalism for the calculation of the collinear splitting functions [65] and the framework of high energy factorization provided by [44], we find in the following that the real contribution to the QCD operator definition of the unpolarized gluon TMD yields precisely the previously derived off-shell TMD splitting kernel. Our current study provides therefore a possibility to recover the so far missing virtual corrections to these off-shell splitting kernels.
The outline of this paper is as follows. In Sec. II we give a precise definition of the goal of this paper in more technical terms, in particular the definition of the gluon TMD PDFs, Sec. III contains a brief review of Lipatov's high energy effective action and presents among other details an extension of the framework of [62] to k T -factorization. In Sec. IV we present the results of our NLO calculation, while Sec. V discusses aspects related to the interplay of CS and BFKL evolution. In Sec. VI we summarize our result and provide an outlook on future research.

II. THE SETUP OF OUR STUDY
The starting point of our study is the TMD factorization of a suitable perturbative process. To be specific we will refer in the following to the transverse momentum distribution of a Higgs boson, as discussed for instance in [64], see also [63]. With p H and m H transverse momentum and mass of the Higgs boson, this factorization is valid for |p H | m H and reads: where y H = 1/2 ln(x A /x B ) is the rapidity of the Higgs boson while x A,B denote the hadron momentum fractions of gluons stemming from hadron A, B respectively and where y c denotes the rapidity which divides soft gluons from hadron A and B; µ is the renormalization point of the cross-section. To be specific we consider scattering of two hadrons with light-like momenta p A and p B which serve to define the light-cone directions which yields the following Sudakov decomposition of a generic four-momentum, and n ± · k T = 0. Here, k T is the embedding of the Euclidean vector k into Minkowski space, so k 2 T = −k 2 . For Eq. (1), the top quark is considered to be integrated out and C t is the corresponding Wilson coefficient; H is the square of the on-shell gluon form factor at time-like momentum transfer q 2 = m 2 H , with infrared divergences subtracted, see [66]. To leading order in perturbation theory they equal one, while the precise NLO expression are not of interest for the following discussion and can be found for instance in [64]. σ 0 (µ) is finally the collinear Born level cross-section for the process gg → H, with v 246 GeV the Higgs vacuum expectation value; α s (µ) denotes finally the strong coupling constant at the renormalization point µ. For an unpolarized hadron, the TMD correlator Γ ij , i, j = 1, 2 can be further decomposed [67] where f g/B (x B , ζ B ; q b , µ) denotes the unpolarized TMD gluon distriubtion and h g/B (x B , ζ B ; q b , µ) the linearly polarized TMD gluon distribution in an unpolarized hadron. In terms of QCD fields the TMD PDF is defined as [63,64] whereS(2y c , σ; µ, ξ) denotes the soft factor and lim σ→∞ n(σ) = n − , with σ → ∞ a suitable regulator whose precise implementation will be given in Eq. (48) below. Gauge links are in general given as a combination of a longitudinal and a transverse gauge link [2], where the transverse gauge link is placed at light-cone infinity. Working in covariant gauge, the gauge field at infinity vanishes and the transverse gauge link therefore equals one. We will therefore in the following not consider the transverse gauge link. The longitudinal gauge link is on the other hand given by where v µ (x) = −it a v a µ (x) denotes the gluonic field and For the soft factor there exists various prescriptions in the literature, see e.g. [2,63,64,[68][69][70]. To keep the discussion as general as possible, we will consider below the most general soft factor introduced in [2,69], S(2y c , σ; µ, ξ) = S (2y c , 2y n ; ξ) S(σ, −2y c , ; ξ)S(σ, 2y n ; ξ) , with where n 1,2 (y 1,2 ) are tilted Wilson lines such that n 1 is placed at rapidity [94] y 1 /2 and n 2 at rapidity −y 2 /2. For a precise definition of the light-cone directions see Eq. (48) further below. The goal of the following sections is to study this gluon TMD in the high energy limit at next-to-leading order. In particular we will discuss the factorization of this TMD PDF into a perturbative coefficient, the BFKL gluon Green's function, and a hadronic impact factor. The latter two will then form the so-called unintegrated gluon density within high energy factorization. Our study is limited to the exchange of 2 reggeized gluons. It is known from various studies, that the gluon TMD also receives corrections due to the exchange of multiple reggeized gluons, which are of importance to take into account corrections due to high gluon densities and their possible saturation, see e.g. [17,24,27,29,71]. While these are very interesting questions -in particular since they can provide modifications of the region of very small transverse momenta due to the emergence of a saturation scale -we do not consider these effects in the following. Instead, great care will be taken to provide a complete discussion of various factorization scales and parameters both due to factorization in the soft-limit and the high energy limit, as well as UV renormalization. The study of these effects is somehow more straightforward, if the observable is restricted to 2 reggeized gluon exchange, which is the reason why we focus on this limit in the following. The obtained results may then be generalized at a later stage to the case of multiple reggeized gluon exchange.
Another motivation for this study is to link the above gluon TMD to the TMD splitting kernels derived in [32,34]. Below we will demonstrate that the TMD gluonto-gluon splitting of [34] arises directly from the real contributions to the 1-loop coefficient. We believe that this is a very interesting result, since it allows us to connect the framework of real TMD splitting kernels to the above operator definitions of TMD PDFs.

III. THE HIGH-ENERGY EFFECTIVE ACTION
Since the current study requires a small but important generalization in comparison to the framework presented in [29], we begin our study with a short review of the high energy effective action and the resulting calculational framework for NLO calculations. Our treatment of high energy factorization is based on Lipatov's high energy effective action [44]. Within this framework, QCD amplitudes are in the high energy limit decomposed into gauge invariant sub-amplitudes which are localized in rapidity space and describe the coupling of quarks (ψ), gluon (v µ ) and ghost (φ) fields to a new degree of freedom, the reggeized gluon field A ± (x). The latter is introduced as a convenient tool to reconstruct the complete QCD amplitudes in the high energy limit out of the sub-amplitudes restricted to small rapidity intervals. Lipatov's effective action is then obtained by adding an induced term S ind. to the QCD action S QCD , where the induced term S ind. describes the coupling of the gluonic field v µ = −it a v a µ (x) to the reggeized gluon field A ± (x) = −it a A a ± (x). High energy factorized amplitudes reveal strong ordering in plus and minus components of momenta which is reflected in the following kinematic constraint obeyed by the reggeized gluon field: Even though the reggeized gluon field is charged under the QCD gauge group SU(N c ), it is invariant under local gauge transformation δA ± = 0. Its kinetic term and the gauge invariant coupling to the QCD gluon field are contained in the induced term For a more in depth discussion of the effective action we refer to the reviews [52,53]. Due to the induced term in Eq. (12), the Feynman rules of the effective action comprise, apart from the usual QCD Feynman rules, the propagator of the reggeized gluon and an infinite number of so-called induced vertices. Vertices and propagators needed for the current study are collected in Fig. 1. Determination of NLO corrections using this effective action approach has been addressed recently in series of publications [46][47][48][49][50]. For a discussion of the analogous high energy effective for flavor exchange [72] at NLO see e.g. [73][74][75].

A. Determination of NLO coefficients
The framework for the determination of NLO corrections has been established in [62] within the determination of the NLO forward Higgs production coefficient in the infinite top mass limit. We will therefore frequently refer to process as an example process in the following, where we further assume that the particles in the fragmentation region of the scattering particles are widely separated in rapidity.
The partonic impact factor of the quark with momentum p b will be later on replaced by the hadronic impact factor, which forms together with the BFKL Green's function the unintegrated gluon density. As in [62], we will consider matrix elements normalized to match corresponding collinear matrix elements for vanishing virtuality of the reggeized gluon state. We therefore have where we average over incoming parton color as well as the color of the reggeized gluon and sum over the color of produced particles; X (n) a denotes the n-particle system produced in the regarding fragmentation region. With we arrive at the following definition of an off-shell partonic cross-section dσ a+ and the corresponding impact factorĥ k T (k): q, a, ± k, c, ν = −iq 2 δ ac (n ± ) ν , Note that this impact factor can in principle be arbitrarily differential, as far as the formulation of high energy factorization is concerned; for a corresponding definition of the other impact factor we refer to [62]. The above expression is subject to so-called rapidity divergences which are understood to be regulated through lower cut-offs on the rapidity of all particles, η i > −ρ/2 with ρ → ∞ and i = 1, . . . , n for n the number of particles produced in the fragmentation region of the initial parton a. For virtual corrections, the regularization is achieved through tilting light-cone directions of the high energy effective action, Below we will also comment on the possibility to regularize rapidity divergences through tilting light-cone direction also in the case of real corrections, see Sec. III D.
As shown through explicit results [46,[48][49][50]62], impact factors contain beyond leading order configurations, which reproduce factorized contributions with internal reggeized gluon exchange. It is therefore necessary to subtract these contributions, see also Fig. 2. To this end one defines the bare one-loop 2-reggeized-gluon Green's function G B (k 1 , k 2 ) as well as the impact factors through the following perturbative expansion Using the following convolution convention we then define the following subtracted bare NLO coefficient, and [95] dσ NLO where we added the super-scripts 'k T ' and 'ugd' to indicate that the impact factor of the particle with momentum p a,b refers to the hard event and the unintegrated gluon distribution respectively. While rapidity divergences cancel for the the above expression, its elements still depend on the regulator. As a next step we therefore define a renormalized Green's function G R through which yields where The transition functions Z ± have a twofold purpose: they both serve to cancel ρ-dependent terms between impact factors and Green's function and allow to define the BFKL kernel. In particular, where

B. Transition function and finite terms
In the following we generalize the treatment given in [62] through including as well the most general finite contribution into our discussion. The need to include finite contribution into this transition factors has been first realized in the determination of the 2-loop gluon Regge trajectory in [48][49][50], where both divergent and finite terms could be simply taken to exponentiate, since one is dealing with 1-reggeized gluon exchange contributions only. A suitable generalization, which both reduces to the exponential ansatz of [50] and obeys Eq. (28), is then given by which is sufficient for a discussion up to NLO accuracy. As we will see in the following, these finite contributions serve a twofold purpose. At first they remove a potential finite contribution in the bare Green's function, where which then yields directly the 1-loop BFKL kernel To keep the treatment as general as possible, the finite contribution is on the other hand then split up into 2 terms: The contributionf ± is used to transfer finite terms contained in the Green's function to the impact factors. We take in the following this function to be identical for both plus and minus directionf ±,(1) (k 1 , k 2 ) =f (1) (k 1 , k 2 ). Indeed, since we are essentially dealing with contributions due to the gluon polarization tensor in the high energy limit, such a symmetric treatment appears to be the appropriate one. One finds The second contributionf ± was absent in the discussion of [62]. It needs to satisfy to 1-loop the following requirement, to ensure absence of finite terms in the Green's function: Note that similar constraints can be imposed on higher order contributions to the functionsf ± to achieve a Green's function without finite contributions at higher orders. Including all contributions, one finally arrives at following expression for the NLO coefficient: where are given in Eq. (33) and Eq. (35) respectively. The functionf (1) (k 1 , k 2 ) as well as the evolution parameters η a,b are on the other hand still undetermined. They are in principle arbitrary, but should be chosen such that both impact factors are free of large logarithms.

C. Scale setting for kT -factorization
The parameters η a,b as well as the functionf are at first arbitrary; the former define through the combination η a − η b the evolution parameter of the 2 reggeized gluon Green's function G R , see also the related discussion in [53,62]. As usually, it is necessary to chose these parameters such that the next-to-leading order corrections to impact factors are under perturbative control and that large logarithms in the center of mass energy are resumed through the 2 reggeized gluon Green's function. Within the k T -factorization setup, one of the impact factors, e.g. the coefficient C R,b (η b , k 2 ) in our example, is to be replaced by the hadronic impact factor, which then builds together with the BFKL Green's function the unintegrated gluon density. Even though the hadronic impact factor is naturally a non-perturbative object, at the very least for transverse momenta |k 2 | 1.5 GeV, it must have an expansion in terms of collinear parton distribution functions and corresponding collinear coefficients, see e.g. [76]. In order to have a complete matching of the resulting expression for the unintegrated gluon density to the collinear gluon distribution in the double-logarithmic limit, it is natural to chose the evolution parameter η a −η b to coincide with the fraction of the hadronic momentum carried on by the gluon x M 2 a /s with M a the invariant mass of the system produced in the fragmentation region [96]. Note that at NLO, an asymmetric scale choice, i.e. choosing the reference scale of the center-ofmass energy √ s to be of the order of a typical scale of one of the impact factors, leads to a modification of the NLO BFKL kernel, see e.g. [14,77]. To repeat this exercise within the context of the high energy effective action, we reconsider Eq. (27), but focuse now on the factorization parameters η a,b and the finite terms introduced above: where M a,b is the invariant mass of the produced final state corresponding to each of the impact factors and C R collects all terms which are independent of both η a,b andf . To equal the evolution parameter of the Green's function with the hadron momentum fraction carried on by the reggeized gluon entering the impact factor C R,a , we set where M, M 0 are so far unspecified reference scale and x 0 = M 0 /M is a parameter of order one, which allows to estimate the scale uncertainty associated with high energy factorization. In the following we chose M to be of the order of M a , i.e. the hard scale. While this is a natural choice for the hard impact factor, it introduces the same scale into the hadronic impact factor, characterized in general by small transverse momenta. We therefore find in the perturbative region of the hadronic impact factor a large collinear logarithm, which at first spoils the convergence of the perturbative expansion. This logarithm can however be absorbed into thef function through settinḡ which eliminates the logarithm in M from the hadronic impact factor. Note that the choice of |k 2 | as the relative scale is somewhat arbitary, and other choices are equally possible, see e.g. [76]. It is interesting to compare this situation to the case where the parameters η a,b are identified with the rapidities of the system produced in the regarding fragmentation region, η a = ln p + a /M a , which allows to verify that the presented treatment agrees -after setting x 0 = 1 -with the one derived in [77], based on a study of ladder diagrams within the Quasi-Multi-Regge-Kinematics in the context of the definition of the NLO inclusive jet vertex. For the NLO BFKL kernel one finally obtains the following contribution which is independent of the parameter x 0 and where K (2) denotes the NLO BFKL kernel if η a,b is identified with the rapidities of the external particles withf (1) = 0. The above expression is in agreement with [77] and [14]. A more detailed discussion of possible choices of the functionf (1) will be presented elsewhere.
Summing up we have the general definition of the unintegrated gluon density where we suppressed the dependence on the invariant mass M b since the latter can in general be expressed in terms of the transverse momentum. The k T -factorization scheme fixes thenf (1) through Eq. (40), while ∆η ab = η a − η b is set to ∆η k T ab ≡ ln 1/x g . The high energy factorized cross-section is then obtained as where 'A'might either denote a parton a, a partonic impact factor convoluted with a parton distribution function of a hadron A or a colorless initial state which allows for a perturbative treatment. Concluding we remark that in [78] a definition of the unintegrated gluon density has been proposed in terms of a operator definition with the high energy gluonic field in light-cone gauge, which requires the inclusion of both so-called 2,3, and 4 body contributions. While an interpretation of such contributions in terms of induced vertices Fig. 1 of the high energy effective action appears to be possible, the precise relation remains unclear.

D. Regularization of rapidity divergences
As pointed out in the above discussion, high energy factorized matrix elements are subject to so-called rapidity divergences. While they cancel at the level of observables after subtraction of factorizing contribution and use of the transition function, intermediate results beyond leading order require a regulator in order to arrive at well-defined matrix elements. While for real production a cut-off on the rapidity of produced particles provides a natural way to regulate such divergences, a consistent regularization is more difficult to achieve in the case of virtual diagrams. As already pointed out in Sec. III A, a suitable way to regulate these divergences in the case of virtual corrections is to tilt the light-cone directions of the high energy effective action away from the light-cone, see Eq. (20). Note that from a formal point of view this a very attractive way of regularizing rapidity divergences, since gauge invariance of the high energy effective action does not depend on the property n ± 2 = 0; tilting light cone directions provides therefore a gauge invariant regulator, similar to dimensional regularization. Nevertheless the current treatment, see e.g. [46,48,62] is somewhat unsatisfactory, since it treats real (cut-off) and virtual (tilting) corrections on somewhat different grounds. At the same time, tilting light-cone directions is also a frequently used regulator for the determination of the 1-loop corrections to TMD PDFs within collinear factorization, see e.g. [2,63] and references therein, which is being use for both real and virtual corrections. It seems therefore natural to regulate rapidity divergences through tilting light-cone directions also in the case of real corrections.
From a technical point, this does not imply any major complications. However the real part of the 1-loop Green's function Eq. (31) would receive a finite correction which in covariant Feynman gauge is related to the square of the induced diagrams (last 2 diagrams in Fig. 3). As can be seen already at the level of diagrams, such contributions arise as well for the corresponding impact factors, which contain an identical diagram once calculating the correction due to the emission of an additional real gluon. As a consequence, it is straightforward to show that the corresponding contributions cancel, once subtracted impact factor and central contribution are being combined. Moreover, such a contribution may be easily absorbed into a generalized version of the functionf (1) , Eq. (35). While including such contributions does not provide any substantial complication, one deals in that case with an entirely spurious contribution, which merely arises due to our choice of our regulator and which has no physical meaning. It seems therefore natural to employ a regulator which avoids such a contribution altogether, at least at 1-loop. The modified regulator is essentially identical to the previously used tilted light-cone vectors, while the tilted elements are taken now to be complex, i.e. we will use in the following As a consequence one has for virtual corrections while real corrections yields |n a,b | 2 = 0, n a · (n b ) * + c.c. = 4(1 + e −2ρ ).
The spurious self-energy like contributions are therefore absent. The only disadvantage of this method is that terms of the form ln(n a,b 2 ) in virtual corrections can give rise to undesired imaginary parts due to space-like n a,b 2 . While at cross-section level such imaginary parts cancel naturally, if one limits oneself to NLO corrections, a consistent treatment of such contribution at amplitude level would require to absorb this imaginary part into the parameter ρ, e.g. through a suitable replacement ρ →ρ = ρ − iπ/2 etc. in the transition functions.
In the following calculation we will meet rapidity divergences which originate both from high energy factorization and the QCD operator definition of the TMD gluon distribution and the corresponding soft factor, which we will consistently regulate through the tilting as described in Eq. (45), while we reserve the use of the regulator ρ → ∞ for rapidity divergences due to high energy factorization. To be specific we define in the following the tilted Wilson lines of the TMD definition Eq. (7) and Eq. (11) as

IV. DETERMINATION OF THE GLUON TMD
The goal of the following section is to determine the gluon TMD Eq. (7) within high energy factorization, i.e. we aim at the determination of the following coefficient C gg * , implicitly defined through Note that the TMD PDFs at first do not depend on the proton momentum fraction x, since high energy factorization requires to integrate over this longitudinal momentum fraction. Such a dependence therefore only arises through a special choice for the parameters η a,b . To allow for a separate discussion of the different contributions of the gluon TMD, we further define where for the moment we define the TMD PDF in 2 + 2 dimensions, since individual expressions are divergent and q − = xp − B . We therefore obtain In the following we evaluate the above gluon TMD for an initial reggeized gluon state with polarization n + at 1-loop. To be precise we consider where the reggeized gluon state r b (k) is defined with the normalization Eq. (17), appropriate for matching to collinear factorization in the limit of vanishing transverse momentum,i.e.
v a µ (ξ) r b while high energy factorization requires to integrate over the minus momentum, Feynman rules for the determination of perturbative corrections are summarized in Fig. 4. With the following convention to denote the perturbative expansion in α s for a generic quantity A(α s ) where A (n) ∼ α n s , we have finally at leading order and therefore and the TMD gluon distributions are up to an overall factor [97] of 1/π at leading order identical to the unintegrated gluon density [17] where the k T -factorization scheme defined in Sec. III C yields expression which are closest to conventional collinear factorization results. Note that the unintegrated gluon density is therefore directly related to the operator definition of the gluon TMD. Moreover, in the dilute limit i.e. considering only 2 reggeized gluon exchange, the unintegrated gluon density is universal [98]. We further stress that the distribution of linearly polarized gluons in an unpolarized hadron is non-zero within high energy factorization already at tree level, in contrast to the result found within collinear factorization [64]. From a technical point of view this is of course easily understood, since the initial gluon carries within high energy factorization already finite k and therefore gives rise to such a distribution.

A. One-loop calculation without soft-factor
To regularize infrared and ultraviolet divergences we use dimensional regularization in d = 4 + 2 . The vir-tual 1-loop correction is provided by the set of diagrams Fig. 5. The last diagram in the second line vanish within dimensional regularization, since it is scale-less. Including the contribution from the complex conjugate amplitude, we obtain the following result To obtain the projections on the distribution for the unpolarized TMD PDF f g and the linearly polarized gluons h g in an unpolarized hadron, we use which amounts to replace the overall tensor structure q i q j /q 2 by one for both the unpolarized and the linearly polarized TMD. In particular, due to the presence of nonzero initial transverse momentum, the virtual correction for both TMD distributions are non-zero and agree with each other. Diagrams for real corrections are depicted in Fig. 6. Parameterizing k − = q − /z and correspondingly l − = q − (1 − z)/z, a straightforward calculation yields where where we used l = k −q. Note that the splitting function corresponding to f g coincides with the real transverse momentum splitting function derived in [34] if we set ρ = ∞ = σ. As demonstrated in [34], this splitting function coincides with the DGLAP splitting function in the limit k → 0, reduces to the real part of the leading order BFKL kernel in the limit z → 0 and yields the real leading order kernel of the CCFM equation in the limit q → k. The present calculation provides on the other hand an opportunity to determine the still missing virtual contribution to this splitting function. We further note -in agreement with the result presented in [64] -that the coefficient of the singularity z → 1 vanishes for the polarized splitting, Eq. (64) in the collinear limit k → 0, after averaging additionally over the azimuthal angle of the incoming gluon. For finite k this singularity is on the other hand present and requires a treatment similar to the case of the unpolarized gluon TMD PDF. With the integrated real corrections -relevant for the current discussion which is based on high energy factorizationwe finally have: where the remaining integrals over z are finite and π is defined in Eq. (32). If inserted into Eq. (49), the convolution integral over the reggeized gluon momentum k gives still rise to an infra-red singularity. A possibility to extract these singularities is the use of a phase space slicing parameter, see e.g. [79][80][81] for an example within the high energy effective action. While this is sufficient to demonstrate finiteness at a formal level, the use of such phase space slicing parameters is in general complicated for numerical studies at NLO accuracy. For the case of collinear NLO calculation, the by now conventional tool to overcome this difficulty are subtraction methods, in particular dipole subtraction [82]. In [62] this has been slightly generalized and applied to the case of divergences which arise due to convolution integrals of transverse momenta. In particular the following decomposition has been proposed: with The expression in the squared brackets on the right-hand side vanish in the limit |l| → 0, and G(k) is a function which parameterizes the transverse momentum dependence of the reggeized gluon state. The function κ(l) is such that the integral on the right hand side of (69) is well-defined, which in practice means that it does not behave worse than ln |l| for |l| → 0 and |l| → ∞. Furthermore, it should be such that the integral in the second line of (68) can be calculated analytically. Note that the factor q 2 /[l 2 + (q + l) 2 ] is needed to achieve convergence in the ultraviolet. In the current setup we merely require the case κ(l) = 1 with and where the corresponding convolution integral can now be defined in d = 2 dimensions.

B. Soft factor, counter terms and renormalization
The above result contains both rapidity divergences due to high energy factorization (ρ → ∞) and the TMD definition (σ → ∞) as well as single and double poles in 1/ . Rapidity divergences due to high energy factorization require the subtraction of high energy factorized contributions to the above correlator as well as the application of the transition function, as given by Eq. (37). Rapidity divergences due to soft radiation require the soft factor Eq. (10) which we did not include so far. With the 1-loop expansion and taking the limit[99] σ, y n → ∞, we obtain finally for the soft function in momentum space at 1-loop where y c is a finite rapidity and takes a role related to a factorization scale, similar to the parameter η a in the case of high energy factorization; in particular it enters directly the defintion of the scales ζ A,B defined in Eq. (2). The 1-loop contribution to the TMD gluon densities due to the soft function is finally given bŷ Note that the 1-loop soft-function, which consists of its real part only within dimensional regularization, agrees with the real part of the 1-loop BFKL kernel. We believe that this coincidence is limited to the 1-loop case and does not hint at a general universality of the rapidity dependence of the soft function. In particular, at the diagrammatic level the soft-function does not give rise to the complete Lipatov vertex, but only to terms related to the induced contributions, see also the discussion in [83] in the context of soft-collinear effective theory for high energy scattering.
As a last step we need to combine virtual (Eq. (59) and real (Eqs. (66)(67)) corrections, with the soft factor, making use of the appropriate projections. Note that the contribution of the soft-factor merely amounts to a replacement of the regulator σ by the factorization parameter y c . We obtain for the 1-loop high energy subtracted and renormalized 1-loop coefficientsC (1)f,h the following result: where While the above expression no longer carries rapidity divergences, it still comes with several poles in 1/ , which are of ultraviolet origin and which require renormalization. The corresponding renormalization constant is identical for both unpolarized and linearly polarized gluons and is obtained as which gives rise to the following anomalous dimension where we used dα s /d ln µ = 2 α s and ζ = (q − ) 2 e 2yc . Note that the above anomalous dimension agrees with the corresponding result obtained within a treatment based on collinear factorization [64]. This is indeed to be expected since it arises due to the renormalization of ultraviolet divergences, which are naturally independent of the non-zero transverse momentum of the initial state gluon. We however stress that the linearly polarized TMD gluon distribution does not give rise to the above anomalous dimension within collinear factorization, since the corresponding distribution vanishes within collinear factorization at tree-level; the 1-loop result is therefore not renormalized. We finally obtained for the renormalized coefficients where we made now the dependence on various parameters y c and η a explicit. The above expressions for the 1-loop cofficients of unpolarized and linearly polarized gluon TMD are one of the main results of this work.

V. EVOLUTION
The coefficients Eqs. (81), (82) depend on three factorization scales and/or parameters: µ (renormalization scale), y c (evolution parameter of the soft function) and η a (evolution parameter of the unintegrated gluon density). In addition we still have the dependence on the functionf (1) , which is also related to the unintegrated gluon density. The dependence on the renormalization scale and the factorization parameter y c gives rise to CSS resummation framework, [4][5][6]. In the treatment established for collinear initial states, see e.g. [2,64,84] it is customary to consider to this end the Fourier-transform of the TMD coefficient to transverse coordinate spacẽ where i = f, h and then to evolve the coefficient in coordinate space, with the evolution operator In the above expression, the CS kernelK CS (b, µ i ) is in general assumed to have both perturbative and nonperturbative contributions; for a detailed discussion see [2,64,84]. The non-perturbative contribution arises due to the inverse Fourier transform, which requires to integrate over large values of b, well into the non-perturbative region. Note that a similar statement applies in principle for the convolution integral of Eq. (97), if the unin-tegrated gluon distribution does not drop-off sufficiently fast for small values of transverse momentum k. The actual evolution takes place in two steps: first one evolves the coefficient at a certain initial renormalization scale µ i from an initial rapidity y c,i -parameterized through ζ B,k = (q − ) 2 e 2y c,k , k = i, f -to a final rapidity y c,f . The second step evolves then the TMD PDF from the initial to the final renormalization scale. While the value of the final renormalization scale is of the order of the hard scale, i.e. the Higgs mass for the current example, the initial renormalization scale µ i must be chosen such that it minimizes the perturbative correction to the TMD coefficients. In collinear calculations it it is naturally taken to be of the order of the transverse momentum q or its inverse conjugate coordinate µ b = 2e −2γ E /|b| with γ E 0.577216 the Euler constant. In the current setup, the optimal choice is far from apparent, since the coefficients depend on multiple scales due to non-zero initial transverse momenta. In particular, the transverse momenta k and q are at least at first not necessarily of the same order of magnitude.

A. A comparison of the kernels of CS and BFKL evolution
Both CS and BFKL evolution describe evolution in rapidity. It is therefore natural to expect that both evolution and their respective kernels have a certain overlap. For the derivation of the CS kernel, we follow closely [2,84], where the kernel of the Collins-Soper evolution equation is defined through the y c dependence of the renormalized soft factor, and is itself subject to the following renormalization group equation, where Γ A cusp is the cusp anomalous dimension in the adjoint representation, see [64,85,86] for higher order terms. At 1-loop one finds While the representation in transverse coordinate space is very useful for a direct solution of the CS-equation through exponentiation of the CS kernel as in Eq. (85), it is also instructive to formulate the CS-equation in momentum space, which allows for a direct comparison with the BFKL equation. In particular, defining in complete analogy to the BFKL treatment a CS Green's function G i (∆y, k 1 , k 2 ), i =BFKL, CS, such that one finds at 1-loop the following simple relation between both kernels: The presence of this factor can be explained as follows.
As it is well known, the virtual correction to the 1-loop BFKL kernel is directly related to the gluon Regge trajectory ω( , k) which in transverse momentum space can be written as The CS equation is on the other hand limited to soft radiation, which implies restriction to transverse momenta |l| |k| and |l − k| |k|. The integrand in the above expression therefore reduces to and the integral vanishes within dimensional regularization through a cancellation of the infra-red and ultra-violet 1/ poles. Removing on the other hand the UV pole through renormalization, one finally ends up with the virtual contribution to the CS kernel. The CS evolution can be therefore understood as the soft approximation to the complete BFKL kernel. Indeed, such an identification is natural, since the CS kernel arises from the rapidity divergence of the soft-function, while the BFKL kernel from the rapidity divergence of partonic cross-sections in the high energy limit √ s → ∞.
The above discussion clearly suggests that the CS and BFKL evolution are closely related to each other, with the CS kernel as the soft limit of the complete BFKL kernel. When considering CS and BFKL evolution for the same quantity, it is therefore necessary to remove the overlap of both evolution equations or alternatively to restrict them to distinct regions in phase space, which are then covered by the regarding evolution equation. In particular the various evolution parameters must obey the following ordering it is needed to separate the phase space covered by BFKL evolution (rapidity evolution of the entire cross-section) and CS evolution (rapidity evolution of soft gluons only). Clearly, soft gluons form a sub-set of the complete crosssection and cannot be evolved separately from the latter in rapidity.

B. kT factorization and alternative schemes
In the following we investigate our result for a specific choice of the evolution parameter of the unintegrated gluon density ,which we fix to coincide with the proton momentum fraction x, i.e. we consider now the equivalent of Eq. (49), but with the following choices for the parameters of the high energy evolution, following the results of Sec. III C ∆η ab = ln x 0 x , η a = ln where M is a still unspecified reference scale. The coefficients take the following form with l = k − q. Even though both η a and the function f (1) depend within this scheme on a certain reference scale M , the dependence on this scale cancels between both contributions, and we remain only with the parameter x 0 which is of order one. The scale M remains there-fore unspecified and can be used to satisfy the ordering condition Eq. (93). A possible and suitable choice is then ζ i B = C · q 2 , with C another constant of order one, which eliminates a potential large logarithm in the coefficients and which specifies eventually M = |q|. The complete resummed gluon TMDs take then the following final form where ζ i B = M 2 and µ i are yet unspecified scales. Reading the above expressions from the right to the left, one first evolves the unintegrated gluon density through BFKL evolution up to the hadron momentum fraction x, with a corresponding factorization uncertainty parametrized through x 0 . The unintegrated gluon distri-bution is then convoluted with the NLO TMD coefficient in transverse momentum k. The resulting expression defines then the gluon TMD with transverse momentum q at a scale ζ i B = M 2 and renormalization point µ i . While M is an arbitary scale, introduced to define the k T -factorization scheme, it is naturally chosen to be of the order of the transverse momentum q ; the same is true for the choice of the renormalization point µ i . In a next step it is therefore needed to evolve this gluon TMD both in ζ (rapidity evolution of soft gluons) and finally in the renormalization scale to its final values, where rapidity evolution of solution gives rise to a convolution in transverse momentum q . The above expressions for CSS evolution (combined evolution in ζ and µ) are the conventional expresses found in the literature, while they are expressed in transverse momentum instead of transverse coordinate space. In particular, the CS Green's function in transverse momentum space is obtained from the frequently used transverse coordinate expression through At 1-loop,K CS is given by Eq. (86). For perturbative higher orders see e.g. Sec. III of [64], where it is needed to include a relative factor of 2 with respect to the convention employed in this paper. Apart from perturbative higher order corrections, one might also consider RG evolution from the scale µ i to a suitable renormalization point of the CS Green's function; finally it is also possible to include a model for non-perturbative effects. If one restricts oneself on the other hand to the leading order kernel Eq. (86), the above integral can be easily evaluated and one finds G LL CS (∆y, q 1 , q 2 , µ) = Γ(1 − α s ∆y) (q 1 − q 2 ) 2 Γ(α s ∆y) with α s = α s C A /π. Note that convergence of the Fourier integral requires α s ∆y < 1 for the above expression. While the kernels of CSS and BFKL evolution in Eq. (97) are well known, our result provides as a new element the perturbative coefficient which connects both evolution equations up to NLO accuracy. In particular, a complete next-to-leading logaritmic resummation of both BFKL and CSS logarithms requires to combine the NLO coefficient Eqs. (95), (96) with the unintegratd gluon distribution evolved with the NLO BFKL kernel [14] as well as CSS evolution with NLO anomalous dimension Eq. (80), and the corresponding expression for the CS kernel, see [64] for a compact summary up to NNLO accuracy of these elements in the gluonic channel. It would be very interesting to compare this result to the low x expansion of exact N 3 LO results for the gluon TMD PDFs [87]. From a technical point of view, this would require to construct a partonic unintegrated gluon distribtion, following Eq. (43), but using NLO quark and gluon impact factors.
While the identification of the evolution parameter according to the k T -scheme, as used above, provides a direct generalization of the collinear result and is been often employed in fits of the unintegrated gluon density, see e.g. [88][89][90][91], it is not necessarily the most adequate to describe rapidity evolution of the system. An alternative form would be to identify η a with the maximal rapidity of the soft gluonic system, η a = y c or with the rapidity of the hard final state, i.e. the Higgs boson η a = y H . While η a = y c eliminates entirely the need for CS evolution, it also ignores the rapidity of the hard event (Higgs boson) in the energy evolution of the TMD PDF. The choice appears therefore to be possilbe, but inadequate. The choice η a = y H evolves the unintegrated gluon density through BFKL evolution up to the rapidity of the hard even, while CS evolution addresses the possible differences y c − y H , where both y c > y H and y c < y H is possible. Introducing furthermore the parameter δy, which allows to address the scale uncertainty associated with high energy factorization, i.e. high energy evolution describes dependence onȳ H = y H + δy and e ±δy is taken to be of order one, one finds where ζ i B = (M 2 H + q 2 )e −2δy , ζ f B = (M 2 H + q 2 )e 2(ȳ H −yc) and y 0 is a parameter of the order of the hadron rapidity.
Note that within this frame, the TMD PDFs no longer depend on the hadron momentum fraction, but rather on rapidity. While this might appear strange at first sight, it is natural from the point of view of high energy factorization where the momentum fraction is -in the case of the k T -factorization scheme -merely an evolution parameter fixed through the kinematics of the final state, while the above rapidity scheme uses a different choice for this evolution parameter.

C. Relation to previous results in the literature
Before we conclude, we briefly discuss the relation of the above results to results in the literature. In principle the TMD gluon distribution has been already study within high energy factorization at NLO in [22], previous studies with a similar scope are [20] and [92,93]. At the level of real corrections, the contributions seem to be identical at the level of Feynman diagrams, leaving aside the absence of the multiple reggeized gluon exchange in the current discussion. The set of virtual diagrams of [22] is on the other hand clearly reduced with respect to the ones considered in this work, Fig. 5. Indeed, the authors of [22] seem to consider only self-energy corrections to the Wilson line in their approach. This difference is direcly related to the fact that [22] makes use of the so-called 'shock-wave picture' for the calculation of next-to-leading order corrections. While this is a frequently used frame for the calculations within the CGC-framework at both LO and NLO, it does not allow to recover the term proportional to β 0 in the anomalous dimension Eq. (80), as already noted by the authors of [22]. We believe that this constitutes an important advantage of the framework of the high energy effective action, since it does not only allow to recover the double-logarithmic contribution to CSS resummation, (i.e. the Sudakov form factor of [22]), but also its single logarithmic terms.
Another point in which we somehow differ with [22], see also the discussion in [20], is the statement that BFKL and CS evolution cover by default distinct regions of phase space. As outlined above, this problem does strictly speaking not occur, if the evolution variable of the gluon density is identified with the hadron momentum fraction. To clarify this point further and to put this into context with the discussion in [22], we consider in the following the combined z → 0 and z → 1 singularities of the real corrections of our 1-loop result. Since the treatment is slighly more involved for the linearly polarized gluon TMD, Eq. (64), we focus in the following on the unpolarized case, Eq. (63). Removing regulators through taken the limits ρ, σ → ∞, and keeping only singular terms, we find that the TMD splitting function reduces to 1 0 dz 1 (q − k) 2P (0)f gg,r (z, qk) At first sight, the poles at z = 0 and z = 1 are therefore indeed well separated. With the rapidity of the produced gluon equal to η l = ln |q−k|z (1−z)q − , the above integral can be however rewritten as where we re-inserted the previously removed regulators ρ, σ as cut-offs on the rapidity integral. As for the real part of the 1-loop BFKL kernel within the high energy effective action, see [46,48,53] for an explicit construction, the above expression is proportional to an integral which extends over the entire range of rapidity. In contrast to the derivation of the BFKL kernel, the above integral is however split up into a 'soft' and a 'hard' part, where the distinction into 'soft' and 'hard' is essentially achieved through the virtual corrections. As a consequence -while we somehow agree with [22] that both pieces are well separated -care is needed to avoid overlap between both contributions. In particular, a gluon with a certain fixed rapidity may be either counted as soft or hard but never as both. At the same time, 'gaps' in rapidity should be avoided for a consistent and correct description.

VI. CONCLUSIONS AND OUTLOOK
In this paper we extended the framework established in [62] for next-to-leading order corrections within Lipatov's high energy effective action to the case where the transition function contains an additional finite contributions. Using this extension we were able to address the special, but important case of impact factors which possess a strong hierarchy with respect to their transverse scales, as it is the case within the k T -factorization setup. The latter is characterized by a impact factor with a hard scale, i.e. the off-shell partonic coefficient, and a hadronic impact factor, characterized by transverse momenta in the non-perturbative domain. The resulting expression have been found to agree with existing results in the literature, which have been established through a study of multi particle production amplitudes in the (Quasi-)Multi-Regge Kinematics [77]. Establishing this formalism at NLO within the high energy effective action is the first key result of this paper.
Another key result is the determination of the nextto-leading order corrections to the gluon TMD PDFs in high energy factorization, making use of the established formalism for the renormalization of matrix elements of reggeized gluon fields. While unsubtracted NLO result is subject to both rapidity divergences due to high energy factorization and rapidity divergences due to definition of the TMD PDF, the subtracted and renormalized coefficient is completely free of such divergencies. In particular we stress that rapidity divergences related to the definition of the TMD PDF can be treated using the soft factor, established within a setup based on collinear factorization. The same observation applies to the treatment of ultraviolet divergences and their renormalization. While this behavior was to be expected and indeed constitutes a necessary requirement, it provides a non-trivial check on the correctness of our result. While we confirm earlier results in the literature which state that unpolarized and linearly polarized gluon TMD agree with each other in the dilute i.e. BFKL limit, we find that both distributions differ at NLO, which is directly related to the non-trivial tensor structure of the real NLO corrections.
As a next step we clarified further the relation between BFKL and CS evolution and clarified that they describe both evolution of the system in rapidity, i.e. BFKL of the cross-section and therefore directly related to the hard final state, while CS evolution rapidity evolution of the soft system. Both evolution equation are therefore not independent and care is needed to avoid over-counting. Unlike previous calculations based on the CGC framework, our study further enabled us to recover the finite term proportional to β 0 in the anomalous dimension of the TMD PDFs.
As an important side result of our study we find that the real NLO contribution of the unpolarized gluon TMD yield precisely the off-shell TMD gluon-to-gluon splitting function, determined in [34].
While [34] determined this splitting function from a diagrammatic approach -essentially requiring simultaneous fulfillment of collinear and high energy limit while imposing gauge invariant production vertices -the current study obtains the same result from the QCD operator definition of gluon TMDs. It therefore establishes an important link between both frameworks, which will be of importance to continue with these efforts. In particular the current study provides a possibility to finally determine the still missing virtual corrections to these splitting kernels, and to formulate corresponding evolution equations.
Apart from these efforts, future studies should investigate phenomenological consequences of the derived results, which now allow to use the complete CSS resummation formulation to resum double and single logarithmic contributions with the methods of the renormalization group, extending previous results in the literature such as [22,27,92,93]. Another direction of research should address the inclusion of high density effects, along the lines of [22], but including the complete treatment of factorization scheme dependence established in this paper, asd well as contributions due to quarks.