Cosmological neutrino mass bounds are becoming increasingly stringent. The latest limit within LambdaCDM from Planck 2018ACT lensingDESI is sum mnu lt 0.072rm eV at 95 CL very close to the minimum possible sum of neutrino masses sum mnu gt 0.06rm 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 LambdaCDM 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 wellknown 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 z0.7 bin which contains a BAO measurement in sim 3sigma tension with Planck expectations. Without this bin and in combination with HiLLiPoP the bound relaxes to sum mnu lt 0.11rm 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.
['Astrophysics - Cosmology and Nongalactic Astrophysics', 'High Energy Physics - Experiment', 'High Energy Physics - Phenomenology']
Living at the Edge A Critical Look at the Cosmological Neutrino Mass Bound
{'Living at the Edge: A Critical Look at the Cosmological Neutrino Mass Bound': "Daniel Naredo-Tuero , 1, ∗ Miguel Escudero , 2, †  \nEnrique Fernandez-Martinez , 1, ‡ Xabier Marcano , 1, § and Vivian Poulin 3, ¶  \n1  \nDepartamento de F'ısica Te'orica and Instituto de F'ısica Te'orica UAM/CSIC, Universidad Aut'onoma de Madrid, Cantoblanco, 28049 Madrid, Spain  \n2 Theoretical Physics Department, CERN, 1211 Geneva 23, Switzerland  \n3 Laboratoire univers et particules de Montpellier (LUPM),  \nCentre national de la recherche scientifique (CNRS) et Universit'e de Montpellier, Place Eug'ene Bataillon, 34095 Montpellier C'edex 05, France  \n(Dated: October 29, 2024)  \nCosmological neutrino mass bounds are becoming increasingly stringent. The latest limit within ΛCDM from Planck 2018+ACT lensing+DESI is ∑ m ν < 0 . 072 eV at 95% CL, very close to the minimum possible sum of neutrino masses ( ∑ m ν > 0 . 06 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 ΛCDM, potentially alleviating these anomalies, can alter the constraints. From a profile likelihood analysis, we derive constraints in agreement at the ∼ 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 contains a BAO measurement in ∼ 3 σ tension with Planck expectations. Without this bin, and in combination with HiLLiPoP, the bound relaxes to ∑ m ν < 0 . 11 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.", 'I. INTRODUCTION: NEUTRINO MASS BOUNDS AS OF MID 2024': "In April 2024, the DESI collaboration presented the most stringent bound on the sum of neutrino masses within the standard cosmological model [1]:  \n∑ m ν < 0 . 072 eV [95% CL] , (1)  \nobtained by combining their new DESI-Y1 baryon acoustic oscillation data (BAO) [2, 3] with Planck [4-6] and ACT data [7, 8] (see also [9, 10] for further updated constraints within ΛCDM 1 ).  \nThe bound in Eq. (1) should be compared with the minimum possible value of the sum of neutrino masses  \ngiven the observed mass squared differences from neutrino oscillation experiments. Taking them from NuFITv5.3 [12] and by taking m ν, lightest → 0, one finds at 5 σ CL:  \n∑ m ν > 0 . 057 eV ≃ 0 . 06 eV [NO] , (2)  \ndepending upon the neutrino mass ordering. At present, there is no strong preference for either ordering from global analyses of neutrino oscillation data, see [12-14].  \n∑ m ν > 0 . 096 eV ≃ 0 . 10 eV [IO] , (3)  \nClearly, the current cosmological limit is very close to the minimal value in normal ordering and already disfavours to some degree the inverted one. Importantly, the limit in Eq. (1) is so close to the minimum physical boundary that statistical statements about the neutrino mass need to be taken with care. However, cosmological limits are typically derived within a Bayesian framework, and as such will depend upon the priors used. While the dependence may be weak when the likelihood largely dominates over the prior, and far away from physical boundaries, it has already been established that neutrino masses are strongly sensitive to the choice of prior [1518], even with DESI [1]. Moreover, and perhaps even more surprising, there is still no hint of a non-zero mass in the posterior probability density, and in fact, it has  \nbeen argued that cosmological data may favor a negative effective neutrino mass [19-21]. While this could be the result of a statistical fluctuation or a systematic effect, it could potentially be the indication of new phenomena in cosmology with groundbreaking implications.  \nGiven the current situation, it is urgent to carefully analyze the origin and behaviour of the present constraints on neutrino masses from cosmology. In this paper, we set as goals to investigate simultaneously: i) the data that have been used and the role of potential statistical anomalies in those data, ii) the statistical methods used to derive the results, and whether unwarranted effects in the Bayesian analysis drive the preference for negative neutrino masses, and iii) the extent upon which they rely on the assumption of the standard ΛCDM model.  \nIn order to tackle these questions, we perform an extensive comparison of cosmological neutrino mass bounds derived both from a Bayesian and a frequentist standpoint. In particular, the use of a frequentist analysis framework can lead to new insight into the sensitivity of the bound on the statistical procedure adopted, and can help in addressing the role of priors in the Bayesian limits. In addition, as we build the likelihood profile of the sum of neutrino masses in light of various datasets, it is possible to extrapolate to the unphysical region and study the potential preference for negative masses, as performed in previous analyses of the neutrino mass in cosmology [22-24]. Our analysis is complementary to recent Bayesian method that rely on the use of an 'effective neutrino mass' to model the effect of a negative neutrino mass [19, 21], and although it comes with its own set of approximations, it by-passes the need for an arbitrary definition which will necessarily miss part of the physical effect that a real model would have.  \nSecond, to address the questions of the robustness of the bound to the choice of data, we perform series of analyses comparing results from the latest BAO data from DESI [2, 3], previous ones from SDSS [24], the compilation of uncalibrated supernova (SN) distances from Pantheon+ [25], and most importantly the latest versions of the Planck likelihoods based on the 2020 PR4 data release [26]. This is key because Planck currently dominates the neutrino mass bounds and because these latest likelihoods not only contain roughly O (10%) more statistical power than the 2018 (PR3) release, but importantly because they have a much better handle of a number of systematic effects (see [26] for details). Importantly, the known lensing anomaly which was present at the 2 . 8 σ level in the Planck 2018 Plik likelihood [4] is only present in the 2020 implementations at the 1 . 7 σ or 0 . 75 σ levels, for the CamSpec [27, 28] and HiLLiPoP [29, 30] likelihoods, respectively. This is well known to have an impact on the inference of neutrino masses [31-33]. Yet, Eq. (1) does not make use of these newer data, and should thus be explicitly checked. In particular, we test the impact of the anomaly on the bounds in two ways: first, by simply updating the data to match the newer release (as  \nin [10, 20]) and second, by explicitly marginalizing over the A lens parameter [4, 34], that scales the amplitude of the lensing power spectrum.  \nFinally, we investigate why the bound in Eq. (1) is so strong, despite the statistical power of DESI-Y1 BAO data being a priori not larger (though competitive) than that of the entire SDSS sample. We believe that this is crucial as the community expects upcoming measurements from DESI and Euclid to actually pin-down the absolute neutrino mass scale in cosmology within the next few years, see e.g. [35-38]. If these collaborations do not report any measurement of the neutrino mass, their analyses will clearly signal a breakdown of the standard cosmological model and possibly even new physics in the neutrino sector. This could point towards decaying neutrinos [39-45], non-standard cosmic neutrino backgrounds [46-51], or neutrinos with a time-varying mass [52-59]. Yet, before making such claims, it is important to note that DESI data are in some (arguably small) 2 σ level tension with Planck. In fact, there are two BAO measurements at z = 0 . 5 and z = 0 . 7 that are driving the tension with Planck (respectively at the 2.8 σ and 2 . 6 σ level), and that seem somewhat at odds with SDSS results as well. This raises the questions of the role of those data points in driving the strong bound in Eq. (1). In addition, when combined with Planck and Pantheon+ (or other SN compilation) the compilation of data seem to favor dynamical dark energy over ΛCDM [1]. Such preference is known to also alter the bound given by Eq. (1). Within our joint Bayesian and frequentist framework, we will thus explicitly check the role of those potential outlier points, as well as the preference for a deviation from ΛCDM in driving the hint of negative neutrino mass, and the strong neutrino mass bounds.  \nThe rest of our study is structured as follows: first, in Section II we briefly review the main cosmological implications of massive neutrinos. We will highlight the crucial role of Planck data as well as why BAO data can significantly tighten the constraints on the neutrino mass. In Section III we outline the various data sets that we will use, and describe the statistical procedure (both frequentist and Bayesian) that we use to analyse the data and derive bounds. In Section IV we present our main results focusing first on analyses with Planck data only, exploring in detail the role of the lensing anomaly; we then perform comprehensive analysis of Planck and BAO data, comparing in particular DESI and SDSS; and finally, check the impact of considering extensions to ΛCDM. In Section V we specifically address the potential cosmological preference for a negative neutrino mass. Finally, in Section VI we draw our conclusions. Additional material supporting our findings, including posterior distributions, correlations and comparisons with previous works, is provided in the Appendices.", 'II. COSMOLOGICAL IMPACT OF THE NEUTRINO MASS': "The existence of the Cosmic Neutrino Background (C ν B) is a key prediction of the Standard Model of cosmology [60]. The C ν B would have formed in the early Universe at temperatures of ∼ 2 MeV [61], and according to the ΛCDM model, we should be living in a Universe filled with a number density of neutrinos of n ν ≃ 56 / cm 3 per helicity state. These neutrinos gravitate today primarily as a result of their mass and in fact their energy density today is not negligible:  \nΩ ν h 2 = ∑ m ν n ν ρ crit ≃ ∑ m ν 93 . 2 eV ≃ 0 . 0012 ∑ m ν 0 . 12 eV . (4)  \nThis is O (1)% of the dark matter energy density.  \nThe implications of the neutrino mass in cosmology have been explored and discussed in depth and length, see [60-65] for reviews. In essence, the cosmological implications of neutrino masses can be understood as follows: 1) Neutrinos are always a relevant component of the energy density of the Universe and therefore contribute to its expansion rate, H ∝ √ ρ . 2) Neutrinos were propagating ultrarelativistically until the Universe cooled down to T ν ≃ m ν / 3 which occurs at z NR ν ≃ 190 m ν / (0 . 1 eV). 3) As a result of their ultrarelativistic speeds, neutrinos are not able to cluster on scales smaller than L ≃ 20 Mpc0 . 1 eV /m ν today. In consequence, and for massive neutrinos which become non-relativistic after recombination m ν ≲ 0 . 6 eV, one of their main cosmological implications is to suppress the amount of structure formation on scales smaller than L ≃ 20 Mpc0 . 1 eV /m ν . This suppression is of course strongly dependent upon the energy density neutrinos represent in the Universe which is directly proportional to the neutrino masses as given in Eq. (4). This discussion clearly highlights that the best way to search for the neutrino mass is arguably by directly observing probes of the Large Scale Structure in the Universe. This is one of the primary goals of the ongoing DESI and Euclid surveys (see [35-37]). Nevertheless, as of today the bound on the neutrino mass is dominated by Planck with the aid of BAO and Supernova measurements to break relevant parameter degeneracies.  \nHow is then Planck sensitive to the neutrino mass if neutrinos with m ν < 0 . 6 eV became non-relativisitic after the CMB was formed? There are two effects that matter [66]: Firstly, as a result of neutrinos traveling at relativistic speeds, they lead to a suppression of the matter power spectrum. This reduces the lensing that CMB photons experience on their way from the last scattering surface until today, and this in turn leads to sharper peaks in the CMB power spectrum. This effect is depicted in Fig. 1. In this figure, several other cosmological parameters have been fixed: the total matter densities of baryons ( ω b ) and dark matter ( ω cdm ), the Thomson optical width to reionization τ , the shape and amplitude of the primordial matter power spectrum ( n s and A s ), as  \nFIG. 1. Impact of a non-zero neutrino mass on the TT power spectrum. Inspired by Figure 26.2 of [67] but showing the Planck error bars taken from the binned PR3 data release a . Note that it is precisely in the range of angular scales where the main impact of neutrino masses appears where the lensing anomaly is present in some Planck likelihood implementations.  \n<!-- image -->  \n- a https://irsa.ipac.caltech.edu/data/Planck/release\\_3/ ancillary-data/  \nwell as the angular scale of the first CMB peak, θ s . From Fig. 1 we can clearly see that the main constraining power for neutrino masses from Planck will come from rather small angular scales. Importantly, this lensing effect is also dependent upon what the values of A s and ω m are.  \nTherefore, in order to make precise inferences on the neutrino mass from CMB observations one needs to have i) a very good understanding of how the CMB is being lensed on small angular scales, ii) control over A s and ω m , and iii) control over Ω m and more generally the expansion history at late-times.  \nSecondly, ∑ m ν impacts the angular diameter distance to the last-scattering surface as neutrinos contribute to the Universe's expansion. It is however possible to exploit the geometrical degeneracy in the CMB by adjusting H 0 to compensate the effect of ∑ m ν on the last scattering surface and ensure the angular diameter distance to the last-scattering surface is left unaffected. This is done in Fig. 1, by keeping θ s fixed, and requires a smaller H 0 . However, in this process Ω m ≡ ω m /h 2 will increase as a result of the necessary decrease in the reduced Hubble parameter h ≡ H 0 / 100 km/s/Mpc. This suggests that an accurate probe of the late-time expansion history can help break the geometrical degeneracy (namely Ω m and/or H 0 ) and further constrain ∑ m ν .  \nIn this context, regarding i), the Planck legacy analysis did report the so-called 'lensing anomaly' which can impact the inferences on the neutrino masses. The anomaly was parametrized introducing an ad hoc variable A lens to change the lensing power relative to its actual physi-  \ncal value, so as to account for possible systematic uncertainties, and Planck data preferred values of A lens > 1. However, the magnitude of the A lens was known to vary depending on the specific likelihood (and specific CMB dataset) used to perform the analysis. While the anomaly is 2 . 8 σ with the official Planck collaboration 2018 TTTEEE likelihood, dubbed plik , it is reduced after the final PR4 data release to the 1.7 σ level with the alternative CamSpec [27, 28] likelihood, and down to 0.75 σ with the HiLLiPoP [29, 30] one.  \nImportantly, the CamSpec and HiLLiPoP likelihoods have been recently updated since the 2018 analysis, in light of a new set of maps produced by the Planck collaboration called NPIPE. The NPIPE maps exploit a number of improvements in the processing of time ordered data to allow for an increase in the signal-to-noise ratio at small scales. They also allow to use a larger sky fraction, and incorporate a better handling of a number of systematic errors thanks to dedicated mock data. This results in a roughly ∼ 10% stronger constraining power on ΛCDM parameters, and importantly, the lensing anomaly seems to be significantly reduced or even absent in those data. In fact, it has been shown that they lead to weaker bounds to ∑ m ν than the Plik 2018 likelihood, see [10, 28, 30].  \nRegarding ii), the best way to obtain better measurements on A s is from large scale CMB polarization measurements which will be provided by LiteBIRD [68] but on a ∼ 10 year time-scale. Importantly, regarding iii), the improvement is happening now, as DESI is taking data and has published already the 1st year data release, while Euclid is on space and will start collecting data in one year as well. In this regard, it is important to highlight that the compatibility between the new DESIY1 data and Planck is at the 2 σ level, and is thus worse than for SDSS. As a result, this (arguably small) tension may impact the neutrino mass bound. This is illustrated in Fig. 2, where we compare constraints under ΛCDM in the Ω m -H 0 r d plane from SDSS, DESI and DESI without the data points at z = 0 . 7. We also superimpose the posteriors from a fit to Planck 2018 (with lensing), highlighting the correlations with ∑ m ν with colored points. As argued in Ref. [1], it is clear that the mismatch in H 0 r d between SDSS and DESI is driving the discrepancy, and a stronger bound to ∑ m ν . However, removing the data points at z = 0 . 7 can significantly shift the mean of the posterior distributions (by about ∼ 0 . 8 σ ) without affecting the error bars. This suggests than an analysis without these data points may lead to significantly weaker bounds, and would argue in favor of a potential statistical fluke driving these bounds.  \nIn summary, it is very important to understand how relevant are possible systematic effects in Planck CMB data on our inferences of the neutrino mass in cosmology, as well as what are the implications of adding different sets of BAO data. This is particularly the case given that the direct combination of Planck + DESI-Y1  \nFIG. 2. Implications of BAO measurements of Ω m and H 0 r d for ∑ m ν inferences. We show the posterior density contours using Planck data (grey and dots), as well as the regions favoured by the full SDSS BAO sample (in blue), DESI-Y1 (in red), and DESI-Y1 without the z = 0 . 7 bin, which contains a 2 . 6 σ outlier (in black dashed). We can clearly see that the z = 0 . 7 leads to a relevant shift on the parameter space with implications for the neutrino mass.  \n<!-- image -->  \nBAO yields a bound ∑ m ν < 0 . 072 eV which is very close to the minimum allowed value from neutrino oscillation experiments ∑ m ν > 0 . 06 eV.", 'III.1. Cosmological Data: CMB, BAO and Supernova': 'In what follows, we will first perform a comprehensive analysis of the bound on neutrino masses coming from considering Planck data alone, to highlight the role of potential anomalies (whether a statistical fluke or a systematic effect) in the data, and how subsequent data releases have affected those bounds. To that end, we consider the following likelihood combinations:  \n- · Planck18-PR3 - We consider the default plik Planck legacy likelihoods for both TT, TE, EE high ℓ spectra as well as the large scale (low ℓ ) EE polarization likelihood SimAll , and also the large scale TT Commander likelihood.\n- · CamSpec22-PR4 -We consider the new Planck CamSpec NPIPE12 7 TTTEEE likelihood for both TT, TE, EE high ℓ spectra [27, 28] as well  \nas the large scale (low ℓ ) EE polarization likelihood SimAll and the large scale TT Commander one.  \n- · HiLLiPoP23-PR4 -We consider the new HiLLiPoP [30] likelihood for both TT, TE, EE high ℓ spectra, and the LoLLiPoP [29, 69] EE, EB and BB low ℓ spectra.\n- · Lensing-PR3 - Unless otherwise stated (no lensing) we also consider the reconstructed gravitational lensing potential power spectra from Planck PR3 [5] 2 . We note, however, that there is a newer PR4 one [6], as well as a complementary one by ACT [7, 8]. Including these data sets will not alter our conclusions and will only improve the bound by ∼ 10% (compare Eq. (1) with Eq. (5a)).  \nSecondly, we will investigate the role of BAO data, and the impact of potential statistical fluctuations within DESI data on the bounds. We include the following data combinations:  \n- · DESI-Y1 - We consider the full set of BAO data as reported in [1, 2].\n- · DESI-Y1-no07 - We consider the full set of BAO data as reported in [1, 2] but without including the two data points at z = 0 . 7. Since this corresponds precisely to the redshift of dark energy domination and this bin contains a data point in ∼ 3 σ tension 3 with Planck predictions, it is interesting to explore its impact on the ∑ m ν bound to understand what would happen if these outliers are not confirmed by future data, as highlighted in Figure 2.\n- · SDSS-full - We consider the full set of BAO data from SDSS as detailed in [24].\n- · DESI-Y1+SDSS - We consider the DESI and SDSS combination as done in [1] that makes use of the BAO measurements from the survey with the largest effective volume in a given redshift bin. Note that this includes the DESI-Y1 BAO measurements at z = 0 . 7 which are in ∼ 3 σ tension with Planck ΛCDM predictions.  \nFinally, we also consider analyses including uncalibrated luminosity distance-redshift measurements from type Ia SN:  \n- · SN-Pantheon -We make use of the Pantheon+ catalog of uncalibrated luminosity distance of type Ia supernovae (SNeIa) in the range 0 . 01 < z < 2 . 3 [25]. We note that there are alternative compilations such as the Union3 [71] and DES-Y5 SNeIa [72] that could be used. We do not anticipate them to strongly impact our conclusions. We leave a dedicated analysis to future work.', 'III.2. Analysis methodology: comparing Bayesian and frequentist framework': "In the cosmology community, it has become standard to perform analyses through a Bayesian framework, as those are typically numerically less expensive than frequentists analyses given the very large number of (cosmological and nuisance) parameters that must be considered. These also have the claimed advantage to incorporate prior knowledge (or lack thereof) in a straightforward manner, as priors are rooted in the definition of the posterior distribution within Bayes theoreom. In addition, there are now advanced tools (based on the notion of Bayesian evidence) to perform model comparison that are routinely used in Cosmology, and can (somewhat) easily tackle the problem known as 'look elsewhere effect' in the frequentist framework, and quantify the abstract notion of 'Occam's razor' that is often put forward. Yet, given that cosmology bounds on the sum of neutrino masses are pushing against the physical boundary, they can be affected by prior effects: the credible intervals built from the Bayesian posteriors become largely influenced by the choice of priors, rather than the data likelihood. As the choice of priors carry a level of arbitrariness, these effects can lead to constraints that are not robust. This is particularly relevant given the current context, and the apparent strength of the cosmological bound.  \nIndeed, over the past years, almost the entirety of each subsequent data release increased progressively the existing bound, with a preferred value at ∑ m ν = 0 and barely any hint for a non-vanishing value, which is surprising and worth investigating. In fact, present constraints have almost exhausted the available parameter space given the lower bound from neutrino oscillations and are already in significant tension with the minimal value implied by an inverted ordering, still allowed by oscillation results. In this context, subtleties when deriving the limits must be taken into consideration. In a Bayesian framework, the results are therefore priordependent (see e.g. Ref. [1]), as the preferred value appears to always be the smallest possible one given the prior (ie., ∑ m ν = 0, 0.06 or 0.1 eV when disregarding the results from oscillations or assuming a normal or inverted ordering, respectively) and thus (artificially) relax the larger the minimal allowed mass is.  \nAs an alternative to the Bayesian approach, and to  \ncompare and contrast the results, we make use of the frequentist framework, relying on a series of χ 2 optimizations at fixed neutrino mass, to build a profile likelihood and derive confidence intervals. The main advantages of the frequentist approach are that i) constraints are insensitive to the specific choice of priors, and ii) the presence of the physical boundary can be accounted for in a statistically consistent manner [73]. Thus, a direct comparison between the two approaches may provide a better handle on the relevance of these subtleties, as well as on the robustness of the cosmological bound.  \nAll power spectra in our study are obtained from the Boltzmann code CLASS [74]. The main physical quantity impacting cosmological observations is the total energy density in non-relativistic neutrinos and as such, for simplicity, we model the neutrinos as fully degenerate with a mass m ν = ∑ m ν / 3, and unless specified, vary it within the prior range m ν ∈ [0 , 1] eV. Cosmological data cannot be sensitive to the mass splittings and this choice can only cause relatively small changes in ∆ χ 2 , see e.g. [23]. We additionally vary the following six cosmological parameters, within large flat priors (when applicable): the angular size of the sound horizon θ s , the physical baryon ω b and dark matter ω cdm densities, the amplitude A s and tilt n s of the primordial power spectrum for scalar modes, the optical depth to reionization τ reio .  \nOn the Bayesian side, we perform a Markov Chain Monte Carlo (MCMC) sampling of the posterior distribution, using the publicly available MontePython code [75, 76]. Cosmological and nuisance parameters are varied according to the 'fast' and 'slow' parameters decomposition [77]. The chains are then marginalized with GetDist [78] in order to extract the bounds. On the frequentist side, instead of posterior marginalization, the relevant procedure is likelihood profiling, which requires the minimization of the χ 2 function for a fixed value of the parameter of interest ( ∑ m ν in our case), varying simultaneously the N -1 remaining cosmological and nuisance parameters. In particular, the various likelihood combinations described in Section III contain of the order of N ∼ 20 -25 parameters making the minimization a highly non-trivial task. In this context, we have opted to perform the numerical minimization with the simulated annealing algorithm implemented in Procoli [79] 4 . Although this procedure is straightforward, the HiLLiPoP likelihood requires significantly more care to successfully converge. We attribute this to the fact that HiLLiPoP is the sole CMB likelihood that is not binned, and is therefore noisier than the others. For this reason, our  \nprofiling is not performed over a regular grid of ∑ m ν . This is, however, not a major cause for concerns, since to properly calibrate the frequentist test statistics we will in any case fit the likelihood profiles. In addition, we can extrapolate the resulting fits to the unphysical negative neutrino mass region, allowing us to asses the potential preference for negative masses recently displayed by some datasets and discussed in [19-21].  \nConsequently, in order to gauge the impact of approaching the physical boundary at ∑ m ν = 0, we will choose two different prescriptions to extract the frequentist bound:  \n- 1. Bounded Likelihood (B.L.): focusing only on the physical region ∑ m ν > 0, we derive the bound via the standard ∆ χ 2 cut assuming Wilks' theorem, see e.g. [83].\n- 2. Feldman-Cousins (F.C.): we fit the χ 2 profile in the physical region to a parabola and extrapolate into the unphysical region so as to find where the true minimum would lie. It is then possible to extract the corrected upper bound from Table X of Ref. [73]. Compared with the previous procedure, this prescription has the advantage of guaranteeing proper coverage of the interval. However, it relies on the extrapolation in order to find the position and depth of the minimum.  \nIn order to assess whether a χ 2 profile is Gaussian and the Feldman-Cousins method can be safely applied, we will perform a parabolic fit of all points that lie below a certain ∆ χ 2 max . By varying this parameter in the interval ∆ χ 2 max ∈ [2 , 4], we can have a measure of how much the best-fit parabola or, equivalently, the F.C. bound, depends on the amount of points that are being fit. Following this method, we will only quote frequentist bounds on the profiles that exhibit stable best-fit parabolas. We find that this is the case for all of the analyses and that the limits derived from the Feldman-Cousins procedure do not appreciatively depend upon ∆ χ 2 max .", 'IV.1. Planck Only Analyses': "Planck Legacy vs Planck 2020 likelihoods - In Fig. 3 we compare the profile likelihoods for the neutrino mass using the latest implementations of the Planck likelihoods. We also show the first Planck 2013 results from [22] in red that already highlighted the potential preference for 'negative' neutrino masses. We show results from the 2018 legacy PR3 (black) as well as the latest implementations of the PR4 data release using CamSpec (green) and HiLLiPoP (blue). As discussed in the introduction, these  \n6  \nFIG. 3. Profile likelihoods for the neutrino masses within ΛCDM for three different versions of the Planck likelihoods: Planck18-PR3 , CamSpec22-PR4 , and HiLLiPoP23-PR4 . For comparison purposes we also show the Planck 2013 results from [22] in red where the potential trend for a best fit in the 'negative' regime was first highlighted. We clearly see that the bound on the neutrino masses changes significantly for each implementation of the likelihood, being HiLLiPoP the one giving the looser constraints. Solid and dashed lines correspond to parabolic fits where the ∆ χ 2 points up to 4 or 2 were used in the fit, respectively.  \n<!-- image -->  \nimplementations feature different levels of lensing discrepancies which are critical for neutrino mass inferences. The lensing anomaly in the Planck18-PR3 one is 2 . 8 σ , in CamSpec22-PR4 it is 1 . 7 σ , while in HiLLiPoP23-PR4 it is only 0 . 75 σ . From Fig. 3, one can see that the strength of the neutrino mass bound directly anti-correlates with the level of the anomaly, and can relax by up to a factor of 2 when going from Planck18-PR3 to HiLLiPoP23-PR4 . One can also clearly notice that the extrapolation to the unphysical region indicates preference for negative neutrino masses for Planck18-PR3 and CamSpec22-PR4 , the two likelihoods that carry some residual lensing anomaly. However, for HiLLiPoP23-PR4 , which has no statistically significant A lens anomaly, the minimum of the χ 2 is consistent with positive (albeit small) neutrino masses.  \nIn Table I we explicitly show the 95% CL bounds on the neutrino mass for the various likelihoods. Given the impact of the lensing anomaly in plik and CamSpec, their preference for a best-fit in the negative mass region is very strong. As such, our samples in the physical region lie very far from the ∆ χ 2 minimum where the Gaussian approximation holds and the extrapolation has extremely large uncertainties as shown in Fig. 3. Therefore, the simple Feldman-Cousins prescription when a Gaussian behaviour is observed cannot be implemented in these two cases. Instead, a full boostraping of the parameter space to calibrate the test statistic and correctly asses at which values of the ∆ χ 2 lay the cuts for the confidence  \nlevels of interest would be needed. This procedure would be extremely computing-expensive and is unfortunately not feasible. On the other hand, for the HiLLiPoP implementation (which does not feature a lensing anomaly), solid frequentits limits can be obtained. A direct comparison shows that this frequentists and Bayesian bounds agree within 20%.  \nThe impact of A lens - To investigate further the potential preference for negative neutrino masses and the role of the lensing anomaly in driving this preference, we show in Fig. 4 the results from analyses that vary in addition the A lens parameter that controls the lensing of temperature and polarization fluctuations ( C ϕϕ L = A lens C ϕϕ L | ΛCDM [34]). This parameter is well known to be correlated with the neutrino mass [4], and this is confirmed through our frequentist analysis, see the right panel of Fig. 15. From Fig. 4 one can see that the χ 2 parabolas become all more or less similar, and that when the A lens parameter is allowed to vary we find no preference for a negative neutrino mass. Therefore, our sampling of the ∆ χ 2 in the physical region is now closer to the minimum and we can reliably extrapolate and derive robust frequentist confidence levels also for Planck18-PR3 and CamSpec22-PR4 . These are reported in Table I. We see that the neutrino mass bounds can be relaxed by up to a factor of ∼ 2 when the A lens parameter is introduced, becoming comparable to the latest laboratory bounds. This is in agreement with the findings of [33]. Importantly, when the A lens  \nparameter is allowed to vary, the profile likelihoods resemble Gaussians and we are able to obtain frequentist limits for ∑ m ν . Direct comparison between Bayesian and frequentist limits shows a 10 -20% agreement depending upon the specific likelihood.  \nThough the A lens parameter is unphysical, this exercise clearly highlights the importance of internal inconsistencies in Planck CMB data for inferences of the neutrino mass in cosmology. Marginalizing over these anomalies, it is remarkable that Planck constrains become only as strong as direct laboratory bounds. However, we stress that this is only part of the cosmological constraining power, as we have ignored the Lensing-PR3 likelihood, which provide additional sensitivity to CMB lensing, as well as BAO and SNe data that are sensitive to the background effects of neutrinos.  \nIn Fig. 5 we show the results but now including the Lensing-PR3 lensing likelihood and with the physical condition A lens = 1. By comparing the results with those without lensing in Fig. 3 we can clearly see that while the impact on the Planck18-PR3 analysis is mild at the level of the 2 σ limit, the extrapolated behaviour to negative neutrino masses shows significantly weaker support for a negative best fit. This is because the lensing likelihood does not feature a lensing anomaly. The shift for the CamSpec22-PR4 and HiLLiPoP23-PR4 cases is substantially more pronounced and as shown in Table I the 95% CL limit improves by a factor of 1.5 when the lensing likelihood is added. We note that for the data set combination HiLLiPoP23-PR4 + Lensing-PR3 the minimization procedure was rather challenging highlighting that a nonglobal minimum of the χ 2 at around ∑ m ν ∼ 0 . 2 eV may be present. In Appendix A we show the Bayesian posterior in Fig. 11 and there is a priori no evidence for multi-modality in it.", 'IV.2. Planck + BAO: A close look at the DESI results': "As discussed in Section II, there are two critical effects of the neutrino mass in the CMB: one is on CMB lensing at small angular scales (high ℓ ), see Fig. 1, and the other is on the angular diameter distance to recombination, which lead to a strong correlation with other cosmological parameters such as Ω m , see Fig. 2. This is why including the Lensing-PR3 likelihood and BAO data is critical, as they can help break degeneracies between ∑ m ν and other cosmological parameters, and allow for an increased sensitivity of the neutrino mass when combined with CMB observations (note that geometric BAO data on their own are not sensitive to the neutrino mass).  \nIn particular, as stated in the introduction, one of the main goals of our study is to understand why the bound from Planck+DESI-Y1 is so strong, and whether it really points to negative neutrino masses.  \nIn order to gauge how unexpected is the result  \n6  \nFIG. 4. Neutrino mass profile likelihoods using the full Planck temperature and polarization data for ΛCDM allowing to vary the unphysical A lens parameter which is strongly correlated with ∑ m ν . We can see that the bounds are significantly relaxed and comparable to the KATRIN laboratory upper limit.  \n<!-- image -->  \n6  \nFIG. 5. Neutrino mass profile likelihoods using the full Planck temperature, polarization and lensing data for ΛCDM. This should be compared with Fig. 3 that does not include lensing.  \n<!-- image -->  \nfrom Planck+DESI-Y1, it is instructive to first estimate the sensitivity of this data combination to neutrino mass, as for instance done in Ref. [84]. To do so, we perform a mock data analysis, making use of the MontePython and the Fake planck realistic likelihood [84], that we combine with our own mock likelihood of DESI-Y1, that simply makes use of the covariance matrix from the data, replacing the measurements with a fiducial model prediction. For the fiducial, we study three configurations: ∑ m ν = 0 , 0 . 06 , 0 . 1  \nTABLE I. Upper limits at 95% CL on the neutrino mass within ΛCDM using various Planck likelihoods. We show the Bayesian limits compared with the two frequentist approaches (B.L. = Bounded Likelihood, and F.C. = Feldman-Cousins) described in Section III.2. We do not quote frequentist F.C. bounds for Planck18-PR3 and CamSpec22-PR4 (no lensing) due to their non-Gaussian behaviour in the physical region, which precludes a reliable extrapolation into the unphysical regime.  \nFIG. 6. Reconstructed posterior distribution of ∑ m ν in a mock analysis of Planck18-PR3 +Lensing+BAO with three configurations: ∑ m ν = 0 , 0 . 06 , 0 . 1 eV. We also show the posterior reconstructed from the real data for comparison.  \n<!-- image -->  \neV, fixing the other ΛCDM parameters as reconstructed from Planck18-PR3 +Lensing+BAO [4]. We run MCMC chains, and show the reconstructed posteriors in Fig. 6. It is clear that the real data (black line) provide significantly stronger constraints than expected, and in fact, our mock data suggests that the combination of Planck+DESI-Y1 has the sensitivity necessary to detect at small significance non-zero neutrino masses, assuming these respect the laboratory lower limits. Barring issue with our mock likelihoods, it is thus clear that there is some 'anomaly' in the real data (whether due to a fluke, or new physics). As already discussed and as we will elaborate in the following with analyses of the real data, the two results that  \nseem to be driving this unexpected constraints beyond the sensitivity are Planck likelihoods through the lensing anomaly and the DESI-Y1 data at z = 0 . 7. Indeed, when we compare instead the posterior for HiLLiPoP, with a significantly reduced lensing anomaly, and remove the DESI outliers, the posterior is much closer to the expected sensitivity as shown by the purple line in Fig. 6. When in presence of such anomalous behaviour from the actual data when compared to the expected sensitivities, an alternative method to the Feldman-Cousins analysis was proposed by Lokhov and Tkachov [85] ensuring also proper coverage when deriving the confidence intervals but conservatively replacing the anomalous results pushing into the unphysical region with the expected sensitivity in absence of a signal. Since we find that the main contributors driving the anomalous results can be identified and that their impact can be almost entirely removed for instance through the HiLLiPoP + DESI-Y1 combination, we prefer not to implement this method as it would hide the different behaviours of the different datasets and cosmological models we want to analize and compare here. However, we do find the comparison with the expected sensitivities presented in Fig. 6, on which the LT method is based, very illustrating. Furthermore, we can compare our expected sensitivities of Planck+DESI-Y1 to the ones obtained in previous forecasts using the full DESI reach [84]. We find a 1 σ sensitivity of ∼ 0 . 06 eV while for Planck+DESI-full Ref. [84] reports a 1 σ sensitivity of 0 . 04 eV. This highlights what is the level of improvement expected from the additional 5 years of BAO data, although an even higher sensitivity can of course be obtained from a full shape analysis of the matter power spectrum, see [35-38].  \nWe now turn to the real data. In Fig. 7 we show the analyses of Planck data 5 combined with: (a) DESI-Y1 BAO results, (b) the full SDSS combination (which have a similar statistical power in the Ω m vs H 0 r d plane as DESI-Y1, see Fig. 2), (c) the DESI/SDSS BAO combination that takes the BAO results at each redshift bin from the survey that has greatest statistical power, and (d) the DESI-Y1 BAO combination but without the z = 0 . 7 bin that contains data which deviate by ∼ 3 σ from ΛCDM predictions. Importantly, in all data set combinations, we consider separately the three latest Planck likelihood implementations.  \nFirst, from panel (a), we confirm that the extrapolation of the results from Planck+DESI-Y1 to the unphysical region seem to favor negative masses. However, in our approach, given uncertainties in the profile represented by the bands, the preference remains at low statistical significance. In fact, we find that it is the combination with HiLLiPoP23-PR4 that seems to lead to the largest preference (reaching at most ∼ 1 σ ). This may appear counter-intuitive given results pre-  \n6  \n6  \nFIG. 7. Profile likelihoods of the neutrino mass within ΛCDM when using Planck+BAO data comparing scenarios with DESIY1 (a), the full SDSS data (b), the combination of DESI-Y1+SDSS (c), and the DESI-Y1 data set but removing the ouliers at z = 0 . 7 (d). We can clearly notice that systematically the HiLLiPoP23 likelihood implementation gives the weakest constraints, and that the two BAO measurements at z = 0 . 7 have a significant impact on both the bound on the neutrino mass as well as the potential preference for a negative best fit.  \n<!-- image -->  \nsented in previous section. We attribute this to the fact that HiLLiPoP23-PR4 has the weakest constraining power on neutrino masses, and thus the statistical power from DESI can be more pronounced. Second, by comparing the results from Planck+DESI-Y1 and Planck+SDSS in panel (b), one can clearly see that the bound on the neutrino mass are significantly weaker in the latter case and that the preference for negative neutrino masses further decrease. This is particularly striking for HiLLiPoP23-PR4 that is now fully compatible with positive masses. However, when comparing the Planck+DESI-Y1 and Planck+DESI-Y1+SDSS (panel (c)) one notices that the profiles are quite similar and only slightly shifted towards more positive values for the latter case. This suggests that the lowz data points from DESI, that are replaced by SDSS is this analysis,  \ndo not play a significant role in the preference for negative neutrino masses. Rather, we can compare the results from Planck+DESI-Y1 and Planck+DESI-Y1no07 (panel (d)) which excludes the outliers at z = 0 . 7. We find that removing these data points loosens the bound and also makes all likelihood combinations to peak at around ∑ m ν ≃ 0, thus removing the potential preference for negative neutrino masses. Although it is obvious that removing two data points from DESI-Y1 is not a valid statistical procedure, this exercise shows that DESIY1 BAO data at z = 0 . 7 have a very significant impact on the neutrino mass bound. We note that the DESIY1+SDSS BAO combination does include these outliers in the combination as the effective volume of DESI is already larger than that of SDSS at z = 0 . 7.  \nFinally, although Fig. 7 shows that the form of the  \nprofiles (in particular the extrapolation) is different for the three Planck implementations, the differences at the level of the 95% CL bound is only at the 10 -20% level, depending upon the specific BAO data used. This is good news because it shows that the potentially large systematic shifts on ∑ m ν when using Planck CMB data alone are not present at the same level when CMB lensing and geometric BAO data are included. Our results are in qualitative agreement with the recent Bayesian analysis of [10].", 'IV.3. Planck + DESI + SN: Impact of the Dark Energy equation of state': 'So far, we have explored constraints to ∑ m ν assuming a flat ΛCDM background, and showed how BAO data can help strengthen the bound by breaking the degeneracy with Ω m . However, in models with more parameters controlling the late-time expansion history, it is expected that additional degeneracies with the neutrino masses will appear. Chief amongst those is the well-known degeneracy with the equation of state of dark energy, w (see e.g. [86].) Given the tentative evidence for a time-varying equation of state of dark energy from DESI when combined with SN data [1], it is relevant to explore how it can impact the bound on the neutrino masses. Following [1], we model the equation of state of dark energy to vary according to the Chevalier-Polarski-Linder parametrization w ( a ) = w 0 +(1 -a ) w a [87, 88], where a is the scale factor, and vary w 0 ∈ [ -3 , 2] and w a ∈ [ -3 , 1].  \nIn Fig. 8 we show the likelihood profile of ∑ m ν built from the combination of Planck, DESI-Y1 BAO and the SN Pantheon sample, in ΛCDM (in black) and in the w 0 w a CDM cosmology (in green). One can see that when the equation of state of dark energy is allowed to vary the bound on the neutrino mass is relaxed, in good agreement with [1]. We find, however, that the 95% bound can be roughly 30% weaker for HiLLiPoP than with Plik or CamSpec. Interestingly we also note that the potential evidence for a negative neutrino mass vanishes.  \nSomething important to highlight is that, while allowing for the equation of state of dark energy to vary relaxes the bound on the neutrino mass, the best fit for the equation of state differs significantly from the cosmological constant value and in particular suggests that w 0 ≳ -1 today. Thus, while the bound may be relaxed, large neutrino masses would require dark energy to behave very differently than a cosmological constant. We note that if one restricts the analysis to constant equation of state, namely w ( a ) = w 0 where only w 0 is allowed to vary, the bound on the neutrino mass remains very similar to that in ΛCDM, see [1].  \nLastly, we investigate whether varying the A lens parameter may remove the preference for negative neutrino masses, despite the inclusion of DESI-Y1 BAO and SN data data. Our results are shown in Fig. 8 in blue. One  \ncan notice that this shifts the best fit to the positive regime and that the bound becomes again a factor of ∼ 2 weaker than when compared to the standard case where A lens = 1. This suggests that, regardless of the behavior of DESI data, it is the lensing anomaly that dominates the preference for negative neutrino masses. Note though, that removing (most of) the constraining power from lensing by including A lens does not remove all the sensitivity to neutrino masses altogether. While the constraints relax, they remain significantly stronger than laboratory ones when BAO and SNIa data are included, in the ball park of ∑ m ν ≲ 0 . 2 -0 . 3 eV depending on which CMB likelihood is used.', 'IV.4. Frequentists vs Bayesian Limits: The impact of statistical choices on the neutrino mass bound': "Until now, we have focused our attention on profile likelihoods, as these allowed us to investigate the preference for negative neutrino masses. In this section we will study how the confidence intervals built from the profile likelihood compare with the Bayesian credible intervals built from the posteriors. While the two statistical approaches address distinct questions and thus do not need to necessarily coincide, we find the comparison useful to understand how dependent is the derived constraint on the choice of the statistical procedure. Moreover, if the likelihood is Gaussian and when assuming flat, noninformative priors on the parameter under study, the two approaches should coincide. Therefore, the comparison allows to asses the role of possible prior effects in the bound given by Eq. (1) as well as deviations from Gaussianity. Notice that the extrapolations to negative masses performed in our work as well as in Ref. [24] and [10] rely on the Gaussianity of the likelihood and posterior respectively. Similarly, the study of Ref. [21] relies on a linear expansion of the dependence of the observables with ∑ m ν which, per Wilks theorem, would imply a Gaussian behaviour. The fact that these analyses do not coincide (see section V for a full comparison) implies that the Gaussian behaviour approximation is violated at some degree. As we will show and discuss below, we find that the differences between the Bayesian and frequentist approaches are at the 10% level, which provides a useful estimation of the size of the uncertainties introduced by these approximations. We refer to Appendix A for the posterior distributions.  \nOur results are summarized in Table II for analyses that combine Planck+BAO data, and in Table III for those that also include SN data from the Pantheon sample. These tables include three estimates of the bound to neutrino masses: the Bayesian limit at 95% CL, those derived using Feldman-Cousins procedure (F.C.), as well as those using the naive bounded maximum likelihood (B.L.) (∆ χ 2 = 3 . 84), all at the same confidence level. Firstly, we generally notice a very good agreement be-  \n6  \n<!-- image -->  \n<!-- image -->  \n2  \n6  \nFIG. 8. Neutrino mass profile likelihoods for Planck+DESI-Y1+Pantheon+ data set combinations. We show ΛCDM in black, varying the equation of state of dark energy in green, and allowing for A lens to vary in blue. In the left panel we show the results for plik, in the middle for CamSpec, and in the right panel for Hillipop. We clearly see a similar behaviour for all of them and the potential preference for a negative best fit to dissapear when the equation of state of dark energy is allowed to vary.  \n<!-- image -->  \nTABLE II. Upper limits at 95% CL on the neutrino mass within ΛCDM using various data set combinations of Planck + BAO data. We show the Bayesian limits compared with the two frequentist approaches (B.L. = Bounded Likelihood, and F.C. = Feldman-Cousins) described in Section III.2. We also report the Gaussian fit for our profile likelihoods.  \ntween the two frequentists approaches, with differences between them at the ≲ 5% level only. This suggests that the fact that the minimum lie beyond the physical region does not significantly affect the bounds to neutrino masses. Secondly, and interestingly, we also notice a very good agreement between the frequentists and Bayesian limits. In fact, we find, that the frequentists limits are in many cases ∼ 10% stronger than the Bayesian ones. For example, considering the data set combination of Plik+DESI, we find at 95% CL:  \n∑ m ν < 0 . 084 eV [Bayesian] , (5a)  \n∑ m ν < 0 . 071 eV [Feldman -Cousins] . (5c)  \n∑ m ν < 0 . 074 eV [Bounded -Likelihood] , (5b)  \nOne clearly sees that the three are very similar, with the frequentist ones being slightly more stringent. While the  \ntwo approaches need not necessarily agree, this could be due to two effects. First, it is possible that there are mild prior effects in the Bayesian analysis, that go in the direction of relaxing the bound. Second, it can be difficult to find the absolute minimum of the χ 2 for each simulated value of ∑ m ν for such a large parameter space. If the simulated annealing methods fails to cool to the absolute minimum, the slightly larger values of the χ 2 would lead to slightly tighter frequentist constraints. Nevertheless and regardless of its origin, this effect is only around the 10% level, and we thus conclude that the constraints are robust to the choice of statistical method up to that level of difference. Let us additionally note that the Bayesian constraint we derive here is slightly different than Eq. (1). This is because the DESI collaboration used more constraining CMB lensing data, combining Planck lensing PR4 with ACT lensing, rather than Planck lensing PR3 as we do here. Nevertheless, we do not expect that using  \n6  \nthis lensing data would change the overall trend.  \nSo far we have included in our analyses values for neutrino masses down to the massless limit, ∑ m ν = 0, but we know from the laboratory that there are physical boundaries at either ∑ m ν = 0 . 06 eV for NO or at ∑ m ν = 0 . 10 eV for IO. To gauge the impact of those experimental lower limits on the cosmological neutrino mass bound, we run dedicated Bayesian analyses restricting the prior to ∑ m ν following either the NO or IO constraints. For the frequentist limit, it is sufficient to consider these boundaries as lower limits in our ∆ χ 2 curves. This procedure yields:  \n∑ m ν < 0 . 121 eV [NO -Bayesian] , (6a)  \n∑ m ν < 0 . 096 eV [NO -Feldman -Cousins] , (6c)  \n∑ m ν < 0 . 106 eV [NO -Bounded -Likelihood] , (6b)  \nand for the inverted ordering case:  \n∑ m ν < 0 . 152 eV [IO -Bayesian] , (7a)  \n∑ m ν < 0 . 127 eV [IO -Feldman -Cousins] . (7c)  \n∑ m ν < 0 . 138 eV [IO -Bounded -Likelihood] , (7b)  \nSince in these scenarios the physical boundary is further away from the best fit of the ∆ χ 2 , the Feldman-Cousins correction becomes more relevant and we observe a larger difference compared to the naive bound one would derive simply assuming the applicability of Wilk's theorem, although it is still within 10%. The difference between the frequentist and Bayesian constraints also increases, with up to 20% difference between the Feldman-Cousins result and the Bayesian posterior. Let us stress that, for this particular dataset, the inverted ordering assumption has a p -value of only 1%.  \nImportantly, we have highlighted before that there are two effects that significantly pull the bound on the neutrino mass in Eq. (1): i) the lensing anomaly present in some of the Planck likelihoods, and ii) the outliers in DESI-Y1 at z = 0 . 7. In this context, to be maximally conservative, one can consider the combination of HiLLiPoP+DESIY1no07 for which there is no lensing anomaly in the Planck likelihood and where the outliers in DESI-Y1 data have been removed. The relevant Bayesian and frequentist limits from HiLLiPoP+DESIY1no07 read:  \n∑ m ν < 0 . 125 eV [Bayesian] , (8a)  \n∑ m ν < 0 . 114 eV [Feldman -Cousins] . (8c)  \n∑ m ν < 0 . 114 eV [Bounded -Likelihood] , (8b)  \nHere we can see again a ∼ 10% agreement between Bayesian and frequentist approaches.  \nWhen considering the physical boundary for normal  \nordering, the constraints read:  \n∑ m ν < 0 . 160 eV [NO -Bayesian] , (9a)  \n∑ m ν < 0 . 125 eV [NO -Feldman -Cousins] . (9c)  \n∑ m ν < 0 . 132 eV [NO -Bounded -Likelihood] , (9b)  \nand for the inverted ordering case:  \n∑ m ν < 0 . 179 eV [IO -Bayesian] , (10a)  \n∑ m ν < 0 . 146 eV [IO -Feldman -Cousins] . (10c)  \n∑ m ν < 0 . 156 eV [IO -Bounded -Likelihood] , (10b)  \nWhere, as before, the FC corrections become more relevant given the larger distance between the best fit and the physical boundary.  \nFinally, we can compare our Bayesian bounds with other recent studies. In particular, our limit for the PlanckPR3+DESI+Pantheon perfectly agrees with the one reported in [9]. Ref. [10] also presented analyses including various versions of the new Planck likelihoods. For the data combination Plik+DESI and HiLLiPoP+DESI (with or without SN), we find bounds that are ∼ 10 -20% looser that those reported in [10]. However, for the case HiLLiPoP+SDSS/DESI we find the same limit as [10]. Given that we agree with Ref. [9] when the very same data is considered, but also with Ref. [10] when a subset of the DESI data set is considered, we conjecture that the differences in the limits may stem from a different implementation of the full DESI likelihood. Our implementation matches the one in the Cobaya public repository [89].", "V. DO COSMOLOGICAL DATA PREFER A 'NEGATIVE' NEUTRINO MASS?": "It has been recently emphasized in [19-21] that cosmological data may prefer 'negative' neutrino masses and that such a preference for unphysical values, if not due to a statistical fluke or an unknown systematic effect, may hint for new physics in cosmology. In their analysis, the Authors of [19, 20] define an 'effective neutrino mass' which, when taking negative values, generates an increase of power in the CMB lensing potential, an effect that is opposite to that of positive neutrino masses. In [21] an alternative definition is adopted which also accounts for the other effects of neutrino masses in cosmological observables (in particular on the angular diameter distance, but not only). In this paper, to gauge the preference for neutrino mass, we have simply extrapolated through a Gaussian fit the behaviour of our profiled ∆ χ 2 . Even though the extrapolation introduces some degree of uncertainty, that we have quantified with the procedure detailed in Section III.2, it provides a complementary way to assess the preference for negative values  \nTABLE III. Upper limits at 95% CL on the neutrino mass using Planck + BAO + SN data within extended ΛCDM models, including a time-varying equation of state of dark energy, as well as varying the A lens parameter. We show the Bayesian limits compared with the two frequentist approaches (B.L. = Bounded Likelihood, and F.C. = Feldman-Cousins) described in Section III.2. We also report the Gaussian fit for our profile likelihoods, obtained from the points below ∆ χ 2 = 4.  \n6  \nFIG. 9. Profile likelihoods for the data set combinations of Planck18-PR3+DESI-Y1 (black), HiLLiPoP23-PR4 (blue), Planck18PR3+DESI-Y1-no07 (red), HiLLiPoP23-PR4+DESI-Y1-no07 (green), and compared with χ 2 eff = -2 log P from [20] (purple) and [21] (orange), which have treatments for 'negative' neutrino masses. By comparing the black and blue curves we can clearly see that the bound on the neutrino mass gets relaxed if the HiLLiPoP likelihood (which does not contain a lensing anomaly) is used. However, it is clear from the extrapolated parabolas that there is still some preference for a negative neutrino mass. This, however, disappears when the DESI BAO data at z = 0 . 7 which contains a ∼ 3 σ outlier is removed (see red and green curves).  \n<!-- image -->  \ndisplayed by the different datasets analyzed without the need for an explicit (arbitrary) modelling of the effect of negative neutrino masses. A similar approach, extrapolating instead a Gaussian fit to the posterior distribution, was adopted in Refs. [24] and [10].  \nOur results agree overall fairly well with all previous works [10, 19-21] considering the different treatments performed in each of them. Our main results and the comparison to previous analyses are shown in Fig. 9  \n(see also App. B). In particular, Ref. [21] argued that the Gaussian extrapolation in [10] underestimates the preference for negative masses of present data. Indeed, Refs. [19-21] find rather more negative best fits through their analyses with their respective 'effective masses' peaking around -0 . 15 eV (see purple and orange lines in Fig. 9). We do find some preference for negative neutrino masses in datasets including both the Planck likelihoods affected by the lensing anomaly (Plik 2018 in Fig. 9 but also CamSpec in previous sections) and the full DESI Y1  \nBAO data (black line). However, this preference is significantly weaker, with a best fit at ∑ m ν ≃ -0 . 05 eV in agreement with the extrapolation of the posterior shown in [10] and the comparison and discussion in [21].  \nMore importantly, Ref. [21] only considered the CamSpec Planck likelihood. We find here that this dataset has a preference for negative neutrino masses, because of the presence of a residual lensing anomaly, and that it can be reduced when using alternative likelihoods. Indeed, we find a preference for negative neutrino masses when Planck 2018 is analyzed on its own. Switching for the HiLLiPoP23-PR4 version of Planck (largely unaffected by the lensing anomaly), there is no preference for negative masses, and the bound significantly relaxes. The preference re-appears for HiLLiPoP23-PR4 when DESI-Y1 BAO data are included because these data have a lower constraining power to ∑ m ν . Nevertheless, Fig. 9 shows that the bound on ∑ m ν relaxes when changing from plik 2018 to HiLLiPoP (black vs blue lines) or when removing the z = 0 . 7 outlier bin in DESI (black vs red lines). Finally, when the z = 0 . 7 bin is removed and the HiLLiPoP implementation are both considered, the preference for negative masses essentially disappears (green line) and the bound is about a factor 2 weaker. Let us note that the authors of [20] did investigate the persistence of their effective negative neutrino mass signal when PR4 likelihoods are used. However, the analyses were not performed for HiLLiPoP23-PR4 but rather for the CamSpec22-PR4 likelihood, where the lensing anomaly is still present, albeit at a reduced level. Ref. [20] does show that the preference for negative masses is relaxed in this scenario at a level that, naively, seems consistent with the reduced presence of the anomaly in CamSpec. In this context, it would be very interesting to explore the potential preference for a 'negative' neutrino mass following the implementation of [21] and [19, 20] but using the HiLLiPoP likelihood and also in the absence of the z = 0 . 7 outliers that we have shown pull the neutrino mass bound significantly. Given our results, and if our hypotheses are correct, the preference for negative neutrino masses would vanish in that case.  \nFrom our results, we conclude that there is no significant nor compelling evidence for negative neutrino masses from cosmology. We do confirm some mild preference when Planck likelihoods affected by the lensing anomaly are adopted, pointing towards a potential residual systematic effect in these implementations. We also find that the z = 0 . 7 bin of DESI, in some tension with Planck data, push towards negative values. If these outliers are confirmed, it will be interesting to see how the preference for negative masses evolves in the future. For now, and in agreement with [10], we conclude there is no compelling evidence for negative masses from present datasets and that present constraints are still perfectly compatible with the results from neutrino oscillations. We believe that a key check to see if cosmological data does prefer effective 'negative' neutrino masses is to do a cosmological analysis using the HiLLiPoP likelihood  \nimplementation.", 'VI. CONCLUSIONS': "Recent cosmological constraints on the sum of neutrino masses ( ∑ m ν ) already disfavour its minimum value if the mass ordering is inverted and allow very little parameter space even for normal ordering [1]. Some analyses combining additional datasets even start to rule out the minimum value of ∑ m ν for a normal ordering [9], seemingly implying tension between cosmological probes and neutrino oscillation experiments. In this context, we have critically assessed the constraints on ∑ m ν from cosmological observables investigating their robustness against different statistical methods and determining which observables mainly dominate the budding tension with oscillations results.  \nTo this end, we have analyzed the constraints on ∑ m ν from different datasets both through the usual Bayesian approach and through complementary frequentist methods. For the latter, given the proximity of the best fit to the physical boundary of the parameter space, we have extrapolated the likelihood profiles into the unphysical region so as to apply the appropriate correction [73]. This also allowed us to explore the preference of some datasets for negative neutrino masses recently reported in [10, 1921].  \nFrom the simulations and results presented in Section IV we conclude that:  \n- · Bayesian and frequentist limits on the neutrino mass agree very well and within 10% precision The agreement between the Bayesian and frequentist bounds on ∑ m ν is very good. Regarding the frequentist constraints, an extrapolation into the unphysical region of negative neutrino masses was performed to apply the appropriate correction [73] to the confidence intervals. Even though this procedure introduces a certain degree of uncertainty, we quantified the range of variation with the number of points for larger values of ∑ m ν and χ 2 included and verified that the impact on the confidence intervals derived is negligible. Moreover, we also compared these with the naive cuts at the values of ∆ χ 2 implied by Wilks theorem without venturing into the unphysical region and found almost identical results. When comparing these frequentist confidence intervals with the Bayesian credible regions with flat priors on ∑ m ν we find good agreement between both approaches but consistently tighter results for the frequentist method, by about ∼ 10%. We thus conclude that the cosmological constraints on ∑ m ν are robust at this level against variations of the statistical method employed to derive them.\n- · The potential preference for 'negative' neu-  \nino mass bounds is strongly dependent upon the Planck likelihood implementation used in the analysis When investigating the constraints derived exclusively from Planck data, we find radically different behaviours between the Plik, CamSpec and HiLLiPoP implementations. Indeed, the preference for unphysical negative masses and hence the strongest constraints on ∑ m ν are obtained for Plik and, to a lesser degree CamSpec, and absent for HiLLiPoP, with a factor ∼ 2 weaker constrain. This tendency for weaker constraints for the CamSpec or HiLLiPoP implementations was also recently confirmed in Ref. [10]. These results seem closely correlated to the presence of the lensing anomaly in these datasets, which is reduced in the CamSpec and almost absent in the HiLLiPoP implementations. To confirm this observation, we performed analyses including the A lens parameter and marginalising over its value. Interestingly, these analyses yielded the weakest constraints for the datasets more strongly affected by the anomaly. In particular, the Plik and CamSpec constraints in combination with other datasets relax by more than a factor 2. We conclude that, when the lensing anomaly is not present in the Planck likelihood either through use of the HiLLiPoP implementation or the inclusion and marginalization over A lens , the combined results of Planck with DESI and Pantheon relax to about ∑ m ν < 0 . 11 eV or even ∑ m ν < 0 . 2 eV, both at 95% CL, and in good agreement with both the normal and inverted orderings currently allowed by neutrino oscillation data.  \n- · DESI-Y1 BAO measurements at z = 0 . 7 pull the neutrino mass bound significantly and also induce some preference for negative values When combining Planck with BAO data from DESI we observe a mild preference for unphysical negative neutrino masses, and hence a more stringent constraint. This mild preference for negative masses even appears for the HiLLiPoP implementation, for which it is absent when analyzed alone or in combination with other datasets. Upon closer examination, this preference seems to be driven mostly the BAO measurements at z = 0 . 7, out of which the angular one is in ∼ 3 σ tension with Planck expectations. Indeed, when removing this bin, the preference for negative masses disappears and the corresponding constraints relax accordingly. Similarly, if SDSS BAO data is used instead of DESI, constraints relax by about a 50% in all cases considered. Thus, given the unexpectedly strong constraint derived when including the full DESI dataset, it will be very interesting to confirm if the present trend is confirmed with higher statistics from upcoming DESI data releases.\n- · The hint of dynamical dark energy relaxes  \nthe neutrino mass bound and removes the preference for 'negative' masses The results presented by the DESI collaboration [1] favour a dynamical equation of state for Dark Energy parametrized through w 0 and w a when extending the ΛCDM. We find that our analyses including this effect also relax significantly the constraints on ∑ m ν by a factor 2 in all cases analyzed, in agreement with the results of [1, 19, 21]. In addition, the preference for negative neutrino mass disappears, suggesting that it may be an artifact of using the wrong cosmological model. More data are necessary to test whether the hint of dynamical dark energy, and the related artificially strong constraints to neutrino mass, are due to a statistical fluke, a systematic error, or a real deviation from ΛCDM.  \nAll in all, we find it is premature to infer significant tension between present cosmological probes and neutrino oscillation data, and that the preference for unphysical negative masses is not yet compelling. Acritical test of this statement would be to have the results for the models including 'effective negative' neutrino masses using the HiLLiPoP likelihood implementation of the Planck data which does not feature a lensing anomaly. In any case, this is a very exciting period as we should be on the brink of a discovery of the absolute neutrino mass from cosmology. This would represent an indirect confirmation of the cosmic neutrino background and set a target for laboratory experiments aiming to kinematically measure the neutrino mass or, even more interestingly, its potential Majorana nature. If in the upcoming years a discovery is not made despite unprecedentedly small statistical and systematic errors, it will be a clear call for a change of paradigm, potentially signaling new physics such as neutrino decays [3945], non-standard thermodynamic histories [46-51], or time varying masses [52-59]. It is clear that upcoming years are of utmost importance for neutrino cosmology, as whether neutrino masses are detected or not, we are on the verge of a major breakthrough.  \nAcknowledgments. EFM acknowledges very interesting discussions with Stefano Gariazzo and thanks him for pointing out relevant references. We also warmly thank Thomas Schwetz for very interesting comments and suggestions for our analysis. This project has received support from the European Union's Horizon 2020 research and innovation programme under the Marie Skglyph[suppress]lodowska-Curie grant agreement No 860881-HIDDeN and No 101086085 - ASYMMETRY, and from the Spanish Research Agency (Agencia Estatal de Investigaci'on) through the Grant IFT Centro de Excelencia Severo Ochoa No CEX2020001007-S and Grant PID2022 137127NB-I00 funded by MCIN/AEI/10.13039/501100011033. We acknowledge support from the HPC-Hydra cluster at IFT and the HTCondor cluster at CERN. XM acknowledges fund-  \nng from the European Union's Horizon Europe Programme under the Marie Skglyph[suppress]lodowska-Curie grant agree-  \n- [1] A. G. Adame et al. (DESI), (2024), arXiv:2404.03002 [astro-ph.CO].\n- [2] A. G. Adame et al. (DESI), (2024), arXiv:2404.03000 [astro-ph.CO].\n- [3] A. G. Adame et al. (DESI), (2024), arXiv:2404.03001 [astro-ph.CO].\n- [4] N. Aghanim et al. (Planck), Astron. Astrophys. 641 , A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO].\n- [5] N. Aghanim et al. (Planck), Astron. Astrophys. 641 , A1 (2020), arXiv:1807.06205 [astro-ph.CO].\n- [6] J. Carron, M. Mirmelstein, and A. Lewis, JCAP 09 , 039 (2022), arXiv:2206.07773 [astro-ph.CO].\n- [7] M. S. Madhavacheril et al. (ACT), Astrophys. J. 962 , 113 (2024), arXiv:2304.05203 [astro-ph.CO].\n- [8] F. J. Qu et al. (ACT), Astrophys. J. 962 , 112 (2024), arXiv:2304.05202 [astro-ph.CO].\n- [9] D. Wang, O. Mena, E. Di Valentino, and S. Gariazzo, (2024), arXiv:2405.03368 [astro-ph.CO].\n- [10] I. J. Allali and A. Notari, (2024), arXiv:2406.14554 [astro-ph.CO].\n- [11] M. Aker et al. , (2024), arXiv:2406.13516 [nucl-ex].\n- [12] I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, T. Schwetz, and A. Zhou, JHEP 09 , 178 (2020), arXiv:2007.14792 [hep-ph].\n- [13] P. F. de Salas, D. V. Forero, S. Gariazzo, P. Mart'ınezMirav'e, O. Mena, C. A. Ternes, M. T'ortola, and J. W. F. Valle, JHEP 02 , 071 (2021), arXiv:2006.11237 [hep-ph].\n- [14] F. Capozzi, E. Di Valentino, E. Lisi, A. Marrone, A. Melchiorri, and A. Palazzo, Phys. Rev. D 104 , 083031 (2021), arXiv:2107.00532 [hep-ph].\n- [15] F. Simpson, R. Jimenez, C. Pena-Garay, and L. Verde, JCAP 06 , 029 (2017), arXiv:1703.03425 [astro-ph.CO].\n- [16] S. Gariazzo, M. Archidiacono, P. F. de Salas, O. Mena, C. A. Ternes, and M. T'ortola, JCAP 03 , 011 (2018), arXiv:1801.04946 [hep-ph].\n- [17] S. Gariazzo et al. , JCAP 10 , 010 (2022), arXiv:2205.02195 [hep-ph].\n- [18] S. Gariazzo, O. Mena, and T. Schwetz, Phys. Dark Univ. 40 , 101226 (2023), arXiv:2302.14159 [hep-ph].\n- [19] N. Craig, D. Green, J. Meyers, and S. Rajendran, (2024), arXiv:2405.00836 [astro-ph.CO].\n- [20] D. Green and J. Meyers, (2024), arXiv:2407.07878 [astroph.CO].\n- [21] W. Elbers, C. S. Frenk, A. Jenkins, B. Li, and S. Pascoli, (2024), arXiv:2407.10965 [astro-ph.CO].\n- [22] P. A. R. Ade et al. (Planck), Astron. Astrophys. 566 , A54 (2014), arXiv:1311.1657 [astro-ph.CO].\n- [23] F. Couchot, S. Henrot-Versill'e, O. Perdereau, S. Plaszczynski, B. Rouill'e d'Orfeuil, M. Spinelli, and M. Tristram, Astron. Astrophys. 606 , A104 (2017), arXiv:1703.10829 [astro-ph.CO].\n- [24] S. Alam et al. (eBOSS), Phys. Rev. D 103 , 083533 (2021), arXiv:2007.08991 [astro-ph.CO].\n- [25] D. Brout et al. , Astrophys. J. 938 , 110 (2022), arXiv:2202.04077 [astro-ph.CO].  \nment no. 101066105-PheNUmenal. The work of DNT was supported by the Spanish MIU through the National Program FPU (grant number FPU20/05333).  \n- [26] Y. Akrami et al. (Planck), Astron. Astrophys. 643 , A42 (2020), arXiv:2007.04997 [astro-ph.CO].\n- [27] G. Efstathiou and S. Gratton, (2019), 10.21105/astro.1910.00483, arXiv:1910.00483 [astro-ph.CO].\n- [28] E. Rosenberg, S. Gratton, and G. Efstathiou, Mon. Not. Roy. Astron. Soc. 517 , 4620 (2022), arXiv:2205.10869 [astro-ph.CO].\n- [29] M. Tristram et al. , Phys. Rev. D 105 , 083524 (2022), arXiv:2112.07961 [astro-ph.CO].\n- [30] M. Tristram et al. , Astron. Astrophys. 682 , A37 (2024), arXiv:2309.10034 [astro-ph.CO].\n- [31] S. Roy Choudhury and S. Hannestad, JCAP 07 , 037 (2020), arXiv:1907.12598 [astro-ph.CO].\n- [32] P. Motloch and W. Hu, Phys. Rev. D 101 , 083515 (2020), arXiv:1912.06601 [astro-ph.CO].\n- [33] E. Di Valentino and A. Melchiorri, Astrophys. J. Lett. 931 , L18 (2022), arXiv:2112.02993 [astro-ph.CO].\n- [34] E. Calabrese, A. Slosar, A. Melchiorri, G. F. Smoot, and O. Zahn, Phys. Rev. D 77 , 123531 (2008), arXiv:0803.2309 [astro-ph].\n- [35] A. Font-Ribera, P. McDonald, N. Mostek, B. A. Reid, H.-J. Seo, and A. Slosar, JCAP 05 , 023 (2014), arXiv:1308.4164 [astro-ph.CO].\n- [36] L. Amendola et al. (Euclid Theory Working Group), Living Rev. Rel. 16 , 6 (2013), arXiv:1206.1225 [astroph.CO].\n- [37] A. Aghamousa et al. (DESI), (2016), arXiv:1611.00036 [astro-ph.IM].\n- [38] A. Chudaykin and M. M. Ivanov, JCAP 11 , 034 (2019), arXiv:1907.06666 [astro-ph.CO].\n- [39] M. Escudero, J. Lopez-Pavon, N. Rius, and S. Sandner, JHEP 12 , 119 (2020), arXiv:2007.04994 [hep-ph].\n- [40] Z. Chacko, A. Dev, P. Du, V. Poulin, and Y. Tsai, JHEP 04 , 020 (2020), arXiv:1909.05275 [hep-ph].\n- [41] M. Escudero and M. Fairbairn, Phys. Rev. D 100 , 103531 (2019), arXiv:1907.05425 [hep-ph].\n- [42] Z. Chacko, A. Dev, P. Du, V. Poulin, and Y. Tsai, Phys. Rev. D 103 , 043519 (2021), arXiv:2002.08401 [astroph.CO].\n- [43] G. Franco Abell'an, Z. Chacko, A. Dev, P. Du, V. Poulin, and Y. Tsai, JHEP 08 , 076 (2022), arXiv:2112.13862 [hep-ph].\n- [44] G. Barenboim, J. Z. Chen, S. Hannestad, I. M. Oldengott, T. Tram, and Y. Y. Y. Wong, JCAP 03 , 087 (2021), arXiv:2011.01502 [astro-ph.CO].\n- [45] J. Z. Chen, I. M. Oldengott, G. Pierobon, and Y. Y. Y. Wong, Eur. Phys. J. C 82 , 640 (2022), arXiv:2203.09075 [hep-ph].\n- [46] Y. Farzan and S. Hannestad, JCAP 02 , 058 (2016), arXiv:1510.02201 [hep-ph].\n- [47] M. Escudero, T. Schwetz, and J. Terol-Calvo, JHEP 02 , 142 (2023), [Addendum: JHEP 06, 119 (2024)], arXiv:2211.01729 [hep-ph].\n- [48] J. Alvey, M. Escudero, and N. Sabti, JCAP 02 , 037 (2022), arXiv:2111.12726 [astro-ph.CO].\n- [49] A. Cuoco, J. Lesgourgues, G. Mangano, and S. Pastor, Phys. Rev. D 71 , 123501 (2005), arXiv:astro-ph/0502465.\n- [50] I. M. Oldengott, G. Barenboim, S. Kahlen, J. Salvado, and D. J. Schwarz, JCAP 04 , 049 (2019), arXiv:1901.04352 [astro-ph.CO].\n- [51] J. J. Renk et al. (GAMBIT Cosmology Workgroup), JCAP 02 , 022 (2021), arXiv:2009.03286 [astro-ph.CO].\n- [52] G. Dvali and L. Funcke, Phys. Rev. D 93 , 113002 (2016), arXiv:1602.03191 [hep-ph].\n- [53] G. Dvali, L. Funcke, and T. Vachaspati, Phys. Rev. Lett. 130 , 091601 (2023), arXiv:2112.02107 [hep-ph].\n- [54] C. S. Lorenz, L. Funcke, E. Calabrese, and S. Hannestad, Phys. Rev. D 99 , 023501 (2019), arXiv:1811.01991 [astroph.CO].\n- [55] C. S. Lorenz, L. Funcke, M. Loffler, and E. Calabrese, Phys. Rev. D 104 , 123518 (2021), arXiv:2102.13618 [astro-ph.CO].\n- [56] I. Esteban and J. Salvado, JCAP 05 , 036 (2021), arXiv:2101.05804 [hep-ph].\n- [57] I. Esteban, O. Mena, and J. Salvado, Phys. Rev. D 106 , 083516 (2022), arXiv:2202.04656 [astro-ph.CO].\n- [58] M. Sen and A. Y. Smirnov, JCAP 01 , 040 (2024), arXiv:2306.15718 [hep-ph].\n- [59] M. Sen and A. Y. Smirnov, (2024), arXiv:2407.02462 [hep-ph].\n- [60] J. Lesgourgues, G. Mangano, G. Miele, and S. Pastor, Neutrino Cosmology (Cambridge University Press, 2013).\n- [61] A. D. Dolgov, Phys. Rept. 370 , 333 (2002), arXiv:hepph/0202122.\n- [62] J. Lesgourgues and S. Pastor, Phys. Rept. 429 , 307 (2006), arXiv:astro-ph/0603494.\n- [63] S. Hannestad, Prog. Part. Nucl. Phys. 65 , 185 (2010), arXiv:1007.0658 [hep-ph].\n- [64] Y. Y. Y. Wong, Ann. Rev. Nucl. Part. Sci. 61 , 69 (2011), arXiv:1111.1436 [astro-ph.CO].\n- [65] M. Lattanzi and M. Gerbino, Front. in Phys. 5 , 70 (2018), arXiv:1712.07109 [astro-ph.CO].\n- [66] M. Archidiacono, T. Brinckmann, J. Lesgourgues, and V. Poulin, JCAP 02 , 052 (2017), arXiv:1610.09852 [astro-ph.CO].\n- [67] J. Lesgourgues and L. Verde (Particle Data Group, Neutrinos in Cosmology Review), PTEP 2022 , 083C01 (2022).\n- [68] E. Allys et al. (LiteBIRD), PTEP 2023 , 042F01 (2023), arXiv:2202.02773 [astro-ph.IM].\n- [69] M. Tristram et al. , Astron. Astrophys. 647 , A128 (2021), arXiv:2010.01139 [astro-ph.CO].\n- [70] G. Efstathiou (2024) arXiv:2406.12106 [astro-ph.CO].\n- [71] D. Rubin et al. , (2023), arXiv:2311.12098 [astro-ph.CO].\n- [72] T. M. C. Abbott et al. (DES), (2024), arXiv:2401.02929 [astro-ph.CO].\n- [73] G. J. Feldman and R. D. Cousins, Phys. Rev. D 57 , 3873 (1998), arXiv:physics/9711021.\n- [74] D. Blas, J. Lesgourgues, and T. Tram, JCAP 07 , 034 (2011), arXiv:1104.2933 [astro-ph.CO].\n- [75] T. Brinckmann and J. Lesgourgues, (2018), arXiv:1804.07261 [astro-ph.CO].\n- [76] B. Audren, J. Lesgourgues, K. Benabed, and S. Prunet, JCAP 1302 , 001 (2013), arXiv:1210.7183 [astro-ph.CO].\n- [77] A. Lewis, Phys. Rev. D 87 , 103529 (2013), arXiv:1304.4473 [astro-ph.CO].\n- [78] A. Lewis, (2019), arXiv:1910.13970 [astro-ph.IM].\n- [79] T. Karwal, Y. Patel, A. Bartlett, V. Poulin, T. L. Smith, and D. N. Pfeffer, (2024), arXiv:2401.14225", '[astro-ph.CO].': "- [80] E. B. Holm, A. Nygaard, J. Dakin, S. Hannestad, and T. Tram, (2023), arXiv:2312.02972 [astro-ph.CO].\n- [81] S. Henrot-Versill'e, O. Perdereau, S. Plaszczynski, B. R. d'Orfeuil, M. Spinelli, and M. Tristram, (2016), arXiv:1607.02964 [astro-ph.CO].\n- [82] F. James and M. Roos, Comput. Phys. Commun. 10 , 343 (1975).\n- [83] W. A. Rolke, A. M. Lopez, and J. Conrad, Nucl. Instrum. Meth. A 551 , 493 (2005), arXiv:physics/0403059.\n- [84] T. Brinckmann, D. C. Hooper, M. Archidiacono, J. Lesgourgues, and T. Sprenger, JCAP 01 , 059 (2019), arXiv:1808.05955 [astro-ph.CO].\n- [85] A. V. Lokhov and F. V. Tkachov, Phys. Part. Nucl. 46 , 347 (2015), arXiv:1403.5429 [physics.data-an].\n- [86] A. Boyle and E. Komatsu, JCAP 03 , 035 (2018), arXiv:1712.01857 [astro-ph.CO].\n- [87] M. Chevallier and D. Polarski, Int. J. Mod. Phys. D 10 , 213 (2001), arXiv:gr-qc/0009008.\n- [88] E. V. Linder, Phys. Rev. Lett. 90 , 091301 (2003), arXiv:astro-ph/0208512.\n- [89] J. Torrado and A. Lewis, JCAP 05 , 057 (2021), arXiv:2005.05290 [astro-ph.IM].", 'Appendix A: Posterior distributions': 'In this Appendix we draw the ∑ m ν posterior distributions for some of our analyses. In particular, we present in Fig. 10 the Bayesian analogue of Fig. 3. Secondly, the Bayesian analyses for the different Planck+BAO combinations are shown in Fig. 12 (see Fig. 7 for their frequentist counterpart). Lastly, the Planck+DESI-Y1+Pantheon+ posteriors are presented in Fig. 13, where we are using the same datasets as in the frequentist analyses of Fig. 8.  \nFor the comparison between the bounds extracted from the Bayesian and frequentist approaches, we refer the reader to Tables I, II and III.  \nP  \nFIG. 10. Posterior distributions for the same data set combinations as in Fig. 3 but analyzed within a Bayesian framework.  \n<!-- image -->  \nP  \nFIG. 11. Posterior distributions for the same data set combinations as in Fig. 5 but analyzed within a Bayesian framework.  \n<!-- image -->', 'Appendix B: Comparison with other Gaussian extrapolations': 'For the sake of assessing the robustness of our results when performing the extrapolation of the likelihood profiles into the unphysical ∑ m ν region, we compare our results with those of other works, such as Refs. [10, 21]. In these works, a similar procedure was performed but within a bayesian framework: the authors fitted to a Gaussian their posterior in the physical region and extrapolated to negative masses.  \nFIG. 12. Posterior distributions for the same data set combinations as in Fig. 7 but analyzed within a Bayesian framework.  \n<!-- image -->  \nFIG. 13. Posterior distributions for the same data set combinations as in Fig. 8 but analyzed within a Bayesian framework.  \n<!-- image -->  \nThe comparison is shown in Fig. 14, where we plot our ∆ χ 2 profiles and their -2 log P for similar datasets. We find an overall agreement in the preference for a best fit in the negative neutrino mass region, suggesting that the  \nconclusions drawn from the extrapolation are not largely dependent on the statistical procedure employed.  \nFIG. 15. Correlations of ∑ m ν with w 0 and A lens . The points are extracted from our profile likelihoods for the Planck+DESIY1+Pantheon+ combination, while the shaded regions correspond to the bayesian credible intervals at 1 and 2 σ .  \n<!-- image -->  \n<!-- image -->  \n6  \nFIG. 14. Comparison between our ∆ χ 2 profiles and their gaussian extrapolation with the corresponding extrapolated -2 log P of Ref. [10] and Ref. [21].', 'Appendix C: Correlations with other parameters': 'It has been shown that extending the standard cosmological model either with the A lens parameter or by adding a dynamical dark energy equation of state w ( a ) = w 0 +(1 -a ) w a , the cosmological bound on neutrino masses can be substantially relaxed and thus reduce the tension with oscillation data. In Fig. 15 we present the correlations of ∑ m ν with w 0 and A lens , as derived from our frequentist and Bayesian analyses by plotting the w 0 (left panel) and A lens (right panel) values of the profile likelihood points as well as the Bayesian contours. Most interestingly, as previously stated, HiLLiPoP23-PR4 yields the closest A lens value to ΛCDM.  \n0.6  \n<!-- image -->  \nw  \n<!-- image -->  \nA'}