QCD evolution of (un)polarized gluon TMDPDFs and the Higgs $q_T$-distribution

We provide the proper definition of all the leading-twist (un)polarized gluon transverse momentum dependent parton distribution functions (TMDPDFs), by considering the Higgs boson transverse momentum distribution in hadron-hadron collisions and deriving the factorization theorem in terms of them. We show that the evolution of all the (un)polarized gluon TMDPDFs is driven by a universal evolution kernel, which can be resummed up to next-to-next-to-leading-logarithmic accuracy. Considering the proper definition of gluon TMDPDFs, we perform an explicit next-to-leading-order calculation of the unpolarized ($f_1^g$), linearly polarized ($h_1^{\perp g}$) and helicity ($g_{1L}^g$) gluon TMDPDFs, and show that, as expected, they are free from rapidity divergences. As a byproduct, we obtain the Wilson coefficients of the refactorization of these TMDPDFs at large transverse momentum. In particular, the coefficient of $g_{1L}^g$, which has never been calculated before, constitutes a new and necessary ingredient for a reliable phenomenological extraction of this quantity, for instance at RHIC or the future AFTER@LHC or Electron-Ion Collider. The coefficients of $f_1^g$ and $h_1^{\perp g}$ have never been calculated in the present formalism, although they could be obtained by carefully collecting and recasting previous results in the new TMD formalism. We apply these results to analyze the contribution of linearly polarized gluons at different scales, relevant, for instance, for the inclusive production of the Higgs boson and the $C$-even pseudoscalar bottomonium state $\eta_{b}$. Applying our resummation scheme we finally provide predictions for the Higgs boson $q_T$-distribution at the LHC.


Introduction
Observables sensitive to the transverse momentum of quarks and gluons inside a hadron have a long theoretical and experimental history. They have proven to be valuable tools to test the QCD dynamics at high-energy colliders, extending the information provided by observables integrated over the intrinsic transverse momenta. At large transverse momentum these observables can be computed in perturbation theory, but if the transverse momentum q T is much smaller than the probe of the hard reaction Q, then large logarithms of their ratio appear and resummation becomes a must in order to obtain reliable results. This issue was already addressed in the eighties by Collins, Soper and Sterman [1].
The main hadronic quantities in observables at q T Q are the transverse momentum dependent functions (TMDs), first considered by Ralston and Soper [2,3] and by Collins

JHEP07(2015)158
and Soper [4,5]. The TMDs represent, generally speaking, the probability of finding a parton inside a hadron with a definite transverse momentum, i.e., TMD parton distribution functions (TMDPDFs); or the probability that a quark or gluon fragments into a hadron with a given transverse momentum (TMDFFs). They play an important role in the rich phenomenology of azimuthal and spin asymmetries (see, e.g., [6,7]).
After the pioneering works, much effort has been devoted to properly describe the polarization of the partons/hadrons, the universality of TMDs and other relevant properties. However, the "naive" (old) definitions introduced in [4,5] and considered in subsequent works, suffer from undesired features preventing them from properly represent physical hadronic quantities, such as uncancelled rapidity divergences. Recently Collins [8] and Echevarria-Idilbi-Scimemi [9,10] have revisited and updated the definition of quark TMDs, making it consistent with a generic factorization theorem and free from the bad features. Having at our disposal the proper definition for such quantities allows us to better deal with physical processes where they appear and from which we want to extract sensible information on the hadron structure. Thus, it is the goal of the present work to extend those efforts to the gluon TMDs, relevant for instance in processes such as Higgs boson and quarkonium production in hadron-hadron collisions.
In this work we pay special attention to three of the eight leading-twist gluon TMD-PDFs, calculate them explicitly at next-to-leading order (NLO) and demonstrate that the rapidity divergences cancel in their proper definitions. On one hand, the distributions of unpolarized (f g 1 ) and linearly polarized (h ⊥g 1 ) gluons inside an unpolarized hadron, and, on the other hand, the gluon helicity TMDPDF (g g 1L ), which represents the distribution of longitudinally polarized gluons inside a longitudinally polarized hadron. The calculation not only supports the definitions introduced in this work, but also allows us to extract valuable perturbative ingredients to resum large logarithms and better control their nonperturbative parts, eventually improving our description of experimental data. We emphasize that the calculation of g g 1L is done for the first time, while for f g 1 and h ⊥g 1 one could combine previous results and then carefully recast them into the new TMD formalism.

JHEP07(2015)158
The evolution of the gluon TMDPDFs, as in the case of quark TMDs [34], turns out to be universal, i.e., the same evolution kernel describes the evolution of any of the leading-twist (un)polarized gluon TMDPDFs. It is interesting to contrast this finding with the evolution of the parton distribution functions (PDFs) and double parton distributions (DPDs) which have vastly different evolution depending on the polarization (see, e.g., [35] for a direct comparison for DPDs). The currently known perturbative ingredients allow us to use the evolution equations to resum the large logarithms up to next-to-next-to-leadinglogarithmic (NNLL) accuracy. Moreover, if we consider the perturbative coefficients of the operator product expansion (OPE) of those TMDs at large transverse momentum, also some parts of them turn out to be universal. Exploiting this feature, we introduce a further step to resum the large logarithms that appear in the OPE coefficients, exponentiating the double logarithms and improving the convergence of the resummation. Thus we provide a general framework to deal with (un)polarized gluon TMDPDFs in different processes and account for their perturbative and non-perturbative contributions.
Drawing attention to the distribution of linearly polarized gluons inside an unpolarized hadron, several works have addressed their role at the LHC (see, e.g., [36][37][38][39][40][41]). In particular, in [41] the authors quantified their contribution in the context of the TMD formalism, both for Higgs boson and C-even scalar quarkonium (χ c0 and χ b0 ) production. In the present work we extend their efforts by implementing the currently known perturbative ingredients to the full extent to perform the resummation at NNLL accuracy, providing more accurate predictions and discussing their uncertainty.
The paper is organized as follows. In section 2 we apply the SCET machinery to derive the factorization theorem for the Higgs q T -distribution in polarized hadron-hadron collisions in terms of well-defined gluon TMDPDFs. In section 3 we discuss the QCD evolution of all the leading-twist gluon TMDPDFs, which turns out to be driven by a universal evolution kernel. Next, in section 4 we address the refactorization of TMDPDFs in terms of collinear functions, which applies when the transverse momentum is in the perturbative domain. In section 5 we consider the gluon helicity TMDPDF (g g 1L ), that accounts for longitudinally polarized gluons inside a longitudinally polarized hadron, and perform a numerical study of the function itself and the impact of evolution. In section 6 we analyze the TMDPDFs that contribute in unpolarized hadron-hadron collisions, i.e., unpolarized and linearly polarized gluons (f g 1 and h ⊥g 1 respectively), and give some estimates of their relative contributions at different scales. Then, in section 7 we study the Higgs boson transverse momentum distribution, paying special attention to the role played by linearly polarized gluons and the non-perturbative effects. Finally, conclusions are drawn in section 8.

