  We present a factorization formalism for inclusive production of heavy
quarkonia of large transverse momentum, $p_T$ at collider energies, including
both leading power (LP) and next-to-leading power (NLP) behavior in $p_T$. We
demonstrate that both LP and NLP contributions can be factorized in terms of
perturbatively calculable short-distance partonic coefficient functions and
universal non-perturbative fragmentation functions, and derive the evolution
equations that are implied by the factorization. We identify projection
operators for all channels of the factorized LP and NLP infrared safe
short-distance partonic hard parts, and corresponding operator definitions of
fragmentation functions. For the NLP, we focus on the contributions involving
the production of a heavy quark pair, a necessary condition for producing a
heavy quarkonium. We evaluate the first non-trivial order of evolution kernels
for all relevant fragmentation functions, and discuss the role of NLP
contributions.
