Constraints on the sum of the neutrino masses in dynamical dark energy models with $w(z) \geq -1$ are tighter than those obtained in $Λ$CDM
Sunny Vagnozzi Suhail Dhawan Martina Gerbino Katherine Freese Ariel Goobar Olga Mena
astro-ph.CO gr-qc hep-ph

LCTP-18-03,NORDITA-2018-004,IFIC/18-02




sunny.vagnozzi@fysik.su.se
The Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm University, SE-106 91 Stockholm, Sweden
The Nordic Institute for Theoretical Physics (NORDITA), Roslagstullsbacken 23, SE-106 91 Stockholm, Sweden



The Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm University, SE-106 91 Stockholm, Sweden



The Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm University, SE-106 91 Stockholm, Sweden



The Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm University, SE-106 91 Stockholm, Sweden
The Nordic Institute for Theoretical Physics (NORDITA), Roslagstullsbacken 23, SE-106 91 Stockholm, Sweden
Leinweber Center for Theoretical Physics, Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA



The Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm University, SE-106 91 Stockholm, Sweden



Instituto de Física Corpuscular (IFIC), Universidad de Valencia-CSIC, E-46980, Valencia, Spain




We explore cosmological constraints on the sum of the three active neutrino masses  in the context of dynamical dark energy (DDE) models with equation of state (EoS) parametrized as a function of redshift z by  w(z)=w_0+w_a z/(1+z), and satisfying w(z)≥-1 for all z. We make use of Cosmic Microwave Background data from the Planck satellite, Baryon Acoustic Oscillations measurements, and Supernovae Ia luminosity distance measurements, and perform a Bayesian analysis. We show that, within these models, the bounds on  do not degrade with respect to those obtained in the  case; in fact the bounds are slightly tighter, despite the enlarged parameter space. We explain our results based on the observation that, for fixed choices of w_0 ,w_a such that w(z)≥-1 (but not w=-1 for all z), the upper limit on  is tighter than the  limit because of the well-known degeneracy between w and M_ν. The Bayesian analysis we have carried out then integrates over the possible values of w_0-w_a such that w(z)≥-1, all of which correspond to tighter limits on M_ν than the  limit. We find a 95% credible interval (C.I.) upper bound of <0.13. This bound can be compared with the 95% C.I. upper bounds of <0.16, obtained within the  model, and <0.41, obtained in a DDE model with arbitrary EoS (which allows values of w < -1). Contrary to the results derived for DDE models with arbitrary EoS, we find that a dark energy component with w(z)≥-1 is unable to alleviate the tension between high-redshift observables and direct measurements of the Hubble constant H_0. Finally, in light of the results of this analysis, we also discuss the implications for DDE models of a possible determination of the neutrino mass ordering by laboratory searches.


Constraints on the sum of the neutrino masses in dynamical dark energy models with w(z) ≥ -1 are tighter than those obtained in 
    Olga Mena
    December 30, 2023
================================================================================================================================





§ INTRODUCTION



The nature of the dark energy (DE) driving the accelerated expansion of the Universe remains one of the greatest open problems in cosmology <cit.>. The most economical DE candidate is a cosmological constant (CC) related to vacuum energy density. The equation of state (EoS) of such a component is w_DE=P_DE/ρ_DE=-1, where P_DE and ρ_DE are the pressure and energy density of the DE respectively. The CC model is however at odds with theoretical expectations of the magnitude of the CC, an issue dubbed the cosmological constant problem <cit.>. An alternative solution to this issue posits the existence of a dynamical dark energy (DDE) component <cit.>, which implies a redshift-dependent equation of state w(z).
 
The value w=-1 plays an important role from the theoretical point of view, as it demarcates two very different physical regimes <cit.>. The energy density of a component with w<-1 increases with time as the Universe expands. More importantly, such a component would violate the dominant energy condition, which imposes the inequality | P |≤ρ <cit.> [We note that this issue can be avoided in models in which the effective dark energy component appears to violate the dominant energy condition while the full theory does not (e.g. <cit.>).]. A component with EoS w<-1 is usually referred to as "phantom energy". It has been shown that a Universe dominated by a phantom DE component would end in a Big Rip: the dissociation of any bound system due to the DE energy density becoming infinite in a finite amount of time <cit.>
 [A notable exception is found in the case when w(z) → -1 asymptotically in the future (that is, where the future geometry is asymptotically de Sitter) <cit.>. 
This occurs for instance in bimetric gravity <cit.> as well as in other modified gravity theories (see e.g. <cit.>) or more complex dark energy scenarios <cit.>.]. Examples of non-phantom dark energy models include many quintessence models <cit.> and Cardassian cosmology <cit.>.

Recent works have studied and forecasted cosmological constraints on the sum of the three active neutrino masses, M_ν = ∑_i m_i (where m_i are the masses of the individual mass eigenstates), within the context of the standard cosmological model, the  model, which fixes w=-1 <cit.>. The exact figures vary slightly depending on the datasets used, but combinations of some of the most recent and reliable datasets are converging towards a robust 95% C.I. upper bound of M_ν≲ 0.15 within the  model.

It is the goal of this paper to reconsider these bounds if one retreats from the restricted value of w=-1 assumed by . In our work we consider dynamical dark energy with monotonic redshift dependence w(z) given by the standard Chevallier-Polarski-Linder (CPL) parametrization in Eq. (<ref>) <cit.>. We use a combination of some of the most recent and robust datasets, which include Cosmic Microwave Background (CMB) measurements from the Planck satellite, Baryon Acoustic Oscillation (BAO) measurements from the SDSS and 6dFGS surveys, and Supernovae Type-Ia (SNeIa) luminosity distance measurements from the JLA catalogue.

One might worry that the neutrino mass bounds could weaken dramatically if the parameter space is enlarged to allow for values of the equation of state other than w=-1. Indeed recent work showed that the cosmological bounds on M_ν are weaker when one enlarges the parameter space to other values of w(z) including phantom values w<-1. In fact there exists a well-known degeneracy between the DE EoS w and the sum of the three active neutrino masses M_ν <cit.>. However, the main result of our paper is that the cosmological bounds on neutrino masses in fact become more restrictive for the case of a DE component with w(z)≥-1 than for the standard  case of w=-1. A comprehensive explanation for this effect will be provided.

From neutrino oscillation data, we know that at least two out of the three neutrino mass eigenstates m_i should be massive, as two different mass splittings between the three active neutrinos are measured. We also know that the smallest mass splitting governs solar neutrino transitions and that it is positive. However, current data are not able to determine the sign of the largest mass splitting, which governs atmospheric neutrino transitions. Therefore, we are left with two possibilities: either the largest mass splitting is positive (normal ordering, NO), or it is negative (inverted ordering, IO). Neutrino oscillation data are currently unable to distinguish among the two possible scenarios. Nevertheless, they impose a lower limit to  of M_ν,min≃0.1  eV within IO and M_ν,min≃0.06  eV within NO <cit.>.

