Updated results on neutrino mass and mass hierarchy from cosmology with Planck 2018 likelihoods

TL;DR

This paper updates cosmological bounds on the sum of neutrino masses using Planck 2018 data, considering different mass hierarchies and extended models, finding mild preference for normal hierarchy but no strong evidence.

Contribution

It provides the latest Bayesian constraints on neutrino mass sum across various hierarchies and extended cosmological models using comprehensive Planck and BAO data.

Findings

95% CL bounds: <0.12 eV (degenerate), <0.15 eV (normal), <0.17 eV (inverted)

Normal hierarchy is mildly preferred over inverted

No strong evidence favoring one hierarchy in extended models

Abstract

In this work we update the bounds on $\sum m_{\nu}$ from latest publicly available cosmological data and likelihoods using Bayesian analysis, while explicitly considering particular neutrino mass hierarchies. In the minimal $\Lambda\textrm{CDM}+\sum m_{\nu}$ model with most recent CMB data from Planck 2018 TT,TE,EE, lowE, and lensing; and BAO data from BOSS DR12, MGS, and 6dFGS, we find that at 95\% C.L. the bounds are: $\sum m_{\nu}<0.12$ eV (degenerate), $\sum m_{\nu}<0.15$ eV (normal), $\sum m_{\nu}<0.17$ eV (inverted). The bounds vary across the different mass orderings due to different priors on $\sum m_{\nu}$. Also, we find that the normal hierarchy is very mildly preferred relative to the inverted, using both minimum $\chi^2$ values and Bayesian Evidence ratios. In this paper we also provide bounds on $\sum m_{\nu}$ considering different hierarchies in various extended cosmological models: $\Lambda\textrm{CDM}+\sum m_{\nu}+r$, $w\textrm{CDM}+\sum m_{\nu}$, $w_0 w_a \textrm{CDM}+\sum m_{\nu}$, $w_0 w_a \textrm{CDM}+\sum m_{\nu}$ with $w(z)\geq -1$, $\Lambda \textrm{CDM} + \sum m_{\nu} + \Omega_k$, and $\Lambda \textrm{CDM} + \sum m_{\nu} + A_{\textrm{Lens}}$. We do not find any strong evidence of normal hierarchy over inverted hierarchy in the extended models either.

Tables

Table 1. Table 1 : Flat priors on the main cosmological parameters constrained in this paper in the analyses with CosmoMC.

Parameter	Prior
${\omega }_c$	[0.005, 0.1]
${\omega }_b$	[0.001, 0.99]
${\Theta }_s$	[0.5, 10]
$\tau$	[0.01, 0.8]
$n_{\text{s}}$	[0.8, 1.2]
$\text{ln}\left[{10}^{10}{A}_{\text{s}}\right]$	[2, 4]
$m_0$ (eV)	[0, 1]
$r$	[0, 1]
$w$	[-3, -0.33]
$w_0$	[-3, -0.33]
$w_a$	[-2, 2]
${\Omega }_k$	[-0.3, 0.3]
$A_{\text{Lens}}$	[0, 3]

Table 2. Table 2 : Constraints on the cosmological parameters in the minimal Λ ​ CDM + ∑ m ν Λ CDM subscript 𝑚 𝜈 \Lambda\textrm{CDM}+\sum m_{\nu} model considering three different hierarchies (degenerate, normal, and inverted) with the Base and Base+SNe datasets, where Base ≡ \equiv Planck 2018 TTTEEE+lowE+lensing+BAO. Here m 0 subscript 𝑚 0 m_{0} is the mass of the lightest neutrino in a particular hierarchy and a freely varying parameter in the model, whereas ∑ m ν subscript 𝑚 𝜈 \sum m_{\nu} , H 0 subscript 𝐻 0 H_{0} , σ 8 subscript 𝜎 8 \sigma_{8} , and S 8 subscript 𝑆 8 S_{8} are derived parameters. Marginalized constraints are given at 1 σ 𝜎 \sigma whereas upper or lower bounds are given at 2 σ 𝜎 \sigma . The χ 2 superscript 𝜒 2 \chi^{2} differences are calculated at best-fit points. Details about models and datasets are given in section 2 .

	Base				Base+SNe
	DH	NH	IH		DH	NH	IH
$\Lambda \text{CDM}+\sum m_{\nu }$
${\omega }_c$	$0.1194\pm 0.0009$	$0.1192\pm 0.0009$	$0.1191\pm 0.0009$		$0.1193\pm 0.0009$	$0.1191\pm 0.0009$	$0.1189\pm 0.0009$
${\omega }_b$	$0.02242\pm 0.00013$	${0.02242}_{-0.00014}^{+0.00013}$	$0.02243\pm 0.00013$		$0.02243\pm 0.00013$	$0.02244\pm 0.00013$	$0.02244\pm 0.00013$
${\Theta }_{\text{s}}$	$1.04100\pm 0.00029$	$1.04100\pm 0.00029$	$1.04100\pm 0.00029$		$1.04102\pm 0.00029$	$1.04103\pm 0.00029$	$1.04103\pm 0.00029$
$\tau$	${0.0554}_{-0.0076}^{+0.0068}$	${0.0569}_{-0.0076}^{+0.0066}$	${0.0585}_{-0.0076}^{+0.0069}$		$0.0556\pm 0.0071$	${0.0573}_{-0.0076}^{+0.0069}$	${0.0588}_{-0.0077}^{+0.0068}$
$n_{\text{s}}$	$0.9666\pm 0.0036$	$0.9668\pm 0.0037$	$0.9671\pm 0.0037$		$0.9669\pm 0.0036$	$0.9673\pm 0.0036$	$0.9675\pm 0.0037$
$\text{ln}\left[{10}^{10}{A}_{\text{s}}\right]$	${3.048}_{-0.015}^{+0.014}$	${3.051}_{-0.015}^{+0.014}$	$3.053\pm 0.015$		$3.046\pm 0.014$	$3.049\pm 0.014$	${3.052}_{-0.015}^{+0.014}$
$m_0$ (eV)
$\sum m_{\nu }$ (eV)
$H_0$ (km/s/Mpc)	${67.81}_{-0.46}^{+0.54}$	${67.50}_{-0.44}^{+0.49}$	$67.22\pm 0.45$		${67.89}_{-0.45}^{+0.52}$	$67.59\pm 0.44$	$67.33\pm 0.43$
${\sigma }_8$	${0.814}_{-0.007}^{+0.010}$	${0.806}_{-0.006}^{+0.009}$	${0.799}_{-0.006}^{+0.008}$		${0.815}_{-0.007}^{+0.010}$	${0.806}_{-0.006}^{+0.008}$	${0.799}_{-0.006}^{+0.008}$
$S_8$	$0.827\pm 0.011$	$0.823\pm 0.011$	$0.820\pm 0.011$		$0.826\pm 0.011$	$0.822\pm 0.011$	$0.818\pm 0.011$
$\Delta {\chi }^2={\chi }^2-{\chi }_{IH}^2$	$-2.89$	$-0.95$	0		$-2.73$	$-1.27$	$0$

Table 3. Table 3 : Constraints on selected cosmological parameters in the extended models considering three different hierarchies (degenerate, normal, and inverted) with the Base and Base+SNe datasets. In the Λ ​ CDM + ∑ m ν + r Λ CDM subscript 𝑚 𝜈 𝑟 \Lambda\textrm{CDM}+\sum m_{\nu}+r model Base data also includes BK15.

	Base				Base+SNe
	DH	NH	IH		DH	NH	IH
