Quantifying the tension between cosmological and terrestrial constraints on neutrino masses

TL;DR

This paper assesses the potential tension between cosmological and terrestrial constraints on neutrino masses, emphasizing the importance of neutrino mass ordering and the precision needed in future measurements to distinguish scenarios.

Contribution

It introduces statistical methods to quantify tensions between cosmological and terrestrial neutrino mass limits, highlighting the impact of neutrino mass ordering assumptions.

Findings

Tension depends on whether normal or inverted ordering is assumed.

Rejecting inverted ordering requires future cosmological measurements with <0.02 eV accuracy.

Current and forecasted data show potential for distinguishing mass orderings with improved precision.

Abstract

The sensitivity of cosmology to the total neutrino mass scale $\Sigma m_\nu$ is approaching the minimal values required by oscillation data. We study quantitatively possible tensions between current and forecasted cosmological and terrestrial neutrino mass limits by applying suitable statistical tests such as Bayesian suspiciousness, parameter goodness-of-fit tests, or a parameter difference test. In particular, the tension will depend on whether the normal or the inverted neutrino mass ordering is assumed. We argue, that it makes sense to reject inverted ordering from the cosmology/oscillation comparison only if data are consistent with normal ordering. Our results indicate that, in order to reject inverted ordering with this argument, an accuracy on the sum of neutrino masses $\sigma ({m_\nu})$ of better than 0.02~eV would be required from future cosmological observations.

Tables

Table 1. Table 1: Tension between cosmological and terrestrial neutrino mass determination assuming different cosmological data sets: current data, a future observation with a value for ∑ m ν subscript 𝑚 𝜈 \sum m_{\nu}  consistent with NO (future NO, Eq. ( 20 )), and a non-observation of ∑ m ν subscript 𝑚 𝜈 \sum m_{\nu}  (future 0, Eq. ( 21 )). In each case we show results for the two parameterizations 3 ​ ν 3 𝜈 3\nu and Σ Σ \Sigma (see section  III.2 ) and the two mass orderings. The table shows the test statistics ln ⁡ S 𝑆 \ln S , Q DMAP subscript 𝑄 DMAP Q_{\rm DMAP} from Eqs. ( 11 ) and ( 14 ) respectively, and the probabilities of the data sets being consistent, p S subscript 𝑝 𝑆 p_{S} , 1 − P ​ ( Q DMAP ) 1 𝑃 subscript 𝑄 DMAP 1-P(Q_{\rm DMAP}) and 1 − P ​ ( Δ ) 1 𝑃 Δ 1-P(\Delta) , corresponding to the suspiciousness test, the parameter goodness-of-fit test, and the parameter shift test, respectively. In the right part of the table we show the Bayesian model dimensionalities according to Eqs. ( 12 ) and ( 13 ), indicating with A 𝐴 A the terrestrial and B 𝐵 B the cosmological data sets. The values for p S subscript 𝑝 𝑆 p_{S} [as well as for 1 − P ​ ( Q DMAP ) 1 𝑃 subscript 𝑄 DMAP 1-P(Q_{\rm DMAP}) ] are calculated with the parameter counting according to Eq. ( 15 ), i.e., for 3 (1) dof for the 3 ​ ν 3 𝜈 3\nu ( Σ Σ \Sigma ) parameterization.

Cosmo data	Model	$\ln S$	$p_S$	$Q_{DMAP}$	$1-P\left(Q_{DMAP}\right)$	$1-P\left(\Delta \right)$	$d_A$	$d_B$	$d_{AB}$	$d$
current	$3\nu$ NO	-1.34	0.13 ($1.52\sigma$)	2.62	0.45 ($0.75\sigma$)	0.011 ($2.66\sigma$)	1.94	0.31	2.28	-0.02
	$3\nu$ IO	-3.30	0.023 ($2.29\sigma$)	5.61	0.13 ($1.50\sigma$)	0.0014 ($3.90\sigma$)	1.98	0.33	3.54	-1.23
	$\Sigma$ NO	-1.23	0.063 ($1.86\sigma$)	2.62	0.098 ($1.65\sigma$)	0.021 ($2.31\sigma$)	0.44	1.24	1.47	0.21
	$\Sigma$ IO	-2.86	0.0095 ($2.59\sigma$)	5.59	0.017 ($2.39\sigma$)	0.0014 ($3.20\sigma$)	0.46	1.21	1.65	0.02
future NO	$3\nu$ NO	0.55	0.59 ($0.53\sigma$)	0.033	1 ($0.00\sigma$)	0.23 ($1.19\sigma$)	1.94	1.81	2.32	1.43
	$3\nu$ IO	-1.75	0.089 ($1.70\sigma$)	3.99	0.26 ($1.12\sigma$)	0.018 ($2.54\sigma$)	1.98	1.71	2.76	0.93
	$\Sigma$ NO	0.2	0.44 ($0.78\sigma$)	0.035	0.84 ($0.21\sigma$)	0.1 ($1.62\sigma$)	0.44	0.93	1.00	0.38
	$\Sigma$ IO	-2.16	0.021 ($2.31\sigma$)	3.98	0.043 ($2.03\sigma$)	0.00099 ($3.29\sigma$)	0.46	0.92	1.60	-0.22