If one performs a Bayesian analysis, then we find that our bounds imply a mild preference for NO due to parameter space volume effects. Indeed, the available mass range is larger for NO because it goes all the way down to M_ν≃0.06 rather than only
down to M_ν≃0.1 for IO. We quantify this preference in terms of probability odds, that we compute following Refs. <cit.>. If future laboratory experiments determine the mass ordering to be inverted, and if we exclude non-standard physics either in the neutrino or the gravitational sector, one could conclude that the current accelerated expansion of the Universe is likely driven by a component with w(z)<-1 within DDE models.

The paper is structured as follows: in Sec. <ref>, we introduce the parametrization adopted for the DDE component, and the conditions imposed on the parameters of the DDE model to satisfy w(z)≥-1; in Sec. <ref>, we present the statistical approach and the dataset employed in this analysis. We discuss the results of this analysis in Sec. <ref> and we finally conclude in Sec. <ref>. For the busy reader who wants to skip to the main results, a summary of the bounds obtained is available in Tab. <ref>, and useful visual representations of the same results are also provided in Fig. <ref> and Fig. <ref>.



§ DYNAMICAL DARK ENERGY PARAMETRIZATIONS



The simplest parametrization of a DDE component is the Chevallier-Polarski-Linder (CPL) parametrization <cit.>. In CPL models, the EoS w is parametrized as a function of redshift z as:


    w_DDE(z)    =    w_0 + w_ az/1+z  ,


where w_0 = w_DDE(z=0) denotes the DE EoS at the present time. This equation corresponds to the first two terms in a Taylor expansion of the EoS in powers of the scale factor a=1/(1+z), around the present time. The truncated expansion of Eq. (<ref>) is appropriate if the DE EoS is sufficiently smooth and does not oscillate in cosmic time [For models where the DE EoS oscillates, different parametrizations are required, such as those proposed in Refs. <cit.> or used in recent observational studies <cit.>. It is beyond the scope of this work to extend our considerations on the neutrino mass bounds to these types of DE models.].

It follows from Eq. (<ref>) that the non-phantom [w(z)≥-1; NP] condition can be satisfied by imposing the following hard priors:


    w_0 ≥ -1   ,    w_0+w_a ≥ -1       ( NP).


The first prior imposes the non-phantom condition at present time (z=0). The second prior imposes the same condition in the far past, since z →∞w_DDE(z) = w_0+w_a. The EoS in Eq. (<ref>) is monotonic. Therefore, it is sufficient to impose the NP condition both at the present time and in the far past for the NP condition to hold throughout the Universe expansion history.

The energy density of a  dark energy component corresponding to Eq. (<ref>) takes the form:


    Ω_DE(z) = Ω_DE,0(1+z)^3(1+w_0+w_a)exp ( -3w_az/1+z )  ,


where Ω_DE,0 is the current dark energy density. The DE component dominates over the other components for 0<z≲ z_m DE, with z_m DE≈ 0.3 the redshift of matter-dark energy equality. In this range of redshifts, the energy density of a non-phantom dynamical dark energy model with w(z)≥-1 is always greater than that of a corresponding CC model with the same Ω_DE,0.

A wide class of smooth non-phantom dynamical dark energy models can be probed if we make use of the EoS given by Eq. (<ref>) and we impose the priors in Eq. (<ref>). We shall refer to this class of models with the acronym NPDDE (non-phantom dynamical dark energy). Let us emphasize that the priors in Eq. (<ref>) are crucial in the derivation of the results we obtain. These priors ensure the stability of the dark energy component and differ from priors considered in previous analyses in the literature <cit.>.



§ DATASETS AND ANALYSIS METHODOLOGY



We compute constraints on the sum of the three active neutrino masses M_ν with a combination of the most recent cosmological datasets. We consider measurements of the cosmic microwave background (CMB) temperature anisotropies (TT) from the Planck 2015 data release <cit.>. We impose a Gaussian prior on the optical depth to reionization of τ = 0.055 ± 0.009 as a proxy for measurements of CMB polarisation at large scales from the upcoming Planck 2018 release. This prior choice is motivated by the 2016 Planck reanalyses of low-resolution maps in polarization from the Planck High Frequency Instrument <cit.>. In addition to CMB measurements, we consider Baryon Acoustic Oscillation (BAO) measurements from the following catalogues: the SDSS-III BOSS DR11 CMASS and LOWZ galaxy samples <cit.>, the DR7 Main Galaxy Sample (MGS) <cit.>, and the 6dFGS survey <cit.>. We also include Supernovae Type-Ia (SNeIa) luminosity distance measurements from the SDSS-II/SNLS3 Joint Light-Curve Analysis (JLA) catalogue <cit.>. We refer to the combination of the CMB TT, τ prior, BAO, and SNeIa datasets as “base”.

We also consider the inclusion of CMB polarization and temperature-polarization spectra (TE,EE) at small scales (ℓ>30) from the Planck 2015 data release <cit.> to the baseline dataset. We refer to the combination of the CMB TT,TE,EE, τ prior, BAO, and SNeIa datasets as “pol”.

The cosmological model is described by the usual six parameters of the ΛCDM model: the baryon and cold dark matter physical energy densities Ω_ b h^2 and Ω_ c h^2, the angular scale of the acoustic horizon at decoupling Θ_ s, the optical depth to reionization τ, as well as the amplitude and tilt of the primordial power spectrum A_ s and n_ s. To this set of parameters, we add the DE EoS parameters w_0 and w_a, and the sum of the three active neutrino masses M_ν. We make use of the publicly available Markov Chain Monte Carlo (MCMC) package  <cit.> to efficiently sample the parameter space.

The treatment of M_ν deserves a further comment. Firstly, we assume three massive degenerate neutrinos, i.e. three massive eigenstates with equal mass M_ν/3. This assumption is a valid approximation of the true neutrino mass spectrum, given the current sensitivity of cosmological data <cit.>.

