  Quarkonium production in hadron collisions at low transverse momentum
$q_\perp \ll M$ with $M$ as the quarkonium mass can be used for probing
transverse momentum dependent (TMD) gluon distributions. For this purpose, one
needs to establish the TMD factorization for the process. We examine the
factorization at the one-loop level for the production of $\eta_c$ or $\eta_b$.
The perturbative coefficient in the factorization is determined at one-loop
accuracy. Comparing the factorization derived at tree level and that beyond the
tree level, a soft factor is, in general, needed to completely cancel soft
divergences. We have also discussed possible complications of TMD factorization
of p-wave quarkonium production.
