  We address quarkonium formation at moderate to large transverse momenta,
where the single-parton collinear fragmentation prevails over the
short-distance emission, directly from the hard sub-scattering, of the
constituent heavy-quark pair. We rely on Non-Relativistic-QCD (NRQCD)
Next-to-Leading Order (NLO) calculations for all parton fragmentation channels
to quarkonia, taken as proxies for initial-scale inputs. Preliminary sets of
Variable-Flavor Number-Scheme (VFNS) fragmentation functions (FFs) are built
via a DGLAP scheme that properly accounts for evolution thresholds. Statistical
errors are assessed via a Monte Carlo (MC), replica-like approach aimed at
catching Missing Higher-Order Uncertainties (MHOUs).