Factorization theorem in terms of well-defined TMDPDFs
Below we derive the factorization theorem for the Higgs q T -distribution in polarized hadronhadron collisions, A(P, S A ) + B(P , S B ) → H(m H , q T ) + X, by performing a set of consecutive matchings between different effective field theories, relevant at each scale:

JHEP07(2015)158
In the first step we integrate out the top quark mass, m t , to build an effective ggH coupling. In the second matching we integrate out the mass of the Higgs boson, m H , and obtain a factorized cross-section in terms of well-defined gluon TMDPDFs, which holds for q T m H . Those gluon TMDPDFs will be expressed in terms of fundamental hadronic matrix elements. Finally, in the region Λ QCD q T m H , we can further refactorize the gluon TMDPDFs in terms of the collinear gluon/quark PDFs, integrating out the large scale q T .
Before discussing the steps in the derivation of the factorization theorem, we introduce the notation used through the paper. A generic vector v µ is decomposed as v µ =n · v n µ 2 + n · vn , with n = (1, 0, 0, 1),n = (1, 0, 0, −1), n 2 =n 2 = 0 and n ·n = 2. We also use v T = |v ⊥ |, so that v 2 ⊥ = −v 2 T < 0. The production of the Higgs boson through gluon-gluon fusion is well approximated by the effective local interaction [42][43][44][45][46] where α s (µ) is the QCD coupling at factorization scale µ, F µν,a the gluon field strength tensor, H is the Higgs field and v ≈ 246 GeV is the Higgs vacuum expectation value. The explicit expressions for the Wilson coefficient C t and its evolution can be found in appendix E. Using the effective lagrangian just introduced, the differential cross section for Higgs production is factorized as where s = (P +P ) 2 . This expression manifests the first step in the matching procedure, where we integrate out the large top quark mass through the perturbative coefficient C t . The n f = 5 effective QCD operator is next matched onto the SCET-q T one: where q 2 = m 2 H and the B ⊥µ n(n) operators, which stand for gauge invariant gluon fields, are given by The collinear and soft Wilson lines are path ordered exponentials Wilson lines with calligraphic typography are in the adjoint representation, i.e., the color generators are given by (t a ) bc = −if abc . In order to guarantee gauge invariance among JHEP07(2015)158 regular and singular gauges, transverse gauge links need to be added (as described in [47,48]). In this work we stick to Feynman gauge for perturbative calculations, and thus transverse gauge links do not play any role. The Wilson matching coefficient C H (−q 2 , µ), which corresponds to the infrared finite part of the gluon form factor calculated in pure dimensional regularization, is given at one-loop by In appendix D we report our explicit NLO calculation, and appendix E gives the evolution and higher order contributions. We anticipate here that for the phenomenological study discussed in section 7 we perform the so-called "π-resummation", which consists on choosing the scale in the above coefficient as µ 2 = −q 2 . In this way the convergence of the hard part is improved. See [49,50] for more details. After some standard algebraic manipulations and a Taylor expansion in order to retain the leading order contribution in q T /m H , the cross-section can be written as T )/s and y is the rapidity of the produced Higgs boson. The Born-level cross section is (2.8) The pure collinear matrix elements and the soft function are defined as Notice that, in order to avoid double counting, one needs to subtract the contribution of soft momentum modes (the so-called "zero-bin" in the SCET nomenclature) from the naively calculated collinear matrix elements, thus obtaining the "pure collinear" matrix elements, denoted by the superscript (0) (see, e.g., [51] for more details). Note that in eq. (2.7) we have applied the SCET machinery to decouple the collinear, anticollinear and soft modes, removing the interactions between them from the Lagrangian [21]. Thus the JHEP07(2015)158 factorization of any operator in SCET is straighforward, as well as the final states |X , which can be written in the factorized form |X = |X n ⊗ |Xn ⊗ |X s , describing the collinear, anticollinear and soft states. As shown in appendices A, B and C by performing an explicit NLO perturbative calculation of unpolarized, linearly polarized and helicity gluon TMDPDFs, the collinear and soft matrix elements defined above are individually ill-defined, since they contain uncancelled rapidity divergences. We stress the fact that this issue is independent of the particular regulator used. Thus, we need to combine them in a certain way to cancel these divergences and obtain well-defined hadronic quantities. Based on the work done in [9,10,34] for the quark case and using η n(n) to label generic parameters that regulate the rapidity divergences present in the (anti-)collinear and soft matrix elements, then the TMDPDFs are defined as where ζ A,B are auxiliary energy scales, the twiddle labels the functions in impact parameter space (IPS) and we have split the soft function in rapidity space as The soft function can be split to all orders in perturbation theory, regardless which particular regulator is used, following the same logic as in [10]. In that work this fundamental property was proven for the soft function relevant for quark TMDs. The proof for the soft function that appears in the gluon TMDs follows analogously, simply changing the color representation from the fundamental to the adjoint. The arbitrariness in the choice of the rapidity cutoff to split the soft function, which is not explicitly shown in eq. (2.11), manifests itself as the appearance of the auxiliary energy scales ζ A and ζ B , which are bound together by ζ A ζ B = q 4 = m 4 H . Resorting to the ∆-regulator for definiteness, the soft function is split as where we have explicitly shown the dependence on the regulator parameters and ζ A = m 2 H /α and ζ B = αm 2 H , with α an arbitrary boost invariant real parameter. We point out that an explicit dependence on q 2 = m 2 H has been added as well in the soft function. This is due to the use of the ∆-regulator, which induces such dependence, in consistency with the fact that the soft function represents the cross-talking between the two collinear sectors. On the contrary, pure collinear matrix elements do not have any remnant information about the opposite collinear sector, and thus cannot depend on q 2 = m 2 H .