Next, we impose a top-hat prior of M_ν≥0  eV. For the purposes of obtaining bounds on neutrino mass from cosmology, we ignore the lower limit M_ν,min≃0.06  eV set by neutrino oscillation experiments [We note that the prior M_ν≥0  eV is in principle improper since it is unconstrained for M_ν→∞ (see <cit.> for details about proper priors in the astronomy literature). In practice, we adopt a cutoff at M_ν,max=3 eV, which makes our prior proper for all intents and purposes. This choice ensures that a large region of the parameter space is sampled, including regions where we already expect the posterior probability to be vanishing from previous experiments (i.e. the region M_ν>>1  eV). These unlikely regions of the parameter space are in fact quickly discarded by the MCMC sampling algorithm and only the region of highest posterior probability density is effectively sampled. For the purposes of our analysis, this is equivalent to taking the limit M_ν,max→∞, making our proper prior de facto a proxy for the improper prior we describe above. Our result is unaffected by any other choice of sufficiently high value of M_ν ,max, as long as the posterior probability density for M_ν>M_ν ,max is known to be vanishingly small (from previous experiments, analytical considerations, or any other argument).]. We believe this choice is appropriate for the purpose of this work, because it ensures a bound on M_ν relying exclusively on cosmological data. For recent works discussing different choices of prior on M_ν, see for instance Refs. <cit.>. 

Moreover, it is fair to say that the only truly a priori information about M_ν is its positivity, i.e. M_ν≥ 0  eV. The fact that M_ν≥ 0.06  eV coming from oscillation experiments is not a priori but, in fact, a posteriori of the oscillation experiments. While the fact that M_ν≥ 0.06  eV can be incorporated as a prior, it is perhaps formally more correct to include it as an external oscillations likelihood [We thank the referee for bringing this argument to our attention in a very clear way.]. When viewed from this perspective, it is absolutely clear how our choice of adopting the prior M_ν≥0  eV is actually independent of the choice of relying exclusively on cosmological data, but rather reflects the only genuine a priori information really present in the problem. The choice of relying exclusively on cosmological data is instead reflected in our choice of not including the oscillations likelihood. Nonetheless, in  Appendix <ref>, we briefly discuss the impact including oscillations data. We find that including the oscillations likelihood (in the approximate, but still appropriate to zeroth order, form we choose) has no impact on the conclusions of our work. We further note that for all intents and purposes, as far as upper limits on M_ν (which are the subject of this work) are concerned, the oscillations likelihood can be to zeroth approximation included as a sharp cutoff at M_ν=0.06  eV, because the uncertainty to which M_ν,min≃0.06  eV is subject is extremely tiny (but an uncertainty is nonetheless present, whereas the physical lower bound M_ν≥ 0  eV is instead subject o no uncertainty).

In any case, we believe that the approach adopted in this work also allows for a consistency check of the underlying cosmological model. Suppose that we assume a certain cosmological model, and then obtain a cosmological bound on  which lies significantly below ∼ 0.06  eV. Such a cosmological bound would indicate that either the cosmological model in question is in tension with results from oscillation experiments or that non-standard neutrino physics is required. For example, models with non-standard neutrino interactions leading to a vanishing neutrino energy density today have been proposed <cit.>. In these cases, the cutoff of the M_ν prior at M_ν=0  eV can be viewed as a phenomenological proxy of the  effect of a lower energy density of neutrinos with respect to the limits imposed by neutrino oscillation measurements. Note that we are implicitly excluding the possibility that such a finding could be a signal for unaccounted systematics in the dataset employed.

Finally, we combine results from cosmology and neutrino oscillation experiments to quantify the preference for one of the two neutrino mass orderings. In this part of the work, we do impose lower bounds on neutrino mass from oscillation experiments, although we previously did not use them in obtaining bounds on neutrino masses.  We follow the Bayesian approach illustrated in <cit.>. We denote by π(I), π(N) the prior probabilities for the normal and inverted ordering respectively. Then, we compute p_O (where O=N, I), the posterior probabilities for each of the two orderings, as follows:



    -0.6 cm p_O = π(O)∫_0^∞ dm_0   L(D | m_0, O)/π(N)∫_0^∞ dm_0   L(D | m_0, N)+π(I)∫_0^∞ dm_0   L(D | m_0, I)  .



In Eq. (<ref>), m_0 is the mass of the lightest neutrino eigenstate and L(D | m_0, O) is the likelihood of cosmological data D. The above Eq. (<ref>) implicitly assumes a cutoff m_0 ,max→∞ in the prior probability for m_0. Moreover, m_0 ,max should take the same value for both orderings and be large enough so as not to cut the prior in a region where the posterior would otherwise be significantly different from zero (so that it is effectively the data through the likelihood, rather than the prior itself, which cuts the region of high m_0, see Footnote (<ref>) for a previous related discussion when considering the formally improper prior on M_ν). In this way, both the numerator and the denominator should formally contain a 1/m_0 ,max normalization factor, which then cancels out when taking the ratio appearing in Eq. (<ref>). We take L(D | m_0, O) from the analysis of cosmological data illustrated in this work. We use the values of π(I), π(N) from the global Bayesian analysis of neutrino oscillation measurements <cit.>. For further details about how to compute L(D | m_0, O) and get to Eq. (<ref>), we refer the reader to the thorough discussions in <cit.>. We convey the results in terms of probability odds of normal versus inverted ordering (p_ NO:p_ IO).

We follow the approach of <cit.> as it is a quick, yet reliable, way to quantify the preference for the normal ordering in different cosmological scenarios. The method used in this work should be kept in mind when one compares the results quoted here with results from other works. Indeed, we remind that alternative approaches can be adopted to quantify the statistical preference for the neutrino mass ordering <cit.>. For the sake of comparison, in Appendix <ref> we report an alternative estimate of the sensitivity to the mass ordering based on the Akaike information criterion (AIC). The specific outcomes of each analysis should be interpreted only in light of the method adopted. 




§ RESULTS


In this section, we present the bounds on the sum of the three active neutrino masses; we provide a thorough physical explanation of the results; we discuss the Bayesian statistical approach we have used, as well as the dependence of our results on this approach; and we conclude by commenting on the implications of our results for the determination of the neutrino mass ordering.



 §.§ Bounds on Neutrino Masses



Table <ref> shows the bounds on the sum of the neutrino masses M_ν for three cases: a) a dark energy component satisfying the dominant energy condition, with equation of state (EoS) w(z)≥-1 throughout the expansion history of the Universe (non-phantom dynamical dark energy, NPDDE); b) the standard cosmological model () with cold dark matter and a cosmological constant where w(z)=-1 is fixed; c) a generic dynamical dark energy (DDE) model with EoS given by the CPL parametrization, Eq. (<ref>), with w_0 and w_a free to vary even within the phantom region where w(z)<-1. We refer to this last model as w_0w_aCDM. Constraints on  are presented for the two different combinations of cosmological datasets, base and pol, described at the beginning of Sec. <ref>.




