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/ψ and ϒ 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/ψ, 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.