Unveiling $ν$ secrets with cosmological data: neutrino masses and mass hierarchy
Sunny Vagnozzi Elena Giusarma Olga Mena Katherine Freese Martina Gerbino Shirley Ho Massimiliano Lattanzi





sunny.vagnozzi@fysik.su.se
The Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm University, SE-106 91 Stockholm, Sweden


egiusarm@andrew.cmu.edu
McWilliams Center for Cosmology, Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA
Lawrence Berkeley National Laboratory (LBNL), Physics Division, Berkeley, CA 94720-8153, USA
Berkeley Center for Cosmological Physics, University of California, Berkeley, CA 94720, USA


omena@ific.uv.es
Instituto de Física Corpuscolar (IFIC), Universidad de Valencia-CSIC, E-46980, Valencia, Spain


The Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm University, SE-106 91 Stockholm, Sweden
Michigan 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


McWilliams Center for Cosmology, Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA
Lawrence Berkeley National Laboratory (LBNL), Physics Division, Berkeley, CA 94720-8153, USA
Berkeley Center for Cosmological Physics, University of California, Berkeley, CA 94720, USA


Dipartimento di Fisica e Scienze della Terra, Università di Ferrara, I-44122 Ferrara, Italy
Istituto Nazionale di Fisica Nucleare (INFN), Sezione di Ferrara, I-44122 Ferrara, Italy




Using some of the latest cosmological datasets publicly available, we derive the strongest bounds in the literature on the sum of the three active neutrino masses, M_ν, within the assumption of a background flat ΛCDM cosmology. In the most conservative scheme, combining Planck cosmic microwave background (CMB) temperature anisotropies and baryon acoustic oscillations (BAO) data, as well as the up-to-date constraint on the optical depth to reionization (τ), the tightest 95% confidence level (C.L.) upper bound we find is M_ν<0.151 eV. The addition of Planck high-ℓ polarization data, which however might still be contaminated by systematics, further tightens the bound to M_ν<0.118 eV. A proper model comparison treatment shows that the two aforementioned combinations disfavor the IH at ∼ 64% C.L. and ∼ 71% C.L. respectively.  In addition, we compare the constraining power of measurements of the full-shape galaxy power spectrum versus the BAO signature, from the BOSS survey. Even though the latest BOSS full shape measurements cover a larger volume and benefit from smaller error bars compared to previous similar measurements, the analysis method commonly adopted results in their constraining power still being less powerful than that of the extracted BAO signal. Our work uses only cosmological data; imposing the constraint M_ν>0.06  eV from oscillations data would raise the quoted upper bounds by O(0.1σ) and would not affect our conclusions.


Unveiling ν secrets with cosmological data: neutrino masses and mass hierarchy
    Massimiliano Lattanzi
    December 30, 2023
==============================================================================





§ INTRODUCTION


The discovery of neutrino oscillations, which resulted in the 2015 Nobel Prize in Physics <cit.>, has robustly established the fact that neutrinos are massive <cit.>. The results from oscillation experiments can therefore be successfully explained assuming that the three neutrino flavour eigenstates (ν_e, ν_μ, ν_τ) are quantum superpositions of three mass eigenstates (ν_1, ν_2, ν_3). In analogy to the quark sector, flavour and mass eigenstates are related via a mixing matrix parametrized by three mixing angles (θ_12, θ_13, θ_23) and a CP-violating phase δ_CP.

Global fits <cit.> to oscillation measurements have determined with unprecedented accuracy five mixing parameters, namely, sin^2θ_12, sin^2θ_13, sin^2θ_23, as well as the two mass-squared splittings governing the solar and the atmospheric transitions. The solar mass-squared splitting is given by Δ m_21^2 ≡ m^2_2 - m^2_1 ≃ 7.6× 10^-5 ^2. Because of matter effects in the Sun, we know that the mass eigenstate with the larger electron neutrino fraction is the one with the smallest mass. We identify the lighter state with “1" and the heavier state (which has a smaller electron neutrino fraction) with “2". Consequently, the solar mass-squared splitting is positive. The atmospheric mass-squared splitting is instead given by |Δ m_31^2| ≡ |m^2_3 - m^2_1|≃ 2.5× 10^-3^2. Since the sign of the largest mass-squared splitting |Δ m_31^2| remains unknown, there are two possibilities for the mass ordering: the normal hierarchy (NH, Δ m_31^2 > 0, with m_1<m_2<m_3) and the inverted hierarchy (IH, Δ m_31^2 < 0, and m_3<m_1<m_2). Other unknowns in the neutrino sector are the presence of CP-violation effects (i.e. the value of δ_CP), the θ_23 octant, the Dirac versus Majorana neutrino nature, and, finally, the absolute neutrino mass scale, see Ref. <cit.> for a recent review on unknowns of the neutrino sector.

Cosmology can address two out of the above five unknowns: the absolute mass scale and the mass ordering. Through background effects, cosmology is to zeroth-order sensitive to the absolute neutrino mass scale, that is, to the quantity:


    M_ν≡ m_ν_1 + m_ν_2 + m_ν_3  ,


where m_ν_i denotes the mass of the ith neutrino mass eigenstate. Indeed, the tightest current bounds on the neutrino mass scale come from cosmological probes, see for instance <cit.>. More subtle perturbation effects make cosmology in principle sensitive to the mass hierarchy as well (see e.g. <cit.> for comprehensive reviews on the impact of nonzero neutrino masses on cosmology), although not with current datasets.

As light massive particles, relic neutrinos are relativistic in the early Universe and contribute to the radiation energy density. However, when they turn non-relativistic at late times, their energy density contributes to the total matter density. Thus, relic neutrinos leave a characteristic imprint on cosmological observables, altering both the background evolution and the spectra of matter perturbations and Cosmic Microwave Background (CMB) anisotropies (see <cit.> as well as the recent <cit.> for a detailed review on massive neutrinos in cosmology, in light of both current and future datasets). The effects of massive neutrinos on cosmological observables will be discussed in detail in Sec. <ref>.

Cosmological probes are primarily sensitive to the sum of the three active neutrino masses M_ν. The exact distribution of the total mass among the three mass eigenstates induces sub-percent effects on the different cosmological observables, which are below the sensitivities of ongoing and near future experiments <cit.>. As a result, cosmological constraints on M_ν are usually obtained by making the assumption of a fully degenerate mass spectrum, with the three neutrinos sharing the total mass [m_ν_i=M_ν/3, with i=1,2,3, which we will later refer to as 3deg, see Eq.(<ref>)]. Strictly speaking, this is a valid approximation as long as the mass of the lightest eigenstate, m_0 ≡ m_1  [m_3] in the case of NH [IH], satisfies:


    m_0 ≫| m_i-m_j |    ,   ∀ i,j=1,2,3.


The approximation might fail in capturing the exact behaviour of massive neutrinos when M_ν∼, where =√(Δ m_21^2)+√(Δ m_31^2)≃0.06  [=√(Δ m_31^2)+√(Δ m_31^2+Δ m_21^2)≃0.1 ] is the minimal mass allowed by oscillation measurements in the NH [IH] scenario <cit.>, see Appendix A for detailed discussions. Furthermore, it has been argued that the ability to reach a robust upper bound on the total neutrino mass below =0.1 would imply having discarded at some statistical significance the inverted hierarchy scenario. In this case, one has to provide a rigorous statistical treatment of the preference for one hierarchy over the other <cit.>. 3deg repeated

We will be presenting results obtained within the approximation of three massive degenerate neutrinos. That is, we consider the following mass scheme, which we refer to as 3deg:


    m_1 = m_2 = m_3 = M_ν/3   (3deg)   ,


This approximation has been adopted by the vast majority of works when M_ν is allowed to vary. This includes the Planck collaboration, which recently obtained M_ν < 0.234   eV at 95% C.L. <cit.> through a combination of temperature and low-ℓ polarization anisotropy measurements, within the assumption of a flat ΛCDM+M_ν cosmology. Physically speaking, this choice is dictated by the observation that the impact of the NH and IH mass splittings on cosmological data is tiny if one compares the 3deg approximation to the corresponding NH and IH models with the same value of the total mass M_ν (see Appendix A for further discussions). For the purpose of comparison with previous work, in Appendix B we briefly discuss other less physical approximations which have been introduced in the recent literature, as well as some of the bounds obtained on M_ν within such approximations.

We present the constraints in light of the most recent cosmological data publicly available. In particular, we make use of i) measurements of the temperature and polarization anisotropies of the CMB as reported by the Planck satellite in the 2015 data release; ii) baryon acoustic oscillations (BAO) measurements from the SDSS-III BOSS data release 11 CMASS and LOWZ samples, and from the Six-degree Field Galaxy Survey (6dFGS) and WiggleZ surveys; iii) measurements of the galaxy power spectrum of the CMASS sample from the SDSS-III BOSS data release 12; iv) local measurements of the Hubble parameter (H_0) from the Hubble Space Telescope; v) the latest measurement of the optical depth to reionization (τ) coming from the analysis of the high-frequency channels of the Planck satellite, and vi) cluster counts from the observation of the thermal Sunyaev-Zeldovich (SZ) effect by the Planck satellite.

In addition to providing bounds on M_ν, we also use these bounds to provide a rigorous statistical treatment of the preference for the NH over the IH. We do so by applying the simple but rigorous method proposed in <cit.>, and evaluate both posterior odds for NH against IH, as well as the C.L. at which current datasets can disfavor the IH.

The paper is organized as follows. In Sec. <ref>, we describe our analysis methodology. In Sec. <ref> we instead provide a careful description of the datasets employed, complemented with a full explanation of the physical effects of massive neutrinos on each of them. We showcase our main results in Sec. <ref>, with Sec. <ref> in particular devoted to an analysis of the relative constraining power of shape power spectrum versus geometrical BAO measurements, whereas in Sec. <ref> we provide a rigorous quantification of the exclusion limits on the inverted hierarchy from current datasets. Finally, we draw our conclusions in Sec. <ref>. 

For the reader who wants to skip to the results: the most important results of this paper can be found in Tabs. <ref>, <ref>, <ref>. The first two of these tables present the most constraining 95% C.L. bounds on the sum of the neutrino masses using a combination of CMB (temperature and polarization), BAO, and other external datasets. The bounds in Tab. <ref> have been obtained using also small-scale CMB polarization data which may be contaminated by systematics, yet we present the results as they are useful for comparing to previous work. Finally Tab. <ref> presents exclusion limits on the Inverted Hierarchy neutrino mass ordering, which is disfavored at about 70% C.L. statistical significance.



§ ANALYSIS METHOD



In the following we shall provide a careful description of the statistical methods employed in order to obtain the bounds on the sum of the three active neutrino masses we show in Sec. <ref>, as well as caveats to our analyses. Furthermore, we provide a brief description of the statistical method adopted to quantify the exclusion limits on the IH from our bounds on M_ν. For more details on the latter, we refer the reader to <cit.> where this method was originally described.



 §.§ Bounds on the total neutrino mass



In our work, we perform standard Bayesian inference (see e.g. <cit.> for recent reviews) to derive constraints on the sum of the three active neutrino masses. That is, given a model described by the parameter vector θ, and a set of data x, we derive posterior probabilities of the parameters given the data, p(θ|x), according to:


    p(θ|x) ∝ L(x|θ)p(θ)   ,


where L(x|θ) is the likelihood function of the data given the model parameters, and p(θ) denotes the data-independent prior. We derive the posteriors using the Markov Chain Monte Carlo (MCMC) sampler  with an efficient sampling method <cit.>. To assess the convergence of the generated chains, we employ the Gelman and Rubin statistics <cit.> R-1, which we require to satisfy R-1<0.01 when the datasets do not include SZ cluster counts, R-1<0.03 otherwise (this choice is dictated by time and resource considerations: runs involving SZ cluster counts are more computationally expensive than those that do not include SZ clusters, to achieve the same convergence). In this way, the contribution from statistical fluctuations is roughly a few percent the limits quoted. [Notice that this is a very conservative requirement, as a convergence of 0.05 is typically more than sufficient for the exploration of the posterior of a parameter whose distribution is unimodal <cit.>.]

We work under the assumption of a background flat ΛCDM Universe, and thus consider the following seven-dimensional parameter vector:


    θ≡{Ω_b h^2, Ω _c h^2, Θ_s, τ, n_s, log(10^10A_s), M_ν}  .


Here, Ω_bh^2 and Ω_ch^2 denote the physical baryon and dark matter energy densities respectively, Θ_s is the ratio of the sound horizon to the angular diameter distance at decoupling, τ indicates the optical depth to reionization, whereas the details of the primordial density fluctuations are encoded in the amplitude (A_s) and the spectral index (n_s) of its power spectrum at the pivot scale k_⋆ = 0.05  h   Mpc ^-1. Finally, the sum of the three neutrino masses is denoted by M_ν. For all these parameters, a uniform prior is assumed unless otherwise specified.

Concerning M_ν, we impose the requirement M_ν≥ 0. Thus, we ignore prior information from oscillation experiments, which, as previously stated, set a lower limit of ∼0.06  [0.10 ] for the NH [IH] mass ordering. If we instead had chosen not to ignore prior information from oscillation experiments, the result would be a slight shift of the center of mass of our posteriors on M_ν towards higher values. As a consequence of these shifts, the 95% C.L. upper limits we report would also be shifted to slightly higher values. Nonetheless, in this way we can obtain an independent upper limit on M_ν from cosmology alone, while at the same time making the least amount of assumptions. It also allows us to remain open to the possibility of cosmological models predicting a vanishing neutrino density today, or models where the effect of neutrino masses on cosmological observables is hidden due to degeneracies with other parameters (see e.g. <cit.>). One can get a feeling for the size of the shifts by comparing our results to those of <cit.>, where a prior M_ν≥ 0.06 was assumed. As we see, the size of the shifts is small, of O(0.1σ). We summarize the priors on cosmological parameters, as well as some of the main nuisance parameters, in Tab. <ref>.

All the bounds on M_ν reported in Sec. <ref> are 95% C.L. upper limits. These bounds depend more or less strongly on our assumption of a background flat ΛCDM model, and would differ if one were to consider extended parameter models, for instance scenarios in which the number of relativistic degrees of freedom N_eff and/or the dark energy equation of state w are allowed to vary, or if the assumption of flatness is relaxed, and so on. For recent related studies considering extensions to the minimal ΛCDM model we refer the reader to e.g. <cit.>, as well as Sec. <ref>. For other recent studies which investigate the effect of systematics or the use of datasets not considered here (e.g. cross-correlations between CMB and large-scale structure) see e.g. <cit.>.





 §.§ Model comparison between mass hierarchies



As we discussed previously, several works have argued that reaching an upper bound on M_ν of order 0.1  eV would imply having discarded the IH at some statistical significance. In order to quantify the exclusion limits on the IH, a proper model comparison treatment, thus rigorously taking into account volume effects, is required. Various methods which allow the estimation of the exclusion limits on the IH have been devised in the recent literature, see e.g. <cit.>. Here, we will briefly describe the simple but rigorous model comparison method which we will use in our work, proposed by Hannestad and Schwetz in <cit.>, and based on previous work in <cit.>. The method allows the quantification of the statistical significance at which the IH can be discarded, given the cosmological bounds on M_ν. We refer the reader to the original paper <cit.> for further details.

Let us again consider the likelihood function L of the data x given a set of cosmological parameters θ, the mass of the lightest neutrino m_0 = m_1  [m_3] for NH [IH], and the discrete parameter H representing the mass hierarchy, with H=N  [I] for NH [IH] respectively: L(x|θ,m_0,H). Then, given the prior(s) on cosmological parameters p(θ), we define the likelihood marginalized over cosmological parameters θ assuming a mass hierarchy H, E_H(m_0), as:


    E_H(m_0) ≡∫ dθ  L(x|θ,m_0,H)p(θ) =  L(x| m_0, H)   .


Imposing an uniform prior m_0 ≥ 0  eV and assuming factorizable priors for the other cosmological parameters it is not hard to show that, as a consequence of Bayes' theorem, the posterior probability of a mass hierarchy H given the data x, p_H ≡ p(H|x), can be obtained as below:


    -0.6 cm p_H = p(H)∫_0^∞ dm_0   E_H(m_0)/p(N)∫_0^∞ dm_0   E_N(m_0)+p(I)∫_0^∞ dm_0   E_I(m_0)  ,


where p(N) and p(I) denote priors on the NH and IH respectively, with p(N)+p(I)=1. The posterior odds of NH against IH are then given by p_N/p_I, whereas the C.L. at which the IH is disfavored, which we refer to as CL_ IH, is given by:


    CL_ IH = 1-p_I   .



The expression in Eq. <ref> is correct as long as the assumed prior on m_0 is uniform, and the priors on the other cosmological parameters are factorizable. Different choices of priors on m_0 will of course lead to a larger or smaller preference for the NH. As an example, <cit.> considered the effect of logarithmic priors, showing that this leads to a strong preference for the NH (see, however, <cit.>).

Another valid possibility, which has not explicitly been considered in the recent literature, is that of performing model comparison between the two neutrino hierarchies by imposing an uniform prior on M_ν instead of m_0. In this case, it is easy to show that the posterior odds for NH against IH, p_N/p_I, is given by (considering for simplicity the case where NH and IH are assigned equal priors):


    p_N/p_I≡∫_0.06  eV^∞dM_ν  E(M_ν)/∫_0.10  eV^∞dM_ν  E(M_ν)  ,


where analogously to Eq. (<ref>), we define the marginal likelihood E(M_ν) as:


    E_H(M_ν) ≡∫ dθ  L(x|θ,M_ν,H)p(θ) =  L(x| M_ν, H)   .