$\Lambda \text{CDM}+\sum m_{\nu }+r$
$r$
$m_0$ (eV)
$\sum m_{\nu }$ (eV)
$H_0$ (km/s/Mpc)	${67.78}_{-0.46}^{+0.52}$	${67.48}_{-0.45}^{+0.48}$	$67.21\pm 0.44$		$67.86\pm 0.48$	${67.58}_{-0.45}^{+0.44}$	${67.30}_{-0.44}^{+0.43}$
${\sigma }_8$	${0.815}_{-0.007}^{+0.010}$	${0.807}_{-0.006}^{+0.008}$	${0.800}_{-0.006}^{+0.008}$		${0.816}_{-0.007}^{+0.010}$	${0.808}_{-0.007}^{+0.008}$	${0.800}_{-0.006}^{+0.007}$
$S_8$	$0.832\pm 0.011$	${0.825}_{-0.010}^{+0.011}$	$0.822\pm 0.011$		$0.828\pm 0.011$	$0.824\pm 0.010$	$0.820\pm 0.010$
$w\text{CDM}+\sum m_{\nu }$
$w$	$-{1.042}_{-0.052}^{+0.072}$	$-{1.068}_{-0.052}^{+0.071}$	$-{1.089}_{-0.055}^{+0.070}$		$-{1.025}_{-0.033}^{+0.037}$	$-{1.037}_{-0.032}^{+0.036}$	$-1.048\pm 0.034$
$m_0$ (eV)
$\sum m_{\nu }$ (eV)
$H_0$ (km/s/Mpc)	${68.67}_{-1.59}^{+1.33}$	${69.01}_{-1.60}^{+1.31}$	${69.24}_{-1.61}^{+1.38}$		$68.33\pm 0.82$	$68.31\pm 0.82$	$68.27\pm 0.82$
${\sigma }_8$	${0.821}_{-0.017}^{+0.016}$	$0.820\pm 0.016$	$0.819\pm 0.016$		${0.818}_{-0.011}^{+0.013}$	$0.814\pm 0.012$	$0.810\pm 0.012$
$S_8$	$0.825\pm 0.011$	${0.822}_{-0.011}^{+0.012}$	$0.822\pm 0.011$		$0.826\pm 0.011$	$0.823\pm 0.011$	$0.821\pm 0.011$
$w_0{w}_a\text{CDM}+\sum m_{\nu }$
$w_0$	$-{0.68}_{-0.14}^{+0.26}$	$-{0.68}_{-0.14}^{+0.26}$	$-{0.68}_{-0.13}^{+0.25}$		$-{0.94}_{-0.09}^{+0.08}$	$-0.94\pm 0.09$	$-0.93\pm 0.09$
$w_a$	$-{1.06}_{-0.79}^{+0.37}$				$-{0.41}_{-0.29}^{+0.46}$	$-{0.49}_{-0.30}^{+0.44}$	$-{0.56}_{-0.32}^{+0.43}$
$m_0$ (eV)
$\sum m_{\nu }$ (eV)
$H_0$ (km/s/Mpc)	${65.70}_{-2.47}^{+1.60}$	${65.78}_{-2.47}^{+1.61}$	${65.80}_{-2.43}^{+1.62}$		$68.28\pm 0.83$	$68.23\pm 0.84$	$68.23\pm 0.82$
${\sigma }_8$	${0.795}_{-0.023}^{+0.018}$	${0.792}_{-0.023}^{+0.017}$	${0.790}_{-0.023}^{+0.018}$		${0.817}_{-0.013}^{+0.015}$	${0.813}_{-0.012}^{+0.014}$	${0.811}_{-0.012}^{+0.013}$
$S_8$	$0.837\pm 0.014$	${0.834}_{-0.013}^{+0.014}$	$0.832\pm 0.013$		$0.827\pm 0.012$	$0.826\pm 0.012$	$0.824\pm 0.012$
$w_0{w}_a\text{CDM}+\sum m_{\nu }$($w\left(z\right)\ge -1)$
$w_0$
$w_a$	${0.009}_{-0.067}^{+0.057}$	${0.007}_{-0.058}^{+0.049}$	${0.007}_{-0.050}^{+0.044}$		${0.028}_{-0.056}^{+0.034}$	${0.022}_{-0.047}^{+0.029}$	${0.020}_{-0.043}^{+0.025}$
$m_0$ (eV)
$\sum m_{\nu }$ (eV)
$H_0$ (km/s/Mpc)	${66.64}_{-0.66}^{+0.97}$	${66.46}_{-0.62}^{+0.88}$	${66.33}_{-0.57}^{+0.83}$		${67.23}_{-0.53}^{+0.63}$	${67.01}_{-0.51}^{+0.57}$	${66.81}_{-0.48}^{+0.54}$
${\sigma }_8$	${0.801}_{-0.010}^{+0.012}$	${0.795}_{-0.009}^{+0.011}$	${0.789}_{-0.008}^{+0.010}$		${0.807}_{-0.008}^{+0.010}$	${0.799}_{-0.008}^{+0.009}$	$0.793\pm 0.008$
$S_8$	$0.826\pm 0.011$	$0.823\pm 0.011$	$0.820\pm 0.011$		$0.824\pm 0.010$	$0.821\pm 0.011$	$0.817\pm 0.011$
$\Lambda \text{CDM}+\sum m_{\nu }$+$A_{\text{Lens}}$
$A_{\text{Lens}}$	${1.100}_{-0.056}^{+0.046}$	${1.107}_{-0.055}^{+0.042}$	${1.116}_{-0.050}^{+0.040}$		${1.098}_{-0.058}^{+0.044}$	${1.106}_{-0.053}^{+0.042}$	${1.116}_{-0.050}^{+0.040}$
$m_0$ (eV)
$\sum m_{\nu }$ (eV)
$H_0$ (km/s/Mpc)	${67.76}_{-0.60}^{+0.71}$	${67.66}_{-0.59}^{+0.66}$	${67.56}_{-0.53}^{+0.63}$		$67.86{\pm }_{-0.57}^{+0.68}$	$67.76\pm 0.59$	${67.66}_{-0.51}^{+0.58}$
${\sigma }_8$	${0.782}_{-0.018}^{+0.029}$	${0.777}_{-0.013}^{+0.025}$	${0.772}_{-0.012}^{+0.022}$		${0.784}_{-0.017}^{+0.028}$	${0.779}_{-0.013}^{+0.024}$	${0.773}_{-0.011}^{+0.021}$
$S_8$	${0.793}_{-0.020}^{+0.023}$	${0.790}_{-0.017}^{+0.022}$	${0.786}_{-0.016}^{+0.020}$		${0.793}_{-0.019}^{+0.023}$	${0.790}_{-0.017}^{+0.021}$	${0.785}_{-0.016}^{+0.019}$
$\Lambda \text{CDM}+\sum m_{\nu }+{\Omega }_{\text{k}}$
${\Omega }_{\text{k}}$	$0.0004\pm 0.0020$	${0.0012}_{-0.0021}^{+0.0020}$	${0.0019}_{-0.0021}^{+0.0019}$		$0.0004\pm 0.0021$	$0.0012\pm 0.0020$	$0.0019\pm 0.0020$
$m_0$ (eV)
$\sum m_{\nu }$ (eV)
$H_0$ (km/s/Mpc)	$67.87\pm 0.67$	$67.75\pm 0.67$	${67.67}_{-0.68}^{+0.67}$		$67.95\pm 0.67$	$67.84\pm 0.66$	$67.77\pm 0.66$
${\sigma }_8$	${0.813}_{-0.008}^{+0.011}$	${0.807}_{-0.008}^{+0.010}$	${0.801}_{-0.008}^{+0.009}$		${0.814}_{-0.009}^{+0.010}$	${0.807}_{-0.008}^{+0.010}$	${0.801}_{-0.008}^{+0.009}$
$S_8$	$0.826\pm 0.011$	$0.823\pm 0.011$	$0.819\pm 0.011$		$0.825\pm 0.011$	$0.822\pm 0.011$	$0.818\pm 0.010$

Table 4. Table 4 : The values of Δ ​ χ 2 = Δ ​ AIC = χ 2 − χ IH 2 Δ superscript 𝜒 2 Δ AIC superscript 𝜒 2 subscript superscript 𝜒 2 IH \Delta\chi^{2}=\Delta\textrm{AIC}=\chi^{2}-\chi^{2}_{\rm IH} for various extended models studied in this paper, with Base and Base+SNe dataset combinations. The χ 2 superscript 𝜒 2 \chi^{2} differences are calculated at best-fit points.

	Base				Base+SNe
	DH	NH	IH		DH	NH	IH
$\Lambda \text{CDM}+\sum m_{\nu }+r$	$-5.32$	$-1.39$	0		$-4.46$	$-3.21$	0
$w\text{CDM}+\sum m_{\nu }$	$-1.66$	$-1.59$	0		$-0.84$	$-0.99$	0
$w_0{w}_a\text{CDM}+\sum m_{\nu }$	$+0.75$	$+2.08$	0		$-0.80$	$-0.41$	0
$w_0{w}_a\text{CDM}+\sum m_{\nu }$ with $w\left(z\right)\ge -1$	$-2.71$	$-0.69$	0		$-4.67$	$-1.31$	0
$\Lambda \text{CDM}+\sum m_{\nu }+A_{\text{Lens}}$	$+0.88$	$+0.20$	0		$+0.13$	$+0.92$	0
$\Lambda \text{CDM}+\sum m_{\nu }+{\Omega }_{\text{k}}$	$-3.56$	$-1.96$	0		$-2.74$	$-1.22$	0

Table 5. Table 5 : The values of Odds(NH:IH) and CL IH and strength of evidence for various cosmological models studied in this paper, with Base dataset.

	${\pi }_{\mathrm{NH}}={\pi }_{\mathrm{IH}}=0.5$				${\pi }_{\mathrm{NH}}=0.9985$,${\pi }_{\mathrm{IH}}=0.0015$
	Odds(NH:IH)	CLIH	Strength of Evidence		Odds(NH:IH)	CLIH	Strength of Evidence
$\Lambda \text{CDM}+\sum m_{\nu }$	1.43:1	58.8%	Inconclusive		950:1	99.895%	Strong
$\Lambda \text{CDM}+\sum m_{\nu }+r$	0.93:1	48.1%	Inconclusive		617:1	99.838%	Strong
$w\text{CDM}+\sum m_{\nu }$	1.3:1	56.6%	Inconclusive		867:1	99.885%	Strong
$w_0{w}_a\text{CDM}+\sum m_{\nu }$	1.24:1	55.4%	Inconclusive		827:1	99.879%	Strong
$w_0{w}_a\text{CDM}+\sum m_{\nu }$ with $w\left(z\right)\ge -1$	2:1	66.6%	Inconclusive		1329:1	99.925%	Strong
$\Lambda \text{CDM}+\sum m_{\nu }+A_{\text{Lens}}$	0.75:1	42.9%	Inconclusive		500:1	99.800%	Strong
$\Lambda \text{CDM}+\sum m_{\nu }+{\Omega }_{\text{k}}$	2.3:1	69.7%	Inconclusive		1530:1	99.935%	Strong

Equations

$$\Delta m_{21}^2=7.39_{-0.20}^{+0.21}\times 10^{-5};\Delta m_{31}^2=2.525_{-0.032}^{+0.033}\times 10^{-3}(\text{NH});\Delta m_{32}^2=-2.512_{-0.032}^{+0.034}\times 10^{-3}(\text{IH}).$$

$$\sum m_{\nu }=m_1+m_2+m_3.$$

$$\sum m_{\nu }=m_0+\sqrt{\Delta m_{21}^2+m_0^2}+\sqrt{\Delta m_{31}^2+m_0^2}(\text{NH}),$$

$$\sum m_{\nu }=m_0+\sqrt{|\Delta m_{32}^2|+m_0^2}+\sqrt{|\Delta m_{32}^2|-\Delta m_{21}^2+m_0^2}(\text{IH}).$$

$$\theta \equiv \left[{\omega }_c,{\omega }_b,{\Theta }_{\text{s}},\tau ,n_{\text{s}},\text{ln}[10^{10}{A}_{\text{s}}],m_0\right].$$

$$w(z)=w_0+w_a(1-a)=w_0+w_a\frac{z}{1+z}.$$

$$w_0\ge -1;w_0+w_a\ge -1$$

$$\text{AIC}={\chi }_{\text{best-fit}}^2+2k$$

$$H(z)=\sqrt{{\omega }_{\gamma }(1+z{)}^4+({\omega }_c+{\omega }_b)(1+z{)}^3+{\Omega }_{DE}(z)h^2+{\rho }_{\nu }(z)h^2/{\rho }_{cr,0}}.$$

$${\Omega }_{DE}(z)={\Omega }_{DE}(0)(1+z{)}^{3(1+w_0+w_a)}\text{exp}\left(-3w_a\frac{z}{1+z}\right).$$

$$H(z)=\sqrt{{\omega }_{\gamma }(1+z{)}^4+({\omega }_c+{\omega }_b)(1+z{)}^3+{\omega }_{\Lambda }+{\omega }_k(1+z{)}^2+{\rho }_{\nu }(z)h^2/{\rho }_{cr,0}}$$

$$p_i\equiv P(M_i|D)=\frac{{\pi }_{i}P(D|M_i)}{P(D)},$$

$$P(D|M_i)=\int d{\theta }_{i}P(D|{\theta }_i,M_i)P({\theta }_i|M_i).$$

$$p_i\equiv P(M_i|D)=\frac{{\pi }_{i}P(D|M_i)}{{\Sigma }_i{\pi }_{i}P(D|M_i)}.$$

$$\mathrm{Odds}(\mathrm{NH}:\mathrm{IH})\equiv \frac{p_{\mathrm{NH}}}{p_{\mathrm{IH}}}$$

$$\mathrm{C}{\mathrm{L}}_{\mathrm{IH}}\equiv 1-p_{\mathrm{IH}}.$$

Full text

1 Introduction

