  I review the current status of quarkonium production theory based on the
non-relativistic QCD factorization formalism (NRQCD). While this theory
describes much of the world's data on J/\psi and \Upsilon production, there are
still outstanding problems, most notably the polarization of quarkonia at large
p_T in hadron colliders. In this talk we will present new tests of NRQCD
involving the distribution of quarkonia within jets. The distribution of
hadrons within jets is determined by nonperturbative functions called
fragmenting jet functions (FJFs). FJFs are convolutions of fragmentation
functions, evolved to the jet energy scale, with perturbatively calculable
matching coefficients. I show how the FJFs for quarkonia can be calculated in
NRQCD in terms of a few NRQCD long-distance matrix elements (LDME), so the
dependence of the cross section on the energy fraction of the heavy quarkonium,
z, is sensitive to the underlying production mechanism, and therefore provides
a new test of NRQCD. The jet energy and z dependence of the cross section can
be used to perform an independent extraction of NRQCD LDME. Finally, I describe
ongoing work building on this result. This includes comparison of analytic
resummed calculations with Monte Carlo simulations for two-jet events in e^+e^-
collisions with B mesons, and three-jet events with J/\psi, as well as the
definition of boost invariant jet substructure variables and calculation of a
boost invariant soft function that are necessary for analytic calculations of
jet cross sections at the Large Hadron Collider.
