  We compute bottom mass ($m_b$) corrections to the transverse momentum ($q_T$)
spectrum of Higgs bosons produced by gluon fusion in the regime $q_T \sim m_b
\ll m_H$ at leading power in $m_b/m_H$ and $q_T/m_H$, where the gluons couple
to the Higgs via a top loop. To this end we calculate the quark mass dependence
of the transverse momentum dependent gluon beam functions (aka gluon TMDPDFs)
at two loops in the framework of SCET. These functions represent the collinear
matrix elements in the factorized gluon-fusion cross section for small $q_T$.
We discuss in detail technical subtleties regarding rapidity regulators and
zero-bin subtractions in the calculation of the virtual corrections present for
massive quarks. Combined with the known soft function for $m_b \neq 0$ our
results allow to determine the resummed Higgs $q_T$ distribution in the
top-induced gluon fusion channel at NNLL$^\prime$ (and eventually N$^3$LL) with
full dependence on $m_b/q_T$. We perform a first phenomenological analysis at
fixed order, where the new corrections to the massless approximation lead to
percent-level effects in the peak region of the Higgs $q_T$ spectrum. Upon
resummation they may thus be relevant for state-of-the-art precision
predictions for the LHC.