It is actually easy to show that in the low-mass region of parameter space currently favoured by cosmological data, i.e. M_ν≲ 0.15  eV, the posterior odds for NH against IH one obtains by choosing a flat prior on M_ν [Eq. (<ref>)] or a flat prior on m_0 [Eq. (<ref>)] are to very good approximation equal. It is also interesting to note that, as is easily seen from Eq. (<ref>), cosmological data will always prefer the normal hierarchy over the inverted hierarchy, simply as a consequence of volume effects: that is, the volume of parameter space available to the normal hierarchy (M_ν>0.06  eV) is greater than that available to the inverted hierarchy (M_ν>0.1  eV). For this reason, the way the prior volume is weighted plays a crucial role in determining the preference for one hierarchy over the other (see discussions in <cit.>).

In our work, we choose to follow the prescription of <cit.> (based on a uniform prior on m_0) and hence apply Eq. (<ref>) to determine the preference for the normal hierarchy over the inverted one from cosmological data.



§ DATASETS AND THEIR SENSITIVITY TO M_Ν



We present below a detailed description of the datasets used in our analyses and their modeling, discussing their sensitivity to the sum of the active neutrino masses. For clarity, all the denominations of the combinations of datasets we consider are summarized in Tab. <ref>. For plots comparing cosmological observables in the presence or absence of massive neutrinos, we refer the reader to <cit.> and especially Fig. 1 of the recent <cit.>.





 §.§ Cosmic Microwave Background


Neutrinos leave an imprint on the CMB (both at the background and at the perturbation level) in, at least, five different ways, extensively explored in the literature <cit.>:



  * By delaying the epoch of matter-radiation equality, massive neutrinos lead to an enhanced early integrated Sachs-Wolfe (EISW) effect <cit.>. This effect is due to the time-variation of gravitational potentials which occurs during the radiation-dominated, but not during the matter-dominated era, and leads to an enhancement of the first acoustic peak in particular. Traditionally this has been the most relevant neutrino mass signature as far as CMB data is concerned.

  * Because of the same delay as above, light (f_ν<0.1) massive neutrinos actually increase the comoving sound horizon at decoupling r_s(z_dec), thus increasing the angular size of the sound horizon at decoupling Θ_s and shifting all the peaks to lower multipoles ℓ's <cit.>.

  * By suppressing the structure growth on small scales due to their large thermal velocities (see further details later in Sec. <ref>), reducing the lensing potential and hence the smearing of the high-ℓ multipoles due to gravitational lensing <cit.>. This is a promising route towards determining both the absolute neutrino mass scale and the neutrino mass hierarchy, see e.g. Ref. <cit.>, because it probes the matter distribution in the linear regime at higher redshift, and because the unlensed background is precisely understood. CMB lensing suffers from systematics as well, although these tend to be of instrumental origin and hence decrease with higher resolution. In fact, a combination of CMB-S4 <cit.> lensing and DESI <cit.> BAO is expected to achieve an uncertainty on M_ν of 0.016 eV <cit.>.

  * Massive neutrinos will also lead to a small change in the diffusion scale, which affects the photon diffusion pattern at high-ℓ multipoles <cit.>, although again this effect is important only for neutrinos which are non-relativistic at decoupling, i.e. for M_ν>0.6  eV.

  * Finally, since the enhancement of the first peak due to the EISW depends, in principle, on the precise epoch of transition to the non-relativistic regime of each neutrino species, that is, on the individual neutrino masses, future CMB-only measurements such as those of <cit.> could, although only in a very optimistic scenario, provide some hints to unravel the neutrino mass ordering <cit.>. Current data instead has no sensitivity to this effect. [The effect is below the  level for all multipoles, hence well beyond the reach of Planck. The effect will be below the reach of ground-based Stage-III experiments such as Advanced ACTPol <cit.>, SPT-3G <cit.>, the Simons Array <cit.> and the Simons Observatory <cit.>. It will most likely be below the reach of ground-based Stage-IV experiments such as CMB-S4 <cit.>, or next-generation satellites such as the proposed LiteBIRD <cit.>, COrE <cit.>, and PIXIE <cit.>.]


Although all the above effects may suggest that the CMB is exquisitely sensitive to the neutrino mass, in practice, the shape of the CMB anisotropy spectra is governed by several parameters, some of which are degenerate among themselves <cit.>. We refer the reader to the dedicated study of Ref. <cit.> (see also <cit.>).

To assess the impact of massive neutrinos on the CMB, all characteristic times, scales, and density ratios governing the shape of the CMB anisotropy spectrum should be kept fixed, i.e. keeping z_eq and the angular diameter distance to last-scattering d_A(z_dec) fixed. This would result in: a decrease in the late integrated Sachs-Wolfe (LISW) effect, which however is poorly constrained owing to the fact that the relevant multipole range is cosmic variance limited; a modest change in the diffusion damping scale for M_ν≳ 0.6 eV; and finally, a Δ C_ℓ/C_ℓ∼ -(M_ν/0.1   eV)% depletion of the amplitude of the C_ℓ's for 20 ≲ℓ≲ 200, due to a smaller EISW effect, which also contains a sub- effect due to the individual neutrino masses, essentially impossible to detect.



  §.§.§ Baseline combinations of datasets used, and their definitions, I.


Measurements of the CMB temperature, polarization, and cross-correlation spectra from the Planck 2015 data release <cit.> are included. We consider a combination of the high-ℓ (30 ≤ℓ≤ 2508) TT likelihood, as well as the low-ℓ (2 ≤ℓ≤ 29) TT likelihood based on the CMB maps recovered with : we refer to this combination as PlanckTT. We furthermore include the Planck polarization data in the low-ℓ (2 ≤ℓ≤ 29) likelihood, referring to it as lowP. Our baseline model, consisting of a combination of PlanckTT and lowP, is referred to as base.

In addition to the above, we also consider the high-ℓ (30 ≤ℓ≤ 1996) EE and TE likelihood, which we refer to as highP. In order to ease the comparison of our results to those previously presented in the literature, we shall add high-ℓ polarization measurements to our baseline model separately, referring to the combination of base and highP as basepol. For the purpose of clarity, we have summarized our nomenclature of datasets and their combinations in Tab. <ref>.

