  We combine soft and collinear matrix elements to define transverse momentum
dependent parton distributions (TMDs) for gluons, free from rapidity
divergences. We establish a factorization theorem at next-to-leading order for
the Higgs boson transverse momentum ($q_T$) spectrum, and use it to derive
evolution equations for gluon TMDs. The evolution for all gluon TMDs is driven
by a universal kernel, i.e. the same for all polarizations. In the region of
intermediate $q_T$ we match the unpolarized, helicity and linearly polarized
gluon distributions onto PDFs. We calculate the resummed Higgs transverse
momentum distribution at NNLL, including the contribution of the linearly
polarized gluons, and investigate the impact of non-perturbative models for the
transverse distribution of gluons inside the proton.