future 0	$3\nu$ NO	-4.58	0.0068 ($2.70\sigma$)	8.78	0.032 ($2.14\sigma$)	0.0016 ($3.66\sigma$)	1.94	0.31	2.55	-0.30
	$3\nu$ IO	-13.04	2.3e-06 ($4.74\sigma$)	24.90	1.6e-05 ($4.31\sigma$)	2.5e-05 ($5.28\sigma$)	1.98	0.32	3.51	-1.21
	$\Sigma$ NO	-4.56	0.0015 ($3.18\sigma$)	8.80	0.0028 ($2.99\sigma$)	8.1e-06 ($4.46\sigma$)	0.44	1.00	1.71	-0.26
	$\Sigma$ IO	-12.68	2.8e-07 ($5.13\sigma$)	24.91	5.4e-07 ($5.01\sigma$)	4.1e-10 ($6.25\sigma$)	0.46	1.03	1.87	-0.37

Equations

$${\Omega }_{\nu }{h}^2=\sum _i{m}_i/(93.12\text{ eV}),$$

$$\sum m_{\nu }\equiv \sum _{i=1}^3{m}_i<0.12\mathrm{eV}(95\%\mathrm{CL}).$$

$$\sum m_{\nu }<0.09\mathrm{eV}(95\%\mathrm{CL}),$$

