Updated Bounds on Sum of Neutrino Masses in Various Cosmological Scenarios

# Updated Bounds on Sum of Neutrino Masses in Various Cosmological Scenarios

## Abstract

We present strong bounds on the sum of three active neutrino masses () using selected cosmological datasets and priors in various cosmological models. We use the following baseline datasets: Cosmic Microwave Background (CMB) temperature data from Planck 2015, Baryon Acoustic Oscillations measurements from SDSS-III BOSS DR12, the newly released Type Ia supernovae (SNe Ia) dataset from Pantheon Sample, and a prior on the optical depth to reionization from 2016 Planck Intermediate results. We constrain cosmological parameters with these datasets with a Bayesian analysis in the background of model with 3 massive active neutrinos. For this minimal model we find a upper bound of 0.152 eV at 95 C.L. Adding the high- polarization data from Planck strengthens this bound to 0.118 eV, which is very close to the minimum required mass of 0.1 eV for inverted hierarchy. This bound is reduced to 0.110 eV when we also vary r, the tensor to scalar ratio ( model), and add an additional dataset, BK14, the latest data released from the Bicep-Keck collaboration (which we add only when is varied). This bound is further reduced to 0.101 eV in a cosmology with non-phantom dynamical dark energy ( model with for all ). Considering the model and adding the BK14 data again, the bound can be even further reduced to 0.093 eV. For the model without any constraint on , the bounds however relax to 0.276 eV. Adding a prior on the Hubble constant ( km/sec/Mpc) from Hubble Space Telescope (HST), the above mentioned bounds further improve to 0.117 eV, 0.091 eV, 0.085 eV, 0.082 eV, 0.078 eV and 0.247 eV respectively. This substantial improvement is mostly driven by a more than 3 tension between Planck 2015 and HST measurements of and should be taken cautiously.

a,b,1]Shouvik Roy Choudhury,\noteCorresponding author. a,b,c]Sandhya Choubey \affiliation[a]Harish-Chandra Research Institute
Training School Complex, Anushaktinagar, Mumbai - 400094, India \affiliation[c]Department of Physics, School of Engineering Sciences, KTH Royal Institute of Technology

## 1 Introduction

Neutrino oscillation experiments have put the existence of neutrino mass on a solid footing [1, 2, 3, 4, 5, 6, 7, 8, 9]. There are three mass eigenstates (, , and ) which are quantum superpositions of the 3 flavour eigenstates (, , and ). The absolute neutrino mass scale, defined as the sum of the mass of the neutrino mass eigenstates, is the quantity,

 ∑mν≡m1+m2+m3, (1)

where is the mass of the neutrino mass eigenstate. Tightest bounds on come from cosmology.

Massive but light active neutrinos behave as radiation in the early universe, when their mass is much higher than the temperature. Their energy density adds a contribution to the total radiation energy density () of the universe, which is conventionally parametrized by , an effective number of species of neutrinos,

 ρr=π215[1+78(411)43Neff]T4γ, (2)

where is the temperature of the photon. In the event of Neutrinos decoupling instantaneously from the QED plasma, we would have gotten 3. But neutrinos decouple in the early universe over many Hubble times (temperature T 10 MeV - 0.1 MeV), which changes by a small amount, with the current best theoretical estimate being 3.046. This result is mostly due to (1) the continuation of decoupling of neutrinos during electron-positron annihilation and (2) QED plasma effects (see [10, 11, 12] for discussions on this topic.) As long as neutrinos can be considered as a radiation species during the photon decoupling (temperature T 0.2 eV), CMB data shows no evidence of bounds on not being compatible with its theoretically predicted value [13]. Any departure of from the theoretical prediction would be due to non-standard effects in the active neutrino sector or to the contribution of other relativistic species like a sterile neutrino. However in this work, we are only interested in bounds . We have not considered variation of . So in all our analyses we have fixed the value of to 3.046.

When neutrinos turn non-relativistic at late times, their energy density adds to the total matter density. Effect of massive neutrinos on cosmology has been widely studied in the literature [14, 15, 16, 17, 18, 19]. For masses much smaller than 1 eV, neutrinos are still very much relativistic at the time of photon decoupling, and their mass cannot affect the evolution of photon perturbations. Consequently, the effect can only appear on the background evolution, and secondary anisotropies like Integrated Sachs-Wolfe (ISW) Effect. These can be partially compensated by varying other free parameters of CDM model, and hence CMB anisotropy alone is not a very useful tool to obtain strong bounds on . As they become non-relativistic, sub-eV neutrinos affect late time evolution of matter perturbations considerably. Due to neutrino free streaming effects, increasing suppression of matter power spectrum in small scales happens with increasing neutrino density fraction with respect to the matter density [18]. Therefore, CMB data accompanied with datasets from Baryon Acoustic Oscillations (BAO) and Large Scale Structure (LSS) measurements can produce very strong upper bounds on .

