¹

1 Introduction↩︎

Neutrino oscillation experiments have provided compelling evidence that neutrinos have non-zero masses, constituting the only experimentally established manifestation to date of physics beyond the Standard Model. Consequently, measuring the absolute neutrino mass has become a central objective in both particle physics and cosmology. On one hand, in particle physics, oscillation data have measured the two mass-squared differences among the three neutrino mass eigenstates \(m_1, m_2, m_3\), with solar neutrino measurements yielding \(\Delta m_{21}^2 \approx 7.4 \times 10^{-5}~\text{eV}^2\) and atmospheric neutrino experiments giving \(|\Delta m_{32}^2| \approx 2.4 \times 10^{-3}~\text{eV}^2\) [1]. However, oscillation experiments cannot determine the absolute neutrino mass, and the unknown sign of \(\Delta m_{32}^2\) implies two possible neutrino mass hierarchies, namely normal hierarchy (NH) with \(m_1< m_2\ll m_3\) and inverted hierarchy (IH) with \(m_3\ll m_1< m_2\). Oscillation data imply a lower bound on the total mass of \(\sum m_\nu \gtrsim 0.06\,{\rm eV}\) for NH and \(\sum m_\nu \gtrsim 0.10\,{\rm eV}\) for IH [2], [3]. Meanwhile, \(\beta\)-decay experiments provide direct laboratory measurements of the neutrino mass; the latest result from KATRIN gives \(\sum m_\nu < 1.35\,{\rm eV}\) [4]. On the other hand, cosmological observations provide a highly sensitive probe of the neutrino masses. Neutrinos behave as relativistic radiation in the early universe and contribute to the matter density at late times, while their free-streaming motion suppresses the growth of small-scale structure. These effects imprint on the cosmic microwave background (CMB), big bang nucleosynthesis (BBN), and the large-scale structure of the universe [5]–[10]. Analyses combining Planck 2018 CMB data and Sloan Digital Sky Survey baryon acoustic oscillations (BAO) data within the \(\Lambda\)-cold dark matter (\(\Lambda\)CDM) model yield an upper limit of \(\sum m_\nu < 0.12\,{\rm eV}\) at 95% confidence level [11], highlighting that cosmological observations currently provide far tighter constraints on neutrino mass than particle physics experiments.

However, cosmological measurements of neutrino mass depend sensitively on the cosmological model, especially on the nature of dark energy (DE) [12]. In our earlier systematic studies, we examined in detail the degeneracy between \(\sum m_\nu\) and the DE equation-of-state (EoS) parameter \(w\) [12]–[14]. We found that when the evolution of \(w\) corresponds to quintessence, namely \(w>-1\) at all times, or to a quintom scenario² in which \(w\) evolves from \(w>-1\) at early times to \(w<-1\) at late times, the upper limit on \(\sum m_\nu\) becomes more stringent compared with \(\Lambda\)CDM. In contrast, when the evolution of \(w\) corresponds to phantom, namely \(w<-1\) at all times, or to a quintom scenario in which \(w\) evolves from \(w<-1\) at early times to \(w>-1\) at late times, the upper limit on \(\sum m_\nu\) is relaxed compared with \(\Lambda\)CDM. These results consistently indicate that once \(w\) evolves from a larger value to a smaller value, the upper limit of \(\sum m_\nu\) becomes tighter than in \(\Lambda\)CDM, while the opposite evolution leads to a looser limit³.

