We address the unphysical energy dependence of quarkonium-hadroproduction cross sections at Next-to-Leading Order (NLO) in $$\alpha _s$$ ${\alpha }_s$ which we attribute to an over-subtraction in the factorisation of the collinear singularities inside the PDFs in the $$\overline{\text {MS}}$$ $\overline{\text{MS}}$ scheme. Such over- or under-subtractions have a limited phenomenological relevance in most of the scattering processes in particle physics. On the contrary, it is particularly harmful for $$P_T$$ $P_T$ -integrated charmonium hadroproduction which renders a wide class of NLO results essentially unusable. Indeed, in such processes, $$\alpha _s$$ ${\alpha }_s$ is not so small, the PDFs are not evolved much and can be rather flat for the corresponding momentum fractions and, finally, some process-dependent NLO pieces are either too small or too large. We propose a scale-fixing criterion which avoids such an over-subtraction. We demonstrate its efficiency for $$\eta _{c,b}$$ ${\eta }_{c,b}$ but also for a fictitious light elementary scalar boson. Having provided stable NLO predictions for $$\eta _{c,b}$$ ${\eta }_{c,b}$ $$P_T$$ $P_T$ -integrated cross sections, $$\sigma ^{\mathrm{NLO}}_{\eta _Q}$$ ${\sigma }_{{\eta }_Q}^{\mathrm{NLO}}$ , and discussed the options to study $$\eta _{b}$$ ${\eta }_b$ hadroproduction, we argue that their measurement at the LHC can help better determine the gluon PDF at low scales and tell whether the local minimum in conventional NLO gluon PDFs around $$x=0.001$$ $x=0.001$ at scales below 2 GeV is physical or not.