Quarkonia production in different high-energy processes has recently been proposed in order to probe gluon transverse-momentum-dependent parton distribution and fragmentation functions (TMDs in general). However, no proper factorization theorems have been derived for the discussed processes, but rather just ansatzs, whose main assumption is the factorization of the two soft mechanisms present in the process: soft-gluon radiation and the formation of the bound state. In this paper it is pointed out that, at low transverse momentum, these mechanisms are entangled and thus encoded in a new kind of non-perturbative hadronic quantities beyond the TMDs: the TMD shape functions. This is illustrated by deriving the factorization theorem for the process pp → η c,b at low transverse momentum. * mgechevarria@pv.infn.it