We report progress on the determination and study of quarkonium production within the fragmentation approximation.
Our analyses address the moderate and large transverse-momentum regime, where the collinear fragmentation of a single parton is expected to dominate over the short-distance production, directly from the hard scattering, of the constituent $(Q \bar{Q})$ system. Parton fragmentation channels to pseudoscalar and vector quarkonia are built on the basis of non-Relativistic QCD next-to-leading computations, which we use to model initial-scale fragmentation inputs. Thus, a preliminary family of Variable-Flavor Number-Scheme (VFNS) fragmentation functions, named NRFF1.0, are constructed through standard DGLAP evolution.
Statistical uncertainties are obtained from a Monte Carlo, replica-like approach embodying missing higher-order uncertainties.

