  We discuss heavy quarkonium production through parton fragmentation,
including a review of arguments for the factorization of high-p_T particles
into fragmentation functions for hadronic initial states. We investigate the
further factorization of fragmentation functions in the NRQCD formalism, and
argue that this requires a modification of NRQCD octet production matrix
elements to include nonabelian phases, which makes them gauge invariant. We
describe the calculation of uncanceled infrared divergences in fragmentation
functions that must be factorized at NNLO, and verify that they are absorbed
into the new, gauge invariant matrix elements.