JHEP07(2015)158
Now, using eq. (2.12) the gluon TMDPDFs are defined as These hadronic quantities are free from rapidity divergences, i.e., they have well-behaved evolution properties and can be extracted from experimental data. As already mentioned, when one performs the perturbative calculations of the collinear matrix elements, depending on the particular regulator that is used, the issue of double counting the soft modes arises. In this case, one needs to subtract, on a diagram-by-diagram basis, the soft limit of each collinear contribution (the "zero-bin"). With the ∆-regulator one can show order by order in perturbation theory that the subtraction of the zero-bin for each collinear matrix element is equivalent to divide it by the soft function: . (2.14) The naively calculated collinear matrix elements, which do not have the label (0) anymore, depend on the hard scale q 2 = m 2 H and the two regulator parameters ∆ ± , in contrast with the pure collinear matrix elements. The latter should depend only on the regulator that belongs to each collinear sector. The naive collinear matrix elements depend as well on the spurious regulator for the other sector. Thus, for the particular case of this regulator we can define the TMDPDFs as Notice that the spurious regulator in the naive collinear matrix elements is now cancelled by dividing them by the proper piece of the soft function, thus recovering the correct regulator dependence in the TMDPDFs as in eq. (2.13). It is worth emphasizing that the proper definition of TMDPDFs in eq. (2.10) is independent of the particular regulator used. We have given their definition using the ∆-regulator, although one could conveniently modify eq. (2.13) in order to use, for instance, the rapidity regulator introduced in [52]. One could also follow the lines of [8], where no regulator is used, i.e., the rapidity divergences among the collinear and soft matrix elements are cancelled by combining (before integration) the integrands of the relevant Feynman diagrams order by order in perturbation theory. Now, we can write the cross-section for the Higgs q T -distribution in terms of welldefined gluon TMDPDFs:

JHEP07(2015)158
This factorized cross-section is valid for q T m H . In the next section we discuss the refactorization of the TMDPDFs in terms of collinear functions when the transverse momentum is a perturbative scale, i.e., when Λ QCD q T m H .
Before moving to the next factorization step, we first need to consider the dependence on the hadron spin and separate the unpolarized (U), longitudinally polarized (L) and transversely polarized (T) situations. In [53] the authors obtained the decomposition of collinear correlators at leading-twist. We emphasize the fact that these correlators suffer from rapidity divergences and thus cannot be considered well-defined hadronic quantities. The decomposition, however, is not directly affected by this issue and follows equivalently. Below, given the proper definition of gluon TMDPDFs in eq. (2.10), we extend the decomposition of [53] and write The functions f g 1 , h ⊥g 1 , g g 1L and g g 1T are T -even, while the rest are T -odd. In sections 6 and 5 we will pay special attention to the functions f g 1 , h ⊥g 1 and g g 1L , calculating them explicitly at NLO to show that they are free from rapidity divergences and to obtain the necessary perturbative ingredients to perform the resummation of large logarithms. These three functions are the only TMDPDFs which are matched onto leading twist collinear matrix elements, i.e., the canonical PDFs.
The Wilson line structure in the operator definition of the TMDPDFs gives rise to calculable process dependence. In the types of processes considered here, where two gluons fuse into a color singlet, the Wilson lines are past pointing. In a process where color would flow into the final state, also future pointing Wilson lines play a role. Functions with different Wilson line structure differ by matrix elements containing intrinsically nonlocal gluonic pole contributions [54][55][56][57]. Depending on the number of such gluonic poles being odd or even, the functions are T-odd or T-even. The functions come with specific processdependent gluonic pole factors that can lead to a breaking of universality, in the simplest cases giving rise to a sign change, such as the Sivers function having a different sign in Drell-Yan and in deep-inelastic scattering (DIS) processes. Other functions, such as h ⊥g 1 need to be written as a linear combination of two or even more functions, with the coefficient in the linear combination depending on the Wilson lines and in turn on the color flow in the process [58,59]. However, for both the gluon-gluon fusion into a color singlet considered here, as well as in the "gluon initiated DIS", exactly the same linear combination contributes to the cross section.

JHEP07(2015)158
Finally, we provide the equivalent of eq. (2.17) in IPS, since as we show next, the evolution of TMDPDFs is done in that space. With the Fourier transform given bỹ where for a generic function f (k T ) we represent Worth noting is that while G µν g/A andG µν g/A are each others Fourier transforms, this does not hold true for the individual gluon TMDs which have factors of k n⊥ (b ⊥ ) in the decomposition (e.g., h ⊥g 1 andh 1T ).