Another neutrino property that future cosmological data might be able to determine is neutrino mass hierarchy. Different global fits [20, 21, 22, 23, 24] to neutrino oscillation experiment results have determined two mass squared splittings with considerable accuracy: eV and eV (1 uncertainties) [25]. These are the mass squared splittings dictating the solar and atmospheric transitions respectively. Since the sign of is not known, it can be either +ve, giving rise to normal hierarchy (NH) with , or -ve, giving rise to inverted hierarchy (IH) with . A lower bound can be put on for each of the two scenarios. For NH, it is around eV, for IH it is around eV. These can be calculated from the values of the mass splittings by noting that

 ∑mν=m0+√Δm221+m20+√|Δm31|2+m20                 \emph{(% NH)}, (3)

and

 ∑mν=m0+√|Δm31|2+m20+√|Δm31|2+Δm221+m20   \emph{(IH)}, (4)

where is the lightest neutrino mass. By convention, for NH, and for IH. Substituting gives the minimum sum of masses for each case.

It is to be noted that current cosmological measurements are primarily sensitive to the sum of the three masses, . If the three active neutrinos have different masses, then they will each become non-relativistic at different temperatures, and will produce different levels of suppression to the small scale matter power spectrum. Therefore, same total mass, but different mass splittings should result in slightly different matter power spectra. However current experiments do not have the sensitivity to distinguish these differences with reasonable significance [25]. For this reason, in this work, we present results which are obtained with the approximation of 3 degenerate neutrino masses (from now on DEG), i.e.,

 m1=m2=m3=∑mν3   (DEG). (5)

This approximation is predominant in literature in analyses where is varied. Planck data combined with others led to a bound of 0.23 eV at 95 C.L. (PlanckTT + lowP + lensing + BAO + JLA + H) [13]. Depending on the used data and variations in the analysis, different analyses [26, 27, 28, 29, 30, 31, 32] obtain 95 C.L. upper bounds from current data approaching the value of 0.1 eV, minimum mass required for IH. These results suggest IH is under pressure from cosmology. However, getting a 95 limit of less than the minimum required mass for IH does not rule out IH. Recent papers [25, 31] have suggested a rigorous but simple statistical method of computing the confidence level at which the hypothesis of IH can be rejected against NH using results from both cosmological and oscillations data. These recent analyses indicate that cosmology does slightly prefer NH compared to IH, but no statistically significant conclusion can be reached yet. In this paper, we, however, do not perform this statistical analysis and concentrate only on obtaining bounds on .

In this work we provide bounds on in the background of five different cosmological models: (1) (2) , where is the tensor to scalar ratio, (3) , where we assume Chevallier-Polarski-Linder (CPL) parametrization for dynamical dark energy instead of a simple cosmological constant, (4) with , i.e., we restrict the parameter space to exclude phantom dark energy (a recent study [33] also has explored this model with datasets different from ours), and (5) with , a model extended with both tensors and non-phantom dynamic dark energy.

We use combinations of the following recent datasets: (1) Cosmic Microwave Background (CMB) temperature, polarization and their cross-correlation data from Planck 2015; (2) Baryon Acoustic Oscillations measurements from SDSS-III BOSS DR12, MGS and 6dFGS; (3) the newly released Type Ia supernovae (SNe Ia) luminosity distance measurements from Pantheon Sample; (4) the latest data released from the BICEP/Keck Collaboration for the BB mode of the CMB spectrum; (5) and also local measurements of the Hubble parameter () from the Hubble Space Telescope; (6) the 2016 measurement of the optical depth to reionization () obtained from the analysis of the data from High Frequency Instrument of the Planck satellite; and (7) the galaxy cluster data from the observation of the Sunyaev-Zel’dovich (SZ) signature from thee 2500 square degree South Pole Telescope Sunyaev Zel’dovich (SPT-SZ) survey. We emphasize here that apart from Planck 2015 and the two redshift priors, the other datasets have not been studied widely in literature for obtaining bounds on and we are the first to use combinations of these above mentioned datasets and priors to obtain the very strong bounds presented in this paper, in the above mentioned cosmological models.

The two low redshift priors (on and ) are particularly important in constraining because of presence of significant degeneracy of these two parameters with [29] in the CMB data. Also to be noted that, future high-resolution CMB polarization measurements [34, 35, 36, 10] might be able to reconstruct the lensing potential very accurately and provide even stronger constraints on . However, we do not include the Planck lensing data, as the lensing potential measurements via reconstruction through the four-point functions of CMB data from Planck 2015 are in tension with the constraints obtained from CMB power spectra [13].

This paper is structured as follows: in Section 2 we describe our analysis method, the varying parameters of various cosmological models analyzed in this paper and the priors on the said parameters. We also briefly describe the Chevallier-Polarski-Linder (CPL) parametrization for dynamical dark energy. In Section 3, we briefly describe the various datasets we have used in this work. In Section 4 we present the results of our analysis. We conclude in Section 5.

## 2 Cosmological Models and Analysis Method

As mentioned in the previous section, in this work we have considered 5 different models of cosmology to obtain bounds on the sum of three active neutrino masses. Below we list the vector of parameters to vary in each of these cosmological models.

For model:

 θ≡[ωc, ωb, Θs, τ, ns, ln[1010As],∑mν]. (6)

