  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\to \eta_{c,b}$ at low transverse momentum.
