  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.