For model:

 θ≡[ωc, ωb, Θs, τ, ns, ln[1010As],∑mν,r]. (7)

For both the models (with or without phantom dark energy) :

 θ≡[ωc, ωb, Θs, τ, ns, ln[1010As],∑mν,w0,wa]. (8)

For the model (non-phantom dark energy with tensors) :

 θ≡[ωc, ωb, Θs, τ, ns, ln[1010As],∑mν,w0,wa,r]. (9)

Here the first 7 parameters for all the models are same. Out of them, the first 6 correspond to the model. Here and are the present-day physical CDM and baryon densities respectively. is the the ratio between the sound horizon and the angular diameter distance at decoupling. is the optical depth to reionization. and are the power-law spectral index and power of the primordial curvature perturbations, respectively, at the pivot scale of . The 7th parameter is of course which is of our biggest concern in this work.

In , along with scalar perturbations we also include tensor perturbations and let another parameter to vary, which is the tensor-to-scalar ratio at the pivot scale of . The choice to study this model is motivated by the results from the latest publicly available dataset from the measurement of the BB mode spectrum of the CMB from Bicep-Keck collaboration, namely BK14, which provides an upper bound to 0.07 at 95% C.L, when combined with Planck 2015 and other datasets [37]; while the bound without BK14 data is much more relaxed at 0.12 [13]. We expect this data to modify the constraints on the neutrino related parameters.

For the two dynamical dark dark energy models, we again only concentrate on scalar perturbations, but the background cosmology with the dark energy equation of state (EoS) is replaced by a varying equation of state with the following parametrization in terms of the redshift :

 w(z)=w0+waz1+z. (10)

This parametrization is famously known as the Chevallier-Polarski-Linder (CPL) parametrization [38, 39]. This just uses the first two terms in a Taylor expansion of the EoS in powers of the scale factor . This parametrization is suitable for describing the past expansion history of the universe, especially at high redshifts [39], but other parametrizations might be needed to describe future evolution [40] since as , diverges.

Notice that corresponds to the dark energy EoS today, whereas refers to the dark energy EoS in the very far past. Between these two times, it is easy to see that is a monotonic function. Therefore, to explore only the non-phantom dark energy region of the parameter space i.e. it is sufficient to apply the following hard priors [33]:

 w0≥−1;       w0+wa≥−1. (11)

We abbreviate the model with for all as the NPDDE model, whereas model without any such restriction on the EoS will be simply called the DDE model. A Universe dominated by a phantom dark energy component () would end in a Big Rip in most cosmological models, where dark energy density becomes infinite in a finite time, resulting in dissociation of any bound state, i.e., the “Big Rip” [41]. Such an universe is unphysical in a sense and hence we study the NPDDE model separately.

In our work, we conduct a Bayesian analysis to derive constraints on . For all the parameters listed in Eq. (6), Eq. (7), and Eq. (9), we impose flat priors in our analysis. The prior ranges are listed on the Table 1. We obtain the posteriors using the Markov Chain Monte Carlo (MCMC) sampler CosmoMC [42] which uses CAMB [43] as the Boltzmann solver and the Gelman and Rubin statistics [44] to quantify the convergence of chains.

## 3 Datasets

Below, we provide a description of the datasets used in our analyses. We have used different combinations of these datasets.

Cosmic Microwave Background: Planck 2015:

Measurements of the CMB temperature, polarization, and temperature-polarization cross-correlation spectra from the publicly available Planck 2015 data release [45] are used. We consider a combination of the high- (30 2508) TT likelihood, as well as the low- (2 29) TT likelihood. We call this combination simply as TT. Along with that, we include the Planck polarization data in the low- (2 29) likelihood, and refer to this as lowP. We also consider the high- (30 1996) EE and TE likelihood. This dataset and TT together are referred to as TTTEEE. The high- polarization data might still be contaminated with residual systematics [13], so bounds obtained without the use of high- polarization can be considered slightly more reliable.

Baryon Acoustic Oscillations (BAO) Measurements and Related Galaxy Cluster data:

In this work, we include BAO measurements obtained from various galaxy surveys. We make use of the SDSS-III BOSS DR12 Consensus sample (as described in [46]; uses LOWZ and CMASS galaxy samples at 0.38, 0.51 and 0.61), the DR7 Main Galaxy Sample (MGS) at [47], and the Six-degree-Field Galaxy Survey (6dFGS) survey at [48]. We refer to this combination as BAO. Here is the effective redshift of the particular survey. In some cases, we have also used the full shape measurements of the correlation function and galaxy power spectrum (refer to [46] for details) from the SDSS-III BOSS DR12. We denote this as FS. The full shape of these measurements reveal additional information other than the BAO signal.

Type Ia Supernovae (SNe Ia) Luminosity Distance Measurements:

We also include Supernovae Type-Ia (SNe Ia) luminosity distance measurements from the Pantheon Sample [49] which consists of data from 279 Pan-STARRS1 (PS1) Medium Deep Survey SNe Ia () and combines it with distance estimates of SNe Ia from SDSS, SNLS, various low-z and HST samples. This combined sample of SNe Ia is largest till date and consists of data from a total of 1048 SNe Ia with redshift range . We denote this data set as PAN hereafter. This dataset replaces the Joint Light-curve Analysis (JLA) SNe Ia sample which consists of 740 spectroscopically confirmed type Ia supernovae [50].