Earth based neutrino oscillation experiments [1, 2, 3, 4, 5, 6, 7, 8] have confirmed that neutrinos are massive, which is the first known departure from the Standard Model of particle physics where neutrinos are considered to be massless. The 3 neutrino flavor states ($\nu_{e}$, $\nu_{\mu}$, $\nu_{\tau}$) are quantum superpositions of the 3 mass eigenstates ($\nu_{i}$, with respective distinct masses, $m_{i}$ for $i=1,2,3$). However, because oscillation experiments use ultra relativistic neutrinos they are only sensitive to the mass-squared splittings ($\Delta m_{ij}^{2}=m_{i}^{2}-m_{j}^{2}$) and not the absolute masses, thus keeping the mass of the lightest neutrino unbounded. On the other hand, while magnitudes of $\Delta m_{21}^{2}$ and $\Delta m_{31}^{2}$ are known to considerable accuracy from the current oscillation data, sign of $\Delta m_{31}^{2}$ is unknown. This leads to two possible hierarchies of neutrino masses: $m_{1}<m_{2}\ll m_{3}$ (normal hierarchy or NH) and $m_{3}\ll m_{1}<m_{2}$ (inverted hierarchy or IH) depending on whether $\Delta m_{31}^{2}$ is positive or negative, respectively. As per the latest NuFit 4.0 [9] global analysis of oscillations data, the values of the mass-squared splittings (in units of eV2) are (limits are given at 1$\sigma$):

[TABLE]

Here, the value for $\Delta m_{21}^{2}$ is applicable to both NH and IH, while the other values are for the particular hierarchies as mentioned in the brackets. It is to be noted that for IH, $\Delta m_{32}^{2}$ is provided (whose value is negative) instead of $\Delta m_{31}^{2}$, but since $\Delta m_{21}^{2}\ll\Delta m_{32}^{2}$, sign of $\Delta m_{31}^{2}$ for IH is also negative (since $\Delta m_{31}^{2}=\Delta m_{32}^{2}+\Delta m_{21}^{2}$). See [10, 11, 12, 13, 14, 15] for results from other global analyses.

A solution to the neutrino mass hierarchy problem may come from cosmology, which currently provides the strongest bounds on the absolute neutrino mass scale, defined as the sum of the three active neutrino masses,

[TABLE]

As far as known physics goes, at temperatures $T\gg{\rm MeV}$, neutrinos remain in equilibrium with the primordial plasma through the standard model weak interactions. At around $T\sim{\rm MeV}$ neutrinos decouple from the plasma and start free streaming. When neutrinos are relativistic ($T\gg m_{\nu}$) they contribute to the radiation energy density. This continues until much later when they turn non-relativistic at temperatures $T\sim m_{\nu}$, and then they contribute to the matter energy density. Effects of massive neutrinos on cosmological observables have been extensively studied in the literature [16, 17, 18, 19, 20, 21, 22] and these effects help in constraining the sum of neutrino masses. If we consider neutrinos with masses $\ll 1$ eV, at the time of photon decoupling they are still relativistic, and their mass has very limited effect on the photon perturbations and their evolution. Hence, for the primary CMB signal, the effect of small neutrino masses can only be seen through the background evolution, and secondary anisotropies like Integrated Sachs-Wolfe (ISW) effect, and these too can be compensated partially by varying other free parameters in the $\Lambda$CDM model. Thus, strong bounds on $\sum m_{\nu}$ cannot be obtained with CMB power spectrum data alone. On the other hand, at late times, neutrinos affect the evolution of matter perturbations to a large extent. Due to the free-streaming effect of neutrinos, i.e. large thermal velocities, neutrinos do not cluster on small length scales, and this causes increasing suppression of small scale matter power spectrum with increasing fraction of neutrino energy density with respect to the total matter density [20]. Thus if we augment CMB anisotropy data with data coming Large Scale Structure (LSS), CMB lensing (which is the weak lensing effect on CMB photons due to LSS) measurements etc, strong bounds on $\sum m_{\nu}$ can be obtained. Even so, currently it is only possible to get an upper bound on $\sum m_{\nu}$ from cosmological data alone.

Let us define the mass of the lightest neutrino mass eigenstate to be $m_{0}$. For normal hierarchy, $m_{0}=m_{1}$, whereas for inverted hierarchy, $m_{0}=m_{3}$. In terms of $m_{0}$, the sum of the neutrino masses can be defined as,

[TABLE]

and

[TABLE]

Putting $m_{0}=0$, one can obtain the minimum neutrino mass sums allowed by the two possible hierarchies, and these are $\sum m_{\nu}=0.05885^{+0.00045}_{-0.00044}$ eV (1$\sigma$)(NH) and $\sum m_{\nu}=0.100^{+0.00070}_{-0.00067}$ eV (1$\sigma$)(IH). Assuming that the normal hierarchy is the true one, the way cosmology can help is by constraining the $\sum m_{\nu}$ below the minimum mass required by inverted hierarchy with reasonable statistical significance, e.g., a determination of $\sum m_{\nu}=0.058\pm 0.008$ eV (1$\sigma$) will exclude inverted hierarchy at 5$\sigma$ and also provide cosmological evidence for non-zero neutrino masses at 7.25 $\sigma$. Currently, the most recent and strongest bounds on $\sum m_{\nu}$ in the minimal $\Lambda$CDM+$\sum m_{\nu}$ model are around $\sum m_{\nu}<0.12$ eV (95% C.L.) [23, 24] with CMB and BAO data, whereas some other studies had reported bounds around $\sum m_{\nu}<0.15$ eV (95%) or better [25, 26, 27, 28, 29, 30, 31, 32, 15, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43]. The strongest bounds are very close to the minimum $\sum m_{\nu}$ required for inverted hierarchy and thus it seems that inverted hierarchy is starting to get under pressure from cosmological data. These bounds are obtained with the assumption that all the three neutrino masses are equal ($m_{i}=\sum m_{\nu}/3$ for $i=1,2,3$), an approximation we denote as degenerate hierarchy (DH). There are also studies covering interesting aspects of measurement of neutrino hierarchy from cosmology [44, 25, 45, 46, 47, 48, 49, 50, 51, 52, 43, 53]. The current cosmological data, however, is not yet sensitive enough to make the distinction between the two hierarchies on a level that can be considered statistically significant, but there is a small preference for normal hierarchy [44, 25]. It is to be noted that the bounds depend on the underlying cosmological model, and any extensions to the minimal $\Lambda\textrm{CDM}+\sum m_{\nu}$ model will usually lead to a more relaxed bound on $\sum m_{\nu}$ [23, 25, 34, 35, 54, 55, 56, 57, 58, 59, 60]. However, it can happen that the neutrino mass bound improves when the extension to $\Lambda$CDM is done in such a way that the allowed parameter space of the new parameters prefers neutrino masses which are smaller than what we get in $\Lambda$CDM. In fact this is the case when one incorporates dynamical dark energy in the cosmological model but restricts the parameter space to non-phantom dark energy only, i.e. neutrino mass bounds in a cosmology with non-phantom dynamical dark energy are stronger than that in $\Lambda\textrm{CDM}$ [23, 41, 61]. Another such model where neutrino mass bounds are stronger than that in $\Lambda\textrm{CDM}$ is holographic dark energy (HDE) [62, 63].

From now on, we shall use the abbreviations DH, NH, and IH for degenerate, normal and inverted hierarchies respectively. While the degenerate hierarchy approximation has been predominant in the literature, and makes sense when neutrino masses are relatively large compared to the square-root of the mass-squared splittings (i.e. $m_{i}\gg\sqrt{\Delta m_{ij}^{2}}$), the cosmological neutrino mass bound is becoming strong enough that it should be replaced by a treatment using either the normal or the inverted hierarchy. Hence, in this paper we have updated the bounds on the $\sum m_{\nu}$ while explicitly considering three different hierarchies (degenerate, normal and inverted), using latest datasets and likelihoods that are publicly available, for the minimal $\Lambda\textrm{CDM}+\sum m_{\nu}$ and some of its extensions. Except in the case of extension with the tensor-to-scalar ratio ($r$) parameter, all the other extensions studied in this paper includes new parameters which are considerably correlated with the sum of neutrino masses in the datasets considered. Details of the models are given in the next section. The neutrino mass bounds are supposed to relax in most of the extended models, and the difference between the upper limits obtained for the three hierarchies is supposed to diminish as the individual masses become much larger than the square root of mass-squared splittings. But still, our motivation to study these extended models is to see whether the latest datasets can make a difference.

To implement the normal and inverted hierarchy, we use the mean values of the mass squared splittings given in eq.(1.1) along with the lightest neutrino mass $m_{0}$ to define $m_{1}$, $m_{2}$, and $m_{3}$, and use $m_{0}$ as a free parameter and $\sum m_{\nu}$ as a derived parameter. We ignore the errors in the measurement of the mass-squared splittings from oscillations data since they are small compared to the mean values and incorporating them would have a very small effect on the bounds on $\sum m_{\nu}$. For degenerate hierarchy, we simply have $\sum m_{\nu}=3m_{0}$.

It is to be understood here that the physical parameter cosmological data is primarily sensitive to $\sum m_{\nu}$ and not the individual neutrino masses. Even if we consider very optimistic future cosmological data, same value of $\sum m_{\nu}$ but different neutrino mass hierarchies lead to changes which are smaller than the expected observational errors and modeling uncertainties, i.e. it will be impossible to differentiate individual neutrino masses even with future data [64]. Thus, while analysing cosmological data, it makes sense to vary $\sum m_{\nu}$ as the cosmological parameter. However, same is not true for other experiments such as neutrino oscillations, kinematic measurements from beta decay etc and thus using the lightest mass $m_{0}$ instead of $\sum m_{\nu}$ allows one to be more general in their approach, and leaves the possibility of incorporating non-cosmological datasets for further analysis [46] open. If one uses $\sum m_{\nu}$ as the cosmological parameter, a good way to incorporate DH, NH, and IH is to use the following non-informative flat priors: $\sum m_{\nu}\geq 0$ (DH), $\sum m_{\nu}\geq 0.06$ (NH), $\sum m_{\nu}\geq 0.10$ eV (IH), with suitable upper limits. Due to eq. (1.3) and (1.4), if one uses a flat prior on $m_{0}$ instead of on $\sum m_{\nu}$, the effective prior on $\sum m_{\nu}$ is not guaranteed to be for the NH and IH case, and requires further inspection. However, as shown in figure 6 of Ref [46], a flat prior on the lightest mass leads to an almost flat prior on $\sum m_{\nu}$, except that the prior probability rises only slightly close to the lowest possible $\sum m_{\nu}$ values in each hierarchical scenario. Thus it is expected that the bounds on $\sum m_{\nu}$ coming from the NH and IH scenarios with a flat prior on $m_{0}$ will be slightly stronger than what we can get by using a flat prior on $\sum m_{\nu}$. As we shall see that the differences in the bounds on $\sum m_{\nu}$ obtained using flat prior on $m_{0}$ are very close to the ones with a flat prior on $\sum m_{\nu}$ and thus using $m_{0}$ as a varying parameter is well-motivated.