For the ΛCDM model, we find 95% C.I. upper bounds of M_ν<0.16 for the base dataset and <0.13 for the pol dataset. When we instead assume the more generic w_0w_aCDM model, that also allows for w(z)<-1, the 95% C.I. upper bound on M_ν is significantly relaxed to <0.41 for the base dataset and M_ν<0.37  eV for the pol dataset. These broader bounds are expected, given the well known degeneracy between M_ν and an arbitrary DDE component.

We now consider a NPDDE model and impose w(z)≥-1 throughout the expansion history. In this case, we find the stringent 95% C.I. upper bounds of <0.13 for the base dataset and M_ν<0.11  eV for the pol dataset. Therefore, we find that the constraints on the sum of the neutrino masses in dynamical dark energy models with w(z)≥-1 are slightly tighter than those obtained in Λ CDM, despite the enlarged parameter space (two extra parameters) in NPDDE models. We note that the upper bounds found within the NPDDE model are also very close to the minimal mass allowed in the inverted ordering scenario, M_ν,min≃0.1.

Figure <ref> depicts the one-dimensional posterior probabilities of M_ν for the w_0w_aCDM generic DDE case (in blue), the  case (in black), and the NPDDE model with w(z)≥-1 (in red). Results for the two dataset combinations employed in this work are shown: solid for base, dashed for pol. For each dataset combination, the significant shift of the upper bounds on M_ν to smaller values is visually clear as one moves from the blue to the red curves. The vertical black dotted-dashed line corresponds to the minimal mass of M_ν,min≈ 0.1  eV allowed by neutrino oscillation data within the inverted ordering.







 §.§ Physical Explanation of Results



We have observed that the bounds on  in the NPDDE model are not weaker – and actually slightly tighter – than those in . Here we provide the physical explanation for these results. The reader can refer to Refs. <cit.> for comprehensive reviews of the effects of massive neutrinos in cosmology [Modifying the assumed expansionary history of the Universe will generically lead to different conclusions concerning M_ν. The reason is that the effect of M_ν on cosmological data is degenerate with other parameters governing the expansionary history, such as the DE EoS w: in other words, the effect on cosmological observables of a change in M_ν (for instance, the resulting change in the distance to last-scattering, further discussed later in this Section) can be compensated by adjusting these other parameters. Therefore, an expansionary history which is different from ΛCDM leads to bounds on M_ν which are different from those obtained assuming ΛCDM. Obviously, the same is true if the expansionary history is restricted to a class of models of which ΛCDM represents a particular case: in our case, ΛCDM represents the particular case of the NPDDE class of models, when w_0=-1 and w_a=0.].

The CMB temperature data accurately constrain the position and amplitude of the first acoustic peak in the CMB power spectrum. These constraints entail a very precise determination of the angular size of the sound horizon at decoupling Θ_ s and of the redshift of matter-radiation equality z_eq. Therefore, any change in the DE sector should be compensated by shifts in the other cosmological parameters such that Θ_ s and z_eq remain approximately fixed.

The angular size of the first peak Θ_ s is defined as the ratio between the sound horizon at decoupling r_s and the angular diameter distance to last scattering D_A. The sound horizon at decoupling r_s is essentially fixed by pre-recombination physics. It is thus unaffected by changes in the dark energy sector, which are only relevant at late times. The angular diameter distance to last scattering D_A is instead sensitive to the late-time evolution of the Universe. Therefore, D_A is affected by the physics of dark energy.

In order to keep Θ_ s unchanged in the NPDDE framework, it is necessary that D_A remains fixed as well. Up to proportionality factors, D_A is given by:


    D_A(z_LS) ∝1/H_0∫_0^z_LSdz/E(z)  ,


where z_LS denotes the redshift of last-scattering. The function E(z) denotes the Hubble parameter at redshift z normalized by its value today. In the NPDDE model it is given by:


    E(z)    =   H(z)/H_0
       ≃   √((Ω_c+Ω_b)(1+z)^3+Ω_DE(z)+Ω_ν(z))  .


In the above equation, Ω_c and Ω_b are the current cold dark matter and baryon energy densities respectively, Ω_DE(z) is the dark energy density given by Eq. (<ref>), and Ω_ν(z) is the neutrino energy density. At late times, after neutrinos become nonrelativistic, Ω_νh^2 ≈ M_ν (1+z)^3 /93.14  eV, where h=H_0/(100  km  s^-1  Mpc^-1). In writing Eq. (<ref>), we are neglecting the contribution of the photon energy density Ω_γ, which is negligible at the redshifts under consideration.

It is easy to show that the normalized expansion rate E(z) at late times is higher in a NPDDE Universe than in a ΛCDM one for fixed values of Ω_c, Ω_b, M_ν and Ω_DE,0, as well shall comment more on below and in Fig. <ref>. The integral in Eq. (<ref>) at fixed M_ν is therefore smaller in the NPDDE case than in the  one. In order to keep D_A fixed, one is left with the option of decreasing both H_0 and M_ν. This option is preferred over the choice where one parameter is decreased by a greater amount while the other parameter is kept fixed, since in the latter case the more sizeable decrease of the first parameter can lead to undesired changes in other regions of the CMB spectra, despite Θ_ s being kept fixed. One could argue that the same effect can be obtained by decreasing Ω_c and/or Ω_b. However, this choice would alter the redshift of matter-radiation equality, which is accurately constrained by the amplitude of the first acoustic peak in the CMB power spectrum. Therefore, it is not the preferred choice. This physical explanation for the shifts in the bounds of M_ν is fully supported by the results of our Monte Carlo analyses. In particular, we have verified that the posterior of Θ_ s is nearly unchanged when moving from the ΛCDM scenario to the NPDDE model.




From the explanation above it follows that in NPDDE models we expect a lower Hubble constant H_0 and/or a lower sum of the neutrino masses M_ν compared to the  case. The shifts in H_0 and  are necessary to keep Θ_ s fixed. Therefore, the very strong anti-correlation (degeneracy) between M_ν and H_0 present in ΛCDM is weakened in NPDDE models.  Note also that in NPDDE models we expect a lower value of σ_8, thus reducing the tension between primary CMB and cluster counts/weak lensing measurements (see e.g. <cit.> for some works examining this tension and possible solutions).

