  We discuss factorization in heavy quarkonium production in high energy
collisions using NRQCD. Infrared divergences at NNLO are not matched by
conventional NRQCD matrix elements. However, we show that gauge invariance and
factorization require that conventional NRQCD production matrix elements be
modified to include Wilson lines or non-abelian gauge links. With this
modification NRQCD factorization for heavy quarkonium production is restored at
NNLO.