Recently, the Dark Energy Spectroscopic Instrument (DESI) collaboration has published its second data release (DR2), incorporating BAO measurements from over 14 million extragalactic sources [18]. The combination of DESI DR2 BAO with CMB and type Ia supernova (SN) data indicates a statistically significant deviation (\(2.8 - 4.2\sigma\)) from the \(\Lambda\)CDM paradigm, favoring a dynamically evolving DE model characterized by an EoS that evolves from \(w < -1\) at early times to \(w > -1\) at late times [18]. These findings have stimulated numerous efforts to develop and test various explanations beyond \(\Lambda\)CDM, including various DE models [19]–[48] and alternative cosmological scenarios [49]–[88]. Furthermore, DESI reports that in the \(\Lambda\)CDM model, the cosmological upper limit on the total neutrino mass is \(\sum m_\nu < 0.064\) eV, which is close to the lower limits from oscillation experiments and in tension with the IH [18]. When considering dynamical DE models, the constraint relaxes to \(\sum m_\nu < 0.16\) eV. This is due to the strong degeneracy between the EoS of DE and neutrino mass, and the DESI results favor the evolution of DE EoS crossing from \(w<-1\) to \(w>-1\), which leads to a larger neutrino mass, consistent with our earlier important studies [12]–[14]. Meanwhile, it is worth noting that the DESI collaboration recently reported approximately a \(2\sigma\) preference for a positive total neutrino mass, \(\sum m_\nu = 0.106^{+0.050}_{-0.069}\,\mathrm{eV}\), within a DE framework mediated by the collapse of stars to cosmologically coupled black holes [89]. These results highlight the crucial role of DE dynamics in cosmological measurements of neutrino mass.

Overall, earlier studies [12]–[14] showed that when DE evolves dynamically with the EoS transitioning from \(w<-1\) at early times to \(w>-1\) at late times, the neutrino mass limit is significantly relaxed. The current DESI data provide strong evidence for such dynamical behavior of DE, naturally favoring a larger neutrino mass. Based on this, we include both key neutrino parameters, \(\sum m_\nu\) and \(N_{\rm eff}\), in our analysis, and combine CMB, DESI BAO, DESY5 SN, and weak lensing data to achieve a high-confidence measurement of \(\sum m_\nu\) at the \(2.7\sigma\) level, accompanied by a detailed and systematic discussion of the impact of dynamical DE on neutrino mass measurement.

2 Methodology and data↩︎

4pt

In cosmology, the total energy density of massive neutrinos is given by [90] \[\rho_\nu(a) = \frac{a^{-4}}{\pi^2} \int \frac{q^2\mathrm{d}q}{\mathrm{e}^{q/T_{\nu0}}+1}\sum_i \sqrt{q^2+m_{i}^2a^2}, \label{eq1}\tag{1}\] where \(q\) is the comoving momentum, \(T_{\nu0}=(4/11)^{1/3}T_\mathrm{cmb}=1.945\,\mathrm{K}\) and \(T_\mathrm{cmb}=2.725\,\mathrm{K}\) represent the present neutrino and CMB temperature, respectively. \(m_i\,(i=1,2,3)\) are the mass of each neutrino species.

In the early universe, neutrinos were relativistic, thereby \(\rho_{\nu}\) is proportional to the photon energy density (\(\rho_{\gamma}\)), given by [90] \[\rho_{\nu}(a) \to N_{\rm eff} \frac{7}{8} \left(\frac{4}{11}\right)^{4/3} \rho_\gamma(a), \label{eq2}\tag{2}\] where \(N_\mathrm{eff}\) is the effective number of relativistic neutrino species and in the Standard Model, \(N_\mathrm{eff}=3.044\) [91]–[93]. Therefore, Eq. (1 ) can be rewritten as \[\rho_\nu(a) = N_{\rm eff} \frac{7}{8} \left(\frac{4}{11}\right)^{4/3} \rho_\gamma(a) \sum_i g\left(\frac{m_i}{T_{\nu0}}a\right), \label{eq3}\tag{3}\] where \[g(y) \equiv \frac{40}{7\pi^4}\int_0^\infty \frac{x^2\sqrt{x^2+y^2}}{\mathrm{e}^x+1} \mathrm{d}x. \label{eq4}\tag{4}\] It can be demonstrated that when \(a\to\infty\) (i.e., late-time universe), Eq. (3 ) tends asymptotically to \[\rho_{\nu}(a) \to \frac{\sum m_{\nu}}{93.14h^2~\rm eV} \rho_\mathrm{crit,0}a^{-3},\] where \(h=H_0/(100~\mathrm{km~s^{-1}~Mpc^{-1}})\) is the dimensionless Hubble constant. \(\rho_\mathrm{crit,0}=3H_0^2/8\pi G\) is the current critical density.

