Abstract
Cosmological neutrino mass bounds are becoming increasingly stringent. The latest limit within $\Lambda$CDM from Planck 2018+ACT lensing+DESI is $\sum m_\nu < 0.072\,{\rm eV}$ at 95% CL, very close to the minimum possible sum of neutrino masses ($\sum m_\nu > 0.06\,{\rm eV}$), hinting at vanishing or even ''negative'' cosmological neutrino masses. In this context, it is urgent to carefully evaluate the origin of these cosmological constraints. In this paper, we investigate the robustness of these results in three ways: i) we check the role of potential anomalies in Planck CMB and DESI BAO data; ii) we compare the results for frequentist and Bayesian techniques, as very close to physical boundaries subtleties in the derivation and interpretation of constraints can arise; iii) we investigate how deviations from $\Lambda$CDM, potentially alleviating these anomalies, can alter the constraints. From a profile likelihood analysis, we derive constraints in agreement at the $\sim 10\%$ level with Bayesian posteriors. We find that the weak preference for negative neutrino masses is mostly present for Planck 18 data, affected by the well-known `lensing anomaly'. It disappears when the new Planck 2020 HiLLiPoP is used, leading to significantly weaker constraints. Additionally, the pull towards negative masses in DESI data stems from the $z=0.7$ bin, which is in $\sim 3\sigma$ tension with Planck expectations. Without these outliers, and in combination with HiLLiPoP, the bound relaxes to $\sum m_\nu < 0.11\,{\rm eV}$ at 95% CL. The recent preference for dynamical dark energy alleviates this tension and further weakens the bound. As we are at the dawn of a neutrino mass discovery from cosmology, it will be very exciting to see if this trend is confirmed by future data.