Galaxy Cluster Data from South Pole Telescope:

In this work, we use data from the SPT-SZ survey [51] which provides data from a sample of 377 clusters (identified at ). SPT-SZ is a survey of 2500 deg of the southern sky conducted with the South Pole Telescope (SPT, [52]). These galaxy clusters are recognized by their Sunyaev-Zel’dovich (SZ) effect [53] signature. We call this dataset as SZ from now on.

Optical Depth to Reionization:

The optical depth is proportional to the electron number density integrated along the line of sight, and thus most of the contribution to it comes from the time when the universe re-ionizes. We impose a bound on this optical depth to reionization, , taken from [54], in which Planck collaboration has identified, modelled and removed previously unexplained systematic effects in the polarization data of the Planck High Frequency Instrument (HFI) on large angular scales (low-) (the data was not made publicly available). It is currently the most recent and reliable measurement of from Planck data. We shall hereafter refer to this prior as . We use as a substitute for low- polarization data, and hence we exclude the lowP data whenever we apply the prior, to avoid any double counting.

Hubble Parameter Measurements:

We use a Gaussian prior of km/sec/Mpc on from the recent 2.4% determination of the local value of the Hubble constant by [55], combining the anchor NGC 4258, Milky Way and LMC Cepheids. We shall refer to this prior as R16. It is to be noted that different datasets prefer different values of and there is no clear consensus. For instance recent strong lensing observations [56] of the H0LiCOW program give a slightly lower value of km/sec/Mpc, whereas another measurement [57] prefers a much lower value of km/sec/Mpc. The recent SDSS DR12 BAO data prefers an even lower value of km/sec/Mpc [46]. We chose the R16 value as it is in 3.4 tension with Planck 2015 measurements [13], whose measured value of is km/sec/Mpc assuming with 3 active neutrinos of total mass fixed at eV. Using the R16 prior we get an idea of how the parameter bounds will change if cosmology has to accommodate such a large value of Hubble constant, say, through some new physics.

B Mode Polarization data of CMB:

For the BB mode spectrum of CMB, we use the latest dataset available from BICEP/Keck collaboration which includes all data (spanning the range: ) taken up to and including 2014 [37]. This dataset is referred to as BK14.

## 4 Results on ∑mν

For clarity, we have presented and explained the results on the sum of three active neutrino masses separately for each model (see Section 2 for a description of models) in different subsections. All the quoted upper bounds are at 95% C.L. The main results are summarized in Tables 29. Details about models and datasets are given in Section 2 and Section 3 respectively.

### 4.1 Results for the ΛCDM+∑mν Model

In this subsection, we present the 2 (95% C.L.) upper bounds on for the model for various combinations of datasets. Upper bounds on are given at 2 (95% C.L.) while marginalized limits for any other parameter mentioned in the text are given at 1 (68% C.L.). We have divided these results in two separate sections for convenience of analyzing and presenting. First in Section 4.1.1 we present results obtained without using any priors on the optical depth to reionization () and Hubble constant () and discuss the effects of different datasets on the bounds. In Section 4.1.2 we summarize the results obtained using the said priors.

#### Results without τ and H0 priors

In Tables 2 and 3 we present the bounds without applying any Gaussian prior to the low redshift parameters and . In Table 2 bounds are obtained without the use of the high- polarization data from Planck 2015, while in Table 3 it is included. Figure 1 and 2 shows 1-D marginalized posterior distributions for and respectively, for various data combinations. As mentioned in Section 1, CMB TT data alone is not particularly sensitive to masses much lower than 1 eV. This is clearly reflected in the results. The TT data alone can only constrain eV at 95% C.L.. Addition of the high- E mode polarization auto-correlation and temperature-polarization cross-correlation data leads to a higher constraining capability, reducing the bound to eV. This phenomenon of mass bounds getting stronger with addition of high- polarization data is seen throughout all the analyses we have done, and corroborates well with previous studies [13, 31].

Addition of the lowP data makes the bounds significantly stronger, i.e., eV for TT+lowP and eV for TTTEEE+lowP. This can be attributed to lowP data being able to partially do away with degeneracies present with and other parameters like and . If we consider TT data only, an increase in reduces the smearing of the damping tail [58, 59], which can be compensated by an increase in . The value of also needs to increase, as the Planck TT data severely constrains the quantity , which leads to a degeneracy between these two parameters; variations approximately following the relation . Effects of and are also not independent in cosmology. The value of determines the overall amplitude of matter power spectrum. Increase in increases the amplitude, whereas an increase suppresses matter power spectrum in small scales. The low- polarization data can in principle break this degeneracy between and , and consequently the three-way degeneracy between , and . This is possible through the appearance of the well known “reionization bump” in the range in the polarization spectra whose amplitude is in the EE spectra and in the TE spectra [60], and the bump cannot be reproduced by varying other parameters, thus breaking the degeneracy. Indeed, while the TT data alone prefers a , the TT+lowP data prefers a much lower ; a smaller value of leading to a stronger upper bound of . Refer to Figure 3 for a visualization of this effect. Similar inference can be made for TTTEEE and TTTEEE+lowP. However this degeneracy breaking is only partial. A very precise measurement of low- polarization is needed to completely break the degeneracy.

