  Within the framework of transverse momentum dependent factorization in
combination with nonrelativistic QCD, we study charmonium and bottomonium
production in hadronic collisions. We focus on quarkonium states with even
charge conjugation, for which the color-singlet production mechanism is
expected to be also dominant in the small transverse momentum region, $q_T^2
\ll 4 M_{c,b}^2$. It is shown that the distributions of linearly polarized
gluons inside unpolarized, longitudinally and transversely polarized protons
contribute to the cross sections for scalar and pseudoscalar quarkonia in a
very distinctive, parity-dependent way, whereas their effects on higher angular
momentum states are strongly suppressed. We derive analytical expressions for
single and double spin asymmetries, which would allow for the direct extraction
of the gluon transverse momentum dependent distributions, mirroring the
phenomenological studies of the Drell-Yan processes aimed at the extraction of
their quark counterparts. By adopting Gaussian models for the gluon TMDs, which
fulfill without saturating everywhere their positivity bounds, we provide
numerical predictions for the transverse single-spin asymmetries. These
observables could be measured at LHCSpin, the fixed target experiment planned
at the LHC.