Another important point is that current cosmological datasets are only able to provide an upper bound to the $\sum m_{\nu}$. This bound changes when one changes the neutrino mass hierarchy model from DH to NH to IH. This happens due to volume effects, i.e. change in the prior on $\sum m_{\nu}$. However, in case future cosmological data is able to detect $\sum m_{\nu}$, it has been shown in [65], that irrespective of the correct neutrino mass model (i.e. NH or IH), using the unrealistic DH approximation (with the prior $\sum m_{\nu}\geq 0$) will actually lead to the recovery of the correct values of $\sum m_{\nu}$ with only a small reconstruction bias due to the assumption of wrong fitting model. A more recent study [64] also made a similar conclusion. Thus, the degenerate hierarchy approximation, while unrealistic, is still useful in cosmology. However, while using DH, we must be wary of any reconstruction bias in the recovered values of $\sum m_{\nu}$ in case it is detected from future data.

We would like to point it out here that the choice for a proper prior on neutrino masses has been a topic of intense debate among the concerned researchers, especially because different prior choices can lead to very different bounds on $\sum m_{\nu}$ and very different inferences on the hierarchy issue. See [44, 46, 25, 49, 50, 51, 66, 47, 48] for discussions on this topic.

Among datasets, for CMB anisotropies we use the most recently released Planck 2018 likelihoods for the data on temperature, E mode polarisation and their cross-correlation. Other than Planck CMB anisotropies, we use latest data from measurements of Planck lensing, CMB B mode, BAO, and SNe Ia luminosity distance.

This paper is structured as follows. In section 2 we provide brief details about the cosmological models and datasets used in this paper. In section 3 we provide and explain the results of our MCMC analyses. In section 4 we also include results from additional analyses using nested sampling and provide a comparison of bayesian evidence between NH and IH, and quantify the evidence against IH closely following [44]. In section 5 we conclude.

2 Methodology: Models and datasets

2.1 Models

In this work we have performed our analyses using a variety of different cosmological models which we shall describe below. Note that when we label models we use the term $\sum m_{\nu}$ to refer to the neutrino mass. This is done in order to conform to the standard labelling used in cosmological parameter analyses. In fact all our runs use a flat prior on $m_{0}$, the mass of the lightest mass state, rather than a flat prior on $\sum m_{\nu}$. In total we have investigated the following 7 sets of cosmological models:

• •

i) The minimal $\Lambda\textrm{CDM}+\sum m_{\nu}$ model: Below we list the vector of varying parameters in this model.

[TABLE]

Here the first six parameters correspond to the $\Lambda\textrm{CDM}$ model. $\omega_{c}=\Omega_{c}h^{2}$ and $\omega_{b}=\Omega_{b}h^{2}$ are the present cold dark matter and baryon energy densities respectively. $\Theta_{\textrm{s}}$ is the ratio between sound horizon $r_{s}$ and angular diameter distance $D_{\textrm{A}}$ at the time of photon decoupling. $\tau$ is the optical depth to re-ionization of the universe at late times. $n_{\textrm{s}}$ and $A_{\textrm{s}}$, on the other hand, relate to early universe cosmology. They are the power-law spectral index and power of the primordial scalar perturbations respectively, evaluated at the pivot scale of $k_{*}=0.05h$ Mpc*-1*. As defined in the previous section, the seventh parameter, $m_{0}$ is the mass of the lightest neutrino and it is the parameter of primary concern in this paper.

• •

ii) The $\Lambda\textrm{CDM}+\sum m_{\nu}+r$ model: Apart from the main 7 parameters, in this model we also include the evolution of tensor perturbations in the analysis along with scalar perturbations, and add an additional free parameter $r$, which is the tensor-to-scalar ratio evaluated at the same pivot scale as $n_{s}$ and $A_{s}$.

• •

iii) The $w\textrm{CDM}+\sum m_{\nu}$ model: Here instead of a cosmological constant with a dark energy equation of state (DE EoS hereafter) fixed at $w(z)=-1$ we opt for a DE EoS $w$ which varies as a free parameter, but does not vary in time (i.e. $w$ can assume different values but the values are fixed throughout the evolution history of the universe). Here $z$ denotes the cosmological redshift ($z=1/a-1$, where $a$ is the scale factor of FRW metric). Here, as well as in all models including non-$\Lambda$ dark energy, we use the PPF prescription [67] for incorporating dark energy perturbations.

• •

iv) The $w_{0}w_{a}\textrm{CDM}+\sum m_{\nu}$ model: In this case we incorporate a dynamically varying DE EoS, i.e. $w(z)$ also varies with time. We parametrize the EoS with the $w_{0}-w_{a}$ approach. This is the Chevallier-Polarski-Linder (CPL) parametrization [68, 69]:

[TABLE]

We hereafter may simply refer to this model as DDE.

• •

v) The $w_{0}w_{a}\textrm{CDM}+\sum m_{\nu}$ model with $w(z)\geq-1$: This model has the same parametrization as the previous model, but the dynamical dark energy is forced to stay in the non-phantom range, i.e. $w(z)\geq-1$. This is achieved by noticing that $w(z)$ in eq. (2.2) is a monotonic function, that at present day $w(z)=w_{0}$ (since $a=1$ is the present value of the scale factor by convention), and that at the very early universe ($a\rightarrow 0$) we had $w(z)\rightarrow w_{0}+w_{a}$. It thus suffices to have the following hard priors applied on the CPL parameters to keep them from crossing the phantom barrier [23, 41, 61]:

[TABLE]

We hereafter may simply refer to this model as NPDDE.

• •

vi) The $\Lambda\textrm{CDM}+\sum m_{\nu}+A_{\textrm{Lens}}$ model: In this extended model we include $A_{\textrm{Lens}}$, which is the scaling of the lensing amplitude. In a particular model, the theoretical prediction for the gravitational potential (which generates the weak lensing of the CMB) corresponds to $A_{\textrm{Lens}}=1$. When $A_{\textrm{Lens}}$ is varied, the weak lensing is uncoupled from the primary anisotropies which cause it, and then scaled by the value of $A_{\textrm{Lens}}$ [70]. $A_{\textrm{Lens}}$ acts as a consistency check parameter. In a particular model, if the data prefers $A_{\textrm{Lens}}>1$, it means it prefers more smoothing of its acoustic peaks in the power spectra (typically caused by lensing) than what theoretically should be. The physical reason for this extra smoothing (if it can’t be accounted for statistical fluctuation in the data) may not be extra lensing but may be any new effect that mimics lensing [24]. There is a well-known $A_{\textrm{lens}}$ tension in the Planck high-$l$ data, but not in its measurements of lensing. In the $\Lambda\textrm{CDM}+A_{\rm Lens}$ model, with Planck 2018 TT+TE+EE+lowE data, one finds that the constraint $A_{\textrm{Lens}}=1.180\pm 0.065$ (68%) [24] is more than 2$\sigma$ away from $A_{\rm Lens}=1$. Solving the $A_{\textrm{Lens}}$ anomaly with a statistical significance of 5$\sigma$ or more requires more precise measurement of the E and B mode polarization of the CMB photons in small scales, which will be possible with future CMB experiments like LiteBIRD [71] (space based); and stage III ground based experiments, provided the accurate measurement of $\tau$ from Planck is included. Optimistic Stage IV CMB experiments can even provide a resolution at 10$\sigma$, also providing significant information on scale dependance of $A_{\textrm{Lens}}$, if any. See [72] for details.

• •

vii) The $\Lambda\textrm{CDM}+\sum m_{\nu}+\Omega_{k}$ model: Here we go away from a flat universe to the one which can be curved. The curvature of the universe is parametrized by $\Omega_{k}$, which is called the curvature energy density, and we allow it to vary freely in this model.

We use the publicly available Markov Chain Monte-Carlo package CosmoMC [73] (which uses the Boltzmann solver CAMB [74]) to perform a Bayesian analysis of cosmological datasets and derive constraints on $\sum m_{\nu}$ and other cosmological parameters. We use the Gelman and Rubin statistics [75] to estimate the convergence of the chains. All our chains had achieved the convergence criterion of $R-1<0.01$. We use flat priors on all the the parameters that are varied in a particular model. The priors are listed in table 1.

We also calculate Bayesian evidences for the purpose of Bayesian model comparison between the NH and IH cases in the above cosmological models using the publicly available package CosmoChord [76], an extension to CosmoMC using PolyChord [77, 78] for nested sampling. We used sufficiently broad priors on the cosmological parameters and 2000 live points for each CosmoChord run to keep the error to the estimated evidence small. The methodology to quantify evidence against IH closely follows [44]. The results are given in section 4.

2.2 Datasets

CMB: Planck 2018. We use the high-$l$ (30 $\leq$ $l$ $\leq$ 2508) and low-$l$ (2 $\leq$ $l$ $\leq$ 29) CMB TT likelihood along with high-$l$ E mode polarization and temperature-polarisation cross correlation likelihood from the recent Planck 2018 public release [79] and we call this combination “TTTEEE”. We use the Planck low-$l$ E mode polarization data, and in the text we mention it as “lowE.”

CMB: Planck 2018 lensing. While the CMB anisotropy power spectra is determined from 2-point correlation functions (TT, TE, EE), the power spectra of the lensing potential is proportional to the 4-point correlation functions such as TTTT, TTEB and so on [24]. We use it in our analyses as an additional CMB probe to ascertain neutrino physics properties, as lensing of CMB photons is produced by the gravitational potential of large scale structure which in turn is affected greatly by the free streaming massive neutrinos. We will refer to this dataset simply as “lensing” from now on.

CMB: BICEP2/Keck array data. While running an MCMC analysis for the $\Lambda\textrm{CDM}+\sum m_{\nu}+r$ model, we also use the latest publicly available data (taken up to and including 2015) on the CMB BB spectra from the BICEP2/Keck collaboration (spanning the range: $20<l<330$) [80]. This dataset is referred to as “BK15” in the paper.

Baryon acoustic oscillations (BAO) Measurements. In this paper we use the latest measurements of the BAO signal from different galaxy surveys: SDSS-III BOSS DR12 (galaxy samples at the effective redshifts of $z_{\textrm{eff}}=$ 0.38, 0.51 and 0.61) [81], the DR7 Main Galaxy Sample (MGS) at $z_{\textrm{eff}}=0.15$ [82], and the Six-degree-Field Galaxy Survey (6dFGS) survey at $z_{\textrm{eff}}=0.106$ [83]). We simply name this combined dataset as “BAO”.

All the CMB and BAO data together constitute our “base” dataset:

Base $\equiv$ Planck 2018 TTTEEE + lowE + lensing + BAO.