Evolution of gluon TMDPDFs
The TMDPDFs defined in eq. (2.10) depend on two scales: the factorization scale µ and the energy scale ζ (related to the rapidity cutoff used to separate the two TMDPDFs). Thus, their evolution kernel is such that it connects these two scales between their initial and final values. Below we derive first the part of the kernel that allows us to evolve the TMDPDFs with respect to µ, and then the one that corresponds to ζ. The evolution of (un)polarized gluon TMDPDFs in terms of the renormalization scale µ is governed by the anomalous dimensions: The renormalization group (RG) equation applied to the factorized cross-section in eq. (2.16) implies the following relation among the different anomalous dimensions:

JHEP07(2015)158
where the anomalous dimension of the coefficients H and C t , γ H and γ t respectively, are given in appendix E. Thus where the non-cusp piece is In the equation above γ g is the non-cusp piece of the anomalous dimension of the hard coefficient C H (see appendix E). It should be mentioned that the splitting of γ H into the two anomalous dimensions γ G given in eq. (3.3) is unique following the restriction ζ A ζ B = m 4 H . The coefficients of the perturbative expansions of Γ cusp and γ nc are known up to three loops and they are collected in appendix E. Now we focus our attention on the evolution in terms of the scale ζ. Following the arguments in [10], one can show that the soft function relevant for gluon TMDs can to all orders be written as with a function R s , depending only on b T and µ, and D g related to the cusp anomalous dimension in the adjoint representation by Given eqs. (2.13) and (3.5), one obtains the following evolution equations in ζ: Notice that the same D g term drives the ζ evolution for all gluon TMDPDFs, since the soft function that enters into their definition and gives the entire ζ evolution is spin-independent. The coefficients of the perturbative expansion of the D g term can be completely obtained from the calculation of the soft function. If we write then the first two coefficients are:

JHEP07(2015)158
The finite coefficients d n (0) cannot be determined by eq. (3.6), but, as already mentioned, by a perturbative calculation of the soft function (or the cross-section in full QCD). The coefficient d 1 (0) can be easily extracted from the NLO calculation of the soft function in appendix A, and one gets d 1 (0) = 0. The coefficient d 2 (0) can be obtained from the soft function relevant for DY or SIDIS processes [60] by using the Casimir scaling, i.e., rescaling it by C A /C F (see also [26]): At small b T the D g term can be calculated perturbatively, but at large b T it has to be modelled and extracted from experimental data. However, we can extend the fact that the soft function is universal and spin-independent to this non-perturbative piece, and it can therefore be used to parametrize the non-perturbative contribution to the evolution of all (un)polarized TMDPDFs. Regardless of how the non-perturbative contribution to the D g term is parametrized, we can perform the evolution of all leading-twist gluon TMDPDFs consistently up to NNLL (given the currently known perturbative ingredients, i.e., γ G and D g ): where the evolution kernelR g is given bỹ Solving analytically the evolution equation of the D g term in the small b T region, where µ b = 2e −γ E /b T is the natural scale of the D g term, and implementing the running of the strong coupling consistently with the resummation order, one obtains (see [60] for quark TMDs) (3.14) In this result we have defined a = α s (µ i )/(4π) and X = aβ 0 L T . The β i and Γ A i coefficients are given in appendix E.
As a final remark, we emphasize the fact that the evolution kernel in eq. (3.12) is valid only in the perturbative region of small b T Λ −1 QCD , since the perturbative expression of D g (even its resummed version) breaks down at large b T [60,61].

JHEP07(2015)158 4 Refactorization of TMDPDFs and resummation of large logarithms
As already anticipated, when the transverse momentum is perturbative, we can perform an operator product expansion (OPE) of the TMDPDFs in terms of collinear functions, integrating out the transverse momentum by means of Wilson coefficients. Depending on the particular TMDPDF considered, the collinear functions that will describe its perturbative small-b T region will be different, and also the relevant Wilson coefficients. However, the part of the Wilson coefficients originating from the evolution, which is universal and spinindependent, will be common for all TMDPDFs. Below we give the general expressions for these OPEs, and in the next sections we explicitly calculate the coefficients for the relevant TMDPDFs in the cases of an unpolarized and longitudinally polarized hadron.
For b T Λ −1 QCD we can refactorize the (renormalized) gluon TMDPDFs of a hadron A in terms of (renormalized) collinear quark/gluon distributions: The convolution refers to momentum fraction x for TMDPDFs that are matched onto twist-2 collinear functions (like f g 1 ), while in the case of TMDPDFs that are matched onto twist-3 functions (like the gluon Sivers function f ⊥g 1T ) it would represent a two-dimensional convolution in the two momentum fractions of the collinear function. In the equation above we have thus represented schematically the OPE of any TMDPDF, whereF g/A stands for any of the functions in eq. (2.19) and f j/A the adequate collinear functions in each case. For example, we could consider the unpolarized gluon TMDPDFf g 1 and match it onto the unpolarized collinear gluon/quark PDFs, as shown in section 6; or we could consider the Sivers functionf and match in onto gluon/quark twist-3 collinear functions [62]. The coefficientsC g/j are different for each TMDPDF.
The natural scale for the coefficientsC g/j is µ ∼ 1/b T ∼ q T , which is the large scale that we integrate out when we perform the OPE. Thus, we can choose to set the resummation scale either in impact parameter space or in momentum space. In the following we discuss these two approaches in more detail.

Resummation in impact parameter space
If we perform the resummation of large logarithms in impact parameter space then the resummed TMDPDF is written as: where ζ 0 ∼ µ 2 b and µ 0 ∼ µ b . The superscript P ert signifies that it is only the perturbative part of the TMDPDFs, 1 valid at small b T 1/Λ QCD . Notice that the functions D g