$$\begin{split}\Delta m^{2}_{21}&=(7.50\pm 0.21)\times 10^{-5}\,{\rm eV}^{2}\,,\\ |\Delta m^{2}_{31}|&=\left\{\begin{array}[]{ll}(2.550\pm 0.025)\times 10^{-3}\,{\rm eV}^{2}&\quad\rm(NO)\\ (2.450\pm 0.025)\times 10^{-3}\,{\rm eV}^{2}&\quad\rm(IO)\end{array}\right.\,,\end{split}$$

$$\begin{split}\sum m_{\nu}>\left\{\begin{array}[]{ll}(0.0591\pm 0.00027)\,{\rm eV}&\quad\rm(NO)\\ (0.0997\pm 0.00051)\,{\rm eV}&\quad\rm(IO)\\ \end{array}\right.\end{split}\,.$$

$${\mathcal{Z}}_D=P(D|M)=\int \text{d}\theta {\mathcal{L}}_D(\theta )\Pi (\theta ),$$

$${\mathcal{P}}_D(\theta )=P(\theta |D,M)=\frac{{\mathcal{L}}_D(\theta )\Pi (\theta )}{{\mathcal{Z}}_D},$$

$$R\equiv \frac{{\mathcal{Z}}_{AB}}{{\mathcal{Z}}_A{\mathcal{Z}}_B}.$$

$$\ln I={\mathcal{D}}_A+{\mathcal{D}}_B-{\mathcal{D}}_{AB},$$

$${\mathcal{D}}_D=\int \text{d}\theta {\mathcal{P}}_D\ln \left(\frac{{\mathcal{P}}_D}{\Pi }\right).$$

$$\ln S\equiv \ln R-\ln I.$$

$$d_D=2\int \text{d}\theta {\mathcal{P}}_D{\left[\ln \left(\frac{{\mathcal{P}}_D}{\Pi }\right)-{\mathcal{D}}_D\right]}^2.$$

$$d=d_A+d_B-d_{AB}.$$

$$\begin{array}{rl}Q\equiv & -2\ln {\mathcal{L}}_{AB}({\hat{\theta }}_{AB})\\ & +2\ln {\mathcal{L}}_A({\hat{\theta }}_A)+2\ln {\mathcal{L}}_B({\hat{\theta }}_B).\end{array}$$

$$n=p_A+p_B-p_{AB},$$

$${\mathcal{P}}_{\Delta }(\Delta \theta )=\int {\mathcal{P}}_A(\theta ){\mathcal{P}}_B(\theta -\Delta \theta )d\theta ,$$

$$\Delta ={\int }_{{\mathcal{P}}_{\Delta }(\Delta \theta )>{\mathcal{P}}_{\Delta }(0)}{\mathcal{P}}_{\Delta }(\Delta \theta )d\Delta \theta ,$$

$$m_{\beta }^2=(0.06\pm 0.32){\text{eV}}^2.$$

$$m_{\beta }^2\approx {\left(\frac{\sum m_{\nu }}{3}\right)}^2.$$

$$\sum m_{\nu }=0.06\pm 0.02\text{eV (future NO)},$$

$$\sum m_{\nu }=0\pm 0.02\text{eV (future 0)}.$$

$$B_{\mathrm{NO},\mathrm{IO}}\equiv \frac{{\mathcal{Z}}_{\mathrm{NO}}}{{\mathcal{Z}}_{\mathrm{IO}}}.$$

Full text

I Introduction

Massive neutrinos affect cosmological observables due to the unique behaviour as dark radiation in early times and as dark matter at late times, see Lesgourgues and Pastor (2006) for a review. Effectively, cosmology is sensitive to the total energy density of relic neutrinos. In the standard scenario, after they become non-relativistic, the correspondence between non-relativistic neutrino energy density and neutrino masses is approximately given by Froustey (2021)

[TABLE]

where $h$ is the reduced Hubble parameter. Constraining the non-relativistic neutrino energy density, therefore, allows us to obtain bounds on the sum of the neutrino masses under the assumption that the above equation holds. Combining information from CMB and BAO observations, the Planck collaboration obtains Aghanim et al. (2020)

[TABLE]

Adding more recent data, even more stringent limits can be obtained. For instance, Ref. Di Valentino et al. (2021) finds

[TABLE]

see also Palanque-Delabrouille et al. (2020); di Valentino et al. (2022). In the near future we expect that the DESI Aghamousa et al. (2016) and/or Euclid Amendola et al. (2018) surveys may provide sensitivities to $\sum m_{\nu}$ of down to $0.02\,\mathrm{eV}$ or beyond, see e.g. Font-Ribera et al. (2014); Basse et al. (2014); Hamann et al. (2012); Carbone et al. (2011); Brinckmann et al. (2019).

On the other hand, neutrino oscillation experiments provide accurate determinations of the two neutrino mass-squared splittings de Salas et al. (2021), see also Refs. Esteban et al. (2020); Capozzi et al. (2021):

[TABLE]

where $\Delta m^{2}_{ij}\equiv m_{i}^{2}-m_{j}^{2}$. The sign of $\Delta m^{2}_{31}$ determines the type of neutrino mass ordering, being positive for normal ordering (NO) and negative for inverted ordering (IO). With the mass-splittings determined, oscillation data provide a lower bound on the sum of the neutrino masses, obtained by assuming that the lightest neutrino mass is zero. From Eq. (4), one finds

[TABLE]

Comparing these results with Eqs. (2) and (3), it is possible to notice that cosmological upper bounds are already comparable to the lower bound for IO, and near future sensitivities will probe the NO region case as well. This may happen in two ways: by measuring a value equal or slightly larger than $\sum m_{\nu}\simeq 0.06$ eV and confirming that at least two neutrinos have a positive mass in agreement with oscillation results, or by strengthening the current upper bounds to the level that both minimal values in Eq. (5) will be disfavored by cosmology. Therefore, it is mandatory to quantify a possible tension between cosmology and oscillation data, which constitutes the main goal of this manuscript. Such a tension would have important implications: the absence of a detection of a finite neutrino mass in cosmology as predicted by Eq. (5) could be a striking signal for a non-standard cosmological model beyond the vanilla $\Lambda$CDM model and/or non-standard neutrino properties, see e.g. Alvey et al. (2022) for a discussion.

Furthermore, the tension between the lower bound on $\sum m_{\nu}$ for IO and the bound from cosmology could be used in principle to disfavour IO compared to NO, see e.g. Hannestad and Schwetz (2016); Gerbino et al. (2017); Vagnozzi et al. (2017); Simpson et al. (2017); Gariazzo and Mena (2019); De Salas et al. (2018); Heavens and Sellentin (2018); Gariazzo et al. (2018); Roy Choudhury and Hannestad (2020); Mahony et al. (2020); Hergt et al. (2021); Jimenez et al. (2022); Gariazzo et al. (2022) for an incomplete list of studies on this topic. Typically, in these papers some kind of Bayesian model comparison between NO and IO is performed, leading to posterior odds for the two models.

In the following we will address this question with a slightly different approach, namely by quantifying the tension between cosmology and oscillation data for the two orderings. We argue that it is meaningful to reject IO from a comparison of cosmology and oscillations only if these two data sets are consistent for NO. In the case when there is tension between cosmology and oscillations for both orderings, a relative comparison of the two models can be misleading. We shall explore this putative tension exploiting both current and future cosmological measurements.

The structure of the manuscript is as follows. In Sec. II we describe the different methods commonly exploited in the literature to quantify tensions between two sets of measurements. Section III contains a description of the methodology for the numerical analyses, the parameterizations employed to describe the parameter space and the data involved in quantifying the tension between cosmological and terrestrial neutrino mass measurements. Section IV presents the results from our analyses, including a mass ordering comparison. Finally, we conclude in Sec. V.

II Tension metrics

In this section we provide a brief review of various metrics used to quantify a tension between different data sets. We follow closely the discussion in Ref. Lemos et al. (2021), where a number of tests is reviewed and applied in the context of the $H_{0}$ tension. We refer the interested reader to Ref. Lemos et al. (2021) for further references and more in depth discussions of the various tests. Additional discussions can be found for instance, in the context of cosmology in Lin and Ishak (2021); Raveri and Hu (2019), in the context of Type Ia Supernova analysis in Amendola et al. (2013), and within a frequentist framework in the context of neutrino oscillations in Maltoni and Schwetz (2003).

To fix the notation, in the following $\mathcal{L}_{D}=P(D|\theta,M)$ denotes the likelihood, which is the probability for the data $D$ given a model $M$ with parameters $\theta$, $\Pi=P(\theta|M)$ is the prior for the parameters,

[TABLE]

is the Bayesian evidence, and

[TABLE]

is the posterior density for the parameters $\theta$ for data $D$. Considering now two data sets $D=A,B$, the question posed here is whether these data are consistent within a given model. In order to quantitatively address this question, the following tests can be used:

• •

Bayesian evidence ratio. Consider the ratio

[TABLE]

The numerator corresponds to the evidence when data sets $A$ and $B$ are described by the same set of parameters $\theta$, whereas in the denominator different parameters may be preferred by the two data sets. Values of $R\gg 1\,(\ll 1)$ would indicate agreement (disagreement) between the two data sets. As discussed in Lemos et al. (2021), $R$ is dependent on the prior volume, and small values of $R$, indicating a possible tension between data sets, can be increased by increasing the prior volume. Therefore, we will not use the Bayesian evidence ratio in our tension analysis below.

• •

Bayesian suspiciousness. This test departs from the Bayesian evidence ratio, but the information ratio $I$ based on the Kullback-Leibler divergence is used to remove the prior dependence. Consider the log-information ratio

[TABLE]

where the Kullback-Leibler divergence is defined as

[TABLE]

Using the log-information ratio we can cancel the prior dependence from the Bayesian evidence ratio $R$ and define the suspiciousness Handley and Lemos (2019a):

[TABLE]

As for $R$, positive values of $\ln S$ indicate agreement among the data sets while negative ones indicate disagreement.

For Gaussian posteriors, the quantity $(d-2\ln S)$ follows a $\chi_{d}^{2}$ distribution, where the number of degrees of freedom can be obtained using the Bayesian model dimensionality defined in Handley and Lemos (2019b):

[TABLE]

In order to compute the significance of the tension between two data sets, the relevant Bayesian dimensionality can be obtained using Handley and Lemos (2019b)

[TABLE]

As we will discuss in the following, the Bayesian model dimensionality $d_{D}$ may have problems when dealing with posteriors that are not Gaussian in the parameters under consideration, or when the prior limits impose a significant cut on the posterior shape. In these cases, we will replace the Bayesian model dimensionality with a more naive counting for the number of degrees of freedom, see below.

• •

Parameter goodness-of-fit tests. This test is based on the idea to evaluate the “cost” of explaining data sets together (i.e., with the same parameter values) as compared to describing them separately (i.e., each data set can chose its own preferred parameter values). Therefore this type of tests is sometimes also called “goodness-of-fit loss” tests. We take as an example two data sets $A$ and $B$, as of interest in this study (generalization to more data sets is straight-forward). Compatibility of the data sets is evaluated using the test statistic

[TABLE]

Here $\hat{\theta}_{D}$ denotes the parameter values which “best” describe data set $D$. In the context of frequentist statistics, $\hat{\theta}_{D}$ is taken as the parameter values maximizing the likelihood $\mathcal{L}_{D}(\theta)$ Maltoni and Schwetz (2003). In this case, the test statistic is denoted as $Q\equiv\chi^{2}_{\rm PG}$ and, by construction, $\chi^{2}_{\rm PG}$ is independent of the prior and any re-parameterization (as long as the number of independent parameters remains the same). In the context of Bayesian analysis, $\hat{\theta}_{D}$ is taken at the parameter values at the “maximum a posteriori” (MAP, the point at which the posterior assumes its maximum value), which in general does depend on the prior choice, see Refs. Lemos et al. (2021); Raveri and Hu (2019), where the corresponding test statistic is denoted by $Q_{\rm DMAP}$ (difference of log-likelihoods at their MAP point). For flat uninformative priors ($\Pi(\theta)=const$) maximum likelihood and maximum posterior points are identical and $\chi^{2}_{\rm PG}=Q_{\rm DMAP}$.

Under certain regularity conditions, $Q$ from Eq. (14) is distributed as a $\chi^{2}_{n}$ distribution, where $n$ is the number of parameters in common to both data sets $A$ and $B$,

[TABLE]

where $p_{D}$ denotes the parameters of data set $D$, see Maltoni and Schwetz (2003); Raveri and Hu (2019) for precise definitions.

• •

Parameter differences. This test measures the distance between posterior distributions for the parameters $\theta$, given two different datasets Raveri and Doux (2021); Raveri et al. (2020). Let us define the difference $\Delta\theta=\theta_{1}-\theta_{2}$, where $\theta_{1}$ and $\theta_{2}$ are two points in the shared parameter space. Assuming, as in our case, that the datasets $A$ and $B$ are independent, the posterior distribution for $\Delta\theta$ can be computed using:

[TABLE]

where the integral runs over the entire parameter space of the shared parameters. The probability that there is a parameter shift between the two posteriors is quantified by the posterior mass above the iso-contour of no shift ($\Delta\theta=0$). This can be obtained by performing the following integral:

[TABLE]

which is symmetric for changes of datasets A$\leftrightarrow$B and gives us the probability $\Delta$. If $\Delta$ is close to zero, no shift is present and the two datasets are in agreement. On the contrary, a probability $\Delta$ close to one indicates a tension between datasets.

For all the tests considered below we will report significance in terms of number of standard deviations by converting probabilities into two-sided Gaussian standard deviations.

III Analysis

III.1 Technical details

One of the objectives of this analysis is to compute the Suspiciousness tests for which the calculation of Bayesian evidences and Kullback-Leibler divergences is required. In order to obtain such quantities, we perform our numerical scans with PolyChord Handley et al. (2015) and analyse the results using anesthetic Handley (2019). Implementations of other tests are taken from the code Tensiometer 111https://github.com/mraveri/tensiometer.. Concerning the implementation of the parameter differences test, we adopt the Tensiometer when considering multi-dimensional parameter spaces, while we directly implement the integral in Eq. (17) when dealing with only one parameter.

Our numerical implementation considers different parameterizations (see later) for the neutrino masses, which are constrained using a set of cosmological and terrestrial observations. In order to reduce the random fluctuations that arise from the initial sampling of the live points in PolyChord, we repeat the nested sampling runs several times for each data combination and neutrino mass parameterization, varying the number of live points each time between 500 and 1500. The quoted results are taken as the mean of the tension metrics applied to each run separately.

III.2 Parameterizations

Our interest below is focused on studying different constraints on neutrino masses. Considering a model with three massive neutrinos, there are several possible ways to describe their mass spectrum that have been adopted in the past in the context of cosmological studies, e.g., Hannestad and Schwetz (2016); Simpson et al. (2017); Schwetz et al. (2017); Gariazzo et al. (2018); Heavens and Sellentin (2018); Jimenez et al. (2022); Gariazzo et al. (2022). Below we will present results for two representative examples, denoted by “$3\nu$” and “$\Sigma$”, respectively:

• •

$\mathbf{3\nu}$-parameterization: We consider the three neutrino masses $m_{A}$, $m_{B}$, $m_{C}$ as independent parameters in the analyses. After sampling the parameters, the masses are ordered from the smallest to the largest and assigned to the mass eigenstates, depending on the considered mass ordering: $m_{1}<m_{2}<m_{3}$ for NO, $m_{3}<m_{1}<m_{2}$ for IO. Similar to Jimenez et al. (2022) (see also Schwetz et al. (2017); Gariazzo et al. (2022)), we impose a Gaussian prior on the logarithm of the three neutrino masses, with the same mean $\mu$ and standard deviation $\sigma$. Hence, neutrino masses are sampled according to a log-normal distribution, without any prior boundaries.222We have tested alternative sampling methods, such as sampling the masses or the logarithms of the masses uniformly within a given range and then apply the lognormal distribution Gariazzo et al. (2022), leading to similar results. The mean $\mu$ and standard deviation $\sigma$ are hyper-parameters in the analysis. We sample them considering a uniform prior on their logarithm, with bounds $5\cdot 10^{-4}<\mu/\text{eV}<0.3$ and $5\cdot 10^{-3}<\sigma/\text{eV}<20$, respectively, and marginalize over them.

• •

$\mathbf{\Sigma}$-parameterization: We describe the neutrino masses by means of their sum $\sum m_{\nu}$ and the two mass splittings $\Delta m^{2}_{21}$ and $|\Delta m^{2}_{31}|$ Heavens and Sellentin (2018). As, for practical purposes, the current and future cosmological probes considered here only depend on $\sum m_{\nu}$, it is possible to marginalize first the likelihood of terrestrial data over the two mass splittings and then perform the combined analysis or the compatibility analysis with just one free parameter ($\sum m_{\nu}$). We verified that this procedure leads to very similar results as performing the entire calculation with three free parameters ($\sum m_{\nu}$, $\Delta m^{2}_{21}$, $|\Delta m^{2}_{31}|$). For definiteness we show here only the results sampling $\sum m_{\nu}$ with a linear prior, since our checks using a logarithmic prior provide very similar results.

Following Gariazzo et al. (2018), we have considered also a range of other parameterizations, e.g., using either $\sum m_{\nu}$ or the lightest neutrino mass and the two mass splittings $\Delta m^{2}_{21}$ and $|\Delta m^{2}_{31}|$; with uniform prior distribution either on the parameters themselves or on their logarithm. We identified our benchmark parameterizations $3\nu$ and $\Sigma$ described above as representative examples, and therefore we restrict the discussion to the two of them.

III.3 Cosmological and terrestrial information on neutrino masses

The aim of this study is to determine the level of tension between cosmological measurements of neutrino masses and terrestrial constraints on the masses and mass splittings. For that purpose, we shall consider the following data constraints:

• •

Neutrino oscillation constraints are simulated by a Gaussian likelihood on the solar and atmospheric mass differences with mean and standard deviations according to Eq. (4). Note that the $\Delta\chi^{2}$ between NO and IO from oscillation data does not affect the tension metrics and is therefore not relevant for our analyses.

• •

For terrestrial neutrino mass measurements we include the result from KATRIN by adopting a Gaussian likelihood with Aker et al. (2022)

[TABLE]

The region of interest corresponds to quasi-degenerate neutrinos and we can use the approximation for the effective mass parameter relevant for KATRIN:

[TABLE]

Effectively, this provides an upper bound on $\sum m_{\nu}$ for the terrestrial data.

The combination of these two data sets is denoted as “terrestrial” in the following. For the cosmological data we consider current data, as well as two possible future scenarios:

• •

For current cosmological data, we consider the full posterior distribution obtained using Planck temperature, polarization and lensing data together with Supernovae Ia luminosity distance measurements and Baryon Acoustic Oscillations plus Redshift Distortions from SDSS IV, which corresponds to a 95% CL upper limit $\sum m_{\nu}<0.09$ eV Di Valentino et al. (2021).

• •

For future cosmological probes, we shall consider a precision of 0.02 eV on the sum of neutrino masses Font-Ribera et al. (2014); Basse et al. (2014); Hamann et al. (2012); Carbone et al. (2011); Brinckmann et al. (2019) and two alternative scenarios: either a value for $\sum m_{\nu}$ corresponding to the minimal value as predicted for the NO, see Eq. (5),

[TABLE]

or a hypothetical non-observation of finite neutrino masses in cosmology,

[TABLE]

Note that the latter case, by construction, is in tension with oscillation data. We will use the statistical tests discussed above to quantify this statement. In both cases, we assume a Gaussian likelihood for $\sum m_{\nu}$.

Figure 1 shows the posteriors for various data sets using the $\Sigma$-parameterization. We observe the top-hat shaped distribution for terrestrial data (red curves), with the lower bound provided by oscillations (its value depending on NO or IO) and the upper bound provided by KATRIN. The interplay with the assumed cosmological data sets is apparent from the figure, and below we are going to quantify possible existing tensions among them.

IV Results

Let us now present the results of our analysis about the consistency or possible tension between cosmology and terrestrial neutrino mass determinations. The numerical results for the suspiciousness, parameter goodness-of-fit, and parameter shift tests are summarized in Tab. 1. We show the results for the corresponding test statistics as well as significances. The compatibility is tested assuming either the current cosmological likelihood, or possible future determinations of $\sum m_{\nu}$, with the two cases future NO and future 0 discussed in section III.3. Furthermore, we check how the results depend on the type of the neutrino mass ordering (normal versus inverted) as well as on the parameterization used for the neutrino masses ($3\nu$ versus $\Sigma$, see section III.2).

IV.1 Suspiciousness and parameter goodness-of-fit tests

We start by discussing the results for the suspiciousness and the parameter goodness-of-fit test, which are illustrated graphically in Fig. 2. In the upper panels we show the corresponding test statistics $(d-2\ln S)$ and $Q_{\rm DMAP}$, see Eqs. (11) and (14). We see from the figure, that these quantities are numerically very similar for the two tests, as well as for the two parameterizations. For the parameter goodness-of-fit test, the quantity $Q_{\rm DMAP}$ is obtained by taking $\hat{\theta}_{D}$ in Eq. (14) as parameter value at the maximum of the posterior. For the $\Sigma$ parameterization, we have only one relevant parameter (namely $\sum m_{\nu}$) for which we take a flat linear prior. Hence, in this case maximum posterior (MAP) and maximum likelihood (MLH) are identical and therefore $Q_{\rm DMAP}\equiv\chi^{2}_{\rm PG}$. For the $3\nu$ parameterization we adopt flat priors in the logarithm of the three neutrino masses, constrained by hyper-parameters (see section III.2). Hence, here in principle, MAP and MLH are not identical. However, we see from Fig. 2 and Tab. 1 that the $Q_{\rm DMAP}$ values for the $3\nu$ and $\Sigma$ parameterizations are very close, and hence we find also for $3\nu$ that the relationship $Q_{\rm DMAP}\approx\chi^{2}_{\rm PG}$ holds to good accuracy.

Under certain regularity conditions the quantities shown in the upper panels of Fig. 2 follow a $\chi^{2}_{d}$ distribution, with $d$ corresponding to the effective number of parameters in common to the two data sets, as defined in Eqs. (13) and (15), respectively. We give the Bayesian model dimensionalities obtained for the various data set combinations in the right part of Tab. 1. We observe that in many cases Eq. (13) leads to negative values for $d$, which do not correspond to a meaningful $\chi^{2}$ definition. As we discuss in the Appendix, this follows from the properties of Bayesian model dimensionality and the specific shape of the posteriors in our application. Therefore, using Bayesian dimensionalities appears not suitable in our case to evaluate the effective number of degrees of freedom. Instead, we are going to use the simple parameter counting from Eq. (15) also in the case of the suspiciousness test, which gives $n=3$ or 1 for the $3\nu$ or $\Sigma$ parameterization, corresponding either to the 3 neutrino masses or to the singe parameter $\sum m_{\nu}$, respectively.333Using a $\chi^{2}_{d}$ distribution for $(d-2\ln S)$ in any case requires regularity conditions, such as Gaussian-shaped posteriors. Large deviations from $d_{D}$ from naive parameter counting signals non-Gaussian posteriors. The probabilities reported in Tab. 1 and lower panels of Fig. 2 have to be interpreted with care, and are understood under the assumption that $(d-2\ln S)$ follows a $\chi^{2}_{n}$ distribution with $n$ given in Eq. (15). We observe from the lower panels of Fig. 2 that, although the test statistics themselves are very similar, the tension quantified by the corresponding significance is somewhat stronger in the $\Sigma$ parameterization, due to the smaller number of dof. This effect is a known property of the parameter goodness-of-fit test: introducing more model parameters reduces the tension Maltoni and Schwetz (2003).

Let us now discuss the physics results. Notice that current cosmological data (blue symbols) show mild tension with terrestrial data for NO at level of $1-2\sigma$ and for IO at the level of $2-3\sigma$. We can neither claim significant tension, nor disfavour IO due to strong tension. Assuming a future determination of $\sum m_{\nu}$ of 0.06 eV according to the minimal NO value (green symbols), both tests show full consistency of the the data sets for NO and disfavour IO at the level of $2\sigma$. This is expected and in agreement with the trivial observation that for a determination according to Eq. (20), $\sum m_{\nu}\geq 0.1$ eV can be excluded at $2\sigma$. Note also, that the tension for IO in this case even slightly decreases with respect to current data, due to the finite mean value for $\sum m_{\nu}$. In particular, the PG test in the $3\nu$ parameterization signals a tension for IO only at the $1\sigma$ level, i.e., no tension. Moving now to the hypothetical case of no-mass detection of future cosmological data (magenta symbols), we see very strong tension for IO (above $4\sigma$), however, also significant tension for NO (between 2 and $3\sigma$). Hence, rejection of IO on the basis of this tension becomes problematic, as also the alternative hypothesis suffers from a non-negligible tension.

IV.2 Parameter shift test

Figure 3 depicts the corresponding results for the parameter shift test. In general we observe a similar pattern as for the suspiciousness and parameter goodness-of-fit tests, and the physics interpretation is similar. However, we notice in all cases that the parameter shift leads to a higher tension. According to the parameter shift test, current data shows tension of $\gtrsim 2\sigma$ ($\gtrsim 3\sigma$) for NO (IO). The non-observation of neutrino mass by future cosmology will lead to a (very) strong tension with oscillation data regardless of the mass ordering. Even for future NO, some tension close to the 2$\sigma$ level still remains for NO in the case of the $\Sigma$ parameterization.

The reason for the relative stronger tensions obtained with the parameter shift test is a Bayesian volume effect. This test, as defined in Eq. (17), measures the relative size of the overlap of the posterior volumes in parameter space of the two models. As an example, we can see from Fig. 1 that even for the future NO case, the overlap volume with the terrestrial posterior is rather small. The result of the parameter shift test depends on the available parameter volume of the data sets, in particular on the upper bound on $\sum m_{\nu}$ from KATRIN: the tension will become stronger (weaker) for a weaker (stronger) upper bound on $\sum m_{\nu}$, just by increasing (decreasing) the terrestrial posterior volume in the region far away from the cosmological posterior volume.444We have checked that this effect is still relatively small if we use the final KATRIN sensitivity instead of the present result, for which results of the parameter shift test are rather similar. Note also that there is no systematic trend when switching from the $3\nu$ to the $\Sigma$ parameterizations: while for current data the tension becomes weaker, for future NO as well as future 0 it becomes stronger (both for NO and IO).

IV.3 Mass ordering comparison

Let us now briefly compare the tension measures presented above to a direct model comparison of NO versus IO. To this aim we consider the so-called Bayes factor, in analogy to the Bayesian evidence ratio from Eq. (8):

[TABLE]

This quantity describes the Bayesian odds in favour of NO, i.e., large values of $B_{\rm NO,IO}$ correspond to a preference for NO. We convert Bayes factors into probabilities by using $P_{\rm NO}=B_{\rm NO,IO}/(1+B_{\rm NO,IO})$ and $P_{\rm IO}=1/(1+B_{\rm NO,IO})$ (given equal initial prior probabilities).

Figure 4 shows the logarithm of the Bayes factor. Here we show only the contribution from the available parameter space volume from the interplay of cosmological and terrestrial data, in order to compare with the tension measures discussed above. We note that here the direct contribution from the $\chi^{2}$ difference between NO and IO from oscillation data alone de Salas et al. (2021); Esteban et al. (2020); Capozzi et al. (2021) is not considered, which may provide additional NO/IO discrimination in the Bayes factor, see Gariazzo et al. (2022) for a recent discussion.

The black, red, and dark-blue symbols in the figure correspond to using the prior-only, terrestrial data alone, and current cosmology without terrestrial data, respectively. None of these cases shows any significant MO preference. Note that the slightly non-zero value for $\ln B_{\rm NO,IO}$ for terrestrial data in the $3\nu$ parameterization is a pure volume effect. By comparing the $\ln B_{\rm NO,IO}$ results for the combination of terrestrial and cosmological data (light-blue, green and magenta), we observe a significant dependence on the parameterization. This is in line with the arguments discussed in Gariazzo et al. (2022); Schwetz et al. (2017), where it is stressed that parameterizations with three independent neutrino masses in general lead to a strong preference for NO compared to other parameterizations. Indeed, from Fig. 4 we see approximately a difference of $1\sigma$ between the two considered parameterizations.

Concerning future cosmological data, a measurement $\sum m_{\nu}=0.06\pm 0.02$ eV would provide a significance of approximately $2-3\sigma$ in favor of NO. Hence, from this argument alone (i.e., without using additional information from oscillation data), a precision such as the one considered here is not sufficient for a decisive determination of the mass ordering. Within the case “future 0”, for which the measurements provide a preferred value $\sum m_{\nu}=0$ eV, the preference for NO is strong, close to the $4\sigma$ level (even for the $\Sigma$ parameterization). This result, however, is a consequence of the stronger rejection of the region at $\sum m_{\nu}>0.1$ eV with respect to the one at $\sum m_{\nu}>0.06$ eV, and does not take into account that also the NO solution suffers from a tension between cosmology and oscillation data, as discussed in the previous sections.

V Conclusions

The neutrino mass sensitivity from cosmological data analyses is entering an exciting phase, approaching the minimal values for $\sum m_{\nu}$ as required by oscillation data, i.e., $\sum m_{\nu}\approx 0.06$ (0.1) eV for NO (IO). In this manuscript we discuss quantitative measures to evaluate a possible tension between cosmology and terrestrial neutrino mass determinations. In particular we have applied the Bayesian suspiciousness test, parameter goodness-of-fit tests and Bayesian parameter differences, and studied implications for current cosmological data or sensitivities to be expected in the near future. In the latter case we assume an accuracy of 0.02 eV on $\sum m_{\nu}$ and consider two possible scenarios, either a mean value of $\sum m_{\nu}=0.06$ eV, i.e. the minimal value predicted for NO, or $\sum m_{\nu}=0$, i.e., a hypothetical non-observation of neutrino mass in cosmology. Our main conclusions can be summarized as follows:

• •

Current data show modest tension between cosmology and terrestrial data, at the level of $1-2\sigma$ for NO and $2-3\sigma$ for IO.

• •

If future cosmology will find a value of $\sum m_{\nu}\approx 0.06$ eV (corresponding to NO with vanishing lightest neutrino mass) the tension for IO will be at the level of $2-3\sigma$. Hence, the assumed accuracy on $\sum m_{\nu}$ of 0.02 eV is not sufficient to exclude IO with decisive, strong significance.

• •

If future cosmological measurements do not find evidence for a non-zero neutrino mass, the tension with terrestrial data will be at the level of $2-3\sigma$ for NO and $\gtrsim 4\sigma$ for IO. Only in this case IO can be disfavoured strongly, however, at the price of having a tension between cosmology and terrestrial data present also for NO.

• •

Bayesian suspiciousness and parameter goodness-of-fit tests give very similar results. In both cases, tension quantified in terms of significances depends on the number of model parameters.

• •

We find that Bayesian model dimensionality is not a useful measure for the relevant degrees of freedom in our case of interest; our results are based on “naive” parameter counting.

• •

Parameter differences in general show stronger tensions and depend on priors and parameterizations in a non-trivial way.

In conclusion, in this work we have emphasized the well-known fact that the neutrino mass-ordering sensitivity in cosmological data analyses emerges from the available parameter space in the interplay between cosmology and neutrino oscillation data. In such a situation, the relative comparison of NO and IO in terms of model-comparison may be misleading. Excluding one of the hypotheses makes sense only if the alternative hypothesis provides a good fit to the data. We have used statistical tests as tension diagnostics between data sets in order to address this point. We argue that IO can only be excluded in a meaningful way if cosmology finds a result for $\sum m_{\nu}$ consistent with the NO prediction. Quantitatively we find that, based on this argument, an accuracy better than 0.02 eV from cosmological observations will be absolutely required in order to reject the inverted mass ordering with decisive significance.

If the tension between cosmological measurements and oscillation results is found also in NO case, or, eventually, a (positive) neutrino mass signal is detected in either beta-decay or neutrinoless double beta decay experiments, but cosmological observations prefer $\sum m_{\nu}=0$, a modification of the neutrino sector or of the dark sector of the theory would be required (see also the discussion in Ref. Gerbino et al. (2022)). In such a case, dark energy could have a phantom nature Vagnozzi et al. (2018), or the dark matter sector could show unexpected interactions, hiding the effect of neutrino masses Stadler et al. (2020). On the other hand, a modification in the neutrino sector could also explain such potential tensions. Some examples are neutrino decay Franco Abellán et al. (2022); Escudero et al. (2020) or conversion Farzan and Hannestad (2016); Escudero et al. (2023) into dark radiation, long-range neutrino interactions Esteban and Salvado (2021) or time-dependent neutrino masses Dvali and Funcke (2016). It is therefore crucial to confront future cosmological mass limits with neutrino oscillation results, as a tool to constrain beyond the standard model interactions in the invisible sector (neutrinos, dark matter and dark energy) of the theory.

Acknowledgements.

This project has received support from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreements No 754496 (FELLINI) and No 860881 (HIDDeN). OM is supported by the MCIN/AEI/10.13039/501100011033 of Spain under grant PID2020-113644GB-I00, by the Generalitat Valenciana of Spain under the grant PROMETEO/2019/083 and by the European Union’s Framework Programme for Research and Innovation Horizon 2020 (2014–2020) under grant H2020-MSCA-ITN-2019/860881-HIDDeN.

Appendix A Properties of Bayesian model dimensionality

For the calculation of the number of degrees of freedom in a Bayesian context, the authors of Handley and Lemos (2019b) suggest to employ the Bayesian model dimensionality $d_{D}$ defined in our Eq. (12). In their article, the authors show that a $1$-dimensional Gaussian distribution corresponds to $d_{D}=1$, while different probability distributions may correspond to $d_{D}$ significantly different from 1, see their figure 3. For instance, a top-hat distribution gives $d_{D}=0$. We provide further examples below. In our specific case, as we can see from the last four columns in table 1, many of the $d_{D}$ values are quite significantly different from the naive parameter counting.

Let us first analyse the column $d_{A}$, which corresponds to terrestrial experiments. In the $3\nu$ parameterization, we have $d_{A}\simeq 2$, as expected from the fact that terrestrial experiments constrain two mass splittings, while the absolute scale of neutrino masses is not significantly constrained. This is particularly clear from $d_{A}\simeq 0.45$ in the $\Sigma$ case, for which the $\sum m_{\nu}$ posterior partially resembles a top-hat distribution (see the red curves in figure 1), that would correspond to $d_{D}=0$ according to Handley and Lemos (2019b).

The columns $d_{B}$ and $d_{AB}$ are more complicated to interpret, but we can understand the results in this way. Let us consider a one-dimensional Gaussian distribution, which gives $d_{D}=1$ for a sufficiently wide prior range including the maximum. For example, if we consider a one-dimensional Gaussian on some parameter $x$, with $\mu=0$ and $\sigma=1$, we would still get $d_{D}=1$ if we consider a prior that cuts in half the distribution (such as, for example $0\leq x\leq 10$). On the contrary, we can obtain $d_{D}<1$ if we consider an asymmetric range that includes $\mu$ (for example $d_{D}=0.67$ for $-1\leq x\leq 10$) or $d_{D}>1$ if the central value falls outside the interval (e.g. $d_{D}=1.36$ for $1\leq x\leq 10$). Notice that the former case mimics the KATRIN likelihood from Eq. (18), since the central value falls inside the allowed prior range and we get $d_{D}<1$, while the latter corresponds to the combination of terrestrial and cosmological data, for which the central value from cosmology ($\sum m_{\nu}=0$) is outside the $\sum m_{\nu}$ range allowed by oscillation experiments ($d_{AB}>1$ in most of the cases). The only case $d_{AB}=1$ is obtained when the preferred cosmological value is at the prior edge (“future NO” scenario, fitted with NO).

We also noticed that switching from a linear to a logarithmic prior on the considered parameter(s) can generate significant differences in the value of $d_{D}$: it is therefore reasonable to have much wider fluctuations in the Bayesian model dimensionality in the $3\nu$ scenario, which considers three log-normal distributions for the mass parameters, than in the $\Sigma$ case, for which we have a simple linear prior and one varying parameter. Similar conclusions can be drawn for the $d_{B}$ column.