  We present a perturbative QCD factorization formalism for the production of
heavy quarkonia of large transverse momentum $p_T$ at collider energies, which
includes both the leading power (LP) and next-to-leading power (NLP)
contributions to the cross section in the $m_Q^2/p_T^2$ expansion for heavy
quark mass $m_Q$. We estimate fragmentation functions in the non-relativistic
QCD formalism, and reproduce the bulk of the large enhancement found in
explicit NLO calculations in the color singlet model. Heavy quarkonia produced
from NLP channels prefer longitudinal polarization.