All the measurements described above are analyzed by means of the publicly available Planck likelihoods <cit.>. [http://www.cosmos.esa.int/web/planck/plawww.cosmos.esa.int/web/planck/pla] When considering a prior on the optical depth to reionization τ we shall only consider the TT likelihood in the multipole range 2 ≤ℓ≤ 29. We do so for avoiding double-counting of information, see Sec. <ref>. Of course, these likelihoods depend also on a number of nuisance parameters, which should be (and are) marginalized over. These nuisance parameters describe, for instance, residual foreground contamination, calibration, and beam-leakage (see Refs. <cit.>).

CMB measurements have been complemented with additional probes which will help breaking the parameter degeneracies discussed. These additional datasets include large-scale structure probes and direct measurements of the Hubble parameter, and will be described in what follows. We make the conservative choice of not including lensing potential measurements, despite measuring M_ν via lensing potential reconstruction is the expected target of the next-generation CMB experiments. This choice is dictated by the observation that lensing potential measurements via reconstruction through the temperature 4-point function are known to be in tension with the lensing amplitude as constrained by the CMB power spectra through the A_lens parameter <cit.> (see also <cit.> for relevant work).



 §.§ Galaxy power spectrum


Once CMB data is used to fix the other cosmological parameters, the galaxy power spectrum could in principle be the most sensitive cosmological probe of massive neutrinos among those exploited here. Sub-eV neutrinos behave as a hot dark matter component with large thermal velocities, clustering only on scales below the neutrino free-streaming wavenumber k_fs <cit.>:


    k_fs≃ 0.018  Ω_m^1/2 (M_ν/1  eV )^1/2 h   Mpc^-1  .


On scales below the free-streaming scale (or, correspondingly, for wavenumbers larger than the free-streaming wavenumber), neutrinos cannot cluster as their thermal velocity exceeds the escape velocity of the gravitational potentials on those scales. Conversely, on scales well above the free-streaming scale, neutrinos behave as cold dark matter after the transition to the non-relativistic regime. Massive neutrinos leave their imprint on the galaxy power spectrum in several different ways:


  * For wavenumbers k>k_fs, the power spectrum in the linear perturbation regime is subject to a scale-independent reduction by a factor of (1-f_ν)^2, where f_ν≡Ω_ν/Ω_m is defined as the ratio of the energy content in neutrinos to that in matter <cit.>.

  * In addition, the power-spectrum for wavenumbers k>k_fs is further subject to a scale-dependent step-like suppression, starting at k_fs and saturating at k ∼ 1  h   Mpc^-1. This suppression is due to the absence of neutrino perturbations in the total matter power spectrum, ultimately due to the fact that neutrinos do not cluster on scales k>k_fs. At k ∼ 1  h   Mpc^-1, the suppression reaches a constant amplitude of Δ P(k)/P(k) ≃ -10f_ν <cit.> (the amplitude of the suppression is independent of redshift, however see the point below).

  * The growth rate of the dark matter perturbations is reduced from δ∝ a to δ∝ a^1-3/5f_ν, due to the absence of gravitational back-reaction effects from free-streaming neutrinos. The redshift dependence of this suppression implies that this effect could be disentangled from that of a similar suppression in the primordial power spectrum by measuring the galaxy power spectrum at several redshifts, which amounts to measuring the time-dependence of the neutrino mass effect <cit.>.

  * On very large scales (10^-3 < k < 10^-2), the matter power spectrum is enhanced by the presence of massive neutrinos <cit.>.

  * As in the case of the EISW effect in the CMB, the step-like suppression in the matter power spectrum carries a non-trivial dependence on the individual neutrino masses, as it depends on the time of the transition to the non-relativistic regime for each neutrino mass eigenstate <cit.> (k_fs∝ m_ν_i^1/2), and thus is in principle extremely sensitive to the neutrino mass hierarchy. However, the effect is very small and very hard to measure, even with the most ambitious next-generation large-scale structure surveys <cit.>. Through the same effect, the lensed CMB as well as the lensing potential power spectrum could also be sensitive to the neutrino mass hierarchy.


Notice that, in principle, once CMB data is used to fix the other cosmological parameters, the galaxy power spectrum could be the most sensitive probe of neutrino masses. In practice, the potential of this dataset is limited by several effects. Galaxy surveys have access to a region of k-space k_min<k<k_max where the step-like suppression effect is neither null nor maximal. The minimum wavenumber accessible is limited both by signal-to-noise ratio and by systematics effects, and is typically of order k ∼ 10^-2 h   Mpc^-1, meaning that the fourth effect outlined above is currently not appreciable. The maximum wavenumber accessible is instead limited by the reliability of the non-linear predictions for the matter power spectrum.

At any given redshift, there exists a non-linear wavenumber, above which the galaxy power spectrum is only useful insofar as one is able to model non-linear effects, redshift space distortions, and the possible scale-dependence of the bias (a factor relating the spatial distribution of galaxies and the underlying dark matter density field <cit.>) correctly. The non-linear wavenumber depends not only on the redshift of the sample but also on other characteristics of the sample itself (e.g. whether the galaxies are more or less massive). At the present time, the non-linear wavenumber is approximately k=0.15  h   Mpc^-1, whereas for the galaxy sample we will consider (DR12 CMASS, at an effective redshift of z=0.57, see footnote 4 for the definition of effective redshift) we will show that wavenumbers smaller than k=0.2  h   Mpc^-1 are safe against large non-linear corrections (see also Fig. <ref>, where the galaxy power spectrum has been evaluated for M_ν=0  eV given that the Coyote emulator adopted <cit.> does not fully implement corrections due to non-zero neutrino masses on small scales, and Ref. <cit.>). [The effective redshift consists of the weighted mean redshift of the galaxies of the sample, with the weights described in <cit.>.]

The issue of the scale-dependent bias is indeed more subtle than it might seem, given that neutrinos themselves induce a scale-dependent bias <cit.>. A parametrization of the galaxy power spectrum in the presence of massive neutrinos in terms of a scale-independent bias and a shot-noise component [see Eq.(<ref>)], which in itself adds two extra nuisance parameters, may not capture all the relevant effects at play. Despite these difficulties, the galaxy power spectrum is still a very useful dataset as it helps breaking some of the degeneracies present with CMB-only data, in particular by improving the determination of Ω_mh^2 and n_s, the latter being slightly degenerate with M_ν. Moreover, as we shall show in this paper, the galaxy power spectrum represents a conservative dataset (see Sec. <ref>).

Nonetheless, a great deal of effort is being invested into determining the scale-dependent bias from cosmological datasets. There are several promising routes towards achieving this, for instance through CMB lensing, galaxy lensing, cross-correlations of the former with galaxy or quasar clustering measurements, or higher order correlators of the former datasets, see e.g. Refs. <cit.>. A sensitivity on M_ν of 0.023 eV has been forecasted from a combination of Planck CMB measurements together with weak lensing shear auto-correlation, galaxy auto-correlation, and galaxy-shear cross-correlation from Euclid <cit.>, after marginalization over the bias, with the figure improving to 0.01 eV after including a weak lensing-selected cluster sample from Euclid <cit.>. Similar results are expected to be achieved for certain configurations of the proposed WFIRST survey <cit.>. It is worth considering that the sensitivity of these datasets would be substantially boosted by determining the scale-dependent bias as discussed above.

A conservative cut-off in wavenumber space, required in order to avoid non-linearities when dealing with galaxy power spectrum data, denies access to the modes where the signature of non-zero M_ν is greatest, i.e. those at high k where the free-streaming suppression effect is most evident. One is then brought to question the usefulness of such data when constraining M_ν. Actually, the real power of P(k) rests in its degeneracy breaking ability, when combined with CMB data. For example, P(k) data is extremely useful as far as the determination of certain cosmological parameters  is concerned (e.g. n_s, which is degenerate with M_ν).

The degeneracy breaking effect of P(k), however, is most evident when in combination with CMB data. As an example, let us consider what is usually referred to as the most significant effect of non-zero M_ν on P(k), that is, a step-like suppression of the small-scale power spectrum. This effect is clearest when one increases M_ν while fixing (Ω_m,h). However, as we discussed in Sec. <ref>, the impact of non-zero M_ν on CMB data is best examined fixing Θ_s. If one adjusts h in order to keep Θ_s fixed, and in addition keeps Ω_bh^2 and Ω_ch^2 fixed, the power spectrum will be suppressed on both large and small scales, i.e. the result will be a global increase in amplitude <cit.>. In other words, this reverses the fourth effect listed above. This is just an example of the degeneracy breaking power of P(k) data in combination with CMB data.

Galaxy clustering measurements are addressed by means of the Sloan Digital Sky Survey III (SDSS-III; <cit.>) Baryon Oscillation Spectroscopic Survey (BOSS; <cit.>) DR12 <cit.>. The SDSS-III BOSS DR12 CMASS sample covers an effective  volume of V_eff≈ 7.4   Gpc^3 <cit.>. It contains 777202 massive galaxies in the range 0.43 < z < 0.7, at an effective redshift z = 0.57 (see footnote 4 for the definition of effective redshift), covering 9376.09 deg^2 over the sky. Here we consider the spherically averaged power spectrum of this sample, as measured by Gil-Marín et al. in <cit.>. We refer to this dataset as P(k). The measured galaxy power spectrum P_meas^g consists of a convolution of the true galaxy power spectrum P_true^g with a window function W(k_i,k_j), which accounts for correlations between the measurements at different scales due to the finite size of the survey geometry:


    P_meas^g(k_i) = ∑_j W(k_i,k_j)P_true^g(k_j)


Thus, at each step of the Monte Carlo, we need to convolve the theoretical galaxy power spectrum P_th at the given point in the parameter space with the window function, before comparing it with the measured galaxy power spectrum and constructing the likelihood.

Following previous works <cit.>, we model the theoretical galaxy power spectrum as:


    P_th = b_HF^2P_HFν^m(k,z)+P_HF^s   ,


where P_HFν^m denotes the matter power spectrum calculated at each step by the Boltzmann solver , corrected for non-linear effects using the  method <cit.>. We make use of the modified version of  designed by <cit.> to improve the treatment of non-linearities in the presence of massive neutrinos. In order to reduce the impact of non-linearities we impose the conservative choice of considering a maximum wavenumber k_max = 0.2  h   Mpc^-1. As we show in Fig. <ref> (for M_ν=0  eV), this region is safe against uncertainties due to non-linear evolution, and is also convenient for comparison with other works which have adopted a similar maximum wavenumber cutoff. The smallest wavenumber we are considering is instead of k_min = 0.03  h   Mpc^-1, and is determined by the control over systematics, which dominate at smaller wavenumbers. The parameters b_HF and P_HF^s denote the scale-independent bias and the shot noise contributions: the former reflects the fact that galaxies are biased tracers of the underlying dark matter distribution, whereas the latter arises from the discrete point-like nature of the galaxies as tracers of the dark matter. We impose flat priors in the range [0.1,10] and [0,10000] respectively for b_HF and P_HF^s.

Although in this simple model the bias and shot noise are assumed to be scale-independent, there is no unique prescription for the form of these quantities. In particular, concerning the bias, several theoretically well-motivated scale-dependent functional forms exist in the literature (such as the Q model of <cit.>, that of <cit.>, or that of <cit.> motivated by local primordial non-Gaussianity). It is beyond the scope of our paper to explore the impact of different bias function choices on the neutrino mass bounds. Instead, we simply note that it is not necessarily true that increasing the number of parameters governing the bias shape may result in broader constraints. Indeed, tighter constraints on M_ν may arise in some of the bias parameterizations with more than one parameter involved, because they might have comparable effects on the power spectrum.





 §.§ Baryon acoustic oscillations


Prior to the recombination epoch, photons and baryons in the early Universe behave as a tightly coupled fluid, whose evolution is determined by the interplay between the gravitational pull of potential wells, and the restoring force due to the large pressure of the radiation component. The resulting pressure waves which set up, before freezing at recombination, imprint a characteristic scale on the late-time matter clustering, in the form of a localized peak in the two-point correlation function, or a series of smeared peaks in the power spectrum. This scale corresponds to the sound horizon at the drag epoch, denoted by r_s(z_drag), where the drag epoch is defined as the time when baryons were released from the Compton drag of photons, see Ref. <cit.>. Then, r_s(z_drag) takes the form:


    r_s(z_drag) = ∫_z_drag^∞ dz  c_s(z)/H(z)  ,


where c_s(z) denotes the sound speed and is given by c_s(z) = c/√(3(1+R)), with R=3ρ_b/4ρ_r being the ratio of the baryon to photon momentum density. Finally, the baryon drag epoch z_drag is defined as the redshift such that the baryon drag optical depth τ_drag is equal to one:


    τ_drag(η_drag) = 4/3Ω_r/Ω_b∫_0^z_drag dz  dη/daσ_T x_e(z)/1+z = 1   ,


where σ_T=6.65 × 10^-29  m^2 denotes the Thomson cross-section and x_e(z) represents the fraction of free electrons.

BAO measurements contain geometrical information in the sense that, as a “standard ruler” of known and measured length, they allow for the determination of the angular diameter distance to the redshift of interest, and hence make it possible to map out the expansion history of the Universe after the last scattering. In addition, they are affected by uncertainties due to the non-linear evolution of the matter density field to a lesser extent than the galaxy power spectrum, making them less prone to systematic effects than the latter. An angle-averaged BAO measurement constrains the quantity D_v(z_eff)/r_s(z_drag), where the dilation scale D_v at the effective redshift of the survey z_eff is a combination of the physical angular diameter distance D_A(z) and the Hubble parameter H(z) (which control the radial and the tangential separations within a given cosmology, respectively):


    D_v(z) =  [ (1+z)^2D_A(z)^2cz/H(z) ]^1/3  .


D_v quantifies the dilation in distances when the fiducial cosmology is modified. The power of the BAO technique resides on its ability of resolving the existing degeneracies present when the CMB data alone is used, in particular in sharpening the determination of Ω_m and of the Hubble parameter H_0, discarding the low values of H_0 allowed by the CMB data. 

Massive neutrinos affect both the low-redshift geometry and the growth of structure, and correspondingly BAO measurements. If we increase M_ν, while keeping Ω_bh^2 and Ω_ch^2 fixed, the expansion rate at early times is increased, although only for M_ν>0.6  eV. Therefore, in order to keep fixed the angular scale of the sound horizon at last scattering Θ_s (which is very well constrained by the CMB acoustic peak structure), it is necessary to decrease Ω_Λ. As Ω_Λ decreases, it is found that H(z) decreases for z < 1 <cit.>. It can be shown that an increase in M_ν has a negligible effect on r_s(z_drag). Hence, we conclude that the main effect of massive neutrinos on BAO measurements is to increase  D_v(z)/r_s(z_drag) and decrease H_0, as M_ν is increased (see <cit.>). It is worth noting that there is no parameter degeneracy which can cancel the effect of a non-zero neutrino mass on BAO data alone, as far as the minimal ΛCDM+M_ν extended model is concerned <cit.>.



  §.§.§ Baseline combinations of datasets used, and their definitions, II.


In this work, we make use of BAO measurements extracted from a number of galaxy surveys. When using BAO measurements in combination with the DR12 CMASS P(k), we consider data from the Six-degree Field Galaxy Survey (6dFGS) <cit.>, the WiggleZ survey <cit.>, and the DR11 LOWZ sample <cit.>, as done in <cit.>. We refer to the combination of these three BAO measurements as BAO. When combining BAO with the base CMB dataset and the DR12 CMASS P(k) measurements, we refer to the combination as basePK. When combining BAO with the basepol CMB dataset and the DR12 CMASS P(k) measurements, we refer to the combination as basepolPK. Recall that we have summarized our nomenclature of datasets (including baseline datasets) and their combinations in Tab. <ref>.

The 6dFGS data consists of a measurement of r_s(z_drag)/D_V(z) at z = 0.106 (as per the discussion above, r_s/D_V decreases as M_ν is increased). The WiggleZ data instead consist of measurements of the acoustic parameter A(z) at three redshifts: z = 0.44, z = 0.6, and z = 0.73, where the acoustic parameter is defined as:


    A(z) = 100D_v(z)√(Ω_mh^2)/cz  .


Given the effect of M_ν on D_v(z), A(z) will increase as M_ν increases. Finally, the DR11 LOWZ data consists of a measurement of D_v(z)/r_s(z_drag) (which increases as M_ν is increased) at z = 0.32. 

Since the BAO feature is measured from the galaxy two-point correlation function, to avoid double counting of information, when considering the base and basepol datasets we do not include the DR11 CMASS BAO measurements, as the DR11 CMASS and DR12 CMASS volumes overlap. However, if we drop the DR12 CMASS power spectrum from our datasets, we are allowed to add DR11 CMASS BAO measurements without this leading to double-counting of information. Therefore, for completeness, we consider this case as well. Namely, we drop the DR12 CMASS power spectrum from our datasets, replacing it with the DR11 CMASS BAO measurement. This consists of a measurement of D_v(z_eff)/r_s(z_drag) at z_eff = 0.57.



  §.§.§ Baseline combinations of datasets used, and their definitions, III.


We refer to the combination of the four BAO measurements (6dFGS, WiggleZ, DR11 LOWZ, DR11 CMASS) as BAOFULL. We instead refer to the combination of the base CMB and the BAOFULL datasets with the nomenclature baseBAO. When high-ℓ polarization CMB data is added to this baseBAO dataset, the combination is referred to as basepolBAO, see Tab. <ref>. The comparison between basePK and baseBAO, as well as between basepolPK and basepolBAO, gives insight into the role played by large-scale structure datasets in constraining neutrino masses. In particular, it allows for an assessment of the relative importance of shape information in the form of the power spectrum against geometrical information in the form of BAO measurements when deriving the neutrino mass bounds. For clarity, all the denominations of the combinations of datasets we consider are summarized in Tab. <ref>.

All the BAO measurements used in this work are tabulated in Tab. <ref>. Note that we do not include BAO measurements from the DR7 main galaxy sample <cit.> or from the cross-correlation of DR11 quasars with the Lyα forest absorption <cit.>, and hence our results are not directly comparable to other existing studies which included these measurements.





 §.§ Hubble parameter measurements


Direct measurements of H_0 are very important when considering bounds on M_ν. With CMB data alone, there exists a strong degeneracy between M_ν and H_0 (see e.g. <cit.>). When M_ν is varied, the distance to last scattering changes as well. Defining ω_b ≡Ω_bh^2, ω_c ≡Ω_ch^2, ω_m ≡Ω_mh^2, ω_r ≡Ω_rh^2, ω_ν≡Ω_νh^2, within a flat Universe, this distance is given by:


    χ = c∫_0^z_decdz/√(ω_r(1+z)^4+ω_m(1+z)^3+  ( 1-ω_m/h^2 ))  ,


where ω_m = ω_c + ω_b + ω_ν. The structure of the CMB acoustic peaks leaves little freedom in varying ω_c and ω_b. Therefore, for what concerns the distance to the last scattering, a change in M_ν can be compensated essentially only by a change in h or, in other words, by a change in H_0. This suggests that M_ν and H_0 are strongly anti-correlated: the effect on the CMB of increasing M_ν can be easily compensated by a decrease in H_0, and vice versa.

In light of the above discussion, we expect a prior on the Hubble parameter to help pinning down the allowed values of M_ν from CMB data. Here, we consider two different priors on the Hubble parameter. The first prior we consider is based on a reanalysis of an older measurement based on the Hubble Space Telescope, the original measurement being H_0 = (73.8 ± 2.4)   km  s^-1 Mpc^-1 <cit.>. The original measurement showed a ∼ 2.4σ tension with the value of H_0 derived from fitting CMB data <cit.>. The reanalysis, conducted by Efstathiou in Ref. <cit.>, used the revised geometric maser distance to NGC4258 of Ref. <cit.> as a distance anchor. This reanalysis obtains a more conservative value of H_0 = (70.6 ± 3.3)   km  s^-1 Mpc^-1, which agrees with the extracted H_0 value from CMB-only within 1σ. We refer to this prior as H070p6.

The second prior we consider is based on the most recent HST 2.4% determination of the Hubble parameter in Ref. <cit.>. This measurement benefits from more than twice the number of Cepheid variables used to calibrate luminosity distances, with respect to the previous analysis <cit.>, as well as from improved determinations of distance anchors. The measured value of the Hubble parameter is H_0 = (73.02 ± 1.79)   km  s^-1 Mpc^-1, which is in tension with the CMB-only H_0 value by 3σ. We refer to the corresponding prior as H073p02. [We do not include here the latest 3.8% determination of H_0 by the H0LiCOW program. The measurement, based on gravitational time delays of three multiply-imaged quasar systems, yields H_0 = 71.9^+2.4_-3.0  km  s^-1 Mpc^-1 <cit.>.]

A consideration is in order at this point. Given the strong degeneracy between M_ν and H_0, we expect the introduction of the two aforementioned priors (especially the H073p02 one) to lead to a tighter bound on M_ν. At the same time, we expect this bound to be less reliable and/or robust. In other words, such a bound would be quite artificial, as it would be driven by a combination of the tension between direct and primary CMB determinations of H_0 and the strong M_ν-H_0 degeneracy. We can therefore expect the fit to degrade when any of the two aforementioned priors is introduced. We nonetheless choose to include these prior for a number of reasons. Firstly, the underlying measurement in <cit.> has attracted significant attention and hence it is worth assessing its impact on bounds on M_ν, subject to the strict caveats we discussed, in light of its potential to break the M_ν-H_0 degeneracy. Next, our results including the H_0 priors will serve as a warning of the danger of adding datasets which are inconsistent between each other.



 §.§ Optical depth to reionization


The first generation of galaxies ended the dark ages of the Universe. These galaxies emitted UV photons which gradually ionized the neutral hydrogen which had rendered the Universe transparent following the epoch of recombination, in a process known as reionization (see e.g. Ref. <cit.> for a review). So far, it is not entirely clear when cosmic reionization took place. Cosmological measurements can constrain the optical depth to reionization τ, which, assuming instantaneous reionization (a very common useful approximation), can be related to the redshift of reionization z_re.

Early CMB measurements of τ from WMAP favored an early-reionization scenario  (z_re = 10.6 ± 1.1 in the instantaneous reionization approximation <cit.>), requiring the presence of sources of reionization at z≳ 10. This result was in tension with observations of Ly-α emitters at z≃ 7 (see e.g. <cit.>), that suggest that reionization ended by z≃ 6. However, the results delivered by the Planck collaboration in the 2015 public data release, using the large-scale (low-ℓ) polarization observations of the Planck Low Frequency Instrument (LFI) <cit.> in combination with Planck temperature and lensing data, indicate that τ = 0.066 ± 0.016 <cit.>, corresponding to a significantly lower value for the redshift of instantaneous reionization: z_re = 8.8^+1.2_-1.1 (see also <cit.> for an assessment of the role of the cleaning procedure on the lower estimate of τ, and <cit.> for an alternative indirect method for measuring large-scale polarization and hence constrain τ using only small-scale and lensing polarization maps), and thus reducing the need for high-redshift sources of reionization <cit.>.

The optical depth to reionization is a crucial quantity when considering constraints on the sum of neutrino masses, the reason being that there exist degeneracies between τ and M_ν (see e.g. <cit.>). If we consider CMB data only (focusing on the TT spectrum), an increase in M_ν, which results in a suppression of structure, reduces the smearing of the damping tail. This effect can be compensated by an increase in τ. Due to the well-known degeneracy between A_s and τ from CMB temperature data (which is sensitive to the combination A_se^-2τ), the value of A_s should also be increased accordingly. However, the value of A_s also determines the overall amplitude of the matter power spectrum, which is furthermore affected by the presence of massive neutrinos, which reduce the small-scale clustering. If, in addition to TT data, low-ℓ polarization measurements are considered, the degeneracy between A_s and τ will be largely alleviated and, consequently, also the multiple ones among the A_s, τ, and  M_ν cosmological parameters.

Recently, the Planck collaboration has identified, modeled, and removed previously unaccounted systematic effects in large angular scale polarization data from the Planck High Frequency Instrument (HFI) <cit.> (see also <cit.>). Using the new HFI low-ℓ polarization likelihood (that has not been made publicly available by the Planck collaboration), the constraints on τ have been considerably improved, with a current determination of τ = 0.055 ± 0.009 <cit.>, entirely consistent with the value inferred from LFI.

In this work, we explore the impact on the constraints on M_ν of adding a prior on τ. Specifically, we impose a Gaussian prior on the optical depth to reionization of τ = 0.055 ± 0.009, consistent with the results reported in <cit.>. We refer to this prior as τ0p055. We expect this prior to tighten our bounds on M_ν. However, a prior on τ is a proxy for low-ℓ polarization spectra (low-ℓ C_ℓ^EE, C_ℓ^BB, and C_ℓ^TE). Therefore, as previously stated, when adding a prior on τ, we remove the low-ℓ polarization data from our datasets, in order to avoid double-counting information, while keeping low-ℓ temperature data.



 §.§ Planck SZ clusters


The evolution with mass and redshift of galaxy clusters offers a unique probe of both the physical matter density, Ω_m, and the present amplitude of density fluctuations, characterized by the root mean squared of the linear overdensity in spheres of radius 8  h^-1Mpc, σ_8, for a review see e.g. <cit.>. Both quantities are of crucial importance when extracting neutrino mass bounds from large-scale structure, due to the neutrino free-streaming nature.

CMB measurements are able to map galaxy clusters via the Sunyaev-Zeldovich (SZ) effect, which consists of an energy boost to the CMB photons, which are inverse Compton re-scattered by hot electrons (see e.g. <cit.>). Therefore, the thermal SZ effect imprints a spectral distortion to CMB photons traveling along the cluster line of sight. The distortion consists of an increase in intensity for frequencies higher than 220 GHz, and a decrease for lower frequencies.

We shall here make us of cluster counts from the latest Planck SZ clusters catalogue, consisting of 439 clusters detected via their SZ signal <cit.>. We refer to the dataset as SZ. The cluster counts function is given by the number of clusters of a certain mass M within a redshift range [z,z+dz], i.e. dN/dz:


    dN/dz|_M > M_min = f_skydV(z)/dz∫_M_min^∞ dM  dn/dM(M,z)   .


The dependence on the underlying cosmological model is encoded in the differential volume dV/dz:


    dV(z)/dz = 4π/H(z)∫_0^z dz'  1/H^2(z')  ,


through the dependence of the Hubble parameter H(z) on the basic cosmological parameters, and further through the dependence of the cluster mass function dn/dM (calculated through N-body simulations) on the parameters Ω_m and σ_8.

The largest source of uncertainty in the interpretation of cluster counts measurements resides in the masses of clusters themselves, which in turn can be inferred by X-ray mass proxies, relying however on the assumption of hydrostatic equilibrium. This assumption can be violated by bulk motion or non-thermal sources of pressure, leading to biases in the derived value of the cluster mass. Further systematics in the X-ray analyses can arise e.g. due to instrument calibration or the temperature structure in the gas. Therefore, it is clear that determinations of cluster masses carry a significant uncertainty, with a typical Δ M/M∼ 10-20%, quantified via the cluster mass bias parameter, 1-b:


    M_X = (1-b)M_500  ,


where M_X denotes the X-ray extracted cluster mass, and M_500 the true halo mass, defined as the total mass within a sphere of radius R_500, R_500 being the radius within which the mean overdensity of the cluster is 500 times the critical density at that redshift.

As the cluster mass bias 1-b is crucial in constraining the values of Ω_m and σ_8, and hence the normalization of the matter power spectrum, it plays an important role when constraining M_ν. We impose an uniform prior on the cluster mass bias in the range [0.1,1.3], as done in Ref. <cit.>, in which it is shown that this choice of 1-b leads to the most stringent bounds on the neutrino mass. There exist as well independent lensing measurements of the cluster mass bias, as those provided by the Weighing the Giants project <cit.>, by the Canadian Cluster Comparison Project <cit.>, and by CMB lensing <cit.> (see also Ref. <cit.>). However, we shall not make use of 1-b priors based on these independent measurements, as the resulting value of σ_8 is in slight tension, at the level of 1-2σ, with primary CMB measurements (however, see <cit.>).

The value of σ_8 indicated by weak lensing measurements is smaller than that derived from CMB-only datasets, favoring therefore quite large values of M_ν, large enough to suppress the small-scale clustering in a significant way. Therefore, we restrict ourselves to the case in which the cluster mass bias is allowed to freely vary between 0.1 and 1.3. It has been shown in <cit.> that this choice leads to robust and unbiased neutrino mass limits. In this way, the addition of the SZ dataset can be considered truly reliable.



§ RESULTS ON M_Ν


We begin here by analyzing the results obtained for the different datasets and their combinations, assessing their robustness. The constraining power of geometrical versus shape large-scale structure datasets will be discussed in Sec. <ref>. In Sec. <ref> we apply the method of <cit.> and described in Sec. <ref> to quantify the exclusion limits on the inverted hierarchy given the bounds on M_ν presented in the following. The 95% C.L. upper bounds on M_ν we obtain are summarized in Tabs. <ref>, <ref>, <ref>, <ref>. The C.L.s at which our most constraining datasets disfavor the Inverted Hierarchy, CL_ IH, obtained through our analysis in Sec. <ref>, are reported in Tab. <ref>.