For the $\Lambda\textrm{CDM}+\sum m_{\nu}+r$ model, the “base” dataset shall also include BK15.

Supernovae luminosity distance measurements. We use the most recent Supernovae Type-Ia (SNe Ia) luminosity distance measurements from the Pantheon Sample [84], which consists of distance information of 1048 SNe Ia ($0.01<z<2.3$), largest till date. Out of the 1048, 279 are from the Pan-STARRS1 (PS1) Medium Deep Survey ($0.03<z<0.68$) and rest of them from SDSS, SNLS, various low-$z$ and HST samples. We call this dataset “SNe”.

To understand how the neutrino mass affects the cosmological observables measured by the above datasets and how we arrive at strong constraints on the same using these datasets, we refer the reader to [22].

3 Results

In this section we provide the results of our analyses on the bounds on neutrino masses considering the three different hierarchies (degenerate, normal, and inverted). In section 3.1 we discuss the results in the minimal $\Lambda\textrm{CDM}+\sum m_{\nu}$ model. The results in the extended models is discussed in section 3.2. Details about models and datasets are given in section 2.1 and 2.2 respectively. All the marginalized limits are given at 68% C.L. (1$\sigma$), whereas cases where only upper or lower bounds are available, bounds are given at 95% C.L. ($2\sigma$). The main results are contained in tables 2-4.

3.1 Results in the minimal $\Lambda\textrm{CDM}+\sum m_{\nu}$ model

In this subsection we provide results of our analyses in the $\Lambda\textrm{CDM}+\sum m_{\nu}$ model, considering three different hierarchies (degenerate, normal, and inverted) for two dataset combinations, namely Base and Base+SNe, where Base $\equiv$ Planck 2018 TTTEEE+lowE+lensing+BAO. The main results are contained in table 2 which provides confidence limits on cosmological parameters for the Base and Base+SNe combinations.

In figure 1 we depict the 1-D posterior distributions of $m_{0}$ (mass of the lightest neutrino in a hierarchy) and $\sum m_{\nu}$ in the $\Lambda\textrm{CDM}+\sum m_{\nu}$ model for Base and Base+SNe datasets considering different hierarchies. With the Base data recover the following 95% bound on the mass sum: $\sum m_{\nu}<0.12$ eV which is same as the bound of $\sum m_{\nu}<0.12$ eV quoted by the Planck 2018 collaboration [80] with the same data. The Planck 2015 result with similar data was $\sum m_{\nu}<0.17$ eV (95%, Planck 2015 TT,TE,EE+lowP+BAO). The main reason this improvement happens with Planck 2018 is because of improved measurement of $\tau$. Planck 2015 TT,TE,EE+lowP likelihoods produced a bound of $\tau=0.079\pm 0.017$ (68%) in the base $\Lambda\textrm{CDM}$ model [85], which improved to $\tau=0.0544^{+0.0070}_{-0.0081}$ (68%) with Planck 2018 TT,TE,EE+lowE [24]. The parameters $\tau$ and $\sum m_{\nu}$ are strongly correlated in the Planck TT data and high-$l$ polarization data [23], and this degeneracy can be broken through better measurement of $\tau$ from the low-$l$ polarization data. Since an increase in $\sum m_{\nu}$ leads to suppression of the matter power spectrum which in turn leads to less gravitational lensing of the CMB photons, we see an increase in $\sum m_{\nu}$ leads to a decrease in the smearing of the CMB acoustic peaks (smearing happens due to lensing). This effect due to increasing neutrino masses, however, can be countered with increasing $\tau$, which exponentially suppresses power in the CMB anisotropy spectra (given other parameters are kept fixed). However, the effect of reionization is also found in the low-$l$ polarization data (TE, EE, BB) in the form of a “reionization bump” whose amplitude is proportional to $\tau^{2}$ in the EE and BB spectra, and to $\tau$ in the TE spectra, and this change can’t be compensated by varying other parameters in the model [86]. Hence the improved removal of systematics in the low-$l$ polarization data with Planck 2018 helps in breaking the degeneracy between $\tau$ and $\sum m_{\nu}$, and this leads to stronger upper bounds on $\sum m_{\nu}$. However, not only $\tau$, Planck 2018 also largely corrects various systematic effects previously present in the high-$l$ polarization spectra of Planck 2015, and that also contributes towards obtaining a better bound on the $\sum m_{\nu}$.

In our Base dataset combination, BAO is another effective tool in constraining neutrino masses. BAO data is useful in breaking the strong anti-correlation between the Hubble constant $H_{0}$ and $\sum m_{\nu}$ present in the Planck data. The comoving distance to the last scattering surface in a flat $\Lambda\textrm{CDM}+\sum m_{\nu}$ universe is defined as: $\chi(z_{dec})=\int^{z_{dec}}_{0}dz/H(z)$, where $z_{dec}$ is the redshift of photon decoupling, and $H(z)=\sqrt{\omega_{\gamma}(1+z)^{4}+(\omega_{c}+\omega_{b})(1+z)^{3}+\omega_{\Lambda}+\rho_{\nu}(z)h^{2}/\rho_{cr,0}}$ (note: $\omega_{i}=\Omega_{i}h^{2}$, and $i\equiv\gamma,c,b,\Lambda$ with $\gamma\equiv$ photons, $c\equiv$ CDM, $b\equiv$ baryons, $\Lambda\equiv$ cosmological constant; $\rho_{\nu}(z)$ is the neutrino energy density at a redshift $z$, $\rho_{cr,0}=3H_{0}^{2}/8\pi G$ is the critical density today). In a flat universe $\Omega_{\Lambda}=1-\Omega_{\gamma}-(\Omega_{c}+\Omega_{b})-\Omega_{\nu}$ and at late times $\rho_{\nu}(z)h^{2}/\rho_{cr,0}=\Omega_{\nu}h^{2}\propto\sum m_{\nu}$. Now, given that early universe physics remains unchanged, $\chi(z_{dec})$ is very well constrained through $\Theta_{s}$ which is the most well constrained parameter by the Planck CMB anisotropies. On the other hand, $\Omega_{\gamma}$ and $(\omega_{c}+\omega_{b})$ are also well constrained by the data. Thus any change in $\chi(z_{dec})$ due to an increase in $H_{0}$ (or $h=H_{0}/100$ km/sec/Mpc) has to be compensated by a decrease in $\sum m_{\nu}$ and vice versa, and hence there is a large anti-correlation between $H_{0}$ and $\sum m_{\nu}$. Thus in the $\Lambda\textrm{CDM}+\sum m_{\nu}$ model, lower values of $H_{0}$ correspond to higher values of $\sum m_{\nu}$ and higher values of $H_{0}$ correspond to lower $\sum m_{\nu}$. BAO data breaks this degeneracy partially by rejecting the low $H_{0}$ values preferred by Planck data, well studied in previous literature [87, 25, 23]. For instance, in $\Lambda\textrm{CDM}+\sum m_{\nu}$, Planck 2015 TT,TE,EE+lowP prefers a value of $H_{0}=66.17^{+1.96}_{-0.81}$ km/sec/Mpc, whereas TT,TE,EE+lowP+BAO prefers $H_{0}=67.67^{+0.54}_{-0.51}$ km/sec/Mpc [23].

Apart from CMB power spectra and BAO, our Base dataset also contains Planck 2018 lensing likelihoods, which as per Planck 2018 collaboration [24] prefers a slightly increased lensing power spectrum amplitude compared to Planck 2015, and thus leads to a slightly tighter neutrino mass constraints, contrary to Planck 2015 lensing likelihoods which used to relax the constraints. In our case, without the Planck 2018 lensing likelihoods, we obtained a bound of $\sum m_{\nu}<0.13$ eV (95%, Planck 2018 TT,TE,EE+lowE+BAO), i.e., Planck 2018 lensing has only a small effect in this $\Lambda\textrm{CDM}+\sum m_{\nu}$ model when used along with Planck 2018 CMB anisotropies, and BAO.

We see that the bounds on $\sum m_{\nu}$ do differ significantly across the different hierarchies. We see that with the Base dataset we get the following bounds: $\sum m_{\nu}<0.15$ eV (NH), $\sum m_{\nu}<0.17$ eV (IH) at 95% C.L. As previously stated, we use the mass of the lightest neutrino, $m_{0}$ as the varying parameter and then use the mass-squared splittings given in eq. 1.1 to determine the other masses in a particular hierarchy, and this implicitly puts a prior on the $\sum m_{\nu}$, i.e., $\sum m_{\nu}\geq 0.06$ eV for normal hierarchy, and $\sum m_{\nu}\geq 0.10$ eV for inverted hierarchy. The reason the bound on $\sum m_{\nu}$ differs significantly across the three hierarchies is possibly the priors imposed on $\sum m_{\nu}$ (see figure 1 for the visualization of the same). We cross-check this by running MCMC chains with the degenerate hierarchy, but with the priors: i) $\sum m_{\nu}\geq 0.06$ eV, and ii) $\sum m_{\nu}\geq 0.10$ eV. We found a 95% upper bound of $\sum m_{\nu}<0.15$ eV in the first case, and in the second case it was $\sum m_{\nu}<0.18$ eV. It seems evident that the priors do play an important role in relaxing the bounds. However it is to be noted that the method we have used in this paper for the normal and inverted hierarchy (i.e. with lightest mass $m_{0}$ and the mass-squared splittings from oscillations experiments) produces bounds which are slightly stronger than the degenerate case with the priors $\sum m_{\nu}\geq 0.06$ eV and $\sum m_{\nu}\geq 0.10$ eV respectively, i.e. the two methods are not completely equivalent. The difference is clearer with Planck 2018 TT,TE,EE+lowE data, with which the NH case leads to a 95% bound of $\sum m_{\nu}<0.29$ eV, whereas the DH case with the $\sum m_{\nu}\geq 0.06$ eV prior gives us a bound of $\sum m_{\nu}<0.32$ eV. For the same data, the IH case produces a 95% bound of $\sum m_{\nu}<0.33$ eV, whereas the DH case with the $\sum m_{\nu}\geq 0.10$ eV prior obtains a bound of $\sum m_{\nu}<0.35$ eV. As discussed in the introduction, in case of NH and IH, with the lightest mass ($m_{0}$) parametrization and a uniform prior on $m_{0}$, the implicit prior on $\sum m_{\nu}$ is not completely flat or uniform (i.e. rises at the low values of $m_{0}$). This non-flat prior causes the low values of $\sum m_{\nu}$ to be more favoured compared to the case of a flat prior on $\sum m_{\nu}$, and hence, this can explain the slightly tighter bounds on $\sum m_{\nu}$ with the $m_{0}$ parametrization.