While and are strongly correlated in the Planck TT, and are strongly anti-correlated. Defining (where with photons, CDM, baryons, and cosmological constant) the comoving distance to the last scattering surface at redshift in a flat universe is given by,

 χ(zdec)=∫zdec0dzH(z)∝∫zdec0dz√ωγ(1+z)4+(ωc+ωb)(1+z)3+ωΛ+ρν(z)h2ρcr,0, (12)

where is the neutrino energy density at a redshift , and is the critical density today. scales differently with redshift, depending on whether neutrinos can be considered as radiation or matter. At late times, when neutrinos become non-relativistic, scales as matter (i.e. ) and depends on . Since in a flat universe, , at late times, the last two terms within the square root in the denominator in Eq. 12 give:

 ωΛ+ρν(z)h2ρcr,0=(1−Ωγ)h2−(ωc+ωb)+ων((1+z)3−1). (13)

Now, is well constrained by CMB accoustic peaks. Since , any change to due to increase in can be compensated by decreasing , i.e., , and hence the anti-correlation.

Addition of the BAO data improves the mass bounds significantly by partially breaking the degeneracy between and . We find that addition of the BAO data to TT + lowP reduces the bound to eV from eV. For the TTTEEE+lowP+BAO case, we get eV, which is also much stronger than the bound without BAO data. One can understand such important changes in bounds by understanding the impact of neutrino masses on the quantity which is measured by BAO using spatial correlation of galaxies. Here is the comoving sound horizon at the end of the baryon drag epoch (the epoch at which baryons decouple from photons, slightly after recombination) and changes in has a small effect on it. On the other hand, the dilation scale at the effective redshift of the survey, is a combination of the angular diameter distance and the Hubble parameter ,

 Dν(z)=[(1+z)2D2A(z)czH(z)]1/3      (c≡speed of light), (14)

and it is the quantity which is affected by most. If is increased while is kept fixed, the expansion rate at early times increases. This requires to decrease to keep fixed, which is very well constrained by the CMB power spectra. Decrease in leads to a increase in , which in turn leads to a decrease in both and . BAO data prefers a higher value of than the CMB spectra, and by rejecting the lower values removes the regions with higher values. See [51] for a detailed discussion on this topic. In our analysis we found that TT+lowP data prefers a value of km/sec/Mpc, whereas TT+lowP+BAO prefers km/sec/Mpc, confirming the above. For TTTEEE+lowP and TTTEEE+lowP+BAO these bounds are km/sec/Mpc and km/sec/Mpc respectively. This effect of BAO data rejecting lower values is evident from Figure 2.

As stated before, the Pantheon Sample (PAN) is the newest dataset available on Supernovae type Ia luminosity distance measurements, replacing its predecessor, the Joint Light-curve Analysis (JLA) sample. Observations of SNe Ia at a range of redshifts ( for the Pantheon Sample) can be used to measure the evolution of luminosity distance as a function of redshift, and thereby determining the evolution of the scale factor [61]. This information can be used to constraint cosmological parameters like dark energy equation of state , and . The PAN dataset also provides substantially stronger mass bound when added to the CMB data, albeit not as strong as BAO data. In particular, TT+lowP+PAN gives a bound of eV, whereas for TTTEEE+lowP+PAN we get eV. The 1 constraints on the Hubble constant are and km/sec/Mpc respectively. These are higher than that of CMB only data but lower than that of CMB+BAO data, which explains the weaker bounds from the Pantheon Sample compared to BAO. On the other hand, inclusion of both BAO and PAN data with CMB produces bounds slightly stronger than CMB+BAO. The bound with TT+lowP+BAO+PAN is eV whereas, for TTTEEE+lowP+BAO+PAN, it is eV, both of which far below the eV bound quoted in [13]. Figure 1 depicts this effect of addition of BAO, PAN and BAO+PAN to the Planck data. Also from Figure 2 we see that the CMB+BAO+PAN combination prefers a slightly higher value of than CMB+BAO. The degeneracy breaking between and due to BAO and PAN can be visualized in Figure 4.

#### Results with τ and H0 priors

In the previous section (4.1.1) we described how lower values of and higher values of help in constraining . Thus precise measurement of these two parameters are instrumental in obtaining meaningful bounds on the sum of neutrino masses. Figures 5 and 6 shows 1-D marginalized posterior distributions for and respectively, for various data combinations. In Tables 4 and 5 we have presented the 95% C.L. bounds on where we have utilized the and R16 priors, along with bounds where we have included the FS and SZ datasets.