Table <ref> shows the results for the more conservative approach when considering CMB data; namely, by neglecting high-ℓ polarization data. The limits obtained when the base dataset is considered are very close to those quoted in Ref. <cit.>, where a three degenerate neutrino spectrum with a lower prior on M_ν of 0.06 eV was assumed, whereas we have taken a lower prior of 0 eV. Our choice is driven by the goal of obtaining independent bounds on M_ν from cosmology alone, making the least amount of assumptions. This different choice of prior is the reason for the (small) discrepancy in our 95% C.L. upper limit on M_ν (0.716 eV) and the limit found in Ref. <cit.> (0.754 eV), and, in general, in all the bounds we shall describe in what follows. That is, these discrepancies are due to differences in the volume of the parameter space explored. When P(k) data are added to the base, CMB-only dataset, the neutrino mass limits are considerably improved, reaching M_ν <0.299 eV at 95% C.L..

The limits reported in Table <ref>, while being consistent with those presented in Ref. <cit.> (obtained with an older BOSS full shape power spectrum measurement, the DR9 CMASS P(k)), are slightly less constraining. We attribute this mild slight loss of constraining power to the fact that the DR12 P(k) appears slightly suppressed on small scales with respect to the DR9 P(k), see Fig. <ref>. This fact, already noticed for previous data releases, can ultimately be attributed to a very slight change in power following an increase in the mean galaxy density over time due to the tiling (observational) strategy of the survey <cit.>. The changes are indeed very small, and the broadband shape of the power spectra for different data releases in fact agree very well within error bars. A small suppression in small-scale power, nonetheless, is expected to favor higher values of M_ν, which help explaining the observed suppression, and this explains the slight difference between our results and those of Ref. <cit.>.

While the addition of external datasets, such as a prior on τ or Planck SZ cluster counts, leads to mild improvements in the constraints on M_ν, the tightest bounds are obtained when considering the H073p02 prior on the Hubble parameter, due to the large existing degeneracy between H_0 and M_ν at the CMB level, and only partly broken via P(k) or BAO measurements. However, as previously discussed, this H073p02 measurement shows a significant tension with CMB estimates of the Hubble parameter. [See e.g. Refs. <cit.> for recent works examining this discrepancy and possible solutions.] Therefore, the 95% C.L. limits on M_ν of  <0.164, <0.140, <0.136 eV for the basePK+H073p02, basePK+H073p02+τ0p055 and basePK+H073p02+τ0p055+SZ cases should be regarded as the most aggressive limits one can obtain when considering a prior on H_0 and neglecting high-ℓ polarization data. Indeed, when using the H070p6 prior, a less constraining limit of M_ν <0.219 eV at 95% C.L. is obtained in the basePK+H070p06 case, value that is closer to the limits obtained when additional measurements (not related to H_0 priors) are added to the basePK data combination.

The tension between the H073p02 measurement and primary CMB determinations of H_0 implies that the very strong bounds obtained using such prior are also the least robust and/or reliable. They are almost entirely driven by the aforementioned tension in combination with the strong M_ν-H_0 degeneracy, and hence are somewhat artificial. We expect in fact the quality of the fit to deteriorate in the presence of 2 inconsistent datasets (that is, CMB spectra and H_0 prior). To quantify the worsening in fit, we compute the Δχ^2 associated to the bestfit, for a given combination of datasets before and after the addition of the H_0 prior. For example, for the basePK dataset combination, we find Δχ^2 ≡χ^2_min(basePK+H073p02)-χ^2_min(basePK)=+5.2, confirming as expected a substantial worsening in fit when the H073p02 prior is added to the basePK dataset. The above observation reinforces the fact that any bound on M_ν obtained using the H073p02 prior should be interpreted with considerable caution, as such bound is most likely artificial.

Table <ref> shows the equivalent to Tab. <ref> but including high-ℓ polarization data. Notice that the limits are considerably tightened. As previously discussed, the tightest bounds are obtained when the H073p02 prior is considered. For instance, we obtain M_ν<0.109 eV at 95% C.L. from the 
basepolPK+H073p02+τ0p055 data combination. We caution once more against the very tight bounds obtained with the H073p02 being most likely artificial. This is confirmed for example by the Δχ^2_min=+6.4 between the basepolPK+H073p02 and basepolPK datasets.



 §.§ Geometric vs shape information


In the following, we shall compare the constraining power of geometrical probes in the form of BAO measurements versus shape probes in the form of power spectrum measurements. For that purpose, we shall replace here the DR12 CMASS P(k) and the BAO datasets by the BAOFULL dataset, which consists of BAO measurements from the BOSS DR11 (both CMASS and LOWZ samples) survey, the 6dFGS survey, and the WiggleZ survey, see Tab. <ref> for more details. The main results of this section are summarized in Tabs. <ref> and <ref>, as well as Figs. <ref> and <ref>.
 
Table <ref> shows the equivalent to the third, fourth, sixth, eighth and ninth rows of Tab. <ref>, but with the shape information from the BOSS DR12 CMASS spectrum replaced by the geometrical BAO information from the BOSS DR11 CMASS measurements. Firstly, we notice that all the geometrical bounds are, in general, much more constraining than the shape bounds, as previously studied and noticed in the literature (see e.g  <cit.>, see also <cit.> for recent studies on the subject). These studies have shown that, within the minimal ΛCDM+M_ν scenario, BAO measurements provide tighter constraints on M_ν than data from the full power spectrum shape. Nevertheless, it is very important to assess whether these previous findings still hold with the improved statistics and accuracy of today's large-scale structure data (see the recent Ref. <cit.> for the expectations from future galaxy surveys).

We confirm that this finding still holds with current data. Therefore, current analyses methods of large-scale structure datasets are such that these are still sensitive to massive neutrinos through background rather than perturbation effects, despite the latter are in principle a much more sensitive probe of the effect of massive neutrinos on cosmological observables. However, as we mentioned earlier, this behaviour could be reverted once we are able to determine the amplitude and scale-dependence of the galaxy bias through CMB lensing, cosmic shear, galaxy clustering measurements, and their cross-correlations (see e.g. <cit.>). 

Moreover, it is also worth reminding that BAO measurements do include non-linear information through the reconstruction procedure, whereas the same information is prevented from being used in the power spectrum measurements due to the cutoff we imposed at k=0.2  h   Mpc^-1. In order to fully exploit the constraining power of shape measurements, improvements in our analyses methods are necessary: in particular, it is necessary to improve our understanding of the non-linear regime of the galaxy power spectrum in the presence of massive neutrinos, as well as further our understanding of the galaxy bias at a theoretical and observational level.



The addition of shape measurements requires at least two additional nuisance parameters, which in our case are represented by the bias and shot noise parameters. These two parameters relate the measured galaxy power spectrum to the underlying matter power spectrum, the latter being what one can predict once cosmological parameters are known. [Moreover, at least another nuisance parameter is required in order to account for systematics in the measured galaxy power spectrum, although the impact of this parameter is almost negligible, as we have checked (see Refs. <cit.>)]. The prescription we adopted relating the galaxy to the matter power spectrum is among the simplest choices. However, it is not necessarily true that more sophisticated choices with more nuisance parameters would further degrade the constraining power of shape measurements, particularly if we were to obtain a handle on the functional form of the scale-dependent bias <cit.>. On the other hand, it remains true that the possibility of benefiting from a large number of modes by increasing the value of k_max (which remains one of the factors limiting the constraining power of shape information compared to geometrical one) would require an exquisite knowledge of non-linear corrections, a topic which is the subject of many recent investigations particularly in the scenario where massive neutrinos are present, see e.g. <cit.>. The conclusion, however, remains that improvements in our current analyses methods, as well as further theoretical and modeling advancements, are necessary to exploit the full constraining power of shape measurements (see also <cit.>).

Finally, we notice that, even without considering the high-ℓ polarization data, we obtain the very constraining bound of M_ν<0.114 eV at 95% C.L. for the baseBAO+H073p02+τ0p055+SZ datasets. We caution again against the artificialness of bounds obtained using the H073p02 prior, as the tension with primary CMB determinations in H_0 leads to a degradation in the quality of fit. Nonetheless, even without considering the H_0 prior, we still obtain a very constraining bound of M_ν<0.151 eV at 95% C.L. In any case, results adopting these dataset combinations contribute to reinforcing the previous (weak) cosmological hints favouring the NH scenario <cit.>.

Table <ref> shows the equivalent to Tab. <ref> but with the high-ℓ polarization dataset included, i.e. adding the highP Planck dataset in the analyses. We note that the results are quite impressive, and it is interesting to explore how far could one currently get in pushing the neutrino mass limits by means of the most aggressive and least conservative datasets. The tightest limits we find are M_ν<0.093 eV at 95% C.L. using the basepolBAO+H073p02+τ0p055+SZ dataset, well below the minimal mass allowed within the IH. Therefore, within the less-conservative approach illustrated here, especially due to the use of the H073p02 prior, there exists a weak preference from present cosmological data for a normal hierarchical neutrino mass scheme. Neglecting the information from the H073p02 prior, which leads to an artificially tight bound as previously explained, the preference turns out to be weaker (M_ν<0.118  eV from the basepolBAO+τ0p055 dataset combination) but still present.



We end with a consideration, stemming from the observation that with our current analyses methods BAO measurements are more constraining than full-shape power spectrum ones. This suggests that, despite uncertainties in the modeling of the galaxy power spectrum due to the unknown absolute scale of the latter (in other words, the size of the bias) and non-linear evolution, the galaxy power spectrum actually represents a conservative dataset given that the bounds on M_ν obtained using the corresponding BAO dataset are considerably tighter.

In the remainder of the Section we will be concerned with providing a proper quantification of the statistical significance at which we can disfavor the IH, performing a simple but rigorous model comparison analysis.



 §.§ Exclusion limits on the inverted hierarchy


Here we apply the method of <cit.> and described in Sec. <ref> to determine the statistical significance at which the inverted hierarchy is disfavored given the bounds on M_ν just obtained. Our results are summarized in Tab. <ref>. In order to quantify the exclusion limits on the inverted hierarchy, we apply Eq. (<ref>) to our most constraining dataset combinations, where the criterion for choosing these datasets will be explained below.

Note that in Eq. (<ref>) we set p(N)=p(I)=0.5. That is, we assign equal priors to NH and IH, which not only is a reasonable choice when considering only cosmological datasets <cit.>, but is also the most uninformative and most conservative choice when there is no prior knowledge about the hierarchies. In any case, the formalism we adopt would allow us to introduce informative prior information on the two hierarchies, i.e. p(N) ≠ p(I) ≠ 0.5. It would in this way be possible to include information from oscillation experiments, which suggest a weak preference for the normal hierarchy due to matter effects (see e.g. <cit.>). Including this weak preference does not significantly affect our results, precisely because the current sensitivity to the neutrino mass hierarchy from both cosmology and oscillation experiments is extremely weak (see also e.g. <cit.>).

We choose to only report the statistical significance at which the IH is discarded for the most constraining dataset combinations, that is, those which disfavor the IH at >70% C.L.: we have checked that threshold for reaching a ≈ 70% C.L. exclusion limit of the IH is reached by datasets combinations which disfavor at 95% C.L. values of M_ν greater than ≈ 0.12  eV. In fact, the most constraining bound within our conservative scheme, obtained through the baseBAO+τ0p055 combination (thus disfavoring datasets which exhibit some tension with CMB or galaxy clustering measurements, for a 95% C.L. upper limit on M_ν of 0.151 eV), falls short of this threshold, and is only able to disfavor the IH at 64% C.L., providing posterior odds for NH versus IH of 1.8:1.

The hierarchy discrimination is improved when small-scale polarization is added to the aforementioned datasets combination, or when the H073p02 prior (and eventually SZ cluster counts) are added to the same datasets combination, leading to a 71% C.L. and 72% C.L. exclusion of the IH respectively. Similar levels of statistical significance for the exclusion of the IH are reached when the datasets combinations basepolPK+H073p02+τ0p055, basepolPK+H073p02+τ0p055+SZ, and basepolBAO+H073p02 are considered, leading to 74% C.L., 71% C.L., and 72% C.L. exclusion of the IH respectively. However, it is worth reminding once more that the latter figures relied on the addition of the H073p02 prior, which leads to less reliable bounds. It is also worth noting that our most constraining datasets combination(s), that is, basepolBAO+H073p02+τ0p055(+SZ), only provide a 77% C.L. exclusion of the IH.

Our findings are totally consistent with those of <cit.> and suggest that an improved sensitivity of cosmological datasets is required in order to robustly disfavor the IH, despite current datasets are already able to substantially reduce the volume of parameter space available within this mass ordering. In fact, it has been argued in <cit.> that a sensitivity of at least ≈ 0.02  eV is required in order to provide a 95% C.L. exclusion of the IH. Incidentally, not only does such a sensitivity seem within the reach of post-2020 experiments <cit.>, but it would also provide a detection of M_ν at a significance of at least 3σ, unless non-trivial late-Universe effects are at play (see e.g. <cit.>).



 §.§ Bounds on M_ν in extended parameter spaces: a brief discussion



Thus far we have explored bounds on M_ν within the assumption of a flat background ΛCDM cosmology. We have used different dataset combinations, and have identified the baseBAO dataset (leading to an upper limit of M_ν<0.186  eV) combination as being the one providing one of the strongest bounds while at the same time being one of the most robust to systematics and tensions between datasets.

However, we expect the bounds on M_ν to degrade if we were to open the parameter space: that is, if we were to vary additional parameters other than the 6 base ΛCDM parameters and M_ν. While there is no substantial indication for the need to extend the base set of parameters of the ΛCDM model (see e.g. <cit.>), one is nonetheless legitimately brought to wonder about the robustness of the obtained bounds against extended parameter spaces.

While a detailed study belongs to a follow-up paper in progress <cit.>, we nonetheless decide to present two examples of bounds on M_ν within minimally extended parameter spaces. That is, we allow in one case the dark energy equation of state w to vary within the range [-3,1] (parameter space denoted by ΛCDM+M_ν+w), and in the other case the curvature energy density Ω_k to vary freely within the range [-0.3,0.3] (parameter space denoted by ΛCDM+M_ν+Ω_k). Both parameters are known to be relatively strongly degenerate with M_ν and hence we can expect our allowing them to vary to lead to less stringent bounds on M_ν. In both cases we consider for simplicity the baseBAO dataset, for the reasons described above: therefore, the corresponding bound within the ΛCDM+M_ν parameter space to which we should compare our results to is M_ν<0.186   eV at 95% C.L., as reported in the first row of Tab. <ref>.

For the ΛCDM+M_ν+w extension, where we leave the dark energy equation of state w free to vary within the range [-3,1], we can expect the bounds on M_ν to broaden due to a well-known degeneracy between M_ν and w <cit.>. Specifically, an increase in M_ν can be compensated by a decrease in w, due to the mutual degeneracy with Ω_m. Our results confirm this expectation. With the baseBAO data combination we find M_ν<0.313  eV at 95% C.L., and w=-1.08^+0.09_-0.08 at 68% C.L., with a correlation coefficient between M_ν and w of -0.56. [The correlation coefficient between two parameters i and j (in this case i=M_ν, j=w)is defined as R=C_ij/√(C_iiC_jj), with C the covariance matrix of cosmological parameters.] The degeneracy between M_ν and w is clearly visible in the triangle plot of Fig. <ref>.



For the ΛCDM+M_ν+Ω_k extension, where we leave the curvature energy density Ω_k free to vary within the range [-0.3,0.3], we can again expect the bounds on M_ν to broaden due to the three-parameter geometric degeneracy between h, Ω_νh^2 and Ω_k <cit.>. For the baseBAO data combination we find M_ν<0.299   eV at 95% C.L., and Ω_k=0.001^+0.003_-0.004 at 68% C.L., with a correlation coefficient between M_ν and Ω_k of 0.60. The degeneracy between M_ν and Ω_k is clearly visible in the triangle plot of Fig. <ref>.



A clarification is in order here: when leaving the dark energy equation of state w and the curvature energy density Ω_k free to vary, it would be extremely useful to add supernovae data, given that these are extremely sensitive to these two quantities. We have however chosen not to do so in order to ease comparison with the bound M_ν<0.186  eV obtained for the same baseBAO combination within the ΛCDM+M_ν parameter space. Moreover, in this way we are able to reach a conservative conclusion concerning the robustness of M_ν bounds to the ΛCDM+M_ν+w and ΛCDM+M_ν+Ω_k parameter spaces, as the addition of supernovae data would lead to tighter bounds than the M_ν<0.313  eV and M_ν<0.299  eV quoted.

Of course, as expected, the bounds on M_ν degrade the moment we consider extended parameter spaces. Given our discussion in Sec. <ref>, this means within the extended parameter spaces considered the preference for one hierarchy over another essentially vanishes. However, the last statement is not necessarily always true: for instance, in certain models of dynamical dark energy with specific functional forms of w(z), the constraints on M_ν can get tighter: an example is the holographic dark energy model, within which bounds on M_ν have been shown to be substantially tighter than within a ΛCDM Universe <cit.>. An interesting thing to note, however, is that within better than 1σ uncertainties (i.e. within ∼ 68% C.L.), both w and Ω_k are compatible with the values to which they are fixed within the minimal ΛCDM+M_ν parameter space, that is, -1 and 0 respectively.



§ CONCLUSIONS


Neutrino oscillation experiments provide information on the two mass splittings governing the solar and atmospheric neutrino transitions, but are unable to measure the total neutrino mass scale, M_ν. The sign of the largest mass splitting, the atmospheric mass gap, remains unknown. The two resulting possibilities are the so-called normal (positive sign) or inverted (negative sign) mass hierarchies. While in the normal hierarchy scheme neutrino oscillation results set the minimum allowed total neutrino mass M_ν to be approximately equal to ∼0.06 eV, in the inverted one this lower limit is ∼0.1 eV.