From table 2 we can see an important trend: in the $\Lambda\textrm{CDM}+\sum m_{\nu}$ model with Base dataset, as we go from degenerate to normal to inverted hierarchy there appears a decrease in the preferred values of $H_{0}$ and $\sigma_{8}$. This is also visualized in figure 2. The reason for the decrease in $H_{0}$ is simply the anti-correlation between $H_{0}$ and $\sum m_{\nu}$ as described above, and the fact that the normal hierarchy prefers neutrino masses which are larger than the degenerate case, and the inverted hierarchy prefers masses which are larger than the normal hierarchy. There is, again, a strong anti-correlation present between $\sigma_{8}$ and $\sum m_{\nu}$, since $\sigma_{8}$ is the amplitude of the linear matter power spectrum at a length scale of $8h^{-1}$MPc, and increasing $\sum m_{\nu}$ increases the neutrino energy density and that increases the suppression of the small scale matter power-spectrum leading to a lower $\sigma_{8}$. Thus, from degenerate to normal to inverted hierarchy, the $H_{0}$-tension between Planck and local measurements increases, whereas the $\sigma_{8}$-tension [88] between Planck and cosmic shear measurements decreases. This trend remains true for the other datasets we have considered in this model (as can be seen from table 2). From figure 2 we also observe that the parameter $S_{8}=\sigma_{8}\sqrt{\Omega_{m}/0.3}$ also goes downward, but not as strongly as $\sigma_{8}$, since increasing neutrino masses also cause $\Omega_{m}$ to increase which compensates for the decreased $\sigma_{8}$, i.e. $S_{8}$ is defined such that its correlation with $\sum m_{\nu}$ is small in this model. Again, while $H_{0}$ and $\sigma_{8}$ both are strongly anti-correlated with $\sum m_{\nu}$, the magnitude of the correlation coefficients ($R_{ij}=C_{ij}/\sqrt{C_{ii}C_{jj}}$ where $i$,$j$ are the two parameters, and $C$ is the covariance matrix of the cosmological parameters) decrease from degenerate to normal to inverted hierarchy because of the implicit priors on $\sum m_{\nu}$, i.e the priors help in partially breaking the degeneracies. We find that $R_{H_{0},\Sigma m_{\nu}}=-$0.55 (DH), $-$0.47 (NH), $-$0.43 (IH); and $R_{\sigma_{8},\Sigma m_{\nu}}=-$0.76 (DH), $-$0.72 (NH), $-$0.67 (IH).

Apart from the results with the Base dataset, in table 2, we provide results with the Base+SNe dataset also. In the absence of a varying dark energy equation of state, the SNe data is able to constrain $\Omega_{m}$ effectively [84], while the Planck CMB data can constrain $\Omega_{m}h^{2}$ well. Thus together they can effectively constrain $H_{0}$, and as found out in [23], the Planck+SNe combination actually prefers $H_{0}$ values which are higher than Planck alone, and thus the SNe data can help in breaking the degeneracy between $H_{0}$ and $\sum m_{\nu}$ partially. BAO data is however much more efficient in breaking the degeneracy and with Base+SNe (note that Base already contains BAO). In the DH case, we find the following 95% bounds on the neutrino mass sum in this $\Lambda\textrm{CDM}+\sum m_{\nu}$ model: $\sum m_{\nu}<0.11$ eV, which is only slightly stronger than the bound without the SNe data.

Akaike information criterion (AIC). To compare the goodness of fit of different models to the same data a popular method is to compute the Akaike information criterion (AIC) [89] for each model, defined as :

[TABLE]

where $k$ is the number of parameters in the model. Comparison with another model is done by computing the difference: $\Delta\textrm{AIC}=\Delta\chi^{2}_{\textrm{best-fit}}+2\Delta k$. Between two models the model with a lower AIC is considered the preferred model. The $2\Delta k$ term penalizes the model with greater number of parameters as it is usually able to provide better fit because of the larger parameter space being available to it. In this work however, we are interested in the comparison of the quality of fits between the neutrino mass hierarchies and since we have the same number of parameters in all the different cases, $2\Delta k=0$, and $\Delta\textrm{AIC}=\Delta\chi^{2}_{\textrm{best-fit}}$. We have provided the $\Delta\textrm{AIC}=\chi^{2}_{\textrm{best-fit}}-\chi^{2}_{\textrm{IH,best-fit}}$ values for each neutrino hierarchy in table 2. We find that cosmological datasets slightly prefer the normal hierarchy over the inverted one. In this $\Lambda\textrm{CDM}+\sum m_{\nu}$ model, with the Base data, we find that the NH is preferred to the IH by a $\Delta\textrm{AIC}=-0.95$, which is a very mild result.

This showcases that the current cosmological data is not sensitive enough to differentiate between the two hierarchies, and is completely consistent with the findings in e.g. [44, 53] which both find that a formal sensitivity to $\sum m_{\nu}$ in the 0.01-0.02 eV range is required to guarantee a distinction between the two hierarchies in case the true value of $\sum m_{\nu}$ is sufficiently less than 0.1 eV (i.e. if the neutrino masses follow NH). This is not possible using current data, but will become possible in the near future using high precision data from e.g. EUCLID (see e.g. [90, 91, 92, 93] for discussions on this topic). However, if the true value of $\sum m_{\nu}>0.1$ eV, we won’t be able to tell one hierarchy from the other even with the said sensitivity. With the other dataset combinations also the $\Delta\textrm{AIC}$ remains mild between NH and IH. The (unphysical) degenerate hierarchy on the other hand, produces a better fit to the data compared to both NH and IH in case of all the dataset combinations we have studied.

3.2 Results in the extended models

In this section we present the results in the extended cosmologies. Table 3 contains the constraints on selected cosmological parameters in the extended models with Base and Base+SNe datasets. Marginalized constraints are given at 1$\sigma$ whereas upper or lower bounds are given at 2$\sigma$. Details about models and datasets are given in section 2. Note that the Base data also contains BK15 when we consider the $\Lambda\textrm{CDM}+\sum m_{\nu}+r$ model.

• •

i) Results in the $\Lambda\textrm{CDM}+\sum m_{\nu}+r$ model: In this model we allow the tensor perturbations to vary as well, along with scalar ones. The CMB B-mode polarization comes from two sources: i) primordial gravitational waves, ii) gravitational weak lensing. With the Base data (which also includes BK15) we get the following 95% bounds: $\sum m_{\nu}<0.11$ eV (DH), $\sum m_{\nu}<0.14$ eV (NH), and $\sum m_{\nu}<0.17$ eV (IH). We do not see any big changes on the mass bounds compared to $\Lambda\textrm{CDM}+\sum m_{\nu}$. The bounds are, however, slightly tighter in the $\Lambda\textrm{CDM}+\sum m_{\nu}+r$ model with the BK15 data than without, most likely due to the additional lensing information encoded in the B mode data BK15, since $r$ and $\sum m_{\nu}$ are not correlated, as can be seen from the left panel of figure 3. The effect of the neutrino mass bounds getting slightly tighter with B mode data was previously noticed in [23, 41, 94] with BK14 data [95], predecessor of BK15.

• •

ii) Results in the $w\textrm{CDM}+\sum m_{\nu}$ model: In this model we go away from the cosmological constant, and consider the DE EoS $w$ as a varying parameter. A well-known degeneracy exists between $w$ and $\sum m_{\nu}$ [54], through their mutual degeneracy with $\Omega_{m}$, which considerably relaxes the bounds on $\sum m_{\nu}$ in this model compared to $\Lambda\textrm{CDM}+\sum m_{\nu}$. The visualization of the anti-correlation between the two parameters is available in the right panel of figure 3 for the Base and Base+SNe data. We find a correlation coefficient of $R_{w,\sum m_{\nu}}=-$0.57 (DH), $-$0.52 (NH), $-$0.47 (IH) with the Base data, whereas $R_{w,\sum m_{\nu}}=-$0.45 (DH), $-$0.35 (NH), $-$0.30 (IH) with Base+SNe. The SNe data when used with the Base data, reduces the magnitude of the anti-correlation with better constraints on $w$. This happens, since in the SNe data $w$ and $\Omega_{m}$ are strongly anti-correlated, whereas in the CMB data $w$ and $\Omega_{m}$ are strongly correlated (for instance see figure 20 of [84]). So including the SNe data with Base leads to a large decrease in the $w$-$\Omega_{m}$ correlation, which in turn leads to a decrease in magnitude of $R_{w,\sum m_{\nu}}$. Also from figure 3, we see that due to the anti-correlation lower values of $w$ prefer higher $\sum m_{\nu}$ and higher values of $w$ prefer lower $\sum m_{\nu}$. The Base+SNe combination constraints $w$ better than Base such that the lowest values of $w$ allowed by Base data is rejected, and this leads to stronger bounds on $\sum m_{\nu}$ with Base+SNe, as can be seen from table 3. The SNe data also rejects high $w$ values, but those regions prefer low $\sum m_{\nu}$ values and hence rejecting the high $w$ region does not help in strengthening the upper bound on $\sum m_{\nu}$.

There is also a strong degeneracy between $w$ and $H_{0}$, that exists since lower values of $w$ correspond to higher present day expansion rate and hence higher $H_{0}$ values. This happens since changing $w$ can change the comoving distance to the last scattering surface $\chi(z_{dec})=\int^{z_{dec}}_{0}dz/H(z)$, which is well-constrained by the CMB data, and thus any change in $\chi(z_{dec})$ needs to be compensated with another parameter. We have,

[TABLE]

Here $\Omega_{DE}(z)=\Omega_{DE}(0)(1+z)^{3(1+w)}$ is the energy density of dark energy with EoS $w$, and at late times when dark energy is a dominant component, decreasing the value of $w$ leads to a decrease in $H(z)$, which can be compensated by increasing either $h$ (or $H_{0}$) or $\sum m_{\nu}$ or both. This is why in the $w\textrm{CDM}+\sum m_{\nu}$ model, with Base data, we find that not only $w$ and $H_{0}$ are anti-correlated ($R_{w,H_{0}}=-0.91$) (DH), the correlation between $H_{0}$ and $\sum m_{\nu}$ is inverted from what it was in $\Lambda\textrm{CDM}+\sum m_{\nu}$, i.e. here $H_{0}$ and $\sum m_{\nu}$ are positively correlated with $R_{H_{0},\sum m_{\nu}}=+0.30$ (DH). The SNe data, when combined with Base, constrains $w$ and $H_{0}$ in a much better way and breaks the degeneracy present between $H_{0}$ and $\sum m_{\nu}$, and we find that $R_{H_{0},\sum m_{\nu}}=-0.013$ (DH) with Base+SNe. This is why, in table 3 with the Base data, as we go from DH to NH to IH, the preference for higher and higher neutrino masses and the positive correlation between $H_{0}$ and $\sum m_{\nu}$ leads to slight increase in the preferred values of $H_{0}$, but with Base+SNe, the almost negligible correlation explains the lack of any big change in the preferred $H_{0}$ values across the different hierarchies. On the other hand, the strong degeneracy between $w$ and $H_{0}$ survives, but becomes slightly weaker at $R_{w,H_{0}}=-0.76$ (DH) with Base+SNe. This leads to the $w\textrm{CDM}+\sum m_{\nu}$ model predominantly preferring values of $w<-1$, i.e. the phantom DE region.