The addition of the Gaussian prior significantly improves the bound by strongly breaking the degeneracy between and , which is depicted in Figure 7 Compared to the bound of eV from TT+BAO, TT+BAO+ yields a bound of eV. This change in mass bound can be attributed to a large change in the preferred value of , mostly driven by the prior on (and albeit preferring a slightly lower value of as depicted in Figure 6). For TT+BAO we have the 1 bound of , whereas for TT+BAO+ we have . Similarly for TTTEEE+BAO, we have eV and , and it improves to eV and for TTTEEE+BAO+. We emphasize here again that this use of the prior is well motivated in the sense that, as Planck Collaboration [54] has mentioned in their paper, (1) it is the most accurate bound we currently have on ; (2) such a small value of also fully agrees with other astrophysical measurements of reionization from high redshift sources. For , our tightest bound (except when we remove the MGS data from BAO) without any prior comes from addition of the PAN data. TTTEEE+BAO+PAN+ gives a bound of eV, whereas without the high- polarization data, we achieved eV. This is one of our main results in this paper, and one of the strongest bounds in literature available presently without the use of any prior.

 ∑mν <0.152 eV (95%) (TT+BAO+PAN+τ0p055), (15a) ∑mν <0.118 eV (95%) (TTTEEE+BAO+PAN+τ0p055). (15b)

A prior on helps to break the degeneracy between and in the Planck data. In Figure 8 we demonstrate the same. Addition of the R16 prior ( km/sec/Mpc) leads to even stronger bounds than BAO data; TT+lowP+R16 yields 0.134 eV at 95% C.L., whereas with TTTEEE+lowP+R16 it is 0.125 eV. A very aggressive bound of eV for is obtained with TTTEEE+BAO+PAN+R16+, while the bound with TT+BAO+PAN+R16+ is a bit relaxed at eV. These might be the most stringent bounds ever reported in literature within the minimal model. One can visualize from Figure 5 that the R16 data prefers neutrinos with lower mass much more, due to the preference of significantly higher values of as shown in Figure 6 and the strong anti-correlation present between and . However, as stated before, we need to be cautious with the interpretation of such tight mass bounds, since they are driven by the large 3.4 tension between Planck and R16 measurements of the Hubble constant and since there seems to be no agreement among datasets on the value of .

We notice the bounds can be strengthened further by removal of the DR7 Main Galaxy Sample (MGS) from the BAO data, as can be seen in Tables 4 and 5. We have denoted the MGS removed dataset simply as BAO – MGS. We find that TT + BAO – MGS + PAN + prefers an km/sec/Mpc which is a bit higher than TT + BAO + PAN + , which prefers km/sec/Mpc. The preference of MGS sample for lower values has been discussed in [47]. The lack of MGS data improves the mass bounds to eV for TT + BAO – MGS + PAN + , and eV for TTTEEE + BAO – MGS + PAN + . Adding the R16 prior, we get eV and eV respectively.

Inclusion of the galaxy cluster data from full spectrum measurements (FS) from the SDSS-III BOSS DR12 either worsened or did not help the bounds, as can be seen in Tables 4 and 5. Previous studies [31, 62, 63] have shown that the constraining power of the BAO measurements is higher than that of the full shape measurements in the minimal model, and we find that still to be true for the latest data. Addition of the galaxy cluster data (SZ) from the SPT-SZ survey also worsened the neutrino mass bounds slightly. As shown in Figure 9, both FS and SZ data prefer a slightly lower value of (the normalization of matter power spectrum on scales of Mpc) and thereby favouring slightly larger values of ; as more suppression of matter power spectrum allows for a larger neutrino mass sum, i.e., and are strongly anti-correlated.

### 4.2 Results for the ΛCDM+r+∑mν Model

In this section we present results in the model in Table 6. For the TT+BAO+PAN+ dataset, we see that in the model 0.161 eV, which is a bit relaxed than the 0.152 eV in the minimal model. This is simply due to added degeneracies in a extended parameter space with an extra parameter, , which is the tensor-to-scalar ratio defined at a pivot scale of Mpc. The TT+BAO+PAN+ combination constrains the tensor-to-scalar ratio at 0.13 (95% C.L.), whereas for TTTEEE+BAO+PAN+, we have 0.12 (95% C.L.). Addition of the BK14 data from BICEP/Keck collaboration, which contains information about the CMB BB spectra, strengthens this bound to 0.07 for both the data combinations. It also strengthens the sum of neutrino mass bounds to 0.133 eV and 0.110 eV for TT+BAO+PAN+BK14+ and TTTEEE+BAO+PAN+BK14+ respectively, which are actually lower than the ones quoted in Eq. 15.

CMB B-mode polarization has two well-known sources [64]. The first part comes from the inflationary gravitational waves (IGW), i.e., tensors, which is supposed to produce a bump peaked around (the so called ’recombination bump’) in the BB-mode CMB spectra due to induction of quadruple anisotropies in the CMB within the last scattering surface. The IGW signature cannot be reproduced by scalar perturbations, and the amplitude of the bump depends on the tensor-to-scalar ratio, . The other part comes from the deflection of CMB photons due to gravitational lensing produced by large scale structure at considerably late times, which converts a small fraction of the E mode power into B mode. This lensing BB spectra peaks at around . The ’reionization bump’ is also expected to be present as in the EE spectra, in the 10 region. However, the BK14 data contains information only in the and cannot constrain through the reionization bump.

