  We present robust, model-marginalized limits on both the total neutrino mass ($\sum m_\nu$) and abundance ($N_{\rm eff}$) to minimize the role of parameterizations, priors and models when extracting neutrino properties from cosmology. The cosmological observations we consider are CMB temperature fluctuation and polarization measurements, Supernovae Ia luminosity distances, BAO observations and determinations of the growth rate parameter from the Data Release 16 of the Sloan Digital Sky Survey IV. The degenerate neutrino mass spectrum (which implies $\sum m_\nu>0$) is weakly (moderately) preferred over the normal and inverted hierarchy possibilities, which imply the priors $\sum m_\nu>0.06$ and $\sum m_\nu>0.1$ eV respectively. Concerning the underlying cosmological model, the $\Lambda$CDM minimal scenario is almost always strongly preferred over the possible extensions explored here. The most constraining $95\%$ CL bound on the total neutrino mass in the $\Lambda$CDM+$\sum m_\nu$ picture is $\sum m_\nu< 0.087$ eV. The parameter $N_{\rm eff}$ is restricted to $3.08\pm 0.17$ ($68\%$ CL) in the $\Lambda$CDM+$N_{\rm eff}$ model. These limits barely change when considering the $\Lambda$CDM+$\sum m_\nu$+$N_{\rm eff}$ scenario. Given the robustness and the strong constraining power of the cosmological measurements employed here, the model-marginalized posteriors obtained considering a large spectra of non-minimal cosmologies are very close to the previous bounds, obtained within the $\Lambda$CDM framework in the degenerate neutrino mass spectrum. Future cosmological measurements may improve the current Bayesian evidence favouring the degenerate neutrino mass spectra, challenging therefore the consistency between cosmological neutrino mass bounds and oscillation neutrino measurements, and potentially suggesting a more complicated cosmological model and/or neutrino sector. 
Model marginalized constraints on neutrino properties from cosmology