Currently, cosmology provides the tightest bounds on the total neutrino mass M_ν, i.e. on the sum of the three active neutrino states. If these cosmological bounds turned out to be robustly and significantly smaller than the minimum allowed in the inverted hierarchy, then one would indeed determine the neutrino mass hierarchy via cosmological measurements. In order to prepare ourselves for the hierarchy extraction, an assessment of the cosmological neutrino mass limits, studying their robustness against different priors and assumptions concerning the neutrino mass distribution among the three neutrino mass eigenstates, is mandatory. Moreover, the development and application of rigorous model comparison methods to assess the preference for one hierarchy over the other is necessary. In this work, we have analyzed some of the most recent publicly available datasets to provide updated constraints on the sum of the three active neutrino masses, M_ν, from cosmology.

One very interesting aspect is whether the information concerning the total neutrino mass from the large-scale structure of the universe in its geometrical form (i.e. via the BAO signature) supersedes that of full-shape measurements of the power spectrum. While previous studies have addressed the question with former galaxy clustering datasets, it is timely to explore the situation with current galaxy catalogs, covering much larger volumes, benefiting from smaller error-bars and also from improved, more accurate descriptions of the mildly non-linear regime in the matter power spectrum.

We find that, despite the latest measurements of the galaxy power spectrum cover a vast volume of our universe, the BAO signature extracted from comparable datasets is still more powerful than the full-shape information, within the minimal ΛCDM+M_ν model studied here. This statement is expected to change within the context of extended cosmological models, such as those with non-zero curvature or a time-dependent dark energy equation of state, and we reserve this study to future work <cit.> (whereas a short discussion on the robustness of the bounds on M_ν within extended parameter spaces is provided in Appendix B).

The reason for the supremacy of BAO measurements over shape information is due to the cutoff in k-space imposed when treating the power spectrum. This cutoff is required to avoid the impact of non-linear evolution. It is worth reminding once more that BAO measurements contain non-linear information wrapped in with the reconstruction procedure. This same non-linear information cannot be used in the power spectrum due to the choice of the conservative cutoff in k-space. Moreover, the need for at least two additional nuisance parameters relating the galaxy power spectrum to the underlying matter power spectrum further degrades the constraining power of the latter. Therefore, the stronger constraints obtained through geometrical rather than shape measurements should not be seen as a limitation of the constraining power of the latter, rather as a limitation of methods currently used to analyze these datasets. A deeper understanding of the non-linear regime of the galaxy power spectrum in the presence of massive neutrinos, as well as further understanding of the galaxy bias at a theoretical and observational level, are required: it is worth noting that a lot of effort is being invested into tackling these issues.

Finally, in this work we have presented the tightest up-to-date neutrino mass constraints among those which can be found in the literature. Neglecting the debated prior on the Hubble constant of H_0 = (73.02 ± 1.79)   km  s^-1 Mpc^-1, the tightest 95% C.L. upper bound we find is M_ν<0.151 eV (assuming a degenerate spectrum), from CMB temperature anisotropies, BAO and τ measurements. Adding Planck high-ℓ polarization data tightens the previous bound to M_ν<0.118 eV. Further improvements are possible if a prior on the Hubble parameter is also added. In this less conservative approach, the 95% C.L. neutrino mass upper limit is brought down to the level of ∼ 0.09 eV, indicating a weak preference for the normal neutrino hierarchy due to volume effects. Our work also suggests that we can identify a restricted set of conservative but robust datasets: this includes CMB temperature data, as well as BAO measurements and galaxy power spectrum data, after adequate corrections for non-linearities. These datasets allow us to identify a robust upper bound of ∼ 0.15  eV on M_ν from cosmological data alone.

In addition to providing updated bounds on the total neutrino mass, we have also performed a simple but robust model comparison analysis, aimed at quantifying the exclusion limits on the inverted hierarchy from current datasets. Our findings indicate that, despite the very stringent upper bounds we have just outlined, current data is not able to conclusively favor the NH over the IH. Within our most conservative scheme, we are able to disfavor the IH with a significance of at most 64% C.L., corresponding to posterior odds of NH over IH of 1.8:1. Even the most constraining and less conservative datasets combinations are able at most to disfavor the IH at 77% C.L., with posterior odds of NH against IH of 3.3:1. This suggests that further improvements in sensitivity, down to the level of 0.02  eV, are required in order for cosmology to conclusively disfavor the IH. Fortunately, it looks like a combination of data from near-future CMB experiments and galaxy surveys should be able to reach this target.

We conclude that our findings, while unable to robustly disfavor the inverted neutrino mass ordering, significantly reduce the volume of parameter space allowed within this mass hierarchy. The more robustly future bounds will be able to disfavor the region of parameter space with M_ν > 0.1  eV, the more the IH will be put under pressure with respect to the NH. In other words future cosmological data, in the absence of a neutrino mass detection, are expected to reinforce the current mild preference for the normal hierarchy mass ordering. On the other hand, if the underlying mass hierarchy is the inverted one, a cosmological detection of the neutrino mass scale could be quick approaching. In any case, we expect neutrino cosmology to remain an active and exciting field of discovery in the upcoming years.



§ APPENDIX A: THE 3DEG APPROXIMATION


Throughout the paper we have presented bounds within the 3deg approximation of a neutrino mass spectrum with three massive degenerate mass eigenstates. The choice was motived, as discussed in Sec. <ref>, by the observations that the NH and IH mass splittings have a tiny effect on cosmological data, when compared to the 3deg approximation with the same value of the total mass M_ν. Here we discuss the conditions under which this approximation is mathematically speaking valid. We also briefly discuss why the 3deg approximation is nonetheless physically accurate given the sensitivity of current data.

Mathematically speaking, the 3deg approximation is valid as long as:


    m_0 ≫| m_i-m_j |    ,   ∀ i,j=1,2,3   ,


where m_0=m_1  [m_3] in the NH [IH] scenario (see Sec. <ref> for the definition of the labeling of the three mass eigenstates). Recall that, according to our convention, m_1<m_2<m_3    [m_3<m_1<m_2] in the NH [IH]. Therefore, the 3deg approximation is strictly speaking valid when the absolute neutrino mass scale is much larger than the individual mass splittings. A good candidate for a figure of merit to quantify the goodness of the 3deg approximation can then be obtained by considering the ratio of any given mass difference, over a quantity proportional to the absolute neutrino mass scale. This leads us to consider the following figure(s) of merit:


    ζ_ij≡3| m_i-m_j |/M_ν  ,


where the indices i,j run over i,j=1,2,3. The figures of merit ζ_ij quantify the goodness of the 3deg approximation. In the case where the 3deg approximation were exact (which, of course, is physically impossible given the non-zero mass-squared splittings), one would have ζ_ij=0. The 3deg approximation, then, can be considered valid from a practical point of view as long as ζ_ij is sufficiently small, where the amount of deviation from ζ_ij=0 one can tolerate defines what is sufficiently small and hence the validity criterion for the 3deg approximation.

In Fig. <ref> we plot our figure(s) of merit ζ_ij, for i,j=1,2 (red) and i,j=1,3 (blue) in Eq. (<ref>) and for the NH (solid) and IH (dashed) scenarios (see the caption for details), against the total neutrino mass M_ν. We plot the same quantities, but this time against the lightest neutrino mass m_0 = m_1  [m_3] for the NH [IH], in Fig. <ref>. As we discussed previously, the 3deg approximation would be exact if ζ_ij=0 (which of course cannot be displayed due to the choice of a logarithmic scale for the y axis).





As we already discussed, the decision of whether or not 3deg is a sensible approximation mathematically speaking depends on the amount of deviation from ζ_ij=0 that can be tolerated. As an example, from Fig. <ref> and Fig. <ref> we see that, considering an indicative value of M_ν≈ 0.15  eV, the value of ζ_13≈ 0.4, indicating a ≈ 40% deviation from the exact 3deg scenario, which can hardly be considered small.

This indicates that, within the remaining allowed region of parameter space, the 3deg approximation is mathematically speaking not valid. It is worth remarking that there is a degree of residual model dependency as this conclusion was reached taking at face value the indicative upper limit on M_ν of ≈ 0.15  eV, which has been derived under the assumption of a flat ΛCDM background. One can generically expect the bounds we obtained to be loosened to some extent if considering extended cosmological scenarios (although this needs not necessarily always be the case).

A different issue is, instead, whether the 3deg approximation is physically appropriate, given the sensitivity of current and near-future experiments. The issue has been discussed extensively in the literature, and in particular in some recent works <cit.>. It has been argued that, if M_ν > 0.1  eV, future cosmological observations, while measuring M_ν with high accuracy, will not be able to discriminate between the NH and the IH. In any case, cosmological measurements in combination with laboratory experiments will in this case (M_ν > 0.1  eV) play a key role in unravelling the hierarchy <cit.>. If M_ν < 0.1  eV, most of the discriminatory power in cosmological data between the NH and the IH is essentially due to volume effects: i.e., the fact that oscillation data force ≃0.1 in the IH, implying that the IH has access to a reduced volume of parameter space with respect to the NH.

Another example of the goodness of the 3deg approximation is provided in <cit.> considering a combination of forecasts for COrE, Euclid, and DESI data. Specifically, <cit.> considered a fiducial mock dataset generated implementing the full NH or IH, and then studied whether fitting the fiducial dataset using the 3deg approximation rather than the “true" NH or IH would lead to substantial biases. The findings suggest that, apart from small O(0.1σ) reconstruction biases (which can be removed for M_ν<0.1  eV), the 3deg approximation is able to recover the fiducial value of M_ν (as long as the free parameter is taken to be consistently either M_ν or m_0). This suggests that even with near-future cosmological data the 3deg approximation will still be sufficiently accurate for the purpose of estimating cosmological parameters, and further validates the goodness of the 3deg approximation in our work.

The conclusion is that current cosmological datasets are sensitive to the total neutrino mass M_ν rather than to the individual masses m_i, implying that the 3deg approximation is sufficiently precise for the purpose of obtaining reliable cosmological neutrino mass bounds for the time being. On the other hand, for future high precision cosmological data, which could benefit from increased sensitivity and could reliably have access to non-linear scales of the matter power spectrum, modelling the mass splittings correctly will matter.

In conclusion, although the 3deg approximation is not, mathematically speaking, valid in the remaining volume of parameter space, it is physically speaking a good approximation given the sensitivity of current datasets. However, quantitative claims about disfavoring the inverted hierarchy have to be drawn with care, making use of rigorous model comparison methods.



§ APPENDIX B: THE 1MASS APPROXIMATION


As argued in a number of works, the ability to robustly reach an upper bound on M_ν of ≈ 0.1  eV translates more or less directly into the ability of excluding the inverted hierarchy at a certain statistical significance, as we quantified in Sec. <ref>. In this case it is desirable to check whether one's conclusions are affected by assumptions on the underlying neutrino mass spectrum. Throughout our paper we have presented bounds on M_ν making the assumption of a spectrum of three massive degenerate neutrinos, denoted 3deg. As we have argued extensively (see e.g. Appendix A), given the sensitivity of current data, this assumption does not to any significant extent influence the resulting bounds. Nonetheless, it is interesting and timely to investigate the dependence of neutrino mass bounds under assumptions of different mass spectra, which was recently partly done in <cit.>.

Here, as in <cit.>, we consider (in addition to the 3deg spectrum) the approximation spectrum featuring a single massive eigenstate carrying the total mass M_ν together with two massless species. We refer to this scheme by the name 1mass:


    m_1 = m_2 = 0   , m_3 = M_ν   (1mass)   .


The motivation for the 1mass choice is twofold: i) it is the usual approximation adopted when performing cosmological analyses with the total neutrino mass fixed to =0.06, in order to mimic the minimal mass scenario in the case of the NH (m_1=0, m_2≪ m_3), and ii) it might provide a better description of the underlying neutrino mass ordering in the  M_ν<0.1 mass region, in which m_1∼ m_2≪ m_3, although a complete assessment goes beyond the scope of our work. The latter is the main motivation for exploring the 1mass approximation further, given the recent weak cosmological hints favoring the NH.

Before proceeding, it is useful to clarify why we have chosen to focus on results within the 3deg scheme. As we discussed, it has been observed that the impact of the NH and IH mass splittings on cosmological data is tiny if one compares the 3deg approximation to the corresponding NH and IH models with the same value of the total mass M_ν. However, this does not necessarily hold when the comparison is made between 3deg and 1mass, because the latter always has two pure dark radiation components (see footnote 8 for a definition of dark radiation) throughout the whole expansion history and, in particular, at the present time (on the other hand, NH and IH can have at most one pure radiation component at present time, a situation which occurs in the minimal mass scenario when m_0 = 0  eV and thus only for one specific point in neutrino mass parameter space) [Dark radiation consists of any weakly or non-interacting extra radiation component of the Universe, see e.g. <cit.> for a review and <cit.> for recent relevant work in connection to neutrino physics. For example, sterile neutrinos may in some models have contributed as dark radiation, see e.g. <cit.>, or possibly thermally produced cosmological axions <cit.>. Dark radiation might also arise in dark sectors with additional relativistic degrees of freedom which decouple from the Standard Model as, for instance, hidden photons (see e.g. <cit.>).]. The extra massless component(s) present in the 1mass case, but not in the NH and IH (1mass features only one extra component compared to the NH and IH if these happen to correspond to the minimal mass scenario where m_0 = 0  eV; if m_0 ≠ 0  eV, 1mass possesses two extra massless components), are known to have a non-negligible impact on cosmological observables, in particular the CMB anisotropy spectra <cit.>.

Let us now discuss how the bounds on M_ν change when passing from the 3deg to the 1mass approximation. We observe that when considering the base dataset combinations, and extensions thereof (i.e. the combinations considered in Tab. <ref>, where we report the 3deg results), the bounds obtained within the 1mass approximation are typically more constraining than the 3deg ones, by about ∼ 2-8%. For example, the 95% C.L. upper bound on M_ν is tightened from 0.716  eV to 0.658  eV for the base combination, from 0.299  eV to 0.293  eV for the base+P(k) combination, and from 0.246  eV to 0.234  eV for the basePK combination. When small-scale polarization data is added (see Tab. <ref> for the 3deg results), we observe a reversal in this behaviour: that is, the bounds obtained within the 1mass approximation are looser than the 3deg ones. For example, the 95% C.L. upper bound on M_ν is loosened from 0.485  eV to 0.619  eV for the basepol combination, from 0.275  eV to 0.300  eV for the basepol+P(k) combination, and from 0.215  eV to 0.228  eV for the basepolPK combination.

Regarding the baseBAO and basepolBAO dataset combinations and extensions thereof (see Tabs. <ref>, <ref> for the 3deg results), no clear trend emerges when passing from the 3deg to the 1mass approximation, although we note that the bounds typically degrade slightly: for example, the 95% C.L. upper bound on M_ν is loosened from 0.186  eV to 0.203  eV for the baseBAO combination, and from 0.153  eV to 0.155  eV for the basepolBAO combination.

We choose not to further investigate the reason behind these tiny but noticeable shifts because, as previously stated, the 1mass distribution is less “physical", owing to the presence of two unphysical dark radiation states. Instead, we report these numbers in the interest of noticing how these shifts suggest that, at present, cosmological measurements are starting to become sensitive (albeit in a very weak manner) to the late-time hot dark matter versus dark radiation distribution among the neutrino mass eigenstates, a conclusion which had already been reached in <cit.>.

One of the reasons underlying the choice of studying the 1mass approximation is that this scheme might represent an useful approximation to the minimal mass scenario in the NH. Of course, the possibility that the underlying neutrino hierarchy is inverted is far from being excluded. This raises the question of whether an analogous scheme, which we refer to as 2mass (already studied in <cit.>), might instead approximate the minimal mass scenario in the IH:


    m_3 = 0   , m_1 = m_2 = M_ν/2    (2mass)   .


Of course, the previously discussed considerations concerning the non-physicality of the 1mass approximation (due to the presence of extra pure radiation components) automatically apply to the 2mass approximation as well. Moreover, we note that bounds on M_ν obtained within the 2mass approximation (which features one pure radiation state) are always intermediate between those of the 3deg (which features no pure radiation state) and the 1mass (which features two pure radiation states) ones (see also e.g. <cit.>). This confirms once more that the discrepancy between bounds within these three different approximations are to be attributed to the impact of the unphysical pure radiation states on cosmological observables, in particular the CMB anisotropy spectra. In conclusion, we remark once more that, while the 3deg approximation is sufficiently accurate given the precision of current data, other approximations which introduce non-physical pure radiation states, such as the 1mass and 2mass ones, are not. Adopting these to obtain bounds on M_ν might instead lead to unphysical shifts in the determination of cosmological parameters, and hence should be avoided.













SV, EG, and MG acknowledge Hector Gil-Marín for very useful discussions. SV and OM thank Antonio Cuesta for valuable correspondence. We are very grateful to the anonymous referee for a detailed and constructive report which enormously helped to improve the quality of our paper. 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 Michigan Center for Theoretical Physics. M.G., S.V. and K.F. acknowledge support by the Vetenskapsrådet (Swedish Research Council) through contract No. 638-2013-8993 and the Oskar Klein Centre for Cosmoparticle Physics. M.L. acknowledges support from ASI through ASI/INAF Agreement 2014-024-R.1 for the Planck LFI Activity of Phase E2. 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. E.G. is supported by NSF grant AST1412966. S.H. acknowledges support by NASA-EUCLID11-0004, NSF AST1517593 and NSF AST1412966.






99

nobel
  “The Nobel Prize in Physics 2015”,
  Nobelprize.org, Nobel Media AB 2014, Web, 18 Oct 2016,
  http://www.nobelprize.org/nobel_prizes/physics/laureates/2015www.nobelprize.org/nobel_prizes/physics/laureates/2015

superk
  Y. Fukuda et al. [Super-Kamiokande Collaboration],
  
  Phys. Rev. Lett.  81 (1998) 1562
  
  [hep-ex/9807003].
  
  
  