JHEP07(2015)158
is universal and spin independent, and thus it is the same for any of the TMDPDFs in eq. (2.19). On the contrary, as already mentioned, the coefficientsC g/j are specific for each TMDPDF, as are the collinear functions f j/A which generate the perturbative tail for each TMDPDF at small b T . So far we have addressed the TMDPDFs in the perturbative region. For large b T we need to model them and extract them from experimental data. To do so, one could implement a smooth cutoff that freezes the perturbative contribution slowly as b T gets larger:F where the cutoff prescription could be, for instance: with n an integer number and b c the parameter that determines the separation between the perturbative and non-perturbative regions. For small b T the perturbative contribution dominates, and gets frozen as we increase the b T , sinceb T → b c for large b T . The nonperturbative modelF N P is constrained to be 1 for b T = 0 and plays an increasingly important role as we increase b T .

Resummation in momentum space
Instead of setting µ 0 ∼ µ b , in this case we keep it in momentum space. In this way we avoid hitting the Landau pole in the strong coupling, and we write the TMDPDF as where ζ 0 = C 2 ζ µ 2 b and µ 0 ∼ q T . We have kept explicitly the dependence on the real parameter C ζ , which will be used later on to test the dependence of the results on the rapidity scale, basically varying it between 1/2 and 2. This is because the way it enters in the final resummed expression for the TMDPDF is subtle, contrary to the scale µ 0 , which can be easily identified where it appears. Now the coefficientsC g/j contain large logarithms ln(µ 0 b T ) which are not minimized by the choice µ 0 ∼ q T (this was the case in the previous subsection, when choosing µ 0 ∼ 1/b T ). However we can further split them by using their RG-equation JHEP07 (2015)158 where P i/j (x/z) are the usual DGLAP splitting kernels, so that double logarithms can be partially exponentiated (see [63] for the quark case): Choosing h Γ(γ) (b T ; µ b ) = 0, the first few coefficients for the perturbative expansions of h Γ(γ) are: After the various steps we have performed, the OPE of gluon TMDPDFs can be rewritten as The functions D g , h Γ and h γ above still contain large logarithms L T that need to be resummed when α s L T is of order 1 (also the coefficientsĨ g/j ). We have already calculated the resummed expression for D R g in the previous section, and following the same procedure we can derive the resummed expressions for the terms h Γ(γ) . Let us consider their evolution equations, eq. (4.8):

JHEP07(2015)158
Notice that we have chosen h Γ(γ) (b T ; µ b ) = 0. By expanding the β-function and re-writing α s (µ b ) in terms of α s (µ) at the proper order, as shown in [60] for the quark TMDPDFs case, we can solve these equations and get and where again a s = α s /(4π) and X = a s β 0 L T . The resummed expressions we have just found are valid only in the perturbative region of small b T . With them, we can finally write the resummed TMDPDFs as: Notice that the functions D R g , h R Γ and h R γ are universal and spin independent, and thus are the same for any of the TMDPDFs in eq. (2.19). On the contrary, the coefficients I g/j are specific for each TMDPDF, as are the collinear functions f j/A which generate the perturbative tail for each TMDPDF at small b T .
Finally, in order to parametrize the non-perturbative contribution at large b T we overlap the perturbative expression in eq. (4.14) with a non-perturbative model and writẽ expression applies, and play an increasingly important role as b T gets larger. Moreover, it should be such that it cancels the contribution that comes from the perturbative expressioñ F Pert g/A in the large b T region, where it does not apply. In simple terms, there is no problem with extending the perturbative expression to the whole b T -space, since the model is used to correct for it in the non-perturbative region. This approach was used in [63] to perform a global fit of Drell-Yan and Z-boson data and extract the unpolarized quark TMDPDFs.
As a final remark, we emphasize the fact that the non-perturbative model in this case is not the same as in the previous subsection. They parametrize the non-perturbative region in a different way, depending on what is the treatment of the perturbative contribution, and thus they can be different.

Gluon helicity TMDPDF
The gluon helicity TMDPDF, g g 1L , represents the distribution of longitudinally polarized gluons inside a longitudinally polarized hadron. In appendix C we perform for the first time an explicit NLO calculation of this quantity and show that, if defined as in eq. (2.10), then the rapidity divergences cancel among the collinear and soft matrix elements. We also perform a NLO calculation on the collinear gluon helicity PDF g j/A (x; µ), which we use to extract the OPE Wilson coefficient of gluon helicity TMDPDF: where the longitudinally polarized collinear quark and gluon PDFs are defined as with ξ n(n) the (anti)collinear fermion field. The result of the coefficient is analogous to the one for the unpolarized gluon TMDPDF, as we will show in next section, apart from the DGLAP splitting kernel (similar result was found in the case of unpolarized and helicity quark TMDPDFs in [64]). It reads where the one-loop DGLAP splitting kernels (collected for all polarizations in [65]) are  Table 1. Perturbative orders in logarithmic resummations, both for for resummations in momentum space and in impact parameter space.
In order to illustrate the QCD evolution of gluon helicity TMDPDF we choose the resummation scale in impact parameter space. We also set µ 2 = ζ = Q 2 , use the evolution kernel in eq. (3.12) and separate the perturbative and non-perturbative contributions in a smooth way as in eq. (4.3). Thus the gluon helicity TMDPDF is given bỹ and we implement a simple non-perturbative model The parameters λ g and λ Q have never been extracted from experimental data, and thus we can only guess their values and give predictions by varying them in a reasonable range. What we know is that λ Q is the same among all (un)polarized gluon TMDPDFs, because it parametrizes the scale-dependent part of the non-perturbative model, which is related to the large-b T tail of the universal D g function. Notice that for simplicity we have neglected any x dependences in the non-perturbative model. The gluon helicity TMDPDF is shown in figure 1 at x = 0.01, Q = 20 GeV, with the non-perturbative parameters b c = 1.5 GeV −1 , λ g = 0.3, λ Q = 0.1. We use the latest available parametrizations of collinear gluon and quark helicity PDFs [66] at NLO for our numerical analysis. The running of the strong coupling is implemented at NNNLO with the MSTW routine [67], with a variable flavor number scheme with m c = 1.4 GeV and m b = 4.75 GeV. The input value for the strong coupling is set to α s (M Z ) = 0.1185, and we impose a lower cutoff for the running scale µ such that it never goes below 1 GeV. The bands come from varying both the rapidity and the resummation scales by a factor of 2 around their default value, and keeping the largest variation for each point in k T . It is clear that the theoretical uncertainty gets reduced as we increase the resummation accuracy by including more perturbative ingredients, as schematically illustrated in table 1. In figure 2 we show the gluon helicity TMDPDF at x = 0.01, for different values of the energy scale and the non-perturbative parameters, at NNLL accuracy. In order to be consistent with the factorization theorem we have cut the curves so that the condition k T Q is fulfilled. As can be seen, the larger the scale, the wider is the distribution. Moreover, the larger the value of the non-perturbative parameters, the smaller the helicity distribution is at low transverse momentum. It is interesting to notice that choosing small and equal non-perturbative parameters, gives a qualitatively different gluon helicity TMDPDF at low scales, as in the upper-left panel.
These results have been obtained using the OPE coefficients that have been calculated for the first time in the present work. They are an important perturbative ingredient that will allow us to better fix the non-perturbative parameters with new measurements, that could be performed at RHIC or at the future AFTER@LHC or EIC with longitudinally polarized hadron beams.