While the bound on is stronger due to BK14, this does not seems to be the main effect in tightening of the mass bounds. We found that the correlation coefficient (defined as , where and are the two parameters being considered and is the covariance matrix of cosmological parameters) between and to be in case of TT+BAO+PAN+, and in case of TT+BAO+PAN+BK14+, which implies that the correlation is very small before addition of BK14, and there is also no big enough change in the correlation with the addition of BK14 dataset to account for the change in mass bound. The main effect might be coming from lensing BB spectra. Quantitatively, the correlation coefficient between and in TT+BAO+PAN+ is , and in TT+BAO+PAN+BK14+ it is . We find that BK14 data prefers a slightly larger value of (see in Figure 10), and due to the strong anti-correlation present between and in the data, the mass bounds improve a bit. Similar inference can be made for the results including the high- polarization from Planck. As before, inclusion of the R16 prior improves the bounds even more. For TT+BAO+PAN+BK14+R16+, we have a bound of 0.107 eV and for TTTEEE+BAO+PAN+BK14+R16+ it is 0.085 eV, both of which are tighter than the corresponding bounds in minimal model without the BK14 data.

### 4.3 Results for the w0waCDM+∑mν Model (DDE)

In this section we present results for the (DDE) model. The mass bounds are presented in Table 7. For the DDE model we let the dark energy parameters vary in both the phantom and non-phantom range. There is a well-known strong degeneracy between the dark energy equation of state, and sum of neutrino masses, [65]. An increase in can be compensated by a decrease in , due to the mutual degeneracy with . This degeneracy leads to a large degradation of the mass bounds, as can be seen from Table 7 and comparing with the results from the model for the same datasets (see Tables 4 and 5). Figures 11 and 12 provide the 1-D marginalized posterior distributions for and respectively. From Figure 11 we can clearly observe that for the same dataset, the DDE model allows much larger values of than . For TT+BAO+PAN+ we obtain a bound of 0.305 eV, whereas for TTTEEE+BAO+PAN+ the bound is slightly tighter at 0.276 eV. The dynamical dark energy model also helps to reduce the tension between Planck 2015 and R16, by allowing higher values of along with a broader distribution. (see Figure 12). Imposition of the R16 prior improves the mass bounds. However, the magnitude of this effect is less than what we saw in . This is because and are also degenerate, i.e., a change in can be compensated by a change in instead of . This decreases the magnitude of correlation coefficient between and . This phenomenon of changing correlation across these two models can be looked upon in Figure 13. Quantitatively, for TT+BAO+PAN+, the correlation coefficient between and changes from in to in . Also, the DDE model and R16 have a much smaller tension than and R16.

The - diagram in Figure 14 shows that for CMB+BAO+PAN+ only a very small region which corresponds to completely non-phantom dark energy is allowed. Rest of the allowed region in the parameter space crosses the phantom barrier ( line) at some point in the evolution of the universe. We also find that the datasets are compatible with a cosmological constant (, ). Imposing the R16 prior leads to shifting of the contours towards the phantom region. Thus, the allowed non-phantom region shrinks even more. A recent study [66] showed that the disfavouring of the non-phantom region even persists in a 12 parameter extended space. Our results are also in agreement with Planck collaboration [13] which reported similar contours for the given combination of similar but older datasets (see Figure 28 in that paper, for the combination of TT+lowP+ext, where ’ext’ implies combination of BAO, JLA and a prior). In the next two sections we present our results on neutrino mass bounds in a cosmology with only non-phantom dynamical dark energy.

-values:

Previous studies [67, 68] reported a improvement in fit with DDE models compared to . We found similar improvement in our analysis. We compare the best-fit values of the and models. We define, , when used for the same dataset. For TT+BAO+PAN+, we find ; for TTTEEE+BAO+PAN+ it is . The is better with the R16 prior. For TT+BAO+PAN+R16+, we find , whereas for TTTEEE+BAO+PAN+R16+ it is .

See also [69, 70, 71, 72, 73, 74] for previous studies on massive neutrinos and dynamic dark energy together.

### 4.4 Results for the w0waCDM+∑mν Model with w(z)≥−1 (Npdde)

While the current data prefers the phantom region of the dark energy parameter space, it is also important to look at the non-phantom side of the things, since phantom dark energy is somewhat unphysical [75]. Dark energy models with a single scalar field cannot cross the phantom barrier () and more general models that permit the crossing require extra degrees of freedom to provide gravitational stability [76]. Field theories allowing phantom dark energy are fraught with one or more of the following problems like unstable vacuum, Lorentz violation, ghosts, superluminal modes, non-locality, or instability to quantum corrections.There, however, have also been theories where the field theory does not have any such issues but other effects like photon-axion conversion or modified gravity leads to an apparent (see [77] for a brief review). Nonetheless, there are theories like quintessence [78, 79] which are non-phantom in nature and it is important to consider situations where we do not allow the phantom crossing.