In Fig. <ref> we show the two-dimensional joint H_0-M_ν posterior for the base dataset. The blue contours are obtained in the  model, the red contours in the NPDDE model where w(z)≥-1, and the grey contours for the more generic w_0w_aCDM model where also w(z)<-1 is allowed. The horizontal dashed line corresponds to M_ν,min≃0.1, the minimal value allowed by neutrino oscillation data in the inverted ordering scenario. The difference between the blue contours () and the red contours (NPDDE model) is compatible with the shifts in H_0 and M_ν required to keep Θ_ s fixed. The green band in Fig. <ref> corresponds to the 68% C.I. for H_0 inferred by direct measurements from the Hubble Space Telescope <cit.>. From Fig. <ref>, it is clear that the tension between direct measurements and cosmological estimates of H_0 is not resolved, and actually worsened, within a NPDDE model. The tension can be partially alleviated by a generic dark energy component (the w_0w_aCDM model) able to access the region w(z)<-1 (grey contours) <cit.> [A notable exception to this statement is however provided by running vacuum models (RVM), motivated by quantum field theoretical considerations <cit.>. Studies have shown that RVMs, which appear to be statistically preferred over ΛCDM, can address the H_0 tension, but do so invoking a non-phantom dark energy component (with the value w=-1 being preferred, see e.g. <cit.> for some recent works). Notice that RVMs point towards values of H_0 which are closer to the CMB inferred value.]. We have checked that similar considerations apply to the corresponding contour plot obtained with the pol dataset.

From Fig. <ref> we also see that the anti-correlation (degeneracy) between M_ν and H_0 is weakened when moving from ΛCDM (blue) to NPDDE models (red). The magnitude of the degeneracy is reflected by the tilt of the main axes of the ellipsoidal M_ν-H_0 contours. The contour in the ΛCDM case is visibly more inclined than the NPDDE one. The weakening of the M_ν-H_0 degeneracy can be rigorously quantified by computing the correlation coefficient between the two parameters. The correlation coefficient between two parameters i and j, R_ij, is defined as R_ij = C_ij/√(C_iiC_jj), where C is the covariance matrix of the cosmological parameters (in our case i=M_ν and j=H_0), estimated from our MCMC runs. For the base [pol] dataset, we find a correlation coefficient of -0.43 [-0.50] in the ΛCDM case, which is lowered to -0.14 [-0.16] in the NPDDE case. Therefore, the correlation between the two parameters is strongly reduced in moving from ΛCDM to NPDDE models.





We shall now demonstrate that the late-time expansion rate E(z) is higher in a Universe with w(z)≥-1 compared to . We shall also identify the redshift range in which this effect is most prominent. We define the following quantity:









    E(z) ≡ ( E(z)|_Λ CDM/E(z)|_ NPDDE )^2|_Ω_m , Ω_DE,0 - 1   ,


where |_Λ CDM and |_ NPDDE indicate that E(z) is evaluated in a ΛCDM Universe or in a Universe with w(z)≥-1 respectively. The notation |_Ω_m , Ω_DE,0 denotes that Ω_m=Ω_c+Ω_b+Ω_ν and Ω_DE,0 are kept fixed when moving from ΛCDM to NPDDE. E(z)=0 therefore corresponds to the ΛCDM case. A negative E(z) instead indicates that the expansion rate normalized by H_0 is higher in the NPDDE model compared to ΛCDM. Note that E(z) is closely related to other diagnostics used in the literature to probe the DE evolution, such as the Om diagnostic <cit.>. In Fig. <ref>, E(z) is plotted for three choices of w_0, w_a. All of the choices satisfy the stability priors imposed by Eq. (<ref>) and ensure that w(z)≥-1.




Figure <ref> clearly shows that E(z) is negative at low redshifts, as expected from the above discussion. ℰ(z) also shows a minimum for z ≈ 0.5 for values of w_0 and w_a that are allowed by cosmological data. The four vertical dashed lines indicate the redshift of the four BAO measurements we consider in this work. The grey shaded band refers to the redshift coverage of the JLA Supernovae Type-Ia sample we consider in this analysis. Thus, we see that the measurements adopted in this work cover the redshift range where the difference between E(z) and E(z)=0 is largest. Therefore, the redshift range of current BAO and SNeIa measurements is ideal to probe the dynamics of non-phantom [w(z)≥-1] dark energy.



 §.§ Comment on the Bayesian statistical approach adopted



Here we comment on the a priori counterintuitive fact that the bounds on M_ν in the NPDDE model are tighter than those in ΛCDM, despite the fact that ΛCDM represents the limiting case of NPDDE when w_0=-1 and w_a=0. Indeed  these tighter bounds are a result of our use of a Bayesian statistical approach <cit.>.

To explain our results, we begin by fixing the parameters w_0 and w_a to specific values not corresponding to  (i.e., w_0 ≠ -1 and w_a ≠ 0), yet still satisfying w(z)≥-1.  Following the explanation of the previous Section, we expect that the bounds on M_ν must become ever tighter as the dark energy model gets farther away from ΛCDM.  We will study specific cases below and find that these expectations are met. Therefore, a Bayesian analysis marginalizing over the range of w_0, w_a values satisfying w(z)≥-1 is expected to obtain a bound on M_ν which is slightly tighter than the  one, as shown by the results in Sec. <ref>

Specifically we considered four test cases: a) w_0=-0.95, w_a=0, b) w_0=-0.95, w_a=0.05, c) w_0=-0.9, w_a=0, and d) w_0=-0.85, w_a=0, and found 95% C.I. upper bounds of a) M_ν<0.13, b) M_ν<0.12, c) M_ν<0.11 and d) M_ν<0.08. Indeed the bounds on neutrino mass are tighter than in the case of standard ΛCDM.
 
The posterior distributions of M_ν are shown in Fig. <ref>, in dashed light blue, dashed purple, dashed yellow, and dashed red for cases a)-d) respectively. The  bound is instead represented by the solid black line. It is visually clear that the bounds for these cases are all tighter than the  one. For pedagogical purposes we have also considered two cases where w_0 ≠ -1 and w_a ≠ 0 are instead fixed to values such that w(z)≥-1 is not satisfied: e) w_0=-1.05, w_a=0, and f) w_0=-1.05, w_a=0.05. The corresponding posterior distributions are shown in dashed dark blue and dashed green in Fig. <ref> and corresponds to 95% C.I. upper bounds of e) M_ν<0.19 and f) M_ν<0.18. As per our expectations, the bounds for cases e) and f) are looser than the  one. To further back up this argument, we show a triangular plot in M_ν-w_0-w_a space in Fig. <ref>, where we compare constraints obtained assuming the w_0w_aCDM model (blue contours) and the NPDDE model (red contours). From there it is clear that restricting the allowed region to the NPDDE parameter space inevitably selects the region of parameter space with very low M_ν, due to the direction of the mutual M_ν-w_0-w_a degeneracies.