Gluon TMDPDFs in an unpolarized hadron
As shown in eq. (2.19) there are two gluon distributions that contribute at leading-twist in the case of an unpolarized hadron: the unpolarized (f g 1 ) and the linearly polarized (h ⊥g 1 ) ones. The latter was introduced in [53] and implemented for the first time in the resummation of gluon-gluon fusion process in impact parameter space in [18,19]. In appendices A and B we perform an explicit NLO calculation of those distributions using their proper definition in eq. (2.13), and we show that they are free from rapidity divergences when the collinear and soft matrix elements are combined properly.
As described in the previous section for the gluon helicity TMDPDF, both f g 1 and h ⊥g  are just the unpolarized collinear gluon/quark PDFs: where the unpolarized collinear PDFs are defined as Note that the TMDPDFs for the unpolarized gluon and the linearly polarized gluon are both matched onto the same PDF, but the first non-zero order of the matching coefficient for the linearly polarized gluon is one order higher in α s than for the unpolarized gluon.

JHEP07(2015)158
In appendices A and B we obtain their matching coefficients at NLO by subtracting the collinear PDFs at the same order. Moreover, in section 4 we have shown that the OPE coefficients for TMDs can be further refactorized, and thus the previous OPEs can be written, setting µ 2 = ζ = Q 2 and using the evolution kernel in eq. (3.12), as The perturbative coefficientsC f,h g/j are given in appendices A and B, the one-loop DGLAP splitting kernel P g/g is given in eq. (5.4) and The contribution of unpolarized and/or linearly polarized gluon distributions in unpolarized hadron-hadron collisions depends on the process under study and has been discussed in several works [36,[38][39][40][41]. In this work we focus on the production of Higgs boson and C-even pseudoscalar bottonium state η b [30], since for the production of P -wave quarkonium states (like χ b0 ) there are arguments that suggest a breaking of the factorization [68]. In the considered cases, Higgs boson and η b production, both unpolarized and linearly polarized distributions play a role, and thus one can investigate their relative contribution to the cross-section. We use our results to quantify the contribution of linearly polarized gluons, considering the following ratio: 2 where the numerator and denominator are the two terms in the factorized cross section which determine the relative contribution from linearly polarized and unpolarized gluons to the cross section, for both Higgs boson and C-even pseudoscalar bottonium production. In order to compute this quantity we will insert the TMDs as in eq. (6.3), choosing ζ 0 ∼ µ 2 0 ∼ µ 2 b and using theb T prescription to separate the perturbative from non-perturbative contributions as in eq. (5.6). The latter will be parametrized as:  similar to the model used previously for the helicity TMDPDF. Notice that the parameter λ Q is the same for both functions, since the evolution is universal among all (un)polarized TMDPDFs, that is, their scale-dependence is the same. Having precise estimates for this ratio will help us predict the measurability of both unpolarized and linearly polarized gluon distributions in a given process (or scale), which is the final goal. Using our resummation scheme and the perturbative ingredients at the highest possible order we provide accurate predictions for this quantity.
In figure 3 we show our results for the ratio R at the relevant scales for the transverse momentum distributions of Higgs boson and η b , all at NNLL accuracy. We used the MSTW08nnlo set [67] and selected different values for the non-perturbative parameters in order to check their impact on the result. The running of the strong coupling is implemented at NNNLO with the MSTW routine, with a variable flavor number scheme with m c = 1.4 GeV and m b = 4.75 GeV. The input value for the strong coupling is set to JHEP07(2015)158 α s (M Z ) = 0.1185, and we impose a lower cutoff for the running scale µ such that it never goes below 0.4 GeV. Comparing our results to the ones presented in [41], we have included the contribution of quark PDFs to the collinear expansion of gluon TMDPDFs (through C f,h g/q in eq. (6.3)) and also higher order perturbative ingredients, performing the resummation consistently at NNLL accuracy. The uncertainty bands also represent an improvement with respect to the results in [41]: they allow us to better quantify what is the effect of non-perturbative contributions relative to the scale uncertainty, and whether experimental data can be used to determine the non-perturbative parameters or distinguish between different models of the non-perturbative input. The bands are obtained by independently varying the scales ζ 0 and µ 0 around their default value by a factor of 2, and plotting the maximum uncertainty for each point in q T .
In order to estimate the impact on the ratio of the different non-perturbative parameters, we have chosen several values in a sensible range and selected some combinations in limiting cases. First, the parameters should be positive, since the gluon distributions are supposed to vanish at large b T . Second, given the values found for similar models in the case of quark TMDPDFs (see, e.g., [69]), we have chosen a maximum value of 1. Finally, since our goal is to estimate the contribution of linearly polarized gluons as compared to unpolarized ones, we have chosen the following limiting cases: (i) λ Q = 0.01, λ f = λ h = 0.01. Small evolution parameter λ Q and similar and small parameters λ f and λ h .
The outcome of the numerical study is clear: the lower the scale the more contribution we have from linearly polarized gluons, although this contribution depends on the value of the non-perturbative parameters, which will have to be fixed by fitting experimental data. At the Higgs boson scale the effect of linearly polarized gluons is small, around 1-9%, making it harder to extract their non-perturbative parameters from experimental data. At lower scales, as in the production of η b , their role is enhanced, from 10% up to 70%, and thus experimental data can better determine them. However, it seems plausible that their non-perturbative parameters could be fixed in the near future by properly combining experimental data for different experiments and at different scales. Thus the framework introduced in this paper, with the proper definition of gluon TMDPDFs and their QCD evolution, will be crucial in order to consistently address different processes in terms of the same hadronic quantities and properly extract their non-perturbative parameters.

Higgs boson q T -distribution
After analyzing the contribution of linearly polarized gluons for η b and Higgs boson production in unpolarized hadron-hadron collisions, we apply our results to provide some predic-

JHEP07(2015)158
tions for the Higgs boson transverse momentum distribution at the LHC. The cross-section for this process can be easily obtained from eq. (2.16) if we consider unpolarized protons: It is well-known that the evolution kernel suppresses the TMDPDFs at large b T , and that this effect is enhanced the larger the relevant hard scale Q is, in this case m H [60] (see also, e.g., the discussion in [70] in the context of the Collins-Soper-Sterman approach). Therefore, the larger the Q the more insensitive is the resummed expression to nonperturbative contributions at large b T . Based on this, we fix the resummation scale in momentum space, µ 0 = Q 0 + q T (with Q 0 = 2 GeV), and writẽ where ζ 0 ∼ µ 2 b . The perturbative coefficientsĨ f (h) are derived from the results in appendices A and B:Ĩ We parametrize the TMDPDFs as in eq. (4.15), which allows us to exploit the perturbative results without using any prescription, like theb T . This procedure was already used in [63] to perform a global fit of Drell-Yan data, to obtain the non-perturbative parameters of unpolarized quark TMDPDFs. Following the same procedure, and leaving some room for small non-perturbative effects, we multiply the TMDPDFs in eq. (7.2) by the following non-perturbative models: where we have neglected the dependence on x and Q for simplicity.

JHEP07(2015)158
The resummation of large logarithms in the cross-section in eq. (7.1) is done by evaluating each perturbative coefficient at its natural scale and then evolving them up to a common scale by using their relevant anomalous dimensions. For the TMDPDFs the resummation was already discussed before and led to eq. (7.2). The natural scale for the coefficient C t is µ t ∼ m t , and its evolution is presented in appendix E. For C H (−q 2 , µ) (remember that H = |C H | 2 ) it was discussed in [50] that the choice µ 2 H ∼ −m 2 H leads to a better convergence of the resummed expression, and thus we apply that procedure in our numerical study, using the anomalous dimensions that appear in appendix E. At NLL accuracy we take the coefficients at LO, and at NNLL accuracy at NLO, combined with the TMDPDFs at the corresponding order given in table 1.
In figure 4 we show the Higgs boson transverse momentum distribution at √ s = 8 TeV, for different values of the non-perturbative parameters and both at NLL and NNLL accuracies. We have chosen β h = β f for simplicity, given that in the previous section we showed that the impact of linearly polarized gluons at the Higgs boson scale is small. The choices β f = β h = 0 and β f = β h = 1 give the most extreme scenarios, where the total non-perturbative contribution of both functions is zero or large. We used the MSTW08nnlo set [67] for the input PDFs. The running of the strong coupling is implemented at NNNLO with the MSTW routine, with a variable flavor number scheme with m c = 1.4 GeV and m b = 4.75 GeV. The input value for the strong coupling is set to α s (M Z ) = 0.1185, and we impose a lower cutoff for the running scale µ such that it never goes below Q 0 = 2 GeV. The bands come from varying both the resummation scale µ 0 and the rapidity scale ζ 0 by a factor of 2 around their default values, which gives a much larger contribution than the variation of the scales µ t and µ H . The bands at NNLL get smaller than the ones at NLL, but there is no overlap between them. This is because we have exponentiated the rapidity scale dependence of the parameter C ζ through the resummed h R γ in eq. (4.13). This exponentiation makes the cross-section at a given resummation order contain some contributions of higher orders, and thus the NLL band, in particular, is smaller. At NNLL this issue is much less relevant, as will be clear below when comparing with the prediction with the resummation in impact parameter space.
If we compare the two panels in figure 4, we see that the impact of the non-perturbative contribution leads to a significant change of the distribution. However the choice β f,h = 1 is rather extreme, since the non-perturbative model in eq. (7.4), which is an exponential function, induces rather large corrections to the perturbative expression in the low b T region, exactly where one would expect it to work better. Thus, given that the non-perturbative parameters should probably be smaller, we conclude that the Higgs boson transverse momentum distribution is not very sensitive to those parameters. The same conclusion was drawn in [26], where a Gaussian model was used to parametrize the non-perturbative contributions, and which led to an almost negligible impact on the distribution. This is easy to understand, since the Gaussian function in the low b T region is closer to 1 than the exponential function and therefore has an even smaller impact.
Let us now turn our attention to figure 5, where we present a similar prediction to the one already discussed, but with the resummation performed in impact parameter space. We use here the same settings for the PDFs and the running of the strong coupling as in figure 4. The relevant resummed expressions for the unpolarized and linearly polarized gluon distributions were given in eq. (6.3). As can be seen, the bands at NLL are now bigger compared to the previous approach, and overlap with the NNLL bands. Again, comparing the two panels in this figure we see that the effect of the explored non-perturbative parameters is rather small. Notice that the NNLL curves within the two approaches, in figures 4 and 5, are compatible, and consistent with the recent results found in [27]. Finally, in figure 6 we show the predictions for the distribution at √ s = 13 TeV, at NNLL accuracy and again for extreme values of the non-perturbative parameters, with both resummation approaches. The cross-section is bigger than at √ s = 8 TeV, but the same conclusions regarding the sensitivity to the non-perturbative parameters apply: the range in parameter variation shown in the left figure is rather large, and it seems unlikely that experimental measurements of the Higgs q T distribution at the LHC will be precise enough to fix the non-perturbative parameters of gluon TMDPDFs, apart from excluding the most vivid parameter values.

Conclusions
Using longitudinally polarized hadron (g g 1L ). Having at our disposal the proper definition of gluon TMDPDFs is crucial in order to consistently analyse different processes where they appear.
From the structure of the factorization theorem derived, we conclude that the evolution of all leading-twist (un)polarized gluon TMDPDFs is universal, i.e., the same evolution kernel can be applied to evolve any of them. Moreover, given the currently known perturbative ingredients we have performed the resummation of large logarithms contained in this evolution kernel up to NNLL accuracy.

JHEP07(2015)158
TMDPDFs are functions that contain perturbatively calculable information when the transverse momentum is in the perturbative domain. In this work we have considered all gluon TMDPDFs and discussed their operator product expansion in terms of collinear functions. The OPE Wilson coefficients depend on the particular distribution but we have shown that part of them is the same for all TMDPDFs. We have furthermore resummed those universal pieces at NNLL accuracy, increasing our control over the perturbative ingredients of the TMDPDFs. Moreover we have derived, for the first time, the NLO Wilson coefficient for the gluon helicity TMDPDF g g 1L , which will allow more accurate phenomenological studies of this quantity in the future, e.g., at RHIC, AFTER@LHC or EIC. We have also derived the OPE Wilson coefficients for f g 1 and h ⊥g 1 in the framework presented in this paper.
Using the obtained results we have performed a numerical study of the contribution of linearly polarized gluons for the productions of η b and Higgs boson in unpolarized hadronhadron collisions. The major conclusion is that the larger the relevant hard scale is, the less sensitive is the observable to their non-perturbative contribution, and therefor harder to extract. Thus one would need to combine low-and high-energy experimental data and properly implement the QCD evolution of gluon TMDPDFs in order to extract it. On the other hand, the fact that at large scales the transverse momentum distributions are less sensitive to the non-perturbative parameters of the TMDPDFs allows us to obtain accurate predictions even if currently there is no information on these parameters.
Finally we have provided some predictions for the Higgs boson transverse momentum distribution at the LHC, both at √ s = 8 TeV and √ s = 13 TeV, using the formalism presented in this paper, i.e., expressing it in terms of well-defined gluon TMDPDFs. We have studied the impact of non-perturbative contributions on the distribution and have shown that the sensitivity to them is very small.
A OPE of f g 1 at NLO In this appendix we present the calculation of the unpolarized gluon TMDPDFs at O(α s ), using dimensional regularization with the MS-scheme (µ 2 → µ 2 e γ E /(4π)) for ultra-violet (UV) divergences and the ∆-regulator [9] for IR and rapidity divergences. We use the Keldysh formalism to perform the calculation (see, e.g., [73,74]). Our first goal is to show explicitly the cancellation of rapidity divergences in the properly defined gluon TMDPDFs in eq. (2.13). On the other hand, we will extract the Wilson matching coefficients of the TMDPDF onto its collinear counterparts, as they appear in eq. (6.1).
With the ∆-regulator, we write the poles of the gluon propagators that involve p orp with a real and positive parameters ∆ ± , and for collinear and soft Wilson lines one has Now, given the fact that the soft and collinear matrix elements must reproduce the soft and collinear limits of full QCD, they need to be regulated consistently, and thus δ ± are related with ∆ ± through the large components of the collinear fields, Note that ∆ ± (and hence δ ± ) are regulator parameters, and are set to zero unless they regulate any divergence.
Let us now proceed with the partonic calculation, using eq. (2.15). If we consider a hadron with definite helicity λ and take into account only the functions f f 1 , h ⊥g 1 and g g
. with the particular result when Λ → 0. We have also used the following relations: where f (x) is any function regular at x → 1. Thus, in IPS, the collinear matrix element for the partonic channel of a gluon splitting into a gluon is The mixed divergences in the result above ( 1 ε UV ln∆ + ) are rapidity divergences, which need to be eliminated by combining it with the soft function as in eq. (2.15) in order to get a well-defined TMDPDF. Now we turn our attention to the soft function. Diagram 8a and its Hermitian conjugate give

JHEP07(2015)158
We are ready now, given eq. (6.1), to extract the matching coefficient of the TMDPDF onto the PDF in the g/g channel: For the g/q channel the unpolarized collinear gluon PDF is given by and thus the matching of the TMDPDF onto the PDF in the g/q channel is

at NLO
The calculation in this appendix follows the same logic as in the previous one, so we limit ourselves to provide the relevant results.
For the distribution of linearly polarized gluons inside an unpolarized hadron at NLO only real diagrams contribute. In the g/g channel we have:

JHEP07(2015)158
The first three coefficients are T F n f , The first three coefficients of the anomalous dimension γ g are [83,84] γ g 0 = 0 ,

JHEP07(2015)158
Finally, the coefficients for the QCD β-function are T F n f , (E.12) 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.