The constraints on are shown in Table 8. We find that the restricting ourselves to only the non-phantom sector yields bounds which are even stronger than the minimal model for the same datasets, even though it is an extended parameter space (also previously confirmed in [33]). For TT + BAO + PAN + , we have eV in the NPDDE model, whereas for model, using the same dataset, we had eV. For TTTEEE + BAO + PAN + , in NPDDE, we have eV, compared to eV for . Adding the R16 prior further reduces the allowable mass region, as we have seen throughout this paper. TT + BAO + R16 + PAN + prefers a eV, and TTTEEE + BAO + PAN + R16 + prefers eV, which is below the minimum sum required by the inverted hierarchy.

However this substantial strengthening of neutrino mass bound in NPDDE model compared to DDE model is not surprising when we consider the degeneracy between and . As depicted in Figure 2 of [65], due to strong anti-correlation between and , higher mass sum values prefer a lower value of and on the other hand higher values of for are dominated by very low mass sum values. In NPDDE, what happens is we remove the phantom region, i.e., the portion of the parameter space which likes larger values of neutrino mass sum. As mentioned before, stronger bounds in a NPDDE model compared to was confirmed in a recent study [33], which also confirmed the phenomenon that as one goes away from the line in the non-phantom region of the parameter space the mass bounds get stronger, whereas in the phantom region going away from the line leads to weaker bounds, by running MCMC with separate fixed values of and . A similar effect is seen in the bounds on . Higher values of prefer a lower , and removal of the phantom region of the parameter space leads to a preference towards lower values of . Consequently, a NPDDE model actually increases the tension between Planck CMB data and R16. The alleviation of tension between Planck and R16 in DDE models comes from the phantom region of the plane. One of the consequences of such strong mass bounds is that, if in future neutrino hierarchy is found to be inverted by experiments, an universe with non-phantom dark energy will be less likely than a cosmological constant or phantom dark energy [33]. The 1-D marginalized posteriors for and for the NPDDE model are shown in Figures 15 and 16 respectively.

### 4.5 Results for the w0waCDM+r+∑mν Model with w(z)≥−1 (Npdde+r)

In this section we report results for the model with . We denote this model as “NPDDE+″. The main motivation behind studying this model was to see if we can further strengthen the mass bounds by adding the tensor-to-scalar ratio as a free parameter and adding the BK14 dataset, as in Section 4.2. We find that it is still possible. Once again, the BK14 data prefers a slightly larger value of , as can be observed from Figure 19, which leads to slightly stronger bounds. The 1-D marginalized posterior distributions for and are given in Figures 17 and 18 respectively. The 95% C.L. bounds on are shown in Table 9. Albeit the fact that we don’t know for sure if we live in a universe with non-phantom dark energy or if the debatable R16 prior should be used, the eV bound for TTTEEE + BAO + PAN + BK14 + R16 + dataset for this NPDDE+ model is possibly the strongest bound on ever reported in literature for any kind of cosmological scenario.

## 5 Discussion and Summary

Neutrino oscillation experiments have confirmed that neutrinos are massive with three distinct species. However, still, certain neutrino properties including the sum of the three neutrino masses () have not been precisely determined. Cosmology can put bounds on and in reality, tightest bounds on are obtained from cosmological data. Massive neutrinos leave distinct imprints in the CMB and can be constrained with CMB data. However since neutrinos with masses 1 eV are relativistic during decoupling of photons, CMB data is not particularly sensitive to low values of . Since massive neutrinos also cause suppression in the matter power spectrum, tighter bounds are obtained with large scale structure data. In this work we have used latest cosmological datasets available and provided very strong bounds on the sum of the masses of three active neutrinos in five different cosmological models: , , (DDE), with (NPDDE), and with (NPDDE+). Among datasets, along with CMB data from Planck 2015, we have used BAO measurements from SDSS-III DR12, MGS and 6dFGS; SNe Ia luminosity distance measurements from Pantheon Sample (PAN); the BK14 data from the BICEP/Keck Collaboration; the galaxy cluster data from the SPT-SZ survey and suitable gaussian priors on (R16) and (). The priors help in breaking the mutual degeneracies of and with present in the Planck data. In the minimal model, we obtained a robust bound of eV at 95% C.L. with the use of TT + BAO + PAN + . Adding the high- polarization data tightens the bound to eV. The use of the prior further improves these bounds to eV and eV respectively, showing a weak preference for normal hierarchy. The low bounds obtained with the R16 prior, km/sec/Mpc are debatable since they are driven by the 3.4 tension between Planck data and R16 over the value of . Currently there seem to be no agreement over datasets on the value of . The R16 prior itself is obtained from combining geometric distance calibrations of Cepheids, each of which separately give constraints on : 72.25 2.51, 72.04 2.67, 76.18 2.37, and 74.50 3.27 km/sec/Mpc [55]. Removing the third constraint (obtained from Milkyway cepheids) can reduce the tension and thereby worsen the bounds. While there is no reason to discard the data from Milkyway cepheids we should be cautious while looking at results obtained with the R16 prior. On the other hand, however, there is a possibility that both Planck and R16 might be correct and the discrepancy has to be explained by some new physics, like say, some dark radiation species which contributes to .

In the dynamical dark energy model (DDE) we find that the degeneracy between the dark energy equation of state,