In the Bayesian statistical approach adopted to obtain the results in Sec. <ref>, w_0 and w_a are not fixed, but rather varied. Subsequently, the uncertainty in w_0 and w_a is integrated out by the process of marginalization, leading to the marginalized posterior on M_ν. Heuristically, this procedure can be viewed as a weighted average over the range of prior possibilities of w_0 and w_a, with weights given by the value of the prior in that particular point of parameter space. For each of these prior possibilities of w_0 and w_a, we have already seen that the corresponding bound on M_ν is tighter than the  bound, see examples a)-d) above as well as the green-shaded region in Fig. <ref>. Therefore, the weighted average of such bounds is expected to be not weaker than the  constraint, as confirmed by our results in Sec. <ref>. This explains the a priori counter-intuitive fact that the upper limits on M_ν for the NPDDE model are slightly tighter than the  limit despite the enlarged parameter space.



 §.§ Implications for the determination of the neutrino mass ordering



Finally, we comment on the implications that the results of this work could have for the determination of the neutrino mass ordering. By integrating the posterior distributions of  for both the base and pol datasets (solid red and dashed red lines in Fig. <ref> respectively) it is straightforward to show that a significant fraction (≳ 90%) of the M_ν posterior probability lies in the range M_ν<0.1  eV. This region is precluded to the inverted mass ordering by neutrino oscillation data.

Should non-cosmological probes such as long-baseline neutrino oscillations experiments (for example T2K <cit.>, NOvA <cit.>, or DUNE <cit.>) establish that the neutrino mass ordering is inverted, the viability of dark energy models with w(z)≥-1 could be jeopardized. This conclusion holds if we exclude exotic physics at play in the cosmological neutrino sector and/or in the gravitational sector. Examples of such exotic models are those with non-standard neutrino interactions predicting a vanishing neutrino energy density today <cit.> or mass-varying neutrinos <cit.>, and models of modified gravity where the bounds on M_ν could be significantly different from those in  <cit.>.

Therefore, we have brought to light a subtle and perhaps unexpected connection between two at first glance seemingly disconnected fields: neutrino oscillation experiments and the nature of dark energy. In the near future, results from the former might be able to shed important light on the latter. It is also worth noticing that our findings could also be very interesting in light of the recently revived Swampland conjectures <cit.> (see also e.g. <cit.>), which suggest that it is not possible to construct metastable de Sitter vacua in a controlled way within string theory. As a corollary, if string theory is the correct high-energy description, the current period of accelerated expansion should be sourced by a quintessence field, e.g. through slowly rolling moduli field which naturally arise in string compactification scenarios. If future long baseline experiments should find the neutrino mass ordering to be inverted, this scenario would naturally be put under pressure, with extremely interesting implications concerning viable high-energy theories.

Finally, we quantify the preference for the normal ordering within NPDDE models in terms of probability odds (p_ NO:p_ IO). We adopt the methodology outlined in Sec. <ref>. For the NPDDE model, where w(z)≥-1, we find that the normal ordering is mildly preferred with posterior odds ∼2:1 for the base dataset and ∼3:1 for the pol dataset.

We compare these figures to those obtained in  the generic w_0w_aCDM model. In this case, we find no preference for any of the two orderings for both the base and pol datasets (posterior odds of ∼1:1). When assuming the standard  cosmological scenario, we find a mild preference for normal ordering of ∼2:1 for both the base and pol dataset combinations.

Finally, in Appendix <ref>, we provide an alternative approach to quantify the preference for the normal ordering. This alternative approach is based on the Akaike Information Criterion (AIC) estimator for the relative quality of statistical models. The findings are qualitatively in agreement with those reported in this section.



§ SUMMARY AND DISCUSSION



A dynamical dark energy (DDE) component driving cosmic acceleration provides an alternative to the cosmological constant. In this work, we have explored cosmological constraints on the sum of the three active neutrino masses M_ν within DDE models. We parametrize the dark energy equation of state (EoS) as a function of redshift z through the usual CPL parametrization w(z)=w_0+w_a z/(1+z). Furthermore, we impose the requirement that  the EoS satisfies w(z)≥-1 throughout the expansion history. We refer to this class of models as non-phantom dynamical dark energy (NPDDE). We employ a combination of CMB, BAO and SNeIa measurements. We denote by base the dataset combination not including CMB polarization data at small scales, and by pol the dataset combination which includes these CMB polarization data.

The conclusions we reach are threefold:


  * We find that the constraints on M_ν assuming a NPDDE model are slightly tighter than those obtained within the standard  scenario. This is the opposite of what is found when a generic DDE model with EoS allowed to enter the region where w(z)<-1 (w_0w_aCDM model) is assumed. More in detail, we find 95% C.I. upper bounds of <0.13 for the base dataset and <0.11 for the pol dataset in a NPDDE model. These figures can be compared to the 95% C.I. upper bounds of <0.16 for the base dataset and <0.13 for the pol dataset in a  model. For the w_0w_aCDM model, we find instead the 95% C.I. upper bounds of <0.41 for the base dataset and <0.37 for the pol dataset. We provide a thorough data-supported physical and statistical explanation of these results. The explanation is based on the effects of massive neutrinos and dark energy on the background cosmological evolution, as well as on the Bayesian statistical method adopted.


  * A DDE component with w(z)≥-1 does not alleviate the tension between cosmological and direct measurements of H_0, contrary to what is found in dark energy models with arbitrary w(z). We find that NPDDE models prefer lower values of H_0 than those inferred by direct measurements. We also show that the well known degeneracy between H_0 and M_ν is reduced within NPDDE models. We provide a thorough explanation of this finding.


  * We combine the results of the cosmological analysis with neutrino oscillation data, and quantify the statistical preference for one of the two neutrino mass orderings over the other. The constraints on  in NPDDE models correspond to probability odds of ∼2:1 in favour of normal ordering with respect to inverted ordering for the base dataset combination, and ∼3:1 for the pol dataset. These odds show a mild preference for normal ordering. If laboratory experiments determine that the neutrino mass ordering is inverted, and if the current cosmic acceleration is caused by a dynamical dark energy component, this component would likely be phantom [w(z)<-1], or at least have to cross the phantom divide at some point during the expansion history. The conclusion holds as long as we exclude non-standard scenarios either in the neutrino sector or in the gravity sector. Therefore, this result brings to light a perhaps unexpected connection between two at first glance seemingly disconnected fields: neutrino oscillation experiments and the nature of dark energy. In the near future, results from the former might be able to shed important light on the latter.



The authors would like to thank Per Andersen, Thejs Brinckmann, Eleonora Di Valentino, Enrique Fernández Martínez, Edvard Mörtsell, Eiichiro Komatsu, Massimiliano Lattanzi, Eric Linder, Matteo Martinelli, Savvas Nesseris, Lorenzo Sebastiani, Zachary Slepian, Alessandra Silvestri, Michael Turner, and Sergio Zerbini for enlightening discussions. This work is based on observations obtained with Planck (http://www.esa.int/Planckwww.esa.int/Planck), an ESA science mission with instruments and contributions directly funded by ESA Member States, NASA, and Canada. We acknowledge use of the Planck Legacy Archive. We also acknowledge the use of computing facilities at NERSC. K.F. acknowledges support from DoE grant DE-SC0007859 at the University of Michigan as well as support from the Leinweber Center for Theoretical Physics. K.F., M.G., and S.V. acknowledge support by the Vetenskapsrådet (Swedish Research Council) through contract No. 638-2013-8993 and the Oskar Klein Centre for Cosmoparticle Physics. S.D. and A.G. acknowledge support by the Vetenskapsrådet, the Swedish Space Board, and the K & A Wallenberg foundation. O.M. is supported by PROMETEO II/2014/050, by the Spanish Grant FPA2014–57816-P of the MINECO, by the MINECO Grant SEV-2014-0398 and by the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreements 690575 and 674896. O.M. would like to thank the Fermilab Theoretical Physics Department for its hospitality.

Note— After our paper appeared on the arXiv, the work <cit.> appeared which considers also neutrino mass bounds within quintessence models. Their work differs from ours in the parametrizations adopted for w(z). Moreover, after our paper appeared on the arXiv, two further works <cit.> confirmed our result that the bounds on M_ν in non-phantom dynamical dark energy models are tighter than those obtained in ΛCDM. All codes, chains, and scripts used to produce the results and plots of this work will be made publicly available at https://github.com/sunnyvagnozzi/NPDDE after acceptance of the paper in a journal.


0.3 cm





§ INCLUDING THE NEUTRINO OSCILLATIONS LIKELIHOOD



Throughout this work, we have imposed a top-hat prior on M_ν of M_ν≥0  eV. That is, we allowed values of M_ν below the minimum value set by oscillation experiments of 0.06  eV. The rationale behind this choice, as we outlined in Sec. <ref>, was threefold:



  * To obtain a bound relying exclusively on cosmological data.

  * To remain open to the possibility of models with non-standard neutrino interactions leading to a neutrino energy density which is either vanishing or lower than the expectation in ΛCDM (e.g. <cit.>): at the level of cosmological data these effects can be phenomenologically captured by considering values of M_ν below the lower bound set by oscillation experiments.

  * To provide an (in)consistency test for DE models where the upper bound on M_ν ends up lying significantly below the lower bound set by oscillation experiments. While this possibility has not been realized in our work due to insufficient sensitivity, it might be realized in the near-future thanks to dramatic improvements in the sensitivity of future CMB, large-scale structure, and supernovae distance measurements data sets.


Moreover, as also explained in Sec. <ref>, the positivity of M_ν is the only genuine a priori information present in the problem, whereas the information that M_ν≥ 0.06  eV is not truly a priori, but rather a posteriori of oscillation experiments. Therefore, the formally correct way of incorporating such information is, in fact, by including the neutrino oscillations likelihood. In addition, as discussed in Sec. <ref>, as far as upper limits on M_ν (which are the subject of this work) are concerned, the oscillations likelihood can be to zeroth approximation included as a sharp cutoff at M_ν=0.06  eV, because the uncertainty to which M_ν,min≃0.06  eV is subject is extremely tiny (but an uncertainty is nonetheless present, whereas the physical lower bound M_ν≥ 0  eV is instead subject o no uncertainty). Let us first discuss this simplified case where the oscillations likelihood is simply included as a sharp cutoff in the M_ν prior, before discussing a more physical, but still simple, way of including the oscillations likelihood.

One might at this point wonder whether or not our results are dependent on the choice of prior: M_ν≥ 0  eV versus M_ν≥ 0.06  eV. In fact, the specific bounds on M_ν within a given model (in this case, w_0w_aCDM, ΛCDM, and NPDDE) are certainly affected by the choice of prior. Nonetheless, it is easy to show that if the prior on M_ν is chosen to be flat even when the lower bound from oscillation experiments is enforced, then as a consequence of Bayes' theorem the key result of our paper remains unchanged. That is, the constraints on M_ν in NPDDE models remain tighter than those obtained in ΛCDM, even when the lower limit of M_ν≥0.06  eV is enforced.

Let us denote by x our data and by θ the set of cosmological parameters excluding M_ν. Let us further denote by L(x| M_ν, θ) our likelihood, and by π(M_ν) and π(θ) the prior distributions on M_ν and θ respectively. Note that we are implicitly assuming that the prior on M_ν can be factorized from the prior on the other cosmological parameters, an assumption which is realized. From Bayes' theorem we know that the posterior distribution of M_ν given the data, p(M_ν|x), is given by the following:


    p(M_ν|x) ∝∫ dθ  L(x| M_ν, θ)π(M_ν)π(θ)   .


Assuming we keep a flat prior on M_ν, the only effect of imposing the lower limit from oscillation experiments is to cut π(M_ν) at 0.06  eV instead of 0  eV. From Eq. (<ref>), we see that the result of this operation would be to shift the posterior of M_ν to higher values: this will affect all quantities computed from the distribution, such as the mean and the 95% C.I. upper bound, both of which would increase, hence leading to broader constraints.

However, in our work we are not interested in the bounds on M_ν per se. The purpose of our work is to examine how the upper limits on M_ν change when moving from ΛCDM to NPDDE. In Fig. <ref>, we showed how the posterior of M_ν obtained assuming the NPDDE model is shifted to lower values compared to the one obtained assuming the ΛCDM model. From Eq. (<ref>), it is easy to see how this fact continues to be true even when the lower limit of 0.06  eV set by oscillation experiments is imposed. Therefore, we expect as a consequence of Bayes' theorem that the upper limits on M_ν in NPDDE models will still be tighter than those obtained in ΛCDM regardless of whether a prior of M_ν≥0  eV or M_ν≥0.06  eV is chosen.

To confirm the above statement explicitly, we recompute the posteriors and upper limits on M_ν obtained in Sec. <ref>, this time imposing the lower limit set by oscillation experiments. We recomputed the bounds only for the ΛCDM and NPDDE models (leaving aside w_0w_aCDM, since it is not important for our conclusions), and only for the base dataset (since the pol dataset leads to identical conclusions). For the ΛCDM case, we find that the 95% C.I. upper bound of 0.16  eV degrades to 0.19  eV when imposing M_ν≥0.06  eV. Similarly, when considering the NPDDE model, the limit degrades from 0.13  eV to 0.17  eV. However, we see that in both cases the upper limit obtained assuming NPDDE is tighter than that obtained assuming ΛCDM, confirming the conclusion we reached previously on the basis of Bayes' theorem.

Recall also that we computed the posterior odds for normal ordering versus inverted ordering using the methodology of <cit.> outlined in Sec. <ref>. We wish to clarify that the posterior odds computed in this way are not affected by the choice of prior on M_ν. The reason is that, in Eq. (<ref>), the likelihood of the cosmological data L(D | m_0,O) has been rewritten in terms of the mass of the lightest neutrino mass eigenstate m_0 rather than the total neutrino mass M_ν. The relation between M_ν involves the squared mass-splittings measured by oscillation experiments. Therefore, the methodology adopted factors in, by construction, the information concerning the lower limits on M_ν set for the normal and inverted orderings. As a result, this methodology automatically ignores the region of the M_ν posterior which lies below 0.06  eV. In fact, it is easy to show that in the low-mass region of M_ν parameter space favoured by data, M_ν≲ 0.15  eV, the posterior odds for normal versus inverted ordering calculated using Eq. (<ref>) are well approximated by the following <cit.>:


    p_N/p_O≈∫_0.06  eV^∞dM_ν p(M_ν|x)/∫_0.10  eV^∞dM_ν p(M_ν|x)  ,


where p(M_ν|x) is the posterior of M_ν. The form of Eq. (<ref>) shows how the information on the lower limits of M_ν for both the normal and inverted ordering enters by construction in the methodology adopted.

So far, we considered the extremely simplified case where the information from oscillation experiments is included as a sharp cutoff in the M_ν prior. In fact, the minimum value allowed by oscillation experiments does not come without uncertainty. Using the best-fit values and 1σ intervals from <cit.> for the two mass-squared splittings assuming the normal ordering, we find M_ν ,min≈ 0.0589 ± 0.0005  eV. Therefore, we can to first approximation treat the oscillations likelihood as being a constant down to 0.0589  eV, and as a truncated Gaussian centered at 0.0589  eV and with width 0.0005  eV below. If we simply included the oscillations likelihood as a sharp cutoff, the result on the M_ν posterior would be a sharp drop at 0.06  eV, which is in a sense “unphysical” because this cutoff is actually induced by the value of M_ν ,min inferred from oscillation experiments with some uncertaintly, and not by a physical lower bound which comes without uncertainty (such as M_ν≥ 0  eV). Nonetheless, since M_ν ,min is known to better than ≈ 0.8% precision, it is perfectly reasonable to expect that the impact on the upper limits on M_ν of using our smeared Gaussian approximation versus a sharp cutoff is going to be negligible at best (which we later confirm).

In Fig. <ref>, we show the posteriors we obtained for the ΛCDM (blue) and NPDDE (red) models, for both the case where the oscillations likelihood is not included (dashed) and the case where it is included (solid). We immediately notice two things. The first is that, as expected, posterior drops very sharply below 0.0589  eV, signalling that the sharp cutoff approximation to the oscillations likelihood is in fact a reasonable approximation, given the extremely tiny uncertainty on M_ν ,min. In fact, we find that we recover the upper limits we computed previously from the sharp cutoff approximation (0.17  eV for the NPDDE model and 0.19  eV for the ΛCDM model). The second observation is that, independently of whether or not the oscillations likelihood is included, the NPDDE posterior is always shifted to lower values of M_ν compared to the ΛCDM one, showing that the main conclusions of our paper are stable against the inclusion of the oscillations likelihood.




To conclude, we summarize the findings of this Appendix. The key conclusion of our work, namely the fact that the upper limits on M_ν are tighter in NPDDE models compared to ΛCDM, persists even when the oscillations likelihood is included. Approximating the oscillations likelihood as a sharp cutoff in M_ν (which we have argued is a reasonable zeroth approximation given the very tiny uncertainty on M_ν ,min), we have shown how this result follows simply from Bayes' theorem. Using a more realistic approximation to the oscillations likelihood (treated as a truncated Gaussian), we show the impact of including this likelihood on the M_ν posterior in Fig. <ref>. Since the methodology adopted to compute the posterior odds of normal versus inverted ordering (see Sec. <ref> and <cit.>) by construction takes into account the lower limits on M_ν coming from oscillation experiments for both normal and inverted ordering, the results obtained in Sec. <ref> are unaffected by whether or not the lower limit of M_ν≥0.06  eV is enforced. Therefore, in NPDDE models (and using the base dataset), the preference for normal versus inverted ordering is ∼ 2:1.



§ ESTIMATING THE PREFERENCE FOR THE NORMAL ORDERING THROUGH THE AKAIKE INFORMATION CRITERION


We complete the analysis of this work by quantifying the preference for one neutrino mass ordering over the other using an alternative statistical method based on the Akaike information criterion (AIC) <cit.>. The AIC is a statistical indicator which estimates the relative statistical quality of different models. We use the AIC to estimate the preference for the normal ordering (NO) against the inverted ordering (IO). For a model with k parameters and log-likelihood ln( L) = -χ^2/2, the AIC is given by:


    AIC = 2k+min(χ^2)   ,


where min(χ^2) denotes the minimum value of the χ^2 for the model. The difference between the AICs of two models, Δ AIC, estimates the relative quality of one model against the other. In particular, the model with the lowest AIC is to be considered statistically preferred.

As in Sec. <ref>, we combine results from cosmology and neutrino oscillation experiments to quantify the preference for the normal ordering. We obtain the posterior probability distribution of M_ν from the cosmological analysis, and interpret this posterior as a likelihood for the cosmological dataset (as done in <cit.>). From oscillation measurements, we take the one-dimensional χ^2 projections for the solar and atmospheric mass splittings computed separately for NO and IO, as provided by NuFIT 3.0 (2016) <cit.>. We then compute the global min(χ^2) for both NO and IO in light of the combination of cosmological and oscillation data. 

The number of parameters k is the same in the two scenarios, therefore ΔAIC=Δmin(χ^2). For NPDDE models, we find ΔAIC=AIC_IO-AIC_NO=3.4 for the base dataset and ΔAIC=3.7 for the pol dataset. These values show a mild preference for the NO model, using the scale provided by <cit.>. In , we find ΔAIC=2.7 for the base dataset and ΔAIC=3.1 for the pol dataset. In this case, the preference for NO is even milder. Finally, in generic DDE models with arbitrary EoS, we find ΔAIC=1.9 for the base dataset and ΔAIC=2.0 for the pol dataset. We interpret these values as a minimal preference for NO. The results obtained using the AIC method are qualitatively in agreement with those obtained performing a Bayesian model comparison and reported in Sec. <ref>. Both indicate that current cosmological data, when interpreted in light of a NPDDE model, show a mild preference for the NO.

JHEP