In a flat Friedmann-Robertson-Walker universe, the dimensionless Hubble parameter can be written as \[\begin{align} E^2(a)&\equiv\frac{H^2(a)}{H_0^2}=(\Omega_\mathrm{b}+\Omega_\mathrm{c})a^{-3}+\Omega_\mathrm{de}f(a)\nonumber\\ &+\Omega_\gamma a^{-4}\biggl[1+N_{\rm eff} \frac{7}{8} \left(\frac{4}{11}\right)^{4/3} \sum_i g\left(\frac{m_i}{T_{\nu0}}a\right)\biggr], \end{align}\] where \(\Omega_{\gamma}\), \(\Omega_\mathrm{b}\), \(\Omega_\mathrm{c}\), and \(\Omega_{\rm de}\) represent the current density parameters of photon, baryon, CDM, and DE, respectively, and \(\Omega_\gamma = 2.469\times10^{-5}/h^2\) for \(T_\mathrm{cmb}=2.725\,\mathrm{K}\). Here, \(f(a)\) is the normalized \(a\)-dependent density of DE, given by \[f(a) = \exp\left( -3 \int_{0}^{\ln a} [1 + w(a')] {\rm d}\ln a' \right),\] where \(w(a)\) is the EoS of DE. In this work, we consider three cosmological models, including the \(\Lambda {\rm CDM}+\sum m_\nu+N_\mathrm{eff}\) model with \(w=-1\), \(w {\rm CDM}+\sum m_\nu+N_\mathrm{eff}\) model with a constant \(w\), and \(w_0w_a {\rm CDM}+\sum m_\nu+N_\mathrm{eff}\) model, which adopts the Chevallier-Polarski-Linder (CPL) parameterization form, \(w(a) = w_0 + w_a(1-a)\) [94], [95].

Table ¿tbl:table2? summarizes the full set of free parameters and the uniform priors adopted in the analysis. We compute the theoretical models using the CAMB code [96], [97] and use the publicly available sampler Cobaya [98] to perform Markov Chain Monte Carlo (MCMC) [99], [100] analysis. We assess the convergence of the MCMC chains using the Gelman-Rubin statistics quantity \(R - 1 < 0.02\) [101] and the MCMC chains are analyzed using the public package GetDist [102].

In our main analysis, we consider a joint analysis of the following datasets: (i) CMB: The CMB likelihoods include the temperature and polarization measurements (TT, TE, and EE spectra) from Planck 2018 [11], [103]–[105], as well as the NPIPE PR4 Planck CMB lensing reconstruction [106] and Data Release 6 of the Atacama Cosmology Telescope [107]. (ii) DESI: The BAO measurements from DESI DR2 are summarized in Table IV of Ref. [18]. (iii) PantheonPlus: The PantheonPlus comprises 1550 spectroscopically confirmed supernovae (SNe) from 18 different surveys, covering \(0.01 < z < 2.26\) [108]. (iv) DESY5: The DESY5 sample comprises 1829 photometrically classified SNe with redshifts in the range \(0.025 < z < 1.3\). (v) DESY1: The DESY1 weak lensing data are based on the analysis of 26 million source galaxies and 6.5 million lens galaxies over a 1321 \(\mathrm{deg}^2\) footprint [109], [110].

3 Results and discussions↩︎

In this section, we present the main parameter constraint results in Table 1 and Figs. 1 and 2. Figure 3 shows the relative variation of the dimensionless Hubble expansion rate.

For the \(\Lambda\mathrm{CDM}+\sum m_\nu+N_{\rm eff}\) model, we obtain an upper limit on the total neutrino mass of \(\sum m_\nu<0.084\;\mathrm{eV}\) using CMB+DESI+DESY5 data, and \(\sum m_\nu<0.085\;\mathrm{eV}\) when DESY1 is included. The inclusion of weak lensing data slightly raises the upper limit on \(\sum m_\nu\), since the lower clustering amplitude (\(S_8\)) inferred from weak lensing data allows for larger neutrino masses. Building on this, when a constant DE EoS is introduced, the \(w\mathrm{CDM}+\sum m_\nu+N_{\rm eff}\) model yields \(\sum m_\nu<0.065\;\mathrm{eV}\) and \(w=-0.955\pm0.022\) using CMB+DESI+DESY5+DESY1 data. This indicates that a constant EoS of DE with \(w>-1\) lowers the neutrino mass upper limit, as illustrated in the upper panel of Fig. 1, which is in agreement with previous studies [13], [52], [111].

When we further consider a dynamically evolving DE, the CMB+DESI+DESY5+DESY1 data yield \(\sum m_\nu = 0.098^{+0.016}_{-0.037}\,\mathrm{eV}\), indicating a \(2.7\sigma\) measurement of a neutrino mass (for further robustness checks regarding this result, see Appendix 5). This high-confidence measurement primarily arises for the following reasons:

• DESI favors a dynamical DE EoS that evolves from \(w< -1\) at early times to \(w> -1\) at late times. This evolving behavior, which leads to a larger neutrino mass, is the primary reason for the high-confidence measurement reported in our analysis, as concluded in our previous systematic study [12]–[14], [52], [112]. Meanwhile, the detection of a non-zero neutrino mass is strongly influenced by evidence for the dynamical evolution of DE. As shown in Fig. 2, \(\sum m_\nu\) is anti-correlated with \(w_a\) and positively correlated with \(w_0\). The CMB+DESI+DESY5+DESY1 data give \(w_0=-0.745^{+0.056}_{-0.063}\) and \(w_a = -0.91^{+0.26}_{-0.22}\), corresponding to evidence for a dynamically evolving DE at the \(\sim 4\sigma\) level, thereby favoring a larger neutrino mass. In contrast, replacing DESY5 with PantheonPlus yields CMB+DESI+PantheonPlus+DESY1 results of \(\sum m_\nu < 0.144\,\mathrm{eV}\), where the evidence for DE dynamics drops to \(\sim2\sigma\) (\(w_0 = -0.845 \pm 0.058\), \(w_a = -0.58^{+0.24}_{-0.21}\)), leading only to an upper limit on the neutrino mass. We further test this conclusion by introducing priors with varying evidences for dynamically evolving DE; see Appendix 6 for a detailed discussion.

• Treating \(N_\mathrm{eff}\) as a free parameter also impacts the neutrino mass measurement. In the \(w_0w_a\mathrm{CDM}+\sum m_\nu+N_{\rm eff}\) model using CMB+DESI+DESY5+DESY1 data, we obtain \(N_{\mathrm{eff}} = 2.46^{+0.60}_{-0.24}\). Although this value is consistent with the particle-physics standard value \(3.044\) at the \(1\sigma\) level, its central value lies below the standard value. One possible physical explanation for a reduced \(N_{\mathrm{eff}}\) is a low reheating temperature in the early universe⁴, which prevents active neutrinos from fully thermalizing to the standard Fermi-Dirac distribution [115], [116]. Some studies have shown that in such scenarios the late-time constraint on \(\sum m_\nu\) can be substantially relaxed [113]. This follows because a reduced neutrino number density, or a colder neutrino population (leading to \(N_{\rm eff}<3.044\)), requires a larger total neutrino mass to produce the same late-time effects on background expansion and structure formation.

• The inclusion of weak lensing data also allows for a larger neutrino mass. Specifically, \(\sum m_\nu\) is anti-correlated with \(S_8\), and the lower \(S_8\) values preferred by the DESY1 weak lensing data also push the neutrino mass higher. Without weak lensing, CMB+DESI+DESY5 only impose an upper limit of \(\sum m_\nu < 0.141\,\mathrm{eV}\), with a multi-peaked one-dimensional (1D) marginalized posterior distribution, as shown in lower panel of Fig. 1. The inclusion of DESY1 suppresses the lower-mass peak of the posterior, thereby providing evidence for a non-zero neutrino mass.

Notably, the DESI collaboration recently found that, under a hypothetical unphysical negative neutrino-mass scenario (i.e., one whose effects on the background expansion and structure formation are opposite to those of a positive-mass neutrino), nearly all cosmological observations preferentially indicate a negative neutrino mass, producing a tension of up to \(3\sigma\) with neutrino-oscillation results. Dynamical DE can effectively relieve this tension [117]. The DESI analysis carried out a systematic investigation of neutrino cosmology using a variety of data combinations (DESI DR1 full-shape, DESI DR2 BAO, CMB, and several SN compilations), a range of DE models, and different neutrino-physics assumptions. In contrast, our analysis concentrates on the case of dynamical DE while treating both \(\sum m_{\nu}\) and \(N_{\mathrm{eff}}\) as free parameters, and additionally includes weak gravitational lensing data to further broaden the neutrino mass, enabling a cosmological measurement of a non-zero, positive neutrino mass. As discussed above, we obtain a detection of non-zero neutrino mass at up to \(2.7\sigma\), which is broadly consistent with neutrino-oscillation experiments and demonstrates the considerable potential of cosmological observations to measure the neutrino mass. While the current analysis is phenomenologically driven, these robust observational preferences underscore the imperative for future theoretical works to construct explicit beyond-\(\Lambda\)CDM models. Investigating scenarios such as those involving interactions between dark matter and neutrinos [118], [119], will be crucial to elucidate the underlying physical mechanisms driving these data anomalies.

Next, we further illustrate the impact of DE on neutrino mass measurements by examining the evolution of the late-time relative expansion rate. The CMB TT power spectrum data tightly constrain the position and amplitude of the first acoustic peak. The determination of the first peak depends precisely on the angular size of the sound horizon at decoupling, \(\Theta_{\rm s}\equiv r_{\rm s}/D_\mathrm{A}(z_\mathrm{LS})\), and the redshift of matter-radiation equality, \(z_{\rm eq}\). Here, the sound horizon at decoupling \(r_{\rm s}\) is essentially fixed by the physics of pre-recombination periods, while DE variations are only relevant at late times; therefore, \(r_{\rm s}\) is unaffected by DE [111]. To keep both the angular diameter distance of last scattering \(D_\mathrm{A}(z_\mathrm{LS})\) and \(z_{\rm eq}\) approximately fixed, any change in DE must be compensated by shifts in other cosmological parameters. \(D_\mathrm{A}(z_\mathrm{LS})\) is given by \[\label{eq8} D_\mathrm{A}(z_\mathrm{LS}) = \frac{c}{H_0(1+z_\mathrm{LS})} \int_{0}^{z_\mathrm{LS}} \frac{1}{E(z)} \,\mathrm{d}z,\tag{5}\] where \(z_{\rm LS}\) is the redshift of last scattering. We then define \[\mathcal{E}(z) = \left. \frac{E(z)_{\Lambda\mathrm{CDM}}}{E(z)_{w\mathrm{CDM}/ w_0w_a\mathrm{CDM}}} \right|_{\Omega_{\rm b},\,\Omega_{\rm c},\,\Omega_{\rm de},\,\sum m_\nu,\,N_{\rm eff}} \;-\;1,\] where the notation \(\left.\right|_{\Omega_{\rm b},\,\Omega_{\rm c},\,\Omega_{\rm de},\,\sum m_\nu,\,N_{\rm eff}}\) indicates that these parameters are held fixed to the best-fit values obtained in the \(\Lambda\mathrm{CDM}+\sum m_\nu+N_{\mathrm{eff}}\) model from the \(\mathrm{CMB}+\mathrm{DESI}+\mathrm{DESY5}+\mathrm{DESY1}\) data, with \(\Omega_{\rm b} = 0.048\), \(\Omega_{\rm c} = 0.253\), \(\Omega_{\rm de} = 0.699\), \(\sum m_\nu = 0.085\) eV (the \(2\sigma\) upper limit), and \(N_\mathrm{eff} = 3.02\).

The redshift evolution of \(\mathcal{E}(z)\) is shown in Fig. 3. It is clear that when the curve lies above zero, DE variations increase the integral in Eq. (5 ), and when it lies below zero, they decrease it. In the \(w\)CDM model, \(w>-1\) yields a relative expansion rate higher than that in \(\Lambda\)CDM, resulting in \(\mathcal{E}(z)<0\) (blue solid line). To keep \(D_{\rm A}(z_{\rm LS})\) unchanged, directly reducing \(\Omega_{\rm b}\), \(\Omega_{\rm c}\), or \(N_{\rm eff}\) is undesirable because each shifts \(z_{\rm eq}\) and would introduce unwanted changes elsewhere in the CMB power spectrum. Instead, lowering \(\sum m_\nu\), whose effect on \(z_{\rm eq}\) is minimal, is the most economical option for compensating the reduced integral. Consequently, introducing a constant EoS of DE with \(w>-1\) tightens the upper limit on \(\sum m_\nu\).

For the \(w_0w_a\)CDM model, current data strongly prefer a crossing DE EoS that evolves from \(w<-1\) at early times to \(w>-1\) at late times. As shown in Fig. 3, this DE dynamics leads to a significant alteration in the late-time relative expansion rate of the universe: it exceeds that of the \(\Lambda\)CDM model for \(0<z\lesssim1\), whereas it falls below it for \(z\gtrsim1\). It should be emphasized that CMB observations constrain the position of the first acoustic peak with extremely high precision, which requires the angular diameter distance to the last scattering surface, \(D_{\rm A}(z_{\rm LS})\), to remain highly stable. Since the integration range for \(1/E(z)\) in Eq. (5 ) extends to \(z_{\rm LS}\approx1090\), the negative deviation in the \(z\gtrsim1\) region dominates the entire geometric integral, leading to a substantial overall increase in the integral value compared to the standard model.

To maintain \(D_{\rm A}(z_{\rm LS})\) constant to fit the CMB data, the model must introduce a corresponding physical compensation effect. If one directly adjusts \(\Omega_{\rm m}\) or \(H_0\) to compensate for the integral change, it would inevitably alter the redshift of matter-radiation equality, \(z_{\rm eq}\), thereby spoiling the spectral structure of the CMB power spectrum. Therefore, increasing the total neutrino mass \(\sum m_\nu\), which has a minimal impact on \(z_{\rm eq}\), becomes the most natural choice to offset this geometric effect of DE dynamics. This is highly consistent with the parameter degeneracy mechanisms revealed in our previous series of studies [12]–[14], [52], [112].

Moreover, such crossing DE dynamics (the joint contribution of \(w_0\) and \(w_a\)) tends to yield a lower \(H_0\) in the inference, which not only suppresses \(N_{\rm eff}\) but also slightly alters the theoretical expectation of \(z_{\rm eq}\). Thus, the dynamically evolving DE in the \(w_0w_a\)CDM model intrinsically drives an increase in \(\sum m_\nu\) and a decrease in \(N_{\rm eff}\), a physical picture that corroborates the latest conclusions of the DESI collaboration [18], [117], [120]. Finally, as the evidence supporting dynamically evolving DE strengthens (i.e., the transition from the pink dashed line to the orange dash-dot line in Fig. 3, reflecting a larger deviation of \(w_0, w_a\) from \(-1, 0\)), the evolutionary deviations in the relative expansion rate become more drastic. This demands a larger \(\sum m_\nu\) and a smaller \(N_{\rm eff}\) to absorb the late-time geometric effects brought about by the variation in the DE EoS, thereby providing high-confidence observational evidence for a non-zero neutrino mass at the physical level.

Finally, it is worth noting that Ref. [64] also reported a \(2\sigma\) measurement of \(\sum m_\nu = 0.190 \pm 0.088\,\mathrm{eV}\) in an extended cosmological model with 12 parameters, including dynamical DE, the running of the scalar spectral index (\(\alpha_s\)), and the scaling of the lensing amplitude (\(A_{\rm lens}\)). Compared with this study, we only consider a dynamical DE scenario including two important neutrino parameters (\(\sum m_\nu\) and \(N_{\rm eff}\)), motivated by the key conclusions of our earlier series of works [12]–[14], [52], [112] and DESI’s preference for dynamical DE with its EoS evolving from \(w<-1\) to \(w>-1\). As anticipated, relative to analyses incorporating other extensions, we achieve a higher statistic significance for the neutrino mass measurement. Furthermore, our result lies below but close to the IH bound and therefore does not support the IH, which is consistent with the conclusions of particle physics experiments [121], [122].

4 Conclusion↩︎

In this paper, we investigate the cosmological implications of simultaneously including \(\sum m_\nu\) and \(N_{\rm eff}\) within the framework of dynamical DE. We further clarify the specific impact of dynamical DE on neutrino mass measurements from the perspectives of parameter degeneracies and the relative expansion rate.

Our joint analysis of DESI DR2, CMB, DESY5, and DESY1 data yields \(\sum m_\nu = 0.098^{+0.016}_{-0.037}\,\mathrm{eV}\), indicating a measurement of a non-zero neutrino mass at the \(2.7\sigma\) level within the \(w_0w_a\mathrm{CDM}+\sum m_\nu+N_{\rm eff}\) model. This high-confidence measurement primarily arises from these factors: (i) The current DESI results favor a dynamical DE EoS evolving from \(w<-1\) to \(w>-1\), which tends to increase the neutrino mass. When DESY5 is replaced with PantheonPlus, the evidence for a positive neutrino mass disappears, as the evidence for dynamical DE becomes weaker. Furthermore, only upper limits are obtained in the \(\Lambda\)CDM and \(w\)CDM frameworks. The overall outcomes highlight the conclusion from our previous systematic studies [12], [13]. (ii) Allowing \(N_{\mathrm{eff}}\) to vary (with a best-fit value below the standard value 3.044) also requires a larger \(\sum m_\nu\) to compensate. (iii) The inclusion of weak lensing data, which prefers a relatively lower \(S_8\), also suppresses the possibility of a lower neutrino mass. Finally, we test the robustness of our results by investigating how DE affects the detection of neutrino mass and alters the evidence for dynamical DE, through the evolution of the late-time relative expansion rate. We find that as the evidence for dynamical DE strengthens, the deviations in the relative expansion rate become more pronounced, requiring a larger \(\sum m_\nu\) to compensate for the DE changes. This serves as the most economical way to offset the integral variation and consequently enhances the evidence for a non-zero neutrino mass.

Our results indicate that this represents a highly compelling cosmological scenario, and that the DESI dataset has inaugurated a new era for measuring neutrino mass through cosmological observations. With the forthcoming, more complete DESI data releases and the inclusion of full-shape modeling of the power spectrum, stronger preferences for dynamical DE are likely to emerge, thereby offering promising prospects for further revealing cosmological evidence of a non-zero neutrino mass.

Acknowledgments↩︎

We thank Sheng-Han Zhou, Jia-Le Ling, and Yi-Min Zhang for their helpful discussions. This work was supported by the National Natural Science Foundation of China (Grants Nos. 12533001, 12575049, and 12473001), the National SKA Program of China (Grants Nos. 2022SKA0110200 and 2022SKA0110203), the China Manned Space Program (Grant No. CMS-CSST-2025-A02), and the National 111 Project (Grant No. B16009).

Conflict of Interest↩︎

The authors declare that they have no conflict of interest.

5 Robustness check of the positive neutrino mass↩︎

In this Appendix, we comprehensively test the robustness of our preference for a positive neutrino mass by examining the stability of the results against prior volume effects, boundary truncations, and dataset variations. In parameter inference, traditional MCMC methods can sometimes be subject to Bayesian prior volume effects, which may cause the peak of the posterior distribution to deviate from the true best-fit point. To rigorously examine the statistical robustness of our \(2.7\sigma\) preference for a positive neutrino mass, we employ the profile likelihood (PL) method, which is independent of priors. We fix \(\sum m_\nu\) at a series of discrete values within the range of \([0, 0.15]~\mathrm{eV}\), and at each fixed point, we compute the global maximum likelihood for the remaining cosmological parameters. As illustrated in Fig. 4, the minimum \(\chi^2\) occurs at \(\sum m_\nu = 0.082~\mathrm{eV}\), which is remarkably close to the central value of \(0.098~\mathrm{eV}\) given by the MCMC analysis. More importantly, the \(\Delta \chi^2\) curve rises significantly towards the low-mass end, reaching \(\Delta \chi^2 \approx 6\) as \(\sum m_\nu \to 0\). Therefore, the current joint datasets possess a genuine physical preference for a positive neutrino mass, rather than being a statistical artifact.

To examine the robustness of our results with respect to the choice of datasets, we conducted a comparative analysis by replacing DESY1 with DESY3. Note that the DESY3 weak lensing data are based on the analysis of approximately 100 million source galaxies and 10 million lens galaxies over a 4143 \(\mathrm{deg}^2\) footprint, providing updated cosmic shear and galaxy clustering observations [123], [124]. Using CMB+DESI+DESY5+DESY3, we obtain \(\sum m_\nu = 0.081^{+0.028}_{-0.071}~\mathrm{eV}\). As illustrated in the upper panel of Fig. 5, the peak of the 1D marginalized posterior distribution of the neutrino mass shifts to the left, and the low-mass peak reappears. This is fully consistent and internally coherent with the conclusions presented in our manuscript, namely that the non-zero neutrino mass is primarily driven by the data preference for dynamical DE. Compared to DESY1, although the joint data with DESY3 also yields a relatively lower \(S_8 = 0.8179\pm 0.0079\), it provides weaker evidence for the dynamical evolution of DE (\(w_0 = -0.765\pm 0.059\), \(w_a = -0.75^{+0.27}_{-0.24}\)), as shown in the lower panel of Fig. 5. Based on the relative expansion rate compensation mechanism discussed in detail earlier, the weakened evidence for the dynamical evolution of DE naturally leads to a corresponding decrease in the required \(\sum m_\nu\) compensation. This robustness test not only confirms the sensitivity of the results to the datasets but also inversely verifies the degeneracy relationship between dynamical DE and the neutrino mass.

It is worth noting that some recent analyses have allowed for an unphysical negative neutrino mass to examine statistical behaviors [117], [125], [126]. Under the standard \(\Lambda\)CDM model, since the peak of the posterior distribution lies in the negative half-axis and deviates from the \(\sum m_{\nu}=0\) boundary, such an extension helps reveal potential data inconsistencies or tensions. However, under the \(w_{0}w_{a}\) CDM model considered in this work, the peak of the posterior distribution is already on the positive half-axis and deviates from the zero boundary (as demonstrated by our PL and MCMC results, where the peak is located between 0.082 and 0.098 eV). Since the posterior distribution itself is no longer primarily affected by the truncation boundary, extending the parameter space into the unphysical negative region would not yield any substantial changes to the physical inference in the positive domain.

6 Impact of dynamical DE evidence on neutrino mass↩︎

In this Appendix, we report and discuss the impact of the evidence for dynamical DE on neutrino mass constraints. We adopt the constraints on the \(w_{0}w_{a}\mathrm{CDM} + \sum m_{\nu} + N_{\rm eff}\) model derived from CMB+DESI+DESY5 data as our baseline. We then impose separately a \(2\sigma\) prior (\(w_{0}=-0.86\pm0.07\), \(w_{a}=-0.50\pm0.25\)) and a \(5\sigma\) prior (\(w_{0}=-0.65\pm0.07\), \(w_{a}=-1.25\pm0.25\)). The resulting 1D marginalized posterior distributions of \(\sum m_{\nu}\) are shown in Fig. 6. As expected, for the bimodal structure of the CMB+DESI+DESY5 posterior, the \(2\sigma\) prior suppresses the high-mass peak, yielding an upper limit of \(\sum m_{\nu} < 0.118\;\mathrm{eV}\). Conversely, the \(5\sigma\) prior suppresses the low-mass peak, leading to \(\sum m_{\nu} = 0.083^{+0.046}_{-0.039}\;\mathrm{eV}\). This further confirms the significant impact of evidence for dynamical DE on neutrino mass measurements and ensures the robustness of our conclusions.

References↩︎