sno
  Q. R. Ahmad et al. [SNO Collaboration],
  
  Phys. Rev. Lett.  89 (2002) 011301
  
  [nucl-ex/0204008].
  
  
  
kamland
  T. Araki et al. [KamLAND Collaboration],
  
  Phys. Rev. Lett.  94 (2005) 081801
  
  [hep-ex/0406035].
  
  
  
minos
  P. Adamson et al. [MINOS Collaboration],
  
  Phys. Rev. Lett.  101 (2008) 131802
  
  [arXiv:0806.2237 [hep-ex]].
  
  
  
dayabay
  F. P. An et al. [Daya Bay Collaboration],
  
  Phys. Rev. Lett.  108 (2012) 171803
  
  [arXiv:1203.1669 [hep-ex]].
  
  
  
reno
  J. K. Ahn et al. [RENO Collaboration],
  
  Phys. Rev. Lett.  108 (2012) 191802
  
  [arXiv:1204.0626 [hep-ex]].
  
  
  
doublechooz
  Y. Abe et al. [Double Chooz Collaboration],
  
  Phys. Rev. D 86 (2012) 052008
  
  [arXiv:1207.6632 [hep-ex]].
  
  
  
t2k
  K. Abe et al. [T2K Collaboration],
  
  Phys. Rev. Lett.  112 (2014) 061802
  
  [arXiv:1311.4750 [hep-ex]].
  
  
  
fit1 
  M. C. Gonzalez-Garcia, M. Maltoni and T. Schwetz,
  
  JHEP 1411, 052 (2014)
  
  [arXiv:1409.5439 [hep-ph]].
  
fit2 
  D. V. Forero, M. Tortola and J. W. F. Valle,
  
  Phys. Rev. D 90, no. 9, 093006 (2014)
  
  [arXiv:1405.7540 [hep-ph]].
  
fit3
  I. Esteban, M. C. González-García, M. Maltoni, I. Martínez-Soler and T. Schwetz,
  
  arXiv:1611.01514 [hep-ph].
  
  
  

fit4
  F. Capozzi, E. Di Valentino, E. Lisi, A. Marrone, A. Melchiorri and A. Palazzo,
  
  Phys. Rev. D 95 (2017) no.9,  096014
  
  [arXiv:1703.04471 [hep-ph]].
  
  
  

fit5
  A. Caldwell, A. Merle, O. Schulz and M. Totzauer,
  
  arXiv:1705.01945 [hep-ph].
  
  
  
Capozzi:2016rtj
  F. Capozzi, E. Lisi, A. Marrone, D. Montanino and A. Palazzo,
  
  Nucl. Phys. B 908 (2016) 218
  
  [arXiv:1601.07777 [hep-ph]].
  
  
  

Giusarma:2014zza
  E. Giusarma, E. Di Valentino, M. Lattanzi, A. Melchiorri and O. Mena,
  
  Phys. Rev. D 90 (2014) no.4,  043507
  
  [arXiv:1403.4852 [astro-ph.CO]].
  
  
  
Palanque-Delabrouille:2015pga
  N. Palanque-Delabrouille et al.,
  
  JCAP 1511 (2015) no.11,  011
  
  [arXiv:1506.05976 [astro-ph.CO]].
  
  
  

DiValentino:2015wba
  E. Di Valentino, E. Giusarma, M. Lattanzi, O. Mena, A. Melchiorri and J. Silk,
  
  Phys. Lett. B 752 (2016) 182
  
  [arXiv:1507.08665 [astro-ph.CO]].
  
  
  
DiValentino:2015sam 
  E. Di Valentino, E. Giusarma, O. Mena, A. Melchiorri and J. Silk,
  
  Phys. Rev. D 93, no. 8, 083527 (2016)
  
  [arXiv:1511.00975 [astro-ph.CO]].
  
  

Cuesta:2015iho
  A. J. Cuesta, V. Niro and L. Verde,
  
  Phys. Dark Univ.  13 (2016) 77
  
  [arXiv:1511.05983 [astro-ph.CO]].
  
  
  
Huang:2015wrx
  Q. G. Huang, K. Wang and S. Wang,
  
  Eur. Phys. J. C 76 (2016) no.9,  489
  
  [arXiv:1512.05899 [astro-ph.CO]].
  
  
  
DiValentino:2016ikp 
  E. Di Valentino, S. Gariazzo, M. Gerbino, E. Giusarma and O. Mena,
  
  Phys. Rev. D 93, no. 8, 083523 (2016)
  
  [arXiv:1601.07557 [astro-ph.CO]].
  
  

Giusarma:2016phn
  E. Giusarma, M. Gerbino, O. Mena, S. Vagnozzi, S. Ho and K. Freese,
  
  Phys. Rev. D 94 (2016) no.8,  083522
  
  [arXiv:1605.04320 [astro-ph.CO]].
  
  
  
Lesgourgues:2006nd
  J. Lesgourgues and S. Pastor,
  
  Phys. Rept.  429 (2006) 307
  
  [astro-ph/0603494].
  
  
  
Wong:2011ip
  Y. Y. Y. Wong,
  
  Ann. Rev. Nucl. Part. Sci.  61 (2011) 69
  
  [arXiv:1111.1436 [astro-ph.CO]].
  
  
  
Lesgourgues:2012uu
  J. Lesgourgues and S. Pastor,
  
  Adv. High Energy Phys.  2012 (2012) 608515
  
  [arXiv:1212.6154 [hep-ph]].
  
  
  
Abazajian:2013oma
  K. N. Abazajian et al. [Topical Conveners: K.N. Abazajian, J.E. Carlstrom, A.T. Lee Collaboration],
  
  Astropart. Phys.  63 (2015) 66
  
  [arXiv:1309.5383 [astro-ph.CO]].
  
  

book
  J. Lesgourgues, G. Mangano, G. Miele and S. Pastor,
  “Neutrino Cosmology,”
  Cambridge, UK: Cambridge University Press (2013) 378 p
  
Lesgourgues:2014zoa
  J. Lesgourgues and S. Pastor,
  
  New J. Phys.  16 (2014) 065002
  
  [arXiv:1404.1740 [hep-ph]].
  
  
  

Archidiacono:2016lnv
  M. Archidiacono, T. Brinckmann, J. Lesgourgues and V. Poulin,
  
  JCAP 1702 (2017) no.02,  052
  
  [arXiv:1610.09852 [astro-ph.CO]].
  
  
  

Lesgourgues:2004ps
  J. Lesgourgues, S. Pastor and L. Perotto,
  
  Phys. Rev. D 70 (2004) 045016
  
  [hep-ph/0403296].
  
  
  

Pritchard:2008wy
  J. R. Pritchard and E. Pierpaoli,
  
  Phys. Rev. D 78 (2008) 065009
  
  [arXiv:0805.1920 [astro-ph]].
  
  


DeBernardis:2009di
  F. De Bernardis, T. D. Kitching, A. Heavens and A. Melchiorri,
  
  Phys. Rev. D 80 (2009) 123509
  
  [arXiv:0907.1917 [astro-ph.CO]].
  
  


Jimenez:2010ev
  R. Jiménez, T. Kitching, C. Peña-Garay and L. Verde,
  
  JCAP 1005 (2010) 035
  
  [arXiv:1003.5918 [astro-ph.CO]].
  
  
  

Wagner:2012sw
  C. Wagner, L. Verde and R. Jiménez,
  
  Astrophys. J.  752 (2012) L31
  
  [arXiv:1203.5342 [astro-ph.CO]].
  
  
  

Hannestad:2016fog
  S. Hannestad and T. Schwetz,
  
  JCAP 1611 (2016) no.11,  035
  
  [arXiv:1606.04691 [astro-ph.CO]].
  
  
  

Xu:2016ddc
  L. Xu and Q. G. Huang,
  
  arXiv:1611.05178 [astro-ph.CO].
  
  

martinanew
  M. Gerbino, M. Lattanzi, O. Mena and K. Freese,
  
  arXiv:1611.07847 [astro-ph.CO].
  
  
blennow
  M. Blennow,
  
  JHEP 1401 (2014) 139
  
  [arXiv:1311.3183 [hep-ph]].
  
  
  

planckcosmological
  P. A. R. Ade et al. [Planck Collaboration],
  
  Astron. Astrophys.  594 (2016) A13
  
  [arXiv:1502.01589 [astro-ph.CO]].
  
  
  

bayes
  R. Trotta,
  
  Contemp. Phys.  49 (2008) 71
  
  [arXiv:0803.4089 [astro-ph]].
  
  
  

trottareview
  R. Trotta,
  
  arXiv:1701.01467 [astro-ph.CO].
  
  
  
cosmomc1 
  A. Lewis and S. Bridle,
  
  Phys. Rev. D 66, 103511 (2002)
  
  [astro-ph/0205436].
  
  
  
cosmomc2 
  A. Lewis,
  
  Phys. Rev. D 87, no. 10, 103529 (2013)
  
  [arXiv:1304.4473 [astro-ph.CO]].
  
  
  
gelmanandrubin 
  S. Brooks and A. Gelman,
  
  J. Comp. Graph. Stat.  7, 434-455 (1998).
  
lewiscosmocoffee
  A. Lewis,
  post on http://cosmocoffee.info/viewtopic.php?t=350&sid=adcdc49e35eb4cbf3c69d8273c9980ebCosmoCoffee
  

Beacom:2004yd
  J. F. Beacom, N. F. Bell and S. Dodelson,
  
  Phys. Rev. Lett.  93 (2004) 121302
  
  [astro-ph/0404585].
  
  
  

bellini
  N. Bellomo, E. Bellini, B. Hu, R. Jiménez, C. Peña-Garay and L. Verde,
  
  arXiv:1612.02598 [astro-ph.CO].
  
  

joudaki
  S. Joudaki,
  
  Phys. Rev. D 87 (2013) 083523
  
  [arXiv:1202.0005 [astro-ph.CO]].
  
  
  

archidiaconogiusarmamelchiorrimena
  M. Archidiacono, E. Giusarma, A. Melchiorri and O. Mena,
  
  Phys. Rev. D 86 (2012) 043509
  
  [arXiv:1206.0109 [astro-ph.CO]].
  
  
  
feeney
  S. M. Feeney, H. V. Peiris and L. Verde,
  
  JCAP 1304 (2013) 036
  
  [arXiv:1302.0014 [astro-ph.CO]].
  
  
  
archidiacono1
  M. Archidiacono, N. Fornengo, C. Giunti, S. Hannestad and A. Melchiorri,
  
  Phys. Rev. D 87 (2013) no.12,  125034
  
  [arXiv:1302.6720 [astro-ph.CO]].
  
  
  
archidiacono2
  M. Archidiacono, S. Hannestad, A. Mirizzi, G. Raffelt and Y. Y. Y. Wong,
  
  JCAP 1310 (2013) 020
  
  [arXiv:1307.0615 [astro-ph.CO]].
  
  
  
mirizzi
  A. Mirizzi, G. Mangano, N. Saviano, E. Borriello, C. Giunti, G. Miele and O. Pisanti,
  
  Phys. Lett. B 726 (2013) 8
  
  [arXiv:1303.5368 [astro-ph.CO]].
  
  
  
verde
  L. Verde, S. M. Feeney, D. J. Mortlock and H. V. Peiris,
  
  JCAP 1309 (2013) 013
  
  [arXiv:1307.2904 [astro-ph.CO]].
  
  
  
gariazzo
  S. Gariazzo, C. Giunti and M. Laveder,
  
  JHEP 1311 (2013) 211
  
  [arXiv:1309.3192 [hep-ph]].
  
  
  
archidiacono3
  M. Archidiacono, N. Fornengo, S. Gariazzo, C. Giunti, S. Hannestad and M. Laveder,
  
  JCAP 1406 (2014) 031
  
  [arXiv:1404.1794 [astro-ph.CO]].
  
  
  
bergstrom
  J. Bergström, M. C. González-García, V. Niro and J. Salvado,
  
  JHEP 1410 (2014) 104
  
  [arXiv:1407.3806 [hep-ph]].
  
  
  
rossi
  G. Rossi, C. Yèche, N. Palanque-Delabrouille and J. Lesgourgues,
  
  Phys. Rev. D 92 (2015) no.6,  063505
  
  [arXiv:1412.6763 [astro-ph.CO]].
  
  
  

Zhang:2015rha
  J. F. Zhang, M. M. Zhao, Y. H. Li and X. Zhang,
  
  JCAP 1504 (2015) 038
  
  [arXiv:1502.04028 [astro-ph.CO]].
  
  
  
divalentino1
  E. Di Valentino, A. Melchiorri and J. Silk,
  
  Phys. Rev. D 92 (2015) no.12,  121302
  
  [arXiv:1507.06646 [astro-ph.CO]].
  
  
  
gerbino
  M. Gerbino, M. Lattanzi and A. Melchiorri,
  
  Phys. Rev. D 93 (2016) no.3,  033001
  
  [arXiv:1507.08614 [hep-ph]].
  
  
  
divalentino2
  E. Di Valentino, E. Giusarma, M. Lattanzi, O. Mena, A. Melchiorri and J. Silk,
  
  Phys. Lett. B 752 (2016) 182
  
  [arXiv:1507.08665 [astro-ph.CO]].
  
  
  
zhang
  X. Zhang,
  
  Phys. Rev. D 93 (2016) no.8,  083011
  
  [arXiv:1511.02651 [astro-ph.CO]].
  
  
  
kitching
  T. D. Kitching, L. Verde, A. F. Heavens and R. Jiménez,
  
  Mon. Not. Roy. Astron. Soc.  459 (2016) no.1,  971
  
  [arXiv:1602.02960 [astro-ph.CO]].
  
  
  

moresco
  M. Moresco, R. Jiménez, L. Verde, A. Cimatti, L. Pozzetti, C. Maraston and D. Thomas,
  
  JCAP 1612 (2016) no.12,  039
  doi:10.1088/1475-7516/2016/12/039
  [arXiv:1604.00183 [astro-ph.CO]].
  
  
  

canac
  N. Canac, G. Aslanyan, K. N. Abazajian, R. Easther and L. C. Price,
  
  JCAP 1609 (2016) no.09,  022
  
  [arXiv:1606.03057 [astro-ph.CO]].
  
  
  

archidiacono4
  M. Archidiacono, S. Gariazzo, C. Giunti, S. Hannestad, R. Hansen, M. Laveder and T. Tram,
  
  JCAP 1608 (2016) no.08,  067
  
  [arXiv:1606.07673 [astro-ph.CO]].
  
  
  

kumar
  S. Kumar and R. C. Nunes,
  
  Phys. Rev. D 94 (2016) 123511
  
  [arXiv:1608.02454 [astro-ph.CO]].
  
  
  
bouchet
  E. Di Valentino and F. R. Bouchet,
  
  JCAP 1610 (2016) no.10,  011
  
  [arXiv:1609.00328 [astro-ph.CO]].
  
  
  

Kumar:2017dnp
  S. Kumar and R. C. Nunes,
  
  arXiv:1702.02143 [astro-ph.CO].
  
  
  

Guo:2017hea
  R. Y. Guo, Y. H. Li, J. F. Zhang and X. Zhang,
  
  JCAP 1705 (2017) no.05,  040
  
  [arXiv:1702.04189 [astro-ph.CO]].
  
  
  

Zhang:2017rbg
  X. Zhang,
  
  Sci. China Phys. Mech. Astron.  60 (2017) no.6,  060431
  
  [arXiv:1703.00651 [astro-ph.CO]].
  
  
  

Li:2017iur
  E. K. Li, H. Zhang, M. Du, Z. H. Zhou and L. Xu,
  
  arXiv:1703.01554 [astro-ph.CO].
  
  
  

Yang:2017amu
  W. Yang, R. C. Nunes, S. Pan and D. F. Mota,
  
  Phys. Rev. D 95 (2017) no.10,  103522
  
  [arXiv:1703.02556 [astro-ph.CO]].
  
  
  

Feng:2017nss
  L. Feng, J. F. Zhang and X. Zhang,
  
  Eur. Phys. J. C 77 (2017) no.6,  418
  doi:10.1140/epjc/s10052-017-4986-3
  [arXiv:1703.04884 [astro-ph.CO]].
  
  
  

Dirian:2017pwp
  Y. Dirian,
  
  arXiv:1704.04075 [astro-ph.CO].
  
  
  

Feng:2017mfs
  L. Feng, J. F. Zhang and X. Zhang,
  
  arXiv:1706.06913 [astro-ph.CO].
  
  

Lorenz:2017fgo
  C. S. Lorenz, E. Calabrese and D. Alonso,
  
  arXiv:1706.00730 [astro-ph.CO].
  
  

Couchot:2017pvz
  F. Couchot, S. Henrot-Versillé, O. Perdereau, S. Plaszczynski, B. Rouillé D. 'orfeuil, M. Spinelli and M. Tristram,
  
  arXiv:1703.10829 [astro-ph.CO].
  
  


Doux:2017tsv
  C. Doux, M. Penna-Lima, S. D. P. Vitenti, J. Tréguer, E. Aubourg and K. Ganga,
  
  arXiv:1706.04583 [astro-ph.CO].
  
  

Simpson:2017qvj
  F. Simpson, R. Jimenez, C. Pena-Garay and L. Verde,
  
  JCAP 1706 (2017) no.06,  029
  
  [arXiv:1703.03425 [astro-ph.CO]].
  
  
  

Schwetz:2017fey
  T. Schwetz, K. Freese, M. Gerbino, E. Giusarma, S. Hannestad, M. Lattanzi, O. Mena and S. Vagnozzi,
  
  arXiv:1703.04585 [astro-ph.CO].
  
  
  
Lewis:2006fu 
  A. Lewis and A. Challinor,
  
  Phys. Rept.  429, 1 (2006)
  
  [astro-ph/0601594].
  
Hall:2012kg
  A. C. Hall and A. Challinor,
  
  Mon. Not. Roy. Astron. Soc.  425 (2012) 1170
  
  [arXiv:1205.6172 [astro-ph.CO]].
  
  
  

Ade:2013zuv
  P. A. R. Ade et al. [Planck Collaboration],
  
  Astron. Astrophys.  571 (2014) A16
  
  [arXiv:1303.5076 [astro-ph.CO]].
  
  
  

s41
  K. N. Abazajian et al. [Topical Conveners: K.N. Abazajian, J.E. Carlstrom, A.T. Lee Collaboration],
  
  Astropart. Phys.  63 (2015) 66
  
  [arXiv:1309.5383 [astro-ph.CO]].
  
  
  

s42
  K. N. Abazajian et al.,
  
  Astropart. Phys.  63 (2015) 55
  
  [arXiv:1309.5381 [astro-ph.CO]].
  
  
  

s43
  K. N. Abazajian et al. [CMB-S4 Collaboration],
  
  arXiv:1610.02743 [astro-ph.CO].
  
  
  

Levi:2013gra
  M. Levi et al. [DESI Collaboration],
  
  arXiv:1308.0847 [astro-ph.CO].
  
  
  

Aghamousa:2016zmz
  A. Aghamousa et al. [DESI Collaboration],
  
  arXiv:1611.00036 [astro-ph.IM].
  
  
  

Aghamousa:2016sne
  A. Aghamousa et al. [DESI Collaboration],
  
  arXiv:1611.00037 [astro-ph.IM].
  
  
  
actpol1
  E. Calabrese et al.,
  
  JCAP 1408 (2014) 010
  
  [arXiv:1406.4794 [astro-ph.CO]].
  
  
  

actpol2
  S. W. Henderson et al.,
  
  J. Low. Temp. Phys.  184 (2016) no.3-4,  772
  
  [arXiv:1510.02809 [astro-ph.IM]].
  
  
  

spt-3g
  B. A. Benson et al. [SPT-3G Collaboration],
  
  Proc. SPIE Int. Soc. Opt. Eng.  9153 (2014) 91531P
  
  [arXiv:1407.2973 [astro-ph.IM]].
  
  
  

sa
  A. Suzuki et al. [POLARBEAR Collaboration],
  
  J. Low. Temp. Phys.  184 (2016) no.3-4,  805
  
  [arXiv:1512.07299 [astro-ph.IM]].
  
  
  
so
  https://simonsobservatory.org/simonsobservatory.org
  

litebird
  T. Matsumura et al.,
  
  J. Low. Temp. Phys.  176 (2014) 733
  
  [arXiv:1311.2847 [astro-ph.IM]].
  
  
  

core1
  F. R. Bouchet et al. [COrE Collaboration],
  
  arXiv:1102.2181 [astro-ph.CO].
  
  
  

core2
  E. Di Valentino et al. [the CORE Collaboration],
  
  arXiv:1612.00021 [astro-ph.CO].
  
  
core3
  F. Finelli et al. [CORE Collaboration],
  
  arXiv:1612.08270 [astro-ph.CO].
  
  

pixie
  A. Kogut et al.,
  
  JCAP 1107 (2011) 025
  
  [arXiv:1105.2044 [astro-ph.CO]].
  
  
  
bond
  G. Efstathiou and J. R. Bond,
  
  Mon. Not. Roy. Astron. Soc.  304 (1999) 75
  
  [astro-ph/9807103].
  
  
  
howlett
  C. Howlett, A. Lewis, A. Hall and A. Challinor,
  
  JCAP 1204 (2012) 027
  
  [arXiv:1201.3654 [astro-ph.CO]].
  
  
  

Gerbino:2016sgw
  M. Gerbino, K. Freese, S. Vagnozzi, M. Lattanzi, O. Mena, E. Giusarma and S. Ho,
  
  Phys. Rev. D 95 (2017) no.4,  043512
  
  [arXiv:1610.08830 [astro-ph.CO]].
  
  
  

planckproducts
  R. Adam et al. [Planck Collaboration],
  
  Astron. Astrophys.  594 (2016) A1
  
  [arXiv:1502.01582 [astro-ph.CO]].
  
  
  

plancklikelihood
  N. Aghanim et al. [Planck Collaboration],
  
  Astron. Astrophys.  594 (2016) A11
  
  [arXiv:1507.02704 [astro-ph.CO]].
  
  
  

Calabrese:2008rt
  E. Calabrese, A. Slosar, A. Melchiorri, G. F. Smoot and O. Zahn,
  
  Phys. Rev. D 77 (2008) 123531
  
  [arXiv:0803.2309 [astro-ph]].
  
  
  

DiValentino:2015ola
  E. Di Valentino, A. Melchiorri and J. Silk,
  
  Phys. Rev. D 92 (2015) no.12,  121302
  
  [arXiv:1507.06646 [astro-ph.CO]].
  
  
  

DiValentino:2015bja
  E. Di Valentino, A. Melchiorri and J. Silk,
  
  Phys. Rev. D 93 (2016) no.2,  023513
  
  [arXiv:1509.07501 [astro-ph.CO]].
  
  
  

Addison:2015wyg
  G. E. Addison, Y. Huang, D. J. Watts, C. L. Bennett, M. Halpern, G. Hinshaw and J. L. Weiland,
  
  Astrophys. J.  818 (2016) no.2,  132
  
  [arXiv:1511.00055 [astro-ph.CO]].
  
  
  

neutrinofootprint
  R. Jiménez, C. Peña-Garay and L. Verde,
  
  Phys. Dark Univ.  15 (2017) 31
  
  [arXiv:1602.08430 [astro-ph.CO]].
  
  
  

kaiserb
  N. Kaiser,
  
  Astrophys. J.  284 (1984) L9.
  
  
  
  
  
Heitmann:2008eq
  K. Heitmann, M. White, C. Wagner, S. Habib and D. Higdon,
  
  Astrophys. J.  715 (2010) 104
  
  [arXiv:0812.1052 [astro-ph]].
  
  
  

Heitmann:2013bra
  K. Heitmann, E. Lawrence, J. Kwan, S. Habib and D. Higdon,
  
  Astrophys. J.  780 (2014) 111
  
  [arXiv:1304.7849 [astro-ph.CO]].
  
  
  

Kwan:2013jva
  J. Kwan, K. Heitmann, S. Habib, N. Padmanabhan, H. Finkel, E. Lawrence, N. Frontiere and A. Pope,
  
  Astrophys. J.  810 (2015) no.1,  35
  
  [arXiv:1311.6444 [astro-ph.CO]].
  
  
  

castorina1
  F. Villaescusa-Navarro, F. Marulli, M. Viel, E. Branchini, E. Castorina, E. Sefusatti and S. Saito,
  
  JCAP 1403 (2014) 011
  
  [arXiv:1311.0866 [astro-ph.CO]].
  
  
  

loverde
  M. LoVerde,
  
  Phys. Rev. D 93 (2016) no.10,  103526
  
  [arXiv:1602.08108 [astro-ph.CO]].
  
  
  

Raccanelli:2017kht
  A. Raccanelli, L. Verde and F. Villaescusa-Navarro,
  
  arXiv:1704.07837 [astro-ph.CO].
  
  

Gaztanaga:2011yi
  E. Gaztañaga, M. Eriksen, M. Crocce, F. Castander, P. Fosalba, P. Martí, R. Miquel and A. Cabré,
  
  Mon. Not. Roy. Astron. Soc.  422 (2012) no.4,  2904
  
  [arXiv:1109.4852 [astro-ph.CO]].
  
  
  

Hand:2013xua
  N. Hand et al.,
  
  Phys. Rev. D 91 (2015) no.6,  062001
  
  [arXiv:1311.6200 [astro-ph.CO]].
  
  


Bianchini:2014dla
  F. Bianchini et al.,
  
  Astrophys. J.  802 (2015) no.1,  64
  
  [arXiv:1410.4502 [astro-ph.CO]].
  
  


Pullen:2015vtb
  A. R. Pullen, S. Alam, S. He and S. Ho,
  
  Mon. Not. Roy. Astron. Soc.  460 (2016) no.4,  4098
  
  [arXiv:1511.04457 [astro-ph.CO]].
  
  
  

Bianchini:2015yly
  F. Bianchini et al.,
  
  Astrophys. J.  825 (2016) no.1,  24
  
  [arXiv:1511.05116 [astro-ph.CO]].
  
  
  

Kirk:2015dpw
  D. Kirk et al. [DES Collaboration],
  
  Mon. Not. Roy. Astron. Soc.  459 (2016) no.1,  21
  
  [arXiv:1512.04535 [astro-ph.CO]].
  
  
  

Pujol:2016lfe
  A. Pujol et al.,
  
  Mon. Not. Roy. Astron. Soc.  462 (2016) no.1,  35
  
  [arXiv:1601.00160 [astro-ph.CO]].
  
  
  

Singh:2016xey
  S. Singh, R. Mandelbaum and J. R. Brownstein,
  
  Mon. Not. Roy. Astron. Soc.  464 (2016) no.2,  2120
  
  arXiv:1606.08841 [astro-ph.CO].
  
  
  

Prat:2016xor
  J. Prat et al. [DES Collaboration],
  
  
  [arXiv:1609.08167 [astro-ph.CO]].
  
  

Laureijs:2011gra
  R. Laureijs et al. [EUCLID Collaboration],
  
  arXiv:1110.3193 [astro-ph.CO].
  
  
  

Carbone:2010ik
  C. Carbone, L. Verde, Y. Wang and A. Cimatti,
  
  JCAP 1103 (2011) 030
  
  [arXiv:1012.2868 [astro-ph.CO]].
  
  


Joudaki:2011nw
  S. Joudaki and M. Kaplinghat,
  
  Phys. Rev. D 86 (2012) 023526
  
  [arXiv:1106.0299 [astro-ph.CO]].
  
  
  

Carbone:2011by
  C. Carbone, C. Fedeli, L. Moscardini and A. Cimatti,
  
  JCAP 1203 (2012) 023
  
  [arXiv:1112.4810 [astro-ph.CO]].
  
  
  

Hamann:2012fe
  J. Hamann, S. Hannestad and Y. Y. Y. Wong,
  
  JCAP 1211 (2012) 052
  
  [arXiv:1209.1043 [astro-ph.CO]].
  
  


Basse:2013zua
  T. Basse, O. E. Bjælde, J. Hamann, S. Hannestad and Y. Y. Y. Wong,
  
  JCAP 1405 (2014) 021
  
  [arXiv:1304.2321 [astro-ph.CO]].
  
  
  

Spergel:2013tha
  D. Spergel et al.,
  
  arXiv:1305.5422 [astro-ph.IM].
  
  
  

Poulin:2016nat
  V. Poulin, P. D. Serpico and J. Lesgourgues,
  
  JCAP 1608 (2016) no.08,  036
  
  [arXiv:1606.02073 [astro-ph.CO]].
  
  
  
sdss
  D. J. Eisenstein et al. [SDSS Collaboration],
  
  Astron. J.  142 (2011) 72
  
  [arXiv:1101.1529 [astro-ph.IM]].
  
  
  
bolton
  A. S. Bolton et al. [Cutler Group, LP Collaboration],
  
  Astron. J.  144 (2012) 144
  
  [arXiv:1207.7326 [astro-ph.CO]].
  
  
  
dawson
  K. S. Dawson et al. [BOSS Collaboration],
  
  Astron. J.  145 (2013) 10
  
  [arXiv:1208.0022 [astro-ph.CO]].
  
  
  
smee
  S. Smee et al.,
  
  Astron. J.  146 (2013) 32
  
  [arXiv:1208.2233 [astro-ph.IM]].
  
  
  
alam
  S. Alam et al. [SDSS-III Collaboration],
  
  Astrophys. J. Suppl.  219 (2015) no.1,  12
  
  [arXiv:1501.00963 [astro-ph.IM]].
  
  
  
12
  S. Alam et al. [BOSS Collaboration],
  
  
  [arXiv:1607.03155 [astro-ph.CO]].
  
  
  
reid
  B. Reid et al.,
  
  Mon. Not. Roy. Astron. Soc.  455 (2016) no.2,  1553
  
  [arXiv:1509.06529 [astro-ph.CO]].
  
  
  

Gil-Marin:2015sqa
  H. Gil-Marín et al.,
  
  Mon. Not. Roy. Astron. Soc.  460 (2016) no.4,  4188
  
  [arXiv:1509.06386 [astro-ph.CO]].
  
  
  

acousticscale
  H. J. Seo et al.,
  
  Astrophys. J.  761 (2012) 13
  
  [arXiv:1201.2172 [astro-ph.CO]].
  
  
  

giusarmadeputterhomena
  E. Giusarma, R. de Putter, S. Ho and O. Mena,
  
  Phys. Rev. D 88 (2013) no.6,  063515
  
  [arXiv:1306.5544 [astro-ph.CO]].
  
  


Lewis:1999bs
  A. Lewis, A. Challinor and A. Lasenby,
  
  Astrophys. J.  538 (2000) 473
  
  [astro-ph/9911177].
  
  
  

halofit1
  R. E. Smith et al. [VIRGO Consortium Collaboration],
  
  Mon. Not. Roy. Astron. Soc.  341 (2003) 1311
  
  [astro-ph/0207664].
  
  


halofit2
  R. Takahashi, M. Sato, T. Nishimichi, A. Taruya and M. Oguri,
  
  Astrophys. J.  761 (2012) 152
  
  [arXiv:1208.2701 [astro-ph.CO]].
  
  
  

birdviel
  S. Bird, M. Viel and M. G. Haehnelt,
  
  Mon. Not. Roy. Astron. Soc.  420 (2012) 2551
  
  [arXiv:1109.4416 [astro-ph.CO]].
  
  
  

bqq
  S. Cole et al. [2dFGRS Collaboration],
  
  Mon. Not. Roy. Astron. Soc.  362 (2005) 505
  
  [astro-ph/0501174].
  
  
    

amendola
  L. Amendola, E. Menegoni, C. Di Porto, M. Corsi and E. Branchini,
  
  arXiv:1502.03994 [astro-ph.CO].
  
  
  

Dalal:2007cu
  N. Dalal, O. Dore, D. Huterer and A. Shirokov,
  
  Phys. Rev. D 77 (2008) 123514
  
  [arXiv:0710.4560 [astro-ph]].
  
  
  

Eisenstein:1997ik
  D. J. Eisenstein and W. Hu,
  
  Astrophys. J.  496 (1998) 605
  
  [astro-ph/9709112].
  
  
  

Hou:2011ec
  Z. Hou, R. Keisler, L. Knox, M. Millea and C. Reichardt,
  
  Phys. Rev. D 87 (2013) 083008
  
  [arXiv:1104.2333 [astro-ph.CO]].
  
  
  

spt
  Z. Hou et al.,
  
  Astrophys. J.  782 (2014) 74
  
  [arXiv:1212.6267 [astro-ph.CO]].
  
  
  
6dfgs
  F. Beutler et al.,
  
  Mon. Not. Roy. Astron. Soc.  416 (2011) 3017
  
  [arXiv:1106.3366 [astro-ph.CO]].
  
  
  
wigglez
  C. Blake et al.,
  
  Mon. Not. Roy. Astron. Soc.  418 (2011) 1707
  
  [arXiv:1108.2635 [astro-ph.CO]].
  
  
  
dr11
  L. Anderson et al. [BOSS Collaboration],
  
  Mon. Not. Roy. Astron. Soc.  441 (2014) no.1,  24
  
  [arXiv:1312.4877 [astro-ph.CO]].
  
  
  
mgs
  A. J. Ross, L. Samushia, C. Howlett, W. J. Percival, A. Burden and M. Manera,
  
  Mon. Not. Roy. Astron. Soc.  449 (2015) no.1,  835
  
  [arXiv:1409.3242 [astro-ph.CO]].
  
  
  
lymana
  A. Font-Ribera et al. [BOSS Collaboration],
  
  JCAP 1405 (2014) 027
  
  [arXiv:1311.1767 [astro-ph.CO]].
  
  
  
giusarmadeputtermena
  E. Giusarma, R. De Putter and O. Mena,
  
  Phys. Rev. D 87 (2013) no.4,  043515
  
  [arXiv:1211.2154 [astro-ph.CO]].
  
  
  

riess2011
  A. G. Riess et al.,
  
  Astrophys. J.  730 (2011) 119
   Erratum: [Astrophys. J.  732 (2011) 129]
  
  [arXiv:1103.2976 [astro-ph.CO]].
  
  
  

revisited
  G. Efstathiou,
  
  Mon. Not. Roy. Astron. Soc.  440 (2014) no.2, 1138
  
  [arXiv:1311.3461 [astro-ph.CO]].
  
  
  

ngc4258
  E. M. L. Humphreys, M. J. Reid, J. M. Moran, L. J. Greenhill and A. L. Argon,
  
  Astrophys. J.  775 (2013) 13
  
  [arXiv:1307.6031 [astro-ph.CO]].
  
  
  

riess2016
  A. G. Riess et al.,
  
  Astrophys. J.  826 (2016) no.1,  56
  
  [arXiv:1604.01424 [astro-ph.CO]].
  
  
  

bonvin
  V. Bonvin et al.,
  
  arXiv:1607.01790 [astro-ph.CO].
  
  
  
inthebeginning
  R. Barkana and A. Loeb,
  
  Phys. Rept.  349 (2001) 125
  
  [astro-ph/0010468].
  
  
  

Hinshaw:2012aka 
  G. Hinshaw et al. [WMAP Collaboration],
  
  Astrophys. J. Suppl.  208, 19 (2013)
  
  [arXiv:1212.5226 [astro-ph.CO]].
  
  


Stark:2010qj
  D. P. Stark, R. S. Ellis, K. Chiu, M. Ouchi and A. Bunker,
  
  Mon. Not. Roy. Astron. Soc.  408 (2010) 1628
  
  [arXiv:1003.5244 [astro-ph.CO]].
  
  
  

Pentericci:2014nia 
  L. Pentericci et al.,
  
  Astrophys. J.  793, no. 2, 113 (2014)
  
  [arXiv:1403.5466 [astro-ph.CO]].
  
  
  

Schenker:2014tda 
  M. A. Schenker, R. S. Ellis, N. P. Konidaris and D. P. Stark,
  
  Astrophys. J.  795, no. 1, 20 (2014)
  
  [arXiv:1404.4632 [astro-ph.CO]].
  
  


Treu:2013ida 
  T. Treu, K. B. Schmidt, M. Trenti, L. D. Bradley and M. Stiavelli,
  
  Astrophys. J.  775, L29 (2013)
  
  [arXiv:1308.5985 [astro-ph.CO]].
  
  
  

Tilvi:2014oia 
  V. Tilvi et al.,
  
  Astrophys. J.  794, no. 1, 5 (2014)
  
  [arXiv:1405.4869 [astro-ph.CO]].
  
  
  

Lattanzi:2016dzq
  M. Lattanzi et al.,
  
  JCAP 1702 (2017) no.02,  041
  
  [arXiv:1611.01123 [astro-ph.CO]].
  
  
  

Meerburg:2017lfh
  P. D. Meerburg, J. Meyers, K. M. Smith and A. van Engelen,
  
  arXiv:1701.06992 [astro-ph.CO].
  
  
mesinger
  A. Mesinger, A. Aykutalp, E. Vanzella, L. Pentericci, A. Ferrara and M. Dijkstra,
  
  Mon. Not. Roy. Astron. Soc.  446 (2015) 566
  
  [arXiv:1406.6373 [astro-ph.CO]].
  
  
  
choudhury
  T. R. Choudhury, E. Puchwein, M. G. Haehnelt and J. S. Bolton,
  
  Mon. Not. Roy. Astron. Soc.  452 (2015) no.1,  261
  
  [arXiv:1412.4790 [astro-ph.CO]].
  
  
  
  
Robertson:2015uda 
  B. E. Robertson, R. S. Ellis, S. R. Furlanetto and J. S. Dunlop,
  
  Astrophys. J.  802, no. 2, L19 (2015)
  
  [arXiv:1502.02024 [astro-ph.CO]].
  
  

  
Bouwens:2015vha 
  R. J. Bouwens, G. D. Illingworth, P. A. Oesch, J. Caruana, B. Holwerda, R. Smit and S. Wilkins,
  
  Astrophys. J.  811, no. 2, 140 (2015)
  
  [arXiv:1503.08228 [astro-ph.CO]].
  
  
  
mitrachoudhuryferrara
  S. Mitra, T. R. Choudhury and A. Ferrara,
  
  Mon. Not. Roy. Astron. Soc.  454 (2015) no.1,  L76
  
  [arXiv:1505.05507 [astro-ph.CO]].
  
  
  

Allison:2015qca
  R. Allison, P. Caucal, E. Calabrese, J. Dunkley and T. Louis,
  
  Phys. Rev. D 92 (2015) no.12,  123535
  
  [arXiv:1509.07471 [astro-ph.CO]].
  
  
  

Liu:2015txa
  A. Liu, J. R. Pritchard, R. Allison, A. R. Parsons, U. Seljak and B. D. Sherwin,
  
  Phys. Rev. D 93 (2016) no.4,  043013
  
  [arXiv:1509.08463 [astro-ph.CO]].
  
  
  

Calabrese:2016eii
  E. Calabrese, D. Alonso and J. Dunkley,
  
  arXiv:1611.10269 [astro-ph.CO].
  
  

plancktau
  N. Aghanim et al. [Planck Collaboration],
  
  arXiv:1605.02985 [astro-ph.CO].
  
  
  

planckreionization
  R. Adam et al. [Planck Collaboration],
  
  Astron. Astrophys.  596 (2016) A108
  
  arXiv:1605.03507 [astro-ph.CO].
  
  
  

clusters
  S. W. Allen, A. E. Evrard and A. B. Mantz,
  
  Ann. Rev. Astron. Astrophys.  49 (2011) 409
  
  [arXiv:1103.4829 [astro-ph.CO]].
  
  
  

sz1
  Y. B. Zeldovich and R. A. Sunyaev,
  
  Astrophys. Space Sci.  4 (1969) 301.
  
  
  
  

sz2
  R. A. Sunyaev and Y. B. Zeldovich,
  
  Astrophys. Space Sci.  7 (1970) 3.
  
  
  
sz3
  R. A. Sunyaev and Y. B. Zeldovich,
  
  Ann. Rev. Astron. Astrophys.  18 (1980) 537.
  
  
  
  

plancksz
  P. A. R. Ade et al. [Planck Collaboration],
  
  Astron. Astrophys.  594 (2016) A27
  
  [arXiv:1502.01598 [astro-ph.CO]].
  
  
  

plancksz1
  P. A. R. Ade et al. [Planck Collaboration],
  
  Astron. Astrophys.  594 (2016) A24
  
  [arXiv:1502.01597 [astro-ph.CO]].
  
  
  
wtg
  A. von der Linden et al.,
  
  Mon. Not. Roy. Astron. Soc.  443 (2014) no.3,  1973
  
  [arXiv:1402.2670 [astro-ph.CO]].
  
  
  
melin
  J. B. Melin and J. G. Bartlett,
  
  Astron. Astrophys.  578 (2015) A21
  
  [arXiv:1408.5633 [astro-ph.CO]].
  
  
  
zaldarriaga
  M. Zaldarriaga and U. Seljak,
  
  Phys. Rev. D 59 (1999) 123507
  
  [astro-ph/9810257].
  
  
  

Kitching:2016zkn
  T. D. Kitching, J. Alsing, A. F. Heavens, R. Jiménez, J. D. McEwen and L. Verde,
  
  arXiv:1611.04954 [astro-ph.CO].
  
  


Dvorkin:2014lea
  C. Dvorkin, M. Wyman, D. H. Rudd and W. Hu,
  
  Phys. Rev. D 90 (2014) no.8,  083503
  
  [arXiv:1403.8049 [astro-ph.CO]].
  
  
  
cuestapc
  A. Cuesta, private communication.
  

h01
  A. Pourtsidou and T. Tram,
  
  Phys. Rev. D 94 (2016) no.4,  043518
  
  [arXiv:1604.04222 [astro-ph.CO]].
  
  
  

h02
  S. Grandis, D. Rapetti, A. Saro, J. J. Mohr and J. P. Dietrich,
  
  Mon. Not. Roy. Astron. Soc.  463 (2016) no.2, 1416,
  
  arXiv:1604.06463 [astro-ph.CO].
  
  


h03
  E. Di Valentino, A. Melchiorri and J. Silk,
  
  Phys. Lett. B 761 (2016) 242
  
  [arXiv:1606.00634 [astro-ph.CO]].
  
  
  

h05
  Q. G. Huang and K. Wang,
  
  Eur. Phys. J. C 76 (2016) no.9,  506
  
  [arXiv:1606.05965 [astro-ph.CO]].
  
  
  

h06
  B. L'Huillier and A. Shafieloo,
  
  JCAP 1701 (2017) no.01,  015
  
  [arXiv:1606.06832 [astro-ph.CO]].
  
  
  

h07
  Y. Chen, S. Kumar and B. Ratra,
  
  Astrophys. J.  835 (2017) 86
  
  [arXiv:1606.07316 [astro-ph.CO]].
  
  
  

h08
  T. Tram, R. Vallance and V. Vennin,
  
  JCAP 1701 (2017) no.01,  046
  
  [arXiv:1606.09199 [astro-ph.CO]].
  
  
  

h09
  J. L. Bernal, L. Verde and A. G. Riess,
  
  JCAP 1610 (2016) no.10,  019
  
  [arXiv:1607.05617 [astro-ph.CO]].
  
  
  

h010
  V. V. Luković, R. D'Agostino and N. Vittorio,
  
  Astron. Astrophys.  595 (2016) A109
  
  [arXiv:1607.05677 [astro-ph.CO]].
  
  

h011
  P. Ko and Y. Tang,
  
  Phys. Lett. B 762 (2016) 462
  
  [arXiv:1608.01083 [hep-ph]].
  
  
  

h012
  T. Karwal and M. Kamionkowski,
  
  Phys. Rev. D 94 (2016) no.10,  103523
  
  [arXiv:1608.01309 [astro-ph.CO]].
  
  
  

h013
  A. E. Romano,
  
  arXiv:1609.04081 [astro-ph.CO].
  
  

h014
  S. Joudaki et al.,
  
  arXiv:1610.04606 [astro-ph.CO].
  
  
  

h015
  A. Shafieloo and D. K. Hazra,
  
  JCAP 1704 (2017) no.04,  012
  
  [arXiv:1610.07402 [astro-ph.CO]].
  
  
  

h016
  W. Cardona, M. Kunz and V. Pettorino,
  
  JCAP 1703 (2017) no.03,  056
  
  [arXiv:1611.06088 [astro-ph.CO]].
  
  
  

h017
  S. Bethapudi and S. Desai,
  
  Eur. Phys. J. Plus 132 (2017) no.2,  78
  
  [arXiv:1701.01789 [astro-ph.CO]].
  
  
  

h018
  I. Odderskov, S. Hannestad and J. Brandbyge,
  
  JCAP 1703 (2017) no.03,  022
  
  [arXiv:1701.05391 [astro-ph.CO]].
  
  
  
hamannetal
  J. Hamann, S. Hannestad, J. Lesgourgues, C. Rampf and Y. Y. Y. Wong,
  
  JCAP 1007 (2010) 022
  
  [arXiv:1003.3999 [astro-ph.CO]].
  
  
  

Zhao:2016ecj
  M. M. Zhao, Y. H. Li, J. F. Zhang and X. Zhang,
  
  Mon. Not. Roy. Astron. Soc.  469 (2017) 1713
  
  [arXiv:1608.01219 [astro-ph.CO]].
  
  
  

Wang:2016tsz
  S. Wang, Y. F. Wang, D. M. Xia and X. Zhang,
  
  Phys. Rev. D 94 (2016) no.8,  083519
  
  [arXiv:1608.00672 [astro-ph.CO]].
  
  
  

brandbyge
  J. Brandbyge, S. Hannestad, T. Haugbølle and Y. Y. Y. Wong,
  
  JCAP 1009 (2010) 014
  
  [arXiv:1004.4105 [astro-ph.CO]].
  
  
  

ichiki
  K. Ichiki and M. Takada,
  
  Phys. Rev. D 85 (2012) 063521
  
  [arXiv:1108.4688 [astro-ph.CO]].
  
  
  

castorina2
  E. Castorina, E. Sefusatti, R. K. Sheth, F. Villaescusa-Navarro and M. Viel,
  
  JCAP 1402 (2014) 049
  
  [arXiv:1311.1212 [astro-ph.CO]].
  
  
  

costanzi
  M. Costanzi, F. Villaescusa-Navarro, M. Viel, J. Q. Xia, S. Borgani, E. Castorina and E. Sefusatti,
  
  JCAP 1312 (2013) 012
  
  [arXiv:1311.1514 [astro-ph.CO]].
  
  
  

castorina3
  E. Castorina, C. Carbone, J. Bel, E. Sefusatti and K. Dolag,
  
  JCAP 1507 (2015) no.07,  043
  
  [arXiv:1505.07148 [astro-ph.CO]].
  
  


Carbone:2016nzj
  C. Carbone, M. Petkova and K. Dolag,
  
  JCAP 1607 (2016) no.07,  034
  
  [arXiv:1605.02024 [astro-ph.CO]].
  
  
  

zennaro
  M. Zennaro, J. Bel, F. Villaescusa-Navarro, C. Carbone, E. Sefusatti and L. Guzzo,
  
  Mon. Not. Roy. Astron. Soc.  466 (2017) no.3,  3244
  
  [arXiv:1605.05283 [astro-ph.CO]].
  
  
  

rizzo
  L. A. Rizzo, F. Villaescusa-Navarro, P. Monaco, E. Munari, S. Borgani, E. Castorina and E. Sefusatti,
  
  JCAP 1701 (2017) no.01,  008
  
  [arXiv:1610.07624 [astro-ph.CO]].
  
  
  

Hand:2017ilm
  N. Hand, U. Seljak, F. Beutler and Z. Vlah,
  
  arXiv:1706.02362 [astro-ph.CO].
  
  
  

Modi:2017wds
  C. Modi, M. White and Z. Vlah,
  
  arXiv:1706.03173 [astro-ph.CO].
  
  

Seljak:2017rmr
  U. Seljak, G. Aslanyan, Y. Feng and C. Modi,
  
  arXiv:1706.06645 [astro-ph.CO].
  
  
jaffe
  J. Errard, S. M. Feeney, H. V. Peiris and A. H. Jaffe,
  
  JCAP 1603 (2016) no.03,  052
  
  [arXiv:1509.06770 [astro-ph.CO]].
  
  
  
inprep 
  S. Vagnozzi et al., in preparation.
  

Raveri:2015maa
  M. Raveri,
  
  Phys. Rev. D 93 (2016) no.4,  043522
  
  [arXiv:1510.00688 [astro-ph.CO]].
  
  
  

Heavens:2017hkr
  A. Heavens, Y. Fantaye, E. Sellentin, H. Eggers, Z. Hosenie, S. Kroon and A. Mootoovaloo,
  
  arXiv:1704.03467 [astro-ph.CO].
  
  
  

Hannestad:2005gj
  S. Hannestad,
  
  Phys. Rev. Lett.  95 (2005) 221301
  
  [astro-ph/0505551].
  
  
  
Archidiacono:2013fha 
  M. Archidiacono, E. Giusarma, S. Hannestad and O. Mena,
  
  Adv. High Energy Phys.  2013, 191047 (2013)
  
  [arXiv:1307.0637 [astro-ph.CO]].
  
  
  

Banerjee:2016suz
  A. Banerjee, B. Jain, N. Dalal and J. Shelton,
  
  arXiv:1612.07126 [astro-ph.CO].
  
  
  

Brust:2017nmv
  C. Brust, Y. Cui and K. Sigurdson,
  
  arXiv:1703.10732 [astro-ph.CO].
  
  

DiValentino:2013qma
  E. Di Valentino, A. Melchiorri and O. Mena,
  
  JCAP 1311 (2013) 018
  
  [arXiv:1304.5981 [astro-ph.CO]].
  
  
  

Roland:2016gli
  S. B. Roland and B. Shakya,
  
  JCAP 1705 (2017) no.05,  027
  
  [arXiv:1609.06739 [hep-ph]].
  
  
  
Melchiorri:2007cd
  A. Melchiorri, O. Mena and A. Slosar,
  
  Phys. Rev. D 76 (2007) 041303
  
  [arXiv:0705.2695 [astro-ph]].
  
  
  
Conlon:2013isa
  J. P. Conlon and M. C. D. Marsh,
  
  JHEP 1310 (2013) 214
  
  [arXiv:1304.1804 [hep-ph]].
  
  
  
Ackerman:mha
  L. Ackerman, M. R. Buckley, S. M. Carroll and M. Kamionkowski,
  
  Phys. Rev. D 79 (2009) 023519
  
  [arXiv:0810.5126 [hep-ph]].
  
  

Kaplan:2009de
  D. E. Kaplan, G. Z. Krnjaic, K. R. Rehermann and C. M. Wells,
  
  JCAP 1005 (2010) 021
  
  [arXiv:0909.0753 [hep-ph]].
  
  
  

Cline:2012is
  J. M. Cline, Z. Liu and W. Xue,
  
  Phys. Rev. D 85 (2012) 101302
  
  [arXiv:1201.4858 [hep-ph]].
  
  
  

CyrRacine:2012fz
  F. Y. Cyr-Racine and K. Sigurdson,
  
  Phys. Rev. D 87 (2013) no.10,  103515
  
  [arXiv:1209.5752 [astro-ph.CO]].
  
  
  

Fan:2013yva
  J. Fan, A. Katz, L. Randall and M. Reece,
  
  Phys. Dark Univ.  2 (2013) 139
  
  [arXiv:1303.1521 [astro-ph.CO]].
  
  
  

Vogel:2013raa
  H. Vogel and J. Redondo,
  
  JCAP 1402 (2014) 029
  
  [arXiv:1311.2600 [hep-ph]].
  
  
  

Petraki:2014uza
  K. Petraki, L. Pearce and A. Kusenko,
  
  JCAP 1407 (2014) 039
  
  [arXiv:1403.1077 [hep-ph]].
  
  
  

Foot:2014uba
  R. Foot and S. Vagnozzi,
  
  Phys. Rev. D 91 (2015) 023512
  
  [arXiv:1409.7174 [hep-ph]].
  
  
  

Foot:2014osa
  R. Foot and S. Vagnozzi,
  
  Phys. Lett. B 748 (2015) 61
  
  [arXiv:1412.0762 [hep-ph]].
  
  


Chacko:2015noa
  Z. Chacko, Y. Cui, S. Hong and T. Okui,
  
  Phys. Rev. D 92 (2015) 055033
  
  [arXiv:1505.04192 [hep-ph]].
  
  
  

Foot:2016wvj
  R. Foot and S. Vagnozzi,
  
  JCAP 1607 (2016) no.07,  013
  
  [arXiv:1602.02467 [astro-ph.CO]].
  
  
  

Boddy:2016bbu
  K. K. Boddy, M. Kaplinghat, A. Kwa and A. H. G. Peter,
  
  Phys. Rev. D 94 (2016) no.12,  123017
  
  [arXiv:1609.03592 [hep-ph]].
  
  