• •

iii) Results in the $w_{0}w_{a}\textrm{CDM}+\sum m_{\nu}$ model: In this model the DE EoS $w(z)=w_{0}+w_{a}\frac{z}{1+z}$ is dynamical in nature, but remains close to the value of $w_{0}+w_{a}$ for high redshift values, and only changes significantly during the late times. Like the previous model, there is degeneracy between $w(z)$ and $H_{0}$, and $w(z)$ and $\sum m_{\nu}$. However now we have two parameters $w_{0}$ and $w_{a}$ instead of $w$, with

[TABLE]

Any change in $\chi(z_{dec})$ due to $H_{0}$ can be readily compensated with $w_{0}$ and $w_{a}$, and hence, in this model the correlation between $H_{0}$ and $\sum m_{\nu}$ is very small even with Base data ($R_{H_{0},\sum m_{\nu}}=+0.08$ (DH)). Again, $w_{0}$ and $w_{a}$ are anti-correlated between themselves, since the change in $\chi(z_{dec})$ due to an increase in $w_{0}$ can be countered with a decrease in $w_{a}$. The anti-correlation can be seen clearly in figure 4. As it can be seen from the figure, neither the Base data nor Base+SNe combination prefers the non-phantom DE region. The data prefers regions where the DE is currently phantom or has been phantom in the past. Since the phantom DE region, $w(z)<-1$ prefers neutrino masses which are larger, the mass bounds are essentially very relaxed in this model, as can be seen from table 3. We have $\sum m_{\nu}<0.25$ eV (DH), 0.26 eV (NH), 0.28 eV (IH) with Base data, whereas with Base+SNe, we get $\sum m_{\nu}<0.27$ eV (DH), 0.28 eV (NH), 0.29 eV (IH). The SNe data produces better constraints on the DE parameters. In figure 4, another important thing to notice is that as we go from DH to NH to IH, the 2D contours shift away from the non-phantom DE region due to the preference for higher and higher masses, and in the IH case, Base data rejects the non-phantom region at more than 2$\sigma$.

We can understand how the individual DE parameters $w_{0}$ and $w_{a}$ are affected by $\sum m_{\nu}$ by looking at figure 5. For the Base data, the relevant correlation coefficients are: $R_{w_{0},\sum m_{\nu}}=-0.004$ (DH), $-0.018$ (NH), $-0.056$ (IH), and $R_{w_{a},\sum m_{\nu}}=-0.27$ (DH), $-0.21$ (NH), $-0.17$ (IH). It is clear that $w_{a}$ is highly anti-correlated with $\sum m_{\nu}$. Since $w(z)=w_{0}+w_{a}z/(1+z)$, given a particular value of $w_{0}$, it is $w_{a}$ which determines the late time dynamics of the dark energy equation of state, and thus is more affected by $\sum m_{\nu}$ than $w_{0}$. Increasing $w_{a}$ has the effect of increasing $\Omega_{\rm DE}(z)$ at late times ($z\ll 1$), as can be seen from eq. 3.3. This can be partially compensated by decreasing $\sum m_{\nu}$ to decrease $\rho_{\nu}(z)$ at late times to keep $\chi(z_{dec})$ fixed. This can be also compensated by decreasing $w_{0}$, and thus $w_{0}$ and $w_{a}$ are anti-correlated, whereas $w_{0}$ and $\sum m_{\nu}$ are weakly correlated in the Base data.

• •

iv) Results in the $w_{0}w_{a}\textrm{CDM}+\sum m_{\nu}$ model with $w(z)\geq-1$ : This is essentially the same model as in the previous case, but the parameter space is restricted to the non-phantom range only, i.e. $w(z)\geq-1$. This parameter space corresponds to dark energy field theories modelled with a single scalar field, like quintessence [96], which cannot cross the phantom barrier (the $w(z)=-1$ line). Due to the degeneracy between $w$ and $\sum m_{\nu}$ (which we have discussed for the last two models) the cosmological data actually prefers smaller and smaller neutrino masses as we go deeper in the non-phantom region in the $w_{0}-w_{a}$ parameter space. Thus the bounds on $\sum m_{\nu}$ obtained in this model are even tighter than the $\Lambda\textrm{CDM}+\sum m_{\nu}$ model. This was first noticed in [61], and subsequently in [23, 41]. From table 3, with the base data we find $\sum m_{\nu}<0.10$ eV (DH), 0.13 eV (NH), 0.16 eV (IH). With the Base+SNe combination, the bounds stand at $\sum m_{\nu}<0.09$ eV (DH), 0.13 eV (NH), 0.16 eV (IH). These are the tightest bounds on $\sum m_{\nu}$ reported in this paper for the datasets considered.

• •

v) Results in the $\Lambda\textrm{CDM}+\sum m_{\nu}+A_{\textrm{Lens}}$ model : As mentioned in section 2.1, the $A_{\textrm{Lens}}$ parameter is used to artificially scale the lensing amplitude predicted by the underlying theoretical model. There is a well-known $A_{\textrm{Lens}}$-issue in the CMB anisotropy data which corresponds to preference for $A_{\textrm{Lens}}$ values which are more than 2$\sigma$ away from the theoretical expectation of $A_{\textrm{Lens}}=1$. For instance Planck 2018 TT,TE,EE+lowE data yields an $A_{\textrm{Lens}}=1.180\pm 0.065$ (68%) in the $\Lambda\textrm{CDM}+A_{\textrm{Lens}}$ model, showing a 2.8$\sigma$ tension [24]. The $\Lambda\textrm{CDM}+A_{\textrm{Lens}}$ also provides significantly better fit to the Planck power spectrum data compare to $\Lambda\textrm{CDM}$. Planck power spectrum constrains $A_{\textrm{Lens}}$ by measuring the smoothing of the CMB acoustic peaks due to the gravitational lensing of the CMB photons from large scale structures. Planck lensing data, however constrains $A_{\textrm{Lens}}$ directly, and adding CMB lensing data to the anisotropy power spectrum usually brings back the $A_{\textrm{Lens}}$ values closer to the theoretically expected result [85]. We find it true in this model as well. With Base data, for the DH case, we have $A_{\textrm{Lens}}=1.100^{+0.100}_{-0.096}$ (95%), which is almost consistent with $A_{\textrm{Lens}}=1$ at 2$\sigma$.

In this $\Lambda\textrm{CDM}+\sum m_{\nu}+A_{\textrm{Lens}}$ model, $\sum m_{\nu}$ and $A_{\textrm{Lens}}$ are strongly correlated, since increasing $\sum m_{\nu}$ has an effect of reducing lensing induced smearing of the acoustic peaks since increasing $\sum m_{\nu}$ caused increasing suppression to the small scale matter power. We find that $R_{A_{\textrm{Lens}},\sum m_{\nu}}=+$0.68 (DH), $+$0.60 (NH), $+$0.55 (IH) with the Base data. The positive correlation between these two parameters can be seen in figure 6. This degeneracy causes the bounds on $\sum m_{\nu}$ to relax considerably compared to the $\Lambda\textrm{CDM}+\sum m_{\nu}$ model, as can be seen from table 3. Compared to a 95% bound of $\sum m_{\nu}<0.12$ eV (DH) with Base data in the minimal $\Lambda\textrm{CDM}+\sum m_{\nu}$ model, here we get a bound of $\sum m_{\nu}<0.29$ eV with the same data. Rather this $A_{\textrm{Lens}}$-issue actually explains why the neutrino mass bounds are so strong in the $\Lambda\textrm{CDM}+\sum m_{\nu}$ model, especially with the degenerate hierarchy, as a higher $\sum m_{\nu}$ would decrease lensing of the CMB photons. From DH to NH to IH, while using same data, the preference for higher $A_{\textrm{Lens}}$ values occurs with higher $\sum m_{\nu}$ values, and the discrepancy with $A_{\textrm{Lens}}=1$ increases.

In the absence of varying dark energy EoS, however, $H_{0}$ and $\sum m_{\nu}$ are strongly correlated in this model. The model also produces lower $\sigma_{8}$ and $S_{8}$ values compared to other models considered in the paper, leading to a decrease in the $\sigma_{8}$ tension present between Planck and cosmic shear experiments like CFHTLenS [97], KiDS-450 [98] etc.

• •

v) Results in the $\Lambda\textrm{CDM}+\sum m_{\nu}+\Omega_{\textrm{k}}$ model : In this model we consider the possibility of a non-flat background geometry of the universe. $\Omega_{\textrm{k}}$ and $\sum m_{\nu}$ are correlated positively and strongly with both Base and Base+SNe combinations, as can be seen from the right panel of figure 6. This causes the bounds on $\sum m_{\nu}$ to degrade compared to the $\Lambda\textrm{CDM}+\sum m_{\nu}$ model. With Base data we have: $\sum m_{\nu}<0.15$ eV (DH), $\sum m_{\nu}<0.17$ eV (DH), $\sum m_{\nu}<0.20$ eV (IH). The origin the degeneracy again lies in the tightly constrained angular diameter distance to the last scattering surface, i.e. $D_{A}(z_{dec})\equiv\textrm{sin}\left(\sqrt{K\chi(z_{\rm dec}}\right)/\sqrt{K}$ (where $K=-\Omega_{k}H_{0}^{2}$ is the curvature), leading to a three way geometric degeneracy between $H_{0}$, $\Omega_{k}$ and $\sum m_{\nu}$ [99] (note that when $\Omega_{k}=0$, $D_{A}(z)=\chi(z)$). In this model, the expression for $H(z)$ is given by,

[TABLE]

.

It is to be noted however, that the Planck data actually prefers values of $\Omega_{\textrm{k}}<0$ [24], since closed universe models produce more lensing amplitude compared to a flat universe. However, inclusion of lensing data usually brings the parameters back closer to a flat universe. BAO data helps in partially breaking the three way degeneracy by constraining $H_{0}$. We find that with Base data, in the DH case $\Omega_{\textrm{k}}=0.0004\pm 0.0020$, which is perfectly consistent with $\Omega_{\textrm{k}}=0$. The correlation coefficients are $R_{\Omega_{k},\sum m_{\nu}}=+0.41$, $R_{H_{0},\sum m_{\nu}}=-0.18$, and $R_{\Omega_{k},H_{0}}=+0.60$ for DH, i.e. the three parameters still remain considerably correlated with each other. From DH to NH to IH, preference for higher neutrino masses leads to preference for higher values of $\Omega_{k}$, although with the Base data $\Omega_{\textrm{k}}=0$ is included within 2$\sigma$ ranges of $\Omega_{k}$, for any of the hierarchies.

Akaike information criterion (AIC). Here we compare the goodness of fit of degenerate, normal and inverted hierarchy scenarios in the extended models that we have studied in this section. As in previous section, we use AIC as a measure of goodness of fit. AIC for the degenerate, normal, and inverted hierarchy are denoted as $\textrm{AIC}_{\textrm{DH}}$, $\textrm{AIC}_{\textrm{NH}}$, and $\textrm{AIC}_{\textrm{IH}}$ respectively. Since the DH, NH and IH cases both have the same number of parameters, $\Delta\textrm{AIC}=\textrm{AIC}-\textrm{AIC}_{\textrm{IH}}=\Delta\chi^{2}$, where $\Delta\chi^{2}$ is the $\chi^{2}$ difference between the $\chi^{2}$ value for a particular hierarchy (DH, NH, or IH) and the IH case, at the best-fit points. The $\Delta\textrm{AIC}$ values for each of the extended models has been listed in table 4, for the Base and Base+SNe datasets. We find that in most scenarios $\Delta\textrm{AIC}<0$, i.e. the fit due to DH/NH is better than IH, whereas in a few cases $\Delta\textrm{AIC}>0$. But in none of the models do we see any statistically significant difference between NH and IH.

4 Bayesian Model Comparison

While in the previous section we had chosen a frequentist measure (i.e. $\Delta\rm AIC$) to compare NH and IH, in this section we take a Bayesian approach. Our method closely follows [44], and that of [100, 101]. Consider that any two models, $M_{1}$ and $M_{2}$, with parameter vectors $\theta_{1}$ and $\theta_{2}$, are being tested against the same data $D$. Then, the posterior probability of a particular model $M_{i}$ ($i\equiv 1,2$) given the data $D$ is,

[TABLE]

where $P(D)$ is the prior probability for the data, and $\pi_{i}\equiv P(M_{i})$ is our prior degree of belief in model $M_{i}$. Since we are dealing with probabilities, if $M_{1}$ and $M_{2}$ are the only two possible models for a particular scenario(for instance, if $M_{1}$ and $M_{2}$ represent NH and IH and we know these are the only two possible mass orderings), we can write $\pi_{1}+\pi_{2}=1$. $P(D|M_{i})$ is the Bayesian evidence or the marginal likelihood, and is calculated by marginalizing product of the likelihood $P(D|\theta_{i},M_{i})$ and the parameter prior $P(\theta_{i}|M_{i})$ over all parameters:

[TABLE]

Again, as we can take the $p_{1}+p_{2}=1$, we can rewrite eq. 4.1 as,

[TABLE]

As indicated before, we shall treat $M_{1}$ and $M_{2}$ as NH and IH, and from now on we shall use the notation $p_{\rm NH}$ and $p_{\rm IH}$ for their posterior probabilities. One can now define the quantity Odds(NH:IH) as

[TABLE]

which is simply the ratio of posterior probabilities [100]. If no apriori knowledge of $\pi_{i}$ exists, it is customary to take all the $\pi_{i}$ to be equal, i.e. $\pi_{\rm NH}=\pi_{\rm IH}=0.5$. In that case, Odds(NH:IH) essentially becomes equal to the Bayes’ factor $B_{\rm NH/IH}$, which is defined as the ratio of evidences in the NH and IH case.

However, in case of neutrino mass hierarchy, one can assign $\pi_{i}$ from the Bayesian analysis of neutrino oscillations data. From the Bayesian analysis of neutrino oscillations data in [52], we see that $ln(B_{\rm NH/IH})=6.5$ (taking only the central value), which points to a posterior $\textrm{Odds(NH:IH)}=665:1$. If we take this as our prior belief before doing any analysis with cosmological data, we get $\pi_{\rm NH}=0.9985,\pi_{\rm IH}=0.0015$. We have used both this and the previous $\pi_{\rm NH}=\pi_{\rm IH}=0.5$ case while calculating the Odds(NH:IH).

Alternatively, following [44, 25], one can also consider the confidence limit at which IH is disfavored, as

[TABLE]

Odds(NH:IH) and CLIH and strength of evidence as per Jeffrey’s scale [102] has been provided in table 5, where we have used only the mean values of the estimated Bayesian evidence (by CosmoChord [76]) to compute the numerical quantities. It can be seen from table 5 that the analyses without the input from neutrino oscillations data on $\pi$ yield only a “Inconclusive” evidence for the normal neutrino mass hierarchy, which agrees with with our small $\Delta\rm AIC$ values from previous section. On the other hand, all the evidences are “Strong” with $\pi_{\rm NH}=0.9985,\pi_{\rm IH}=0.0015$, i.e. the preference for normal mass hierarchy is mostly driven by the oscillations data.

5 Discussion and conclusions

Presently the upper bounds on the sum of 3 active neutrino masses, $\sum m_{\nu}$ from analyses of cosmological data in the backdrop of $\Lambda\textrm{CDM}+\sum m_{\nu}$ model are bordering on the minimum sum of neutrino masses required by the inverted hierarchy, which is around 0.1 eV. However, these analyses are usually done with the assumption of 3 degenerate neutrino masses, but terrestrial neutrino oscillation experiments have confirmed that the three neutrino masses are not equal. Thus in this paper we update the bounds on $\sum m_{\nu}$ from latest publicly available cosmological data while explicitly considering particular neutrino mass hierarchies using the results on mass-squared splittings from a global analysis of neutrino oscillations data, NuFit 4.0 [9]. For implementing the normal and inverted hierarchy scenarios, we use the mass of the lightest mass eigenstate, denoted with $m_{0}$, as the varying parameter, and use the mass-squared splittings from NuFit 4.0 to determine the other masses in a particular hierarchy, and thus use the total mass sum, i.e. $\sum m_{\nu}$, as a derived parameter. This approach puts some implicit priors on the neutrino mass sum: $\sum m_{\nu}\geq 0.06$ eV for the normal hierarchy (NH) case, $\sum m_{\nu}\geq 0.10$ eV for the inverted hierarchy (IH) case.

In the minimal $\Lambda\textrm{CDM}+\sum m_{\nu}$ model with Planck 2018 TT,TE,EE, lowE, lensing, and the latest BAO data from various galaxy surveys, we find that at 95% C.L. $\sum m_{\nu}<0.12$ eV in the case of degenerate mass approximation. We call this dataset “Base.” This is similar to the bound of $\sum m_{\nu}<0.12$ eV quoted by the Planck 2018 collaboration using the same data, in the same model. We also find that in the same model we have $\sum m_{\nu}<0.15$ eV in case of NH and $\sum m_{\nu}<0.17$ eV in case of IH; i.e., the bounds vary significantly across the different mass orderings. The main reason for these differences in the bounds is the implicit priors on $\sum m_{\nu}$ when we assume a particular hierarchy. The case of degenerate neutrino masses with a prior $\sum m_{\nu}\geq 0.06$ eV or $\sum m_{\nu}\geq 0.10$ eV produces almost the same bounds as the NH or IH case with the lightest mass ($m_{0}$) parametrization we have used in the paper. The NH and IH case produce bounds which are only slightly stronger than the degenerate case with priors of $\sum m_{\nu}\geq 0.06$ eV and $\sum m_{\nu}\geq 0.10$ eV respectively, which can be attributed to the fact that with the lightest mass ($m_{0}$) parametrization, the implicit prior on $\sum m_{\nu}$ is not completely flat (i.e. rises at the low values of $m_{0}$). Also, we find that the normal hierarchy is very mildly preferred to the inverted: $\Delta\chi^{2}\equiv\chi^{2}_{\textrm{NH}}-\chi^{2}_{\textrm{IH}}=-0.95$ (best-fit). We also study this model against another dataset combination: Base+SNe. For Base and Base+SNe datasets, the $\chi^{2}$ differences between the NH and IH case remain statistically mild, but the bounds on $\sum m_{\nu}$ across the three different mass orderings vary considerably.

In this paper, we also provide bounds on $\sum m_{\nu}$ considering different hierarchies in various extended models: $\Lambda\textrm{CDM}+\sum m_{\nu}+r$, $w\textrm{CDM}+\sum m_{\nu}$, $w_{0}w_{a}\textrm{CDM}+\sum m_{\nu}$, $w_{0}w_{a}\textrm{CDM}+\sum m_{\nu}$ model with $w(z)\geq-1$, $\Lambda\textrm{CDM}+\sum m_{\nu}+A_{\textrm{Lens}}$, and $\Lambda\textrm{CDM}+\sum m_{\nu}+\Omega_{k}$. Here $r$ is the tensor-to-scalar ratio, $w$ is the dark energy equation of state (DE EoS) parameter, $w_{0}$ and $w_{a}$ again parametrize the DE EoS but in a dynamical manner: $w(z)=w_{0}+w_{a}z/(1+z)$ (CPL parametrization), $A_{\textrm{Lens}}$ is the parameter for artificial scaling of the lensing amplitude, and $\Omega_{k}$ is the curvature energy density. We used the Base and Base+SNe dataset combinations to constrain these models.

Consistent with other studies (see e.g. [103] for a very recent example) we found that in some cases the formal bound $\sum m_{\nu}$ could be loosened by up to factor of two. However, in none of the extended models could we find any statistically significant difference in the quality of fit to the data between NH and IH, i.e. the current cosmological data is not sufficiently strong to demarcate the two hierarchies. Apart from checking $\chi^{2}$ differences between NH and IH, we also confirmed the same using computation of Bayesian evidence. The results are provided in section 4. Starting with equal apriori beliefs on NH and IH scenarios, we find no conclusive evidence for the normal hierarchy over inverted hierarchy in any of the cosmological models with the Base dataset. This finding is consistent with previous work showing that if the actual value of $\sum m_{\nu}$ is sufficiently less than 0.1 eV (i.e. the neutrino masses follow NH) a formal sensitivity to $\sum m_{\nu}$ of 0.01-0.02 eV is required to guarantee a conclusive distinction between the two hierarchies [44, 53]. However, if the actual value of $\sum m_{\nu}>$ 0.1 eV, then we will not be able to determine the hierarchy, due to the fact that even the most optimistic future cosmological measurements won’t be able to differentiate among individual neutrino masses [64]. Hence, for a direct detection of the neutrino mass hierarchy, we have to depend on future terrestrial neutrino oscillation experiments like DUNE [104] and JUNO [105].

Acknowledgements

The authors thank the anonymous referee immensely, for the suggestions on improving the manuscript. for SRC thanks the computing facilities at HRI (http://www.hri.res.in/cluster/) and at Aarhus University (http://www.cscaa.dk/), and the warm hospitality at Aarhus University during the completion of the project. SRC would also like to thank the Department of Atomic Energy (DAE) Neutrino Project of HRI. STH is supported by a grant from the Villum Foundation.