Cosmological constraints on the neutrino mass including systematic uncertainties

TL;DR

This paper derives a robust upper limit on the sum of neutrino masses using cosmological data, carefully accounting for systematic uncertainties like foreground modeling and reionisation, within the ΛCDM framework.

Contribution

It provides the first comprehensive analysis of systematic uncertainties affecting neutrino mass constraints from cosmology, leading to a more reliable upper bound.

Findings

Upper limit on neutrino mass: < 0.17 eV at 95% CL.

Foreground systematics contribute about 0.01 eV to the uncertainty.

Current data cannot distinguish between normal and inverted neutrino mass hierarchies.

Abstract

When combining cosmological and oscillations results to constrain the neutrino sector, the question of the propagation of systematic uncertainties is often raised. We address this issue in the context of the derivation of an upper bound on the sum of the neutrino masses ($\Sigma m_\nu$) with recent cosmological data. This work is performed within the ${{\mathrm{\Lambda{CDM}}}}$ model extended to $\Sigma m_\nu$, for which we advocate the use of three mass-degenerate neutrinos. We focus on the study of systematic uncertainties linked to the foregrounds modelling in CMB data analysis, and on the impact of the present knowledge of the reionisation optical depth. This is done through the use of different likelihoods built from Planck data. Limits on $\Sigma m_\nu$ are derived with various combinations of data, including the latest Baryon Acoustic Oscillations (BAO) and Type Ia Supernovae (SN) results. We also discuss the impact of the preference for current CMB data for amplitudes of the gravitational lensing distortions higher than expected within the ${{\mathrm{\Lambda{CDM}}}}$ model, and add the Planck CMB lensing. We then derive a robust upper limit: $\Sigma m_\nu< 0.17\hbox{ eV at }95\% \hbox{CL}$, including 0.01 eV of foreground systematics. We also discuss the neutrino mass repartition and show that today's data do not allow one to disentangle normal from inverted hierarchy. The impact on the other cosmological parameters is also reported, for different assumptions on the neutrino mass repartition, and different high and low multipole CMB likelihoods.

Figures

Tables

Table 1. Table 1: Summary of data and likelihoods with their corresponding acronyms. All are ready to use in the CAMEL  software. Plik , Commander , lowTEB are available through the Planck  PLA.

Acronym	Description
hlpTT	high-$\ell$ HiLLiPOP temperature Planck  likelihood (cf. Sect 2.1)
hlpTTps	high-$\ell$ HiLLiPOP temperature Planck likelihood with an astrophysical model of point sources
PlikTT	public high-$\ell$ temperature Planck likelihood
TT	refers to the temperature CMB data
TE	refers to the TE CMB correlations
ALL	refers to the combination of temperature and polarisation CMB data (incl. TT and TE)
Comm	Commander  low-$\ell$ temperature Planck public likelihood (cf. Sect 2.2)
lowTEB	pixel-based temperature and polarisation low-$\ell$ Planck public likelihood (cf. Sect 2.2)
${\tau }_{\mathrm{reio}}$	auxiliary constraint on ${\tau }_{\mathrm{reio}}$ from Planck reionisation measurement with Lollipop (cf. Sect 2.2)
VHL	very high-$\ell$ data (cf. Sect 3.4)
BAO	latest DR12 BAO data (cf. Sect 2.5)
SNIa	JLA supernovae compilation (cf. Sect 2.6)

Table 2. Table 2: 95 % percent \% CL upper limits on Σ ​ m ν Σ subscript m 𝜈 \mathrm{\Sigma m_{\nu}}  in ν ​ Λ ​ CDM ​ ( 3 ​ ν ) 𝜈 Λ CDM 3 𝜈 \mathrm{\nu\Lambda{CDM}(3\nu)}  (i.e. with A L = 1 subscript A L 1 \mathrm{A_{L}}=1 ) and results on A L subscript A L \mathrm{A_{L}}   (68 % percent \% CL) in the Λ ​ CDM ​ ( 3 ​ ν ) Λ CDM 3 𝜈 \mathrm{\Lambda{CDM}(3\nu)} + A L subscript A L \mathrm{A_{L}}  model (i.e. with Σ ​ m ν = 0.06 ​ eV Σ subscript m 𝜈 0.06 eV \mathrm{\Sigma m_{\nu}}=0.06\thinspace\text{eV} ) obtained when combining the Planck  TT+lowTEB+BAO+ SNI a.

PlanckTT+lowTEB	$\Sigma m_{\nu }$	$A_L$
BAO+SNIa	limit (eV)
hlpTT	0.18	1.16$\pm$0.09
hlpTTps	0.20	1.14$\pm$0.08
PlikTT	0.17	1.19$\pm$0.09

Table 3. Table 3: 95 % percent \% CL upper limits on Σ ​ m ν Σ subscript m 𝜈 \mathrm{\Sigma m_{\nu}}  in ν ​ Λ ​ CDM ​ ( 3 ​ ν ) 𝜈 Λ CDM 3 𝜈 \mathrm{\nu\Lambda{CDM}(3\nu)}  obtained with hlpTE +BAO+ SNI a in combination with lowTEB, or an auxiliary constraint on τ reio subscript 𝜏 reio \tau_{\rm reio} and Commander .

PlanckTE+low-$\ell$	$\Sigma m_{\nu }$
+BAO+SNIa	limit (eV)
hlpTE+lowTEB	0.20
$\text{hlpTE}\text{+}{\tau }_{\mathrm{reio}}+$Commander	0.19

Table 4. Table 4: Results on Σ ​ m ν Σ subscript m 𝜈 \mathrm{\Sigma m_{\nu}}  (95 % percent \% CL upper limits) and A L subscript A L \mathrm{A_{L}}  (68% CL) obtained from a combined fit in the ν ​ Λ ​ CDM ​ ( 3 ​ ν ) 𝜈 Λ CDM 3 𝜈 \mathrm{\nu\Lambda{CDM}(3\nu)} + A L subscript A L \mathrm{A_{L}}  model with Planck  TT+lowTEB+BAO+ SNI a.

PlanckTT+lowTEB	($\Sigma m_{\nu }$[eV],$A_L$)
BAO+SNIa
hlpTT	(0.39, $1.22\pm 0.12$)
hlpTTps	(0.34, $1.18\pm 0.10$)
PlikTT	(0.40, $1.28\pm 0.12$)

Table 5. Table 5: 95 % percent \% CL upper limits on Σ ​ m ν Σ subscript m 𝜈 \mathrm{\Sigma m_{\nu}}  in ν ​ Λ ​ CDM ​ ( 3 ​ ν ) 𝜈 Λ CDM 3 𝜈 \mathrm{\nu\Lambda{CDM}(3\nu)}  (i.e. with A L = 1 subscript A L 1 \mathrm{A_{L}}=1 ) and results on A L subscript A L \mathrm{A_{L}}   (68 % percent \% CL) in the Λ ​ CDM ​ ( 3 ​ ν ) Λ CDM 3 𝜈 \mathrm{\Lambda{CDM}(3\nu)} + A L subscript A L \mathrm{A_{L}}  model (i.e. with Σ ​ m ν = 0.06 ​ eV Σ subscript m 𝜈 0.06 eV \mathrm{\Sigma m_{\nu}}=0.06\thinspace\text{eV} ) obtained when combining Planck  TT+lowTEB+BAO+ SNI a+lensing.

PlanckTT+lowTEB	$\Sigma m_{\nu }$	$A_L$
BAO+SNIa+lensing	limits (eV)
hlpTT	0.21	1.06 $\pm$ 0.05
hlpTTps	0.21	1.06 $\pm$ 0.06
PlikTT	0.23	1.05 $\pm$ 0.06

Table 6. Table 6: 95 % percent \% CL upper limits on Σ ​ m ν Σ subscript m 𝜈 \mathrm{\Sigma m_{\nu}}  obtained with hlpTT +lowTEB+BAO+ SNI a+lensing for different neutrino mass repartition: three degenerate masses, normal hierarchy (NH), inverse hierarchy (IH) and one massive plus two massless neutrinos.

$\nu$ mass setting	$\Sigma m_{\nu }$ limits (eV)
$\left[3\nu \right]\Lambda \mathrm{CDM}\left(3\nu \right)$	0.21
$\left[3\nu \text{NH}\right]$	0.21
$\left[3\nu \text{IH}\right]$	0.21
$\left[1\nu \right]\Lambda \mathrm{CDM}\left(1\nu \right)$	0.23

Table 7. Table 7: 95 % percent \% CL upper limits on Σ ​ m ν Σ subscript m 𝜈 \mathrm{\Sigma m_{\nu}}  in ν ​ Λ ​ CDM ​ ( 3 ​ ν ) 𝜈 Λ CDM 3 𝜈 \mathrm{\nu\Lambda{CDM}(3\nu)}  obtained when combining PlikALL , hlpALL  or hlpALLps   with lowTEB+ SNI a+BAO. ALL refers to the combination TT+TE+EE.

PlanckALL+lowTEB	$\Sigma m_{\nu }$
+SNIa+BAO	limits (eV)
hlpALL	0.18
hlpALLps	0.18
PlikALL	0.15

Table 8. Table 8: 95 % percent \% CL upper limits on Σ ​ m ν Σ subscript m 𝜈 \mathrm{\Sigma m_{\nu}}  in ν ​ Λ ​ CDM ​ ( 3 ​ ν ) 𝜈 Λ CDM 3 𝜈 \mathrm{\nu\Lambda{CDM}(3\nu)}  obtained when combining PlikALL , hlpALL  or hlpALLps   with SNI a+BAO+lensing, using lowTEB for the low- ℓ ℓ \ell (second column) and for the combination of an auxiliary constraint on τ reio subscript 𝜏 reio \tau_{\rm reio} + Commander (third column) . ALL refers to the combination TT+TE+EE.

PlanckALL	+lowTEB	+${\tau }_{\mathrm{reio}}$
+SNIa+BAO+lensing		+Commander
hlpALL	0.20	0.16
hlpALLps	0.21	0.17
PlikALL	0.19	0.17

Table 9. Table 9: Values of the χ 2 = − 2 ​ log ⁡ ℒ superscript 𝜒 2 2 ℒ \chi^{2}=-2\log\cal{L} profiled out over all the other (cosmological and nuisance) parameters as a function of Σ ​ m ν Σ subscript m 𝜈 \mathrm{\Sigma m_{\nu}}  for the hlpALLps +BAO+ SNI a+lensing combination, using the auxiliary constraint on τ reio subscript 𝜏 reio \tau_{\rm reio}  combined with Commander  at low ℓ ℓ \ell . They correspond to the red dots/plain line of Fig.  11 .

$\Sigma m_{\nu }$( eV)	${\chi }^2=-2\log ℒ$
0.01	28613.76
0.02	28613.75
0.04	28613.86
0.06	28614.08
0.08	28614.44
0.10	28614.93
0.12	28615.57
0.14	28616.35
0.16	28617.29
0.18	28618.38
0.20	28619.58
0.25	28623.26
0.30	28627.86
0.35	28633.30
0.40	28639.60

Equations

$$\Delta m_{21}^2=m_2^2-m_1^2$$

$$\Delta m^2=m_3^2-(m_1^2+m_2^2)/2.$$

$$=$$

$$A_{\mathrm{S}Z}=A_{\mathrm{k}SZ}+1.6A_{\mathrm{t}SZ}=9.5\pm 3\text{leavevmode}\mu K^2,$$

$${\tau }_{\mathrm{reio}}=0.058\pm 0.012,$$

$$\mathrm{\Sigma }{\mathrm{m}}_{\nu }<0.17[incl.0.01(foreground\; syst.)]eV\; at\; 95\%\; CL.$$

Full text

Introduction

In the last decade, cosmology has entered a precision era, confirming the six parameters $\mathrm{\Lambda{CDM}}$ concordance model with unprecedented accuracy. This allows us to open the parameters’ space, and to confront the corresponding extensions with data. In the following, we explore the neutrino sector. We only deal with three standard neutrinos species (Schael et al. 2006), and focus on the extension to the sum of the neutrino masses ($\mathrm{\Sigma m_{\nu}}$). Moreover, the neutrino mass splitting scenario has been set up to match the neutrino oscillation results. A three mass-degenerate neutrinos model is advocated for and used throughout this study. It must be noted that the assumptions on the neutrino mass scenario have already been shown to be of particular importance for the derivations of cosmological results (for example in Marulli et al. 2011).

Recent works (for instance Alam et al. 2016; Sherwin et al. 2016; Giusarma et al. 2016; Y che et al. 2017; Vagnozzi et al. 2017) on the derivation of upper bounds on $\mathrm{\Sigma m_{\nu}}$  usually take the Cosmological Microwave Background (CMB) likelihoods as granted. Furthermore, no uncertainty from the analysis of this cosmological probe is propagated until the final results. In this paper, we investigate the systematic uncertainties linked to the modelling of foreground residuals in the Planck CMB likelihood implementations.

To address this issue, the most accurate method would have been to make use of full end to end simulations, including an exhaustive description of the foregrounds. This is not possible given the actual knowledge of the foreground’s physical properties. Instead, we propose a comparison of the results derived from different likelihoods built from the Planck 2015 data release, and based on different foreground assumptions. Namely the public Plik and the HiLLiPOP likelihoods are examined for the high-$\ell$ part. We also investigate the impact of our current knowledge on the reionisation optical depth ($\tau_{\rm reio}$). For the low-$\ell$ part, the lowTEB likelihood is compared to the combination of the Commander likelihood with an auxiliary constraint on the $\tau_{\rm reio}$ parameter, derived from the last Planck 2016 measurements (Planck Collaboration Int. XLVII 2016).

The differences of the impact of the foreground modellings are twofold: On one hand they show up as slight deviations on the $\mathrm{\Sigma m_{\nu}}$ bounds inferred from the different likelihoods, and, on the other hand, they manifest themselves in the form of different values of the amplitude of the gravitational lensing distortions ($\mathrm{A_{L}}$). Indeed, fitting for $\mathrm{A_{L}}$ represents a direct test of the accuracy and robustness of the likelihood with respect to the $\mathrm{\Lambda{CDM}}$ model (Couchot et al. 2017). We also address this point, and discuss how it is linked to $\mathrm{\Sigma m_{\nu}}$.

Derivations of systematic uncertainties on $\mathrm{\Sigma m_{\nu}}$ are performed for different combinations of cosmological data: The Planck temperature and polarisation likelihoods, the latest BAO data from Boss DR12, and SNIa, as well as the direct measurement of the lensing distortion field power spectrum from Planck.

We also address the question of the sensitivity of the combination of those datasets to the neutrinos mass hierarchy.

We start with a description of the standard cosmology, the impact of massive neutrinos, and their mass repartition, as well as the profile likelihood method (Sect. 1). In Sect. 2, we describe the likelihoods and datasets. Turning to the $\mathrm{\Sigma m_{\nu}}$ constraints, we first focus on the results obtained with CMB temperature data for different likelihoods at intermediate multipoles. We investigate different choices for the low-$\ell$ likelihoods, and examine the pros and cons of the use of high-angular-resolution datasets. In Sect. 4, we derive the $\mathrm{\Sigma m_{\nu}}$ constraints obtained when combining CMB temperature, BAO and SNIa data, and check the robustness of the results with respect to the high-$\ell$ likelihoods. The choices for the low-$\ell$ parts are compared. A cross-check of the results is performed using the temperature-polarisation TE correlations. Then, the impact of the observed tension on $\mathrm{A_{L}}$ is further discussed, followed by the combination of the data with the CMB lensing. The neutrino mass hierarchy question is addressed in this context. In Sect. 5.1, we discuss the (TT+TE+EE) combination with BAO and SN data, with and without CMB lensing. Finally, we derive the cosmological parameters and illustrate their variations depending on the assumptions on the neutrino mass repartition, the low-$\ell$ likelihoods, and the fact that we release or do not release $\mathrm{\Sigma m_{\nu}}$ in the fits.

1 Phenomenology and methodology

This section discusses the standard cosmology and the role of neutrinos in the Universe’s thermal history. We then briefly review the current constraints coming from the observation of the neutrino oscillations phenomenon, and discuss the mass hierarchy. A definition of the $\mathrm{\Lambda{CDM}}$ models considered for this paper is given. The statistical methodology based on profile likelihoods is also presented.

1.1 Standard Cosmology

The “standard” cosmological model describes the evolution of a homogeneous and isotropic Universe, the geometry of which is given by the Friedman-Robertson-Walker metric, following General Relativity. In this framework, the theory reduces to the well-known Friedman equations. The Universe is assumed to be filled with several components, of different nature and evolution (matter, radiation,…). Their inhomogeneities are accounted for as small perturbations of the metric. In the $\mathrm{\Lambda{CDM}}$ model, the Universe’s geometry is assumed to be euclidean (no curvature) and its constituents are dominated today by a cosmological constant ($\Lambda$), associated with dark energy, and cold dark matter; it also includes radiation, baryonic matter and three neutrinos. Density anisotropies are assumed to result from the evolution of primordial power spectra, and only purely adiabatic scalar modes are assumed.

The minimal $\mathrm{\Lambda{CDM}}$ model is described with only six parameters. Two of them describe the primordial scalar mode power spectrum: the amplitude ($\mathrm{A_{\rm s}}$), and the spectral index ($n_{\rm s}$). Two other parameters represent the reduced energy densities today: $\omega_{\rm b}=\Omega_{\rm b}h^{2}$, for the baryon, and $\omega_{\rm c}=\Omega_{\rm c}h^{2}$ for the cold dark matter. The last two parameters are the angular size of the sound horizon at decoupling, $\theta_{S}$, and the reionisation optical depth ($\tau_{\rm reio}$). In this chosen parameterisation, $H_{0}$ is derived in a non-trivial way from the above parameters. In addition, the sum of the neutrino masses is usually fixed to $\mathrm{\Sigma m_{\nu}}=0.06$  eV based on oscillation constraints (Forero et al. 2012, 2014; Capozzi et al. 2016): This is discussed in Sect. 1.3.

Departures from the $\mathrm{\Lambda{CDM}}$ model assumptions are often studied by extending its parameter space and testing it against the data, for instance, through the inclusion of $\Omega_{k}$ for non-euclidean geometry, $N_{\hbox{\tiny{eff}}}$ for the number of effective relativistic species, or $Y_{p}$ for the primordial mass fraction of ${}^{4}\hbox{He}$ during BBN. In addition to those physics-related parameters, a phenomenological parameter, $\mathrm{A_{L}}$, has been introduced (Calabrese et al. 2008) to scale the deflection power spectrum which is used to lens the primordial CMB power spectra. This parameter permits to size the (dis-)agreement of the data with the $\mathrm{\Lambda{CDM}}$ lensing distortion predictions. Testing that its value, inferred from data, is compatible with one is a thorough consistency check (we refer to e.g. Calabrese et al. 2008; Planck Collaboration XIII 2016; Couchot et al. 2017). In this work, we use the $\mathrm{A_{L}}$ consistency check in the context of the constraints on $\mathrm{\Sigma m_{\nu}}$. In practice, it means that we check the value of $\mathrm{A_{L}}$ (using $\mathrm{\Lambda{CDM}}$+$\mathrm{A_{L}}$ model) for each dataset on which we then report a $\mathrm{\Sigma m_{\nu}}$ limit (using $\nu\mathrm{\Lambda{CDM}}$ model, i.e. with $\mathrm{A_{L}}=1$).

1.2 Neutrinos in cosmology

One of the generic features of the standard hot big bang model is the existence of a relic cosmic neutrino background. In parallel, the observation of the neutrino oscillation phenomena requires that those particles are massive, and establishes the existence of flavour mixed-mass eigenstates (cf. Sect. 1.3) (Pontecorvo 1957; Maki et al. 1962). As far as cosmology is concerned, depending on the mass of the lightest neutrino (Bilenky et al. 2001), this implies that there are at least two non-relativistic species today. Massive neutrinos therefore impact the energy densities of the Universe and its evolution.

Initially neutrinos are coupled to the primeval plasma. As the Universe cools down, they decouple from the rest of the plasma at a temperature up to a few MeV depending on their flavour (Dolgov 2002). This decoupling is fairly well approximated as an instantaneous process (Kolb & Turner 1994; Dodelson 2003). Given the fact that, with today’s observational constraints, neutrinos can be considered as relativistic at recombination (Lesgourgues & Pastor 2006). In addition, for $\mathrm{m_{\nu}}$ in the range from $10^{-3}$ to 1  eV, they should be counted as radiation at the matter-radiation equality redshift, $z_{\hbox{\tiny{eq}}}$, and as non-relativistic matter today (Lesgourgues & Pastor 2014; Lesgourgues et al. 2013), which is measured through $\Omega_{m}$. $\mathrm{\Sigma m_{\nu}}$ is therefore correlated to both $z_{\hbox{\tiny{eq}}}$ and $\Omega_{m}$.

The induced modified background evolution is reflected in the relative position and amplitude of the peaks of the CMB power spectra (through $z_{\hbox{\tiny{eq}}}$). It also affects the CMB anisotropies power spectra at intermediate or high multipole ($\ell\gtrsim 200$) as potential shifts of the power spectrum due to a change in the angular distance of the sound horizon at decoupling. Finally it also leaves an imprint on the slope of the low-$\ell$ tail due to the late Integrated Sachs Wolfe (ISW) effect. An additional effect of massive neutrinos comes from the fact that they affect the photon temperature through the early ISW effect. As a result a reduction of the CMB temperature power spectrum below $\lesssim 500$ is observed.

On the matter power spectrum side, two effects are induced by the massive neutrinos. In the early Universe, they free-stream out of potential wells, damping matter perturbations on scales smaller than the horizon at the non-relativistic transition. This results in a suppression of the P(k) at large k which also depends on the individual masses repartition (Hu et al. 1998; Lesgourgues & Pastor 2006). At late time, the non-relativistic neutrino masses modify the matter density, which tends to slow down the clustering.

CMB anisotropies are lensed by large-scale structures (LSS). Measuring CMB gravitational lensing therefore provides a constraint on the matter power spectrum on scales where the effects of massive neutrinos are small but still sizeable (Kaplinghat et al. 2003; Lesgourgues et al. 2006).

1.3 Neutrino mass hierarchy

As stated above, we have to choose a neutrino mass splitting scenario to define the $\mathrm{\Lambda{CDM}}$ model. In general, CMB data analyseis that aim at measuring cosmological parameters not related to the neutrino sector (including Planck papers, e.g. Planck Collaboration XIII 2016) are done assuming two massless neutrinos and one massive neutrino, while fixing $\mathrm{\Sigma m_{\nu}}=0.06$  eV.

For the work of this paper, our choice is motivated considering neutrino oscillation data. More precisely, we use the differences of squared neutrino masses deduced from the best fit values of the global 3$\nu$ oscillation analysis based on the work of Capozzi et al. (2016):

[TABLE]

where the two usual scenarios are considered: The normal (NH) and the inverted hierarchy (IH), for which the lightest neutrino is the one of the first and third generation respectively.

Individual masses can be computed numerically under the above assumptions, for each mass hierarchy, as a function of $\mathrm{\Sigma m_{\nu}}$ , as highlighted in Figure 1 (see also Lesgourgues & Pastor (2014)). In each hierarchy, Equations 1-3 impose a lower bound on $\mathrm{\Sigma m_{\nu}}$  , corresponding to the case where the lightest mass is strictly null (numerically, $\sim$0.059 and $\sim$0.099 eV for NH and IH, respectively); also shown in Figure 1 as vertical dashed lines.

Those results show that, given the oscillation constraints, neutrino masses are nearly degenerate for $\mathrm{\Sigma m_{\nu}}\gtrsim 0.25$  eV. Moreover, given the current cosmological probes (essentially CMB and BAO data), we observe almost no difference in $\mathrm{\Sigma m_{\nu}}$ constraints when comparing results obtained with one of the two hierarchies with the case with three mass-degenerate neutrinos for which the mass repartition is such that each neutrino carries $\mathrm{\Sigma m_{\nu}}/3.$ (we refer to Sect. 4.5 and Giusarma et al. 2016; Vagnozzi et al. 2017; Schwetz et al. 2017). Indeed, as shown in Palanque-Delabrouille et al. (2015), the difference is less than 0.3% in the 3D linear matter power spectrum and is reduced even to less than 0.05% when considering the 1D flux power spectrum  (see also Agarwal & Feldman 2011). This justifies the simplifying choice of the three mass-degenerate neutrinos scenario, which is used in this paper.

In Sect. 4.5, we show that this is not equivalent to the configuration where the total mass is entirely given to one massive neutrino with the two other neutrinos being massless.

1.4 Constraints on $\mathrm{\Sigma m_{\nu}}$ and degeneracies

The inference from CMB data of a limit on $\mathrm{\Sigma m_{\nu}}$in the $\mathrm{\Lambda{CDM}}$ framework is not trivial because of degeneracies between parameters. Indeed, the impact of $\mathrm{\Sigma m_{\nu}}$ on the CMB temperature power spectrum is partly degenerated with that of some of the six other parameters.

In particular, the impact of neutrino masses on the angular-diameter distance to last scattering surface is degenerated with $\Omega_{\Lambda}$ (and consequently with the derived parameters $H_{0}$ and $\sigma_{8}$) in flat models and with $\Omega_{k}$ otherwise (Hou et al. 2014). Late-time geometric measurements help in reducing this geometric degeneracy. Indeed, at fixed $\theta_{S}$, the BAO distance parameter $D_{V}(z)$ increases with increasing neutrino mass while the Hubble parameter decreases.

Another example is the correlation of $\mathrm{\Sigma m_{\nu}}$ with $\mathrm{A_{\rm s}}$ (Allison et al. 2015). As explained in Sect 1.2, $\mathrm{\Sigma m_{\nu}}$ can impact the amplitude of the matter power spectrum and thus is directly correlated to $\mathrm{A_{\rm s}}$ and consequently with $\tau_{\rm reio}$ through the amplitude of the first acoustic peak (which scales like $\mathrm{A_{\rm s}}e^{-2\tau_{\rm reio}}$). The constraint on $\mathrm{\Sigma m_{\nu}}$ therefore depends on the low-$\ell$ polarisation likelihood, which drives the constraints on $\tau_{\rm reio}$. The addition of lensing distortions, the amplitude of which is proportional to $\mathrm{A_{\rm s}}$, helps to break this degeneracy.

Moreover, the suppression of the small-scale power in LSS due to massive neutrinos, which imprints on the CMB lensing spectra, can be compensated for by an increase of the cold-dark-matter density, shifting the matter-radiation equality to early times (Hall & Challinor 2012; Pan et al. 2014). This induces an anti-correlation between $\mathrm{\Sigma m_{\nu}}$ and $\Omega_{\rm cdm}$ when using CMB observable. On the contrary, both parameters similarly affect the angular diameter distance so that BAO can help to break this degeneracy.

1.5 Cosmological model

As discussed in the previous sections, the neutrino mass repartition can have significant impact on the constraints for $\mathrm{\Sigma m_{\nu}}$. By $\mathrm{\Lambda{CDM}(1\nu)}$, we refer to the definition used in Planck Collaboration XIII (2016); it assumes two massless and one massive neutrinos.

However, in the following, we adopt a scenario with three mass-degenerate neutrinos, that is, where the neutrino generations equally share the mass ($\mathrm{\Sigma m_{\nu}}$/3). We note that this is also the model adopted in Planck Collaboration XIII (2016) when $\mathrm{\Sigma m_{\nu}}$ constraints have been extracted. We also stick to this scenario when fixing $\mathrm{\Sigma m_{\nu}}$ to 0.06  eV and we refer to it as $\mathrm{\Lambda{CDM}(3\nu)}$.

The notations $\mathrm{\nu\Lambda{CDM}(1\nu)}$ (resp. $\mathrm{\nu\Lambda{CDM}(3\nu)}$) will be used to differentiate the case where we open the parameters’ space to $\mathrm{\Sigma m_{\nu}}$ from the $\mathrm{\Lambda{CDM}(1\nu)}$ (resp. $\mathrm{\Lambda{CDM}(3\nu)}$) case.

To derive the values for the observables from the cosmological model, we make use of the CLASS Boltzmann solver (Blas et al. 2011). Within this software, the non-linear effects on the matter power spectrum evolution can be included using the halofit model recalibrated as proposed in Takahashi et al. (2012) and extended to massive neutrinos as described in Bird et al. (2012). Our baseline setup for the $\mathrm{\Sigma m_{\nu}}$ studies is to use CLASS, including non-linear effects, tuned to a high-precision setting.

In order to compare order of magnitudes in the non-linear effects propagation, we have also used CAMB (Lewis et al. 2000), in which both the Takahashi and the Mead (Mead et al. 2016) models are made available.

1.6 Profile likelihoods

The results described below were obtained from profile likelihood analyses performed with the CAMEL software111camel.in2p3.fr(Henrot-Versill et al. 2016). As described in Planck Collaboration Int. XVI (2014), this method aims at measuring a parameter $\theta$ through the maximisation of the likelihood function $\cal{L}(\theta,\mu)$, where $\mu$ is the full set of cosmological and nuisance parameters excluding $\theta$. For different, fixed $\theta_{i}$ values, a multidimensional minimisation of the $\chi^{2}(\theta_{i},\mu)=-2\ln{\cal{L}}(\theta_{i},\mu)$ function is performed. The absolute minimum, $\chi^{2}_{min}$, of the resulting $\chi^{2}_{min}(\theta_{i})$ curve is by construction the (invariant) global minimum of the problem, that is, the ‘best fit’. From the $\chi^{2}_{min}(\theta_{i})-\chi^{2}_{min}$ curve, the so-called profile likelihood, one can derive an estimate of $\theta$ and its associated uncertainty (James 2007). All minimisations have been performed using the MINUIT software (James 1994). In the $\mathrm{\Sigma m_{\nu}}$ studies presented below, $95\%$ CL upper limits are derived following the Gaussian prescription proposed by Feldman & Cousins (1998) (hereafter denoted F.C.), as described in  Planck Collaboration Int. XVI (2014).

Unless otherwise explicitly stated, we use the frequentist methodology throughout this paper. A comparison with the Bayesian approach has already been presented in Planck Collaboration Int. XVI (2014) and Planck Collaboration XI (2016), showing that results do not depend on the chosen statistical method for the $\mathrm{\Lambda{CDM}}$ model, as well as for $\nu\mathrm{\Lambda{CDM}}$.

2 Likelihoods and datasets

In this Section, we detail the likelihoods that are used hereafter for the derivation of the results on $\mathrm{\Sigma m_{\nu}}$. They are summarised in Table 1 together with their related acronyms.

2.1 Planck high-$\ell$ likelihoods

In order to assess the impact of foreground residuals modelling on the $\mathrm{\Sigma m_{\nu}}$ constraints, we make use of different Planck high-$\ell$ likelihoods (HiLLiPOP and Plik). They both use a Gaussian approximation of the likelihood based on cross-spectra between half-mission maps at the three lowest frequencies (100, 143 and 217 GHz) of Planck-HFI, but rely on different assumptions for modelling foreground residuals. Comparing the results on $\mathrm{\Sigma m_{\nu}}$ obtained with both of these likelihoods is a way to assess a systematic uncertainty on the foreground residuals modelling.

Plik is the public Planck likelihood. It is described in detail in Planck Collaboration XI (2016). It uses empirically motivated power spectrum templates to model residual contamination of foregrounds (including dust, CIB, tSZ, kSZ, SZxCIB and point sources) in the cross-spectra. The foreground residuals in HiLLiPOP are directly derived from Planck measurements (Couchot et al. 2017): This is the main difference between HiLLiPOP and Plik. For $\rm{\Lambda CDM}$ cosmology, both likelihoods have been compared in Planck Collaboration XI (2016).

In any of the Planck  high-$\ell$ likelihoods, the residual amplitudes of the foregrounds are compatible with expectations, with only a mild tension on unresolved point-source amplitudes coming essentially from the 100  GHz frequency (we refer to Sect. 4.3 in Planck Collaboration XI 2016). In order to assess the impact of the point-source modelling on the parameter reconstructions (and in particular $\mathrm{\Sigma m_{\nu}}$), we use two variants of the HiLLiPOP likelihood. The first one, labelled hlpTTps, makes use of a physical model with two unresolved point-source components, corresponding to the radio and IR frequency domains, with fixed frequency scaling factors and number counts tuned on data (Couchot et al. 2017). The second one, labelled hlpTT, uses one free amplitude for unresolved point-sources per cross-frequency leading to six free parameters (as used in Couchot et al. 2017), in a similar way as what is done in Plik. This allows one to alleviate the tension on the point-source amplitudes. Both hlpTTps and hlpTT lead to very similar results in the $\mathrm{\Lambda{CDM}(1\nu)}$ model, with a lower level of correlation between parameters for the former. Comparing results obtained with hlpTTps and hlpTT is therefore useful for assessing their robustness with respect to the unresolved point-source tension.

Both HiLLiPOP and alsoPlik include polarisation information using the EE and TE angular cross-power spectra. Unless otherwise explicitly stated, only the temperature (TT) part is considered in the following.

Together with auxiliary constraints on nuisance parameters (such as the relative and absolute calibration) associated to each likelihood, we can also add a Gaussian constraint to the SZ template amplitudes as suggested in Planck Collaboration XI (2016). This constraint is based on a joint analysis of the Planck-2013 data with those from ACT and SPT (see Sect. 2.3) and reads:

[TABLE]

when normalized at $\ell=3000$. The role of this additional constraint is also discussed in the following.

2.2 Low-$\ell$

At low-$\ell$, two options are investigated to study the impact of one choice or another on the $\mathrm{\Sigma m_{\nu}}$ limit determination:

• •

lowTEBA pixel-based likelihood that relies on the Planck low-frequency instrument 70 GHz maps for polarisation and on a component-separated map using all Planck frequencies for temperature (Commander).

• •

A combination of a temperature-only likelihood, Commander (Planck Collaboration XI 2016), based on a component-separated map using all Planck frequencies, and a Gaussian auxiliary constraint on the reionisation optical depth,

[TABLE]

derived from the last Planck results of the reionisation optical depth (Planck Collaboration Int. XLVII 2016) Lollipop likelihood (Mangilli et al. 2015).

2.3 High-resolution CMB data

High resolution CMB data, namely the ACT, SPT_high, and SPT_low datasets are also used in this work. They are later quoted “VHL” (Very high-$\ell$) when combined altogether. The ACT data are those presented in Das et al. (2014). They correspond to cross power spectra between the 148 and 220 GHz channels built from observations performed on two different sky areas (an equatorial strip of about $300\ \mathrm{deg}^{2}$ and a southern strip of $292\ \mathrm{deg}^{2}$ for the 2008 season, and about $100\ \mathrm{deg}^{2}$ otherwise) and during several seasons (between 2007 and 2010), for multipoles between 1000 and 10000 (for 148$\times$148) and 1500 to 10000 otherwise. For SPT, two distinct datasets are examined. The higher $\ell$ part, dubbed SPT_high, implements the results, described in Reichardt et al. (2012), from the observations of $800\ \mathrm{deg}^{2}$ at 95, 150, and 220 GHz of the SPT-SZ survey. The cross-spectra cover the $\ell$ range between 2000 and 10000. As in Couchot et al. (2017), we prefer not to consider the more recent data from George et al. (2015) because the calibration, based on the Planck 2013 release, leads to a 1% offset with respect to the last Planck data. We also add the Story et al. (2012) dataset, dubbed SPT_low, consisting of a 150 GHz power spectrum, which ranges from $\ell=650$ to 3000, resulting from the analysis of observations of a field of $2540\ \mathrm{deg}^{2}$. Both SPT datasets have an overlap in terms of sky coverage and frequency. We have however checked that this did not bias the results by, for example, removing the 150x150 GHz part from the SPT_high likelihood, as was done in Couchot et al. (2017).

2.4 Planck CMB Lensing

The full sky CMB temperature and polarisation distributions are inhomogeneously affected by gravitational lensing due to large-scale structures. This is reflected in additional correlations between large and small scales, and, in particular, in a smoothing of the power spectra in TT, TE, and EE. From the reconstruction of the four-point correlation functions (Hu & Okamoto 2002), one can reconstruct the power spectrum of the lensing potential $C_{\ell}^{\phi\phi}$ of the lensing potential $\phi$. In the following we make use of the corresponding 2015 temperature lensing likelihood estimated by Planck (Planck Collaboration XV 2016).

2.5 Baryon acoustic oscillations

In Sect. 4, information from the late-time evolution of the Universe geometry are also included. The more accurate and robust constraints on this epoch come from the BAO scale evolution. They bring cosmological parameter constraints that are highly complementary with those extracted from CMB, as their degeneracy directions are different.

BAO generated by acoustic waves in the primordial fluid can be accurately estimated from the two-point correlation function of galaxy surveys. In this work, we use the acoustic-scale distance ratio $D_{\hbox{\tiny{V}}}(z)/r_{\mathrm{drag}}$ measurements from the 6dF Galaxy Survey at $z=0.1$ (Beutler et al. 2014). At higher redshift, we included the BOSS DR12 BAO measurements that recently have been made available (Alam et al. 2016). They consist in constraints on $(D_{M}(z),H(z),f\sigma_{8}(z))$ in three redshift bins, which encompass both BOSS-LowZ and BOSS-CMASS DR11 results. Thanks to the addition of the results on $f\sigma_{8}(z)$ the constraints on $\mathrm{\Sigma m_{\nu}}$ are significantly reduced with respect to previous BAO measurements (Alam et al. 2016). The combination of those measurements is labelled “BAO” in the following. We note that this is an update of the BAO data with respect those used in Planck Collaboration XIII (2016).

2.6 Type Ia supernovae

SNIa also constitute a powerful cosmological probe. The study of the evolution of their apparent magnitude with redshift played a major role in the discovery of late-time acceleration of the Universe. We include the JLA compilation (Betoule et al. 2014), which spans a wide redshift range (from 0.01 to 1.2), while compiling up-to-date photometric data. This is further referenced to as “SNIa” in the following.

3 CMB temperature results

3.1 Orders of magnitude

The differences between the expected $C_{\ell}$ spectra for $\mathrm{\Sigma m_{\nu}}$=0.3 eV and $\mathrm{\Sigma m_{\nu}}$=0.06 eV in the $\mathrm{\Lambda{CDM}(3\nu)}$ model are shown on Fig. 2 in black on the upper panel without considering any non-linearities. The shaded area indicates the CMB spectrum divided by a factor $10^{3}$. The size of the effect of increasing $\mathrm{\Sigma m_{\nu}}$ up to 0.3 eV, except at the first peak, is of the order of $\simeq 3\mu K^{2}$. More interesting is the bottom part of this Figure (with the same color-code) where this difference is divided by the uncertainties estimated on the hlpTT spectra. It shows that a sensitivity of few percent of a $\sigma$ over all the $\ell$ range has to be achieved in order to fit for a 0.3 eV neutrino mass (the example taken here).

The extreme case of the differences between linear and non-linear models of the CMB temperature power spectrum are also illustrated for $\mathrm{\Sigma m_{\nu}}$=0.3 eV: For CLASS, in orange, corresponding to Bird et al. (2012), and for CAMB; where two models are compared, Mead in red and Takahashi  in green (cf. Sect. 1). The plots show that the non-linear effects are of the order of 1$\mu$K and correspond to, at most, $\simeq 1\%$ of a $\sigma$. The difference between those estimations gives a hint towards the theoretical uncertainty associated to the propagation of non-linear effects. In addition to this, it must be kept in mind that when constraining extensions of $\mathrm{\Lambda{CDM}}$ models, all the cosmological parameters are correlated, such that those very small effects have to be disentangled from any other (more or less degenerated) parameter’s configuration.

To conclude, the effect one tries to fit on temperature power spectra to extract information on $\mathrm{\Sigma m_{\nu}}$ is very tiny, and spreads over the whole multipole range. It therefore requires one to master the underlying model used to build the CMB likelihood function to a very high accuracy.

3.2 $\mathrm{\nu\Lambda{CDM}(3\nu)}$

The profile likelihood results on $\mathrm{\Sigma m_{\nu}}$  derived from the 2013 Planck temperature power spectra have been compared with those obtained with a Bayesian analysis in Planck Collaboration Int. XVI (2014) in the $\mathrm{\nu\Lambda{CDM}(1\nu)}$ model. It was then shown that the profile likelihood shape was non-parabolic. We recover the same results with the 2015 data in the $\mathrm{\nu\Lambda{CDM}(3\nu)}$ model: This is illustrated for different high-$\ell$ likelihoods combined with lowTEB on Fig. 3.

Fig. 3 illustrates that the behaviour of the $\Delta\chi^{2}$ as a function of $\mathrm{\Sigma m_{\nu}}$ is almost independent of the choice of the likelihood. Still, the spread of the profile likelihoods gives an indication of the systematic uncertainties linked to this choice. For such particular shapes of the profile likelihood, one cannot simply use the Gaussian confidence level intervals detailed in Feldman & Cousins (1998): One should rely on extensive simulations to properly build the corresponding Neyman construction (Neyman 1937), and apply the FC ordering principle; this is beyond the scope of this work. We do not therefore quote any limit for non-parabolic profile likelihood.

The use of the $\mathrm{A_{\text{SZ}}}$  constraint (cf. Eq. 4) does improve the constraint on $\mathrm{\Sigma m_{\nu}}$This is further discussed in Sect. 3.4, together with the impact of the combination of the VHL data.

3.3 Impact of low-$\ell$ likelihoods

On Fig. 4 are shown several $\mathrm{\Sigma m_{\nu}}$ profile likelihoods corresponding to different choices for the low-$\ell$ likelihoods, while keeping hlpTT for the high-$\ell$ part. They all present the same shape which, as previously, prevents us from extracting upper bounds.

The results obtained when combining hlpTT with lowTEB (in blue) are very close to those obtained with a $\tau_{\rm reio}$ auxiliary constraint+Commander (in green), showing that with those datasets, the results do not significantly depend on the choice of the low-$\ell$ polarisation likelihood. The same conclusion can be derived from the comparison of the results obtained using hlpTT+$\tau_{\rm reio}$ auxiliary constraint (in red).

However, the difference between these two sets of profile likelihoods highlights the impact of Commander. A possible origin of this difference lies in the fact that when adding Commanderin $\mathrm{\Lambda{CDM}(3\nu)}$+$\mathrm{A_{L}}$, one reconstructs a higher $\mathrm{A_{L}}$ value. Indeed, with hlpTT+$\tau_{\rm reio}$, we get $\mathrm{A_{L}}=1.16\pm 0.11$, while we find $\mathrm{A_{L}}=1.20\pm 0.10$ for hlpTT+$\tau_{\rm reio}$+Commander, that is, a higher value with a similar uncertainty. This higher tension with regards to the $\mathrm{\Lambda{CDM}}$ model (which assumes $\mathrm{A_{L}}$=1) artificially leads to a tighter constraint on $\mathrm{\Sigma m_{\nu}}$ (we refer also to Sect. 4.4).

3.4 Impact of VHL data

It was suggested in Planck Collaboration XI (2016) to add a constraint on the SZ amplitudes to mimic the impact of VHL data, and we have shown in Fig. 3 that the use of such a constraint does tighten slightly the constraints on $\mathrm{\Sigma m_{\nu}}$.

In this Section, we try to go one step further by actually using the VHL data themselves to further constrain the foreground residuals amplitudes in the $\mathrm{\nu\Lambda{CDM}(3\nu)}$ case, using the same procedure as the one described in Couchot et al. (2017).

Fig. 5 shows the $\mathrm{\Sigma m_{\nu}}$ profile likelihoods obtained when combining hlpTT+lowTEB with VHL data in green: An apparent $\Delta\chi^{2}$ minimum shows up, around $\mathrm{\Sigma m_{\nu}}\sim 0.7,$  eV with a $\Delta\chi^{2}$ decrease with regards to $\mathrm{\Sigma m_{\nu}}=0$ of the order of two units. This is quite different from the Planck only $\mathrm{\Sigma m_{\nu}}$ profile likelihoods previously studied, even when the $\mathrm{A_{\text{SZ}}}$ constraint has been added (cf. Sect. 3.2). In the $\mathrm{\nu\Lambda{CDM}(1\nu)}$ model, we have checked that the shape of the profile is about the same but for the minimum, which is around $\mathrm{\Sigma m_{\nu}}$=0.4  eV, close to the results obtained by Di Valentino et al. (2013); Hou et al. (2014).

To investigate this particular behaviour, we must stress that, for the combination of Planck with VHL data, one needs to compute the CMB power spectra up to $\ell\simeq\ 5000$. We therefore need to control the foreground residuals modelling, the datasets intercalibration uncertainties, and the uncertainties on non-linear effects models over a very broad range of angular scales.

To tackle the issue of the foreground modelling, several settings have been studied. They are represented on Fig. 5. The blue profile likelihood is built while fixing all the foreground amplitude nuisance parameters to their mean expectation values. It can be compared with two other profile likelihoods (in cyan and in red), built when fitting only the SZ and the CIB templates amplitudes, respectively (these foregrounds dominate at the higher end of the $\ell$ range). The observed variations, regarding both the $\chi^{2}$ rise at low $\mathrm{\Sigma m_{\nu}}$  and the $\mathrm{\Sigma m_{\nu}}$ value at the minimum, with respect to the default case (in blue), show that our combination of Planck and VHL datasets is too sensitive to the foreground residuals modellings to be reliable for the derivation of a limit on $\mathrm{\Sigma m_{\nu}}$. This may also come from the fact that the modelling we have used for the full sky Planck surveys is not accurate enough for the VHL small patches of the sky.

We have also investigated the impact of the uncertainties on the modelisation of non-linear effects. The mean values of the cosmological parameters, derived from the best fits of the hlpTT+lowTEB+VHL for $\mathrm{\Sigma m_{\nu}}$=0.06  eV and for 0.7  eV, were used to compute the temperature $C_{\ell}$ spectra. We have observed that the difference between these spectra was of the same order of magnitude as the difference of spectra expected from two non-linear models for $\mathrm{\Sigma m_{\nu}}$=0.06  eV (namely between Takahashi and Mead cf. Sec. 1.5). As such a difference leads to a variation of up to 2 $\chi^{2}$ units, we could expect that the uncertainties on non-linear models would lead to similar $\chi^{2}$ differences222Still, the proper propagation of the uncertainties of non-linear effects is beyond the scope of this work.. In addition, it must be noted that this difference is also of the order of magnitude of the relative calibration between the different VHL datasets and Planck.

For all those reasons, we have chosen not to include the VHL datasets in the following (we refer also to Addison et al. (2016) for the tensions between VHL datasets and Planck). The potential impact of the uncertainties on non-linear models becomes negligible when one only considers CMB spectra up to $\ell\simeq 2500$ (e.g. for Planck-only data).

4 Adding BAO and SNIa data

As noted in Sect. 1.4, the main degeneracy when using CMB data to constrain flat $\nu\mathrm{\Lambda{CDM}}$ models, is between $\mathrm{\Sigma m_{\nu}}$ and $\Omega_{\Lambda}$ which are both related to the angular-diameter distance to the last scattering surface. This translates into a degeneracy between $\mathrm{\Sigma m_{\nu}}$ and the derived parameters $\sigma_{8}$ and $H_{0}$ as illustrated in Fig. 6. The effect of neutrino free-streaming on structure formation favours lower $\sigma_{8}$ values at large $\mathrm{\Sigma m_{\nu}}$, which in addition require one to lower $H_{0}$. Adding BAO and SNIa data breaks this relation, and substantially tightens the constraint on $\mathrm{\Sigma m_{\nu}}$. In this section, we analyse the combination of Planck CMB data with DR12 BAO and SNIa data (as described in Sect. 1).

4.1 hlpTT, hlpTTps, and PlikTT  comparison

Fig. 7 compares the three Planck likelihoods when they are combined with lowTEB, BAO and SNIa. The impressive improvement with respect to the Planck only results (Fig. 3) can be measured, for example, by the comparison of the range of $\mathrm{\Sigma m_{\nu}}$ values for which the $\Delta\chi^{2}$ is below 3. As expected, those results illustrate that most of the constraint on $\mathrm{\Sigma m_{\nu}}$ does not come from CMB-only data (at decoupling neutrinos act essentially as radiation) but from the combination with late-time probes (where they contribute as matter). In addition, for this combination of probes, the likelihood profiles take on a standard parabolic shape: The derived upper bounds on $\mathrm{\Sigma m_{\nu}}$, using the F.C. prescription, are summarised in Table 2. We also quote the $\mathrm{A_{L}}$ values obtained using the same datasets for the $\mathrm{\Lambda{CDM}}$+$\mathrm{A_{L}}$ model (fixing $\mathrm{\Sigma m_{\nu}}=0.06\thinspace\text{eV}$). We note that they differ from one by roughly 2$\sigma$. The impact on the $\mathrm{\Sigma m_{\nu}}$ limit is discussed in Sect. 4.4.

The profiles of the different high-$\ell$ likelihoods are very similar, giving confidence in the final results that can be derived from their comparison. The spread between the curves reflects the remaining systematic uncertainty linked to the choice of the underlying foreground modelling. We have checked that, for hlpTT and hlpTTps, removing the foreground nuisance parameter auxiliary constraints does not impact the results: This provides an additional proof that the model and the data are in very good agreement. The information added by the $\mathrm{A_{\text{SZ}}}$ constraint is of no use in this particular combination of data within the $\mathrm{\nu\Lambda{CDM}(3\nu)}$ model. The systematic uncertainty on the $\mathrm{\Sigma m_{\nu}}$ limit due the foreground modelling, deduced from this comparison, is therefore estimated to be of the order of 0.03  eV for this particular data combination.

As expected, the main improvement with respect to the Planck only case comes from the addition of the BAO dataset: the contribution on the $\mathrm{\Sigma m_{\nu}}$ limit of the addition of SNIa is of the order of $\simeq 0.01\leavevmode\nobreak\ \thinspace\text{eV}$.

4.2 Impact of low-$\ell$ likelihoods

While in the previous Section we focused on the estimation of the remaining systematic uncertainties linked to the choice of the high-$\ell$ likelihood, a comparison of the low-$\ell$ parts is now performed. We already discussed in Sect. 3.3 the impact of this choice on the results derived from CMB data only; this comparison focuses on the combination of BAO and SNIa data.

The results are summarised in Fig. 8. For the two HiLLiPOP likelihoods, tightening the constraints on $\tau_{\rm reio}$ with the use of $\tau_{\rm reio}$+Commander  in place of lowTEBresults in a limit of 0.15 eV (resp. 0.16 eV) for hlpTTps (resp. hlpTT) and amounts to a few $10^{-2}$ eV decrease compared to the lowTEB case. This decrease is a direct consequence of both the ($\mathrm{\Sigma m_{\nu}}$,$\tau_{\rm reio}$) correlation (Allison et al. 2015), and the smaller value of the reionisation optical depth constraint from $\sim$ 0.07 to 0.058 (Planck Collaboration Int. XLVII 2016).

4.3 Cross-check with TE

As pointed out in Galli et al. (2014) and Couchot et al. (2017), CMB temperature-polarisation cross-correlations (TE) give competitive constraints on $\mathrm{\Lambda{CDM}}$ parameters. The leading advantage of using only these data is that one depends very weakly on foreground residuals and therefore uncertainty linked to the model parametrisation is reduced. In practice, only one foreground nuisance parameter is required: The amplitude of the polarized dust. Nevertheless, the S/N being lower than in the TT case for Planck, a likelihood based on TE spectra is not competitive when constraining extensions to the six $\mathrm{\Lambda{CDM}}$ parameters. Indeed an estimation of the TE-only constraint on $\mathrm{\Sigma m_{\nu}}$ would lead to a limit higher than 1 eV. However, as soon as BAO data are added, one obtains a constraint competitive with TT as shown in Fig. 9. As in the TT case, all profile likelihoods are nicely parabolic, and the corresponding limits are summarised in Table 3.

As for temperature-only data, adding the SNIa data improves only very marginally the results up to 0.01  eV. Tests of the dependencies on the low-$\ell$ likelihoods have also been performed and an example is given in Table 3. As a final result, we obtain $\mathrm{\Sigma m_{\nu}}$¡0.20  eV  at 95% CL as strong as in the TT case, showing that the loss in signal over noise of TE (statistical uncertainty) is balanced by improved control of foreground modelling (systematic uncertainty).

4.4 $\mathrm{A_{L}}$ and $\mathrm{\Sigma m_{\nu}}$

4.4.1 $\mathrm{\nu\Lambda{CDM}(3\nu)}$ model

As previously stated, CMB data tend to favour a value of $\mathrm{A_{L}}$ greater than one. In the combination of Planck high-$\ell$ likelihood with lowTEB, BAO and SNIa, the $\mathrm{A_{L}}$ values estimated in the $\mathrm{\Lambda{CDM}(3\nu)}$+$\mathrm{A_{L}}$ model, are summarised in the third column of Table 2. As expected they are almost identical to the ones obtained with CMB data only.

The fact that $\mathrm{A_{L}}$ is not fully compatible with the $\mathrm{\Lambda{CDM}}$ model, has to be taken into account when stating final statements on $\mathrm{\Sigma m_{\nu}}$ since, otherwise, the results are not obtained within a coherent model: On one side we fix $\mathrm{A_{L}}$ to one by working within a $\nu\mathrm{\Lambda{CDM}}$ model while the data are, at least, $\simeq 2\sigma$ away from this value, and on the other side, fixing $\mathrm{A_{L}}=1$ results, artificially, in a tighter constraint on $\mathrm{\Sigma m_{\nu}}$. This last effect can be seen, for example, in Table 2, for which the higher the $\mathrm{A_{L}}$ value, the tighter the constraint on $\mathrm{\Sigma m_{\nu}}$.

There are two ways to propagate this effect on the $\mathrm{\Sigma m_{\nu}}$ limit determination. The first is to open up the parameter space to $\mathrm{\nu\Lambda{CDM}(3\nu)}$+$\mathrm{A_{L}}$ (as it is done in the Sect. 4.4.2). The second is to better constrain the lensing sector by considering the Planck lensing likelihood and then to fit only for the $\mathrm{\Sigma m_{\nu}}$ extension using the $\mathrm{\nu\Lambda{CDM}(3\nu)}$ model, fixing $\mathrm{A_{L}}=1$ (cf. Sect. 4.4.3).

4.4.2 The $\mathrm{\nu\Lambda{CDM}(3\nu)}$+$\mathrm{A_{L}}$ model

In this Section, we open the $\mathrm{\nu\Lambda{CDM}(3\nu)}$ parameter space to $\mathrm{A_{L}}$ for the combination of Planck high-$\ell$ likelihoods with lowTEB+BAO+SNIa.

The limits derived from the corresponding profile likelihoods are summarised in Table 4. The increase of the limits with respect to those of Table 2 results from two effects. First of all we open up the parameter space, propagating the uncertainty on $\mathrm{A_{L}}$ on the $\mathrm{\Sigma m_{\nu}}$ determination. The second effect is linked to the fact that, as already stated, the CMB data tend to favour a higher $\mathrm{A_{L}}$ value than expected within a $\mathrm{\Lambda{CDM}}$ model. We have observed that this effect propagates as an increase of the baryon energy density, a slight decrease of the cold dark matter energy density, and this shows up, with a fixed geometry, as a higher neutrino energy density. Those two combined effects drive the limit to high values of $\mathrm{\Sigma m_{\nu}}$ when fitting for both $\mathrm{\Sigma m_{\nu}}$ and $\mathrm{A_{L}}$.

4.4.3 Combining with CMB lensing

Another way of tackling the $\mathrm{A_{L}}$ problem is to add the lensing Planck likelihood to the combination (see Sect. 2.4). This allows us to obtain a lower $\mathrm{A_{L}}$ value, as shown in the third column of Table 5 in the $\mathrm{\Lambda{CDM}(3\nu)}$+$\mathrm{A_{L}}$ model. With this combination, the $\mathrm{A_{L}}$ value extracted from the data is fully compatible with the $\mathrm{\Lambda{CDM}}$ model, allowing us to derive a limit on $\mathrm{\Sigma m_{\nu}}$ together with a coherent $\mathrm{A_{L}}$ value.

As expected, in the $\mathrm{\Lambda{CDM}(3\nu)}$ model, the $\mathrm{\Sigma m_{\nu}}$ limits are therefore pushed toward higher values than what has been presented in Table 2: This is exemplified by the second column of Table 5.

4.5 Constraint on the neutrino mass hierarchy

As explained in Sect. 1.4, the neutrino mass repartition leaves a very small signature on the CMB and matter power spectra. In this section, we test whether or not the combination of modern cosmological data is sensitive to it.

We compare the results obtained with four configurations of neutrino mass settings. The first one corresponds to one massive and two massless neutrinos as in $\mathrm{\nu\Lambda{CDM}(1\nu)}$ and is labelled $[1\nu]$. The second one is built under the assumption of three mass-degenerate neutrinos as in $\mathrm{\nu\Lambda{CDM}(3\nu)}$ and is denoted $[3\nu]$. We also discuss the normal hierarchy [3$\nu$ NH] (resp. inverted hierarchy [3$\nu$ IH]) derived from Eq. 1 and Eq. 2 (resp. Eq. 3).

In contrast with the rest of this paper, we did not subtract, in this Section, the minimum of the $\chi^{2}$ to plot the profile likelihoods. This allows us to assess the $\chi^{2}$ difference between the various neutrino configurations. In Fig. 10, we show the results obtained using the combination hlpTT+lowTEB+BAO+SNIa+lensing. The 95% CL upper limits derived from these profile likelihoods are reported in Table 6.

The observed difference between [1$\nu$] and [$3\nu$] illustrates the impact of the choice of the number of massive neutrinos on the derived constraint on $\mathrm{\Sigma m_{\nu}}$. More important is the comparison of the profile likelihoods built for the different hierarchy scenarios. The fact that they are indistinguishable (both in shape and in absolute $\chi^{2}$ values), and, even more, that they are almost identical to the one of the three degenerated masses, shows that there is, with modern data, no hint of a preference for the data towards one scenario or another, for this particular data combination (we refer also to the latest discussion in Schwetz et al. 2017).

5 Adding CMB polarisation

In the previous Section, we derived limits on $\mathrm{\Sigma m_{\nu}}$  from various high-$\ell$ Planck temperature likelihoods combined with BAO and SNIa. All those results were cross-checked with the almost foreground-free TE Planck spectra. In this Section, we combine the temperature and polarisation CMB data from Planck together with BAO, SNIa. As done previously, the CMB lensing is then also added in the combination to address the $\mathrm{A_{L}}$ tension. We then show the final results of this paper on the $\mathrm{\Sigma m_{\nu}}$ determination.

5.1 Combination of TT, TE, EE BAO and SNIa

The 95% CL upper limits on $\mathrm{\Sigma m_{\nu}}$ corresponding to the full TT+TE+EE likelihoods (labelled ALL), combined with BAO, SNIa and lowTEB are summarised in Table 7.

They are very close to the temperature-only upper limit of Table 2, showing that the use of the polarisation information in addition to the temperature does not add much information. They are also very close, showing the consistency of the results with respect to the high-$\ell$ Planck likelihoods when BAO and SNIa are included.

However, for this data combination, we are still left with a 2$\sigma$ tension on $\mathrm{A_{L}}$ (the $\mathrm{A_{L}}$ values are almost the same as the ones of the TT combination of Table 2). The fact that the results from PlikALL are lower than those of HiLLiPOP is linked to the fact that the $\mathrm{A_{L}}$ value of Plik is higher than the one of HiLLiPOP. We will come back to this point in the next Section.

5.2 Combining with CMB lensing

As done in Sect. 4.4.3, we now add to the data combination, the lensing Planck likelihood (see Sect. 2.4). The corresponding profile likelihoods are shown in Fig. 11, and the results are given in Table 8 for $\mathrm{\nu\Lambda{CDM}(3\nu)}$ (i.e. with $\mathrm{A_{L}}=1$). To compare with Table 7, the $\mathrm{\Sigma m_{\nu}}$ limits are higher when lowTEB is used at low-$\ell$, but more robust with respect to the $\mathrm{A_{L}}$ issue thanks to the use of the lensing data. For the ALL case, in the $\mathrm{\Lambda{CDM}(3\nu)}$+$\mathrm{A_{L}}$ model we end up with a value of $\mathrm{A_{L}}$ compatible with one and very comparable with those of Table 5. The limits on $\mathrm{\Sigma m_{\nu}}$ are therefore not artificially lowered by an overly high $\mathrm{A_{L}}$ value. Even though we end up with upper limits that are pushed toward higher bounds if compared to those obtained without the lensing data, we insist on the fact that this data combination is compatible with the $\mathrm{\Lambda{CDM}}$ model.

When making use of the latest $\tau_{\rm reio}$ measurement, we almost recover the results of Table 7. We use the differences between the upper limits obtained with the three Planck likelihoods of Table 8 (last column) to estimate a systematic error coming from the foreground modelling of 0.01 eV.

Table 9 provides the $\chi^{2}=-2\log\cal{L}$ values as a function of $\mathrm{\Sigma m_{\nu}}$, where the likelihood ($\cal{L}$) has been profiled out over the nuisance and cosmological parameters. It corresponds to the combination of hlpALLps+BAO+SNIa+lensing, using the auxiliary constraint on $\tau_{\rm reio}$ combined with Commander at low-$\ell$. This dataset is chosen for the final limits derivation since it corresponds to the most up-to-date results on $\tau_{\rm reio}$. Table 9 can be used for neutrino global fits.

5.3 Cosmological parameters: $\mathrm{\Lambda{CDM}}$ versus $\nu\mathrm{\Lambda{CDM}}$

We compare the $\mathrm{\Lambda{CDM}}$ cosmological parameters and their error bars derived with the profile likelihood method using various combinations of CMB temperature+polarisation high-$\ell$ and low-$\ell$ likelihoods, with the CMB lensing likelihood from Planck, BAO and SNIa datasets.

More precisely, this comparison is done:

when $\mathrm{\Sigma m_{\nu}}$ is, or not, a free parameter, 2. 2.

using different foreground-modelling choices (via the different high-$\ell$ likelihoods), 3. 3.

switching from the publicly available lowTEB low-$\ell$ likelihood to the combination of an auxiliary constraint on $\tau_{\rm reio}$ with Commander, to size the impact of a tighter constraint on $\tau_{\rm reio}$, 4. 4.

between the neutrino mass settings of the $\mathrm{\Lambda{CDM}(1\nu)}$ and $\mathrm{\Lambda{CDM}(3\nu)}$ models.

These results are summarised in Fig. 12. They are very similar to the Planck 2015 results (Planck Collaboration XIII 2016) even though we are using here a new version of the BAO data (DR12). As stated in Sect. 1.6, we have checked that they do not depend on the chosen statistical approach (Bayesian or Frequentist), either for the $\mathrm{\Lambda{CDM}}$ or for the $\nu\mathrm{\Lambda{CDM}}$ model.

The values and uncertainties of the cosmological parameters in the $\mathrm{\nu\Lambda{CDM}(3\nu)}$ model (in red) are similar to those obtained in $\mathrm{\Lambda{CDM}(3\nu)}$  (in blue), but are marginally shifted and with slightly larger $68\%$ CL uncertainties. This is true with lowTEB  (as seen from the hlpALL results, circles) as well as with an auxiliary constraint on $\tau_{\rm reio}$ with Commander  for both hlpALL and hlpALLps (shown with squares). The increase of the uncertainties is related to the addition of $\mathrm{\Sigma m_{\nu}}$ in the fit. The small shifts of the mean values are nearly the same for all the tested cases. This could be the result of a best fit value of $\mathrm{\Sigma m_{\nu}}$ slightly different from 0.06  eVassumed in the $\mathrm{\Lambda{CDM}(3\nu)}$ model.

Switching from lowTEB (plain line on Fig. 12) to an auxiliary constraint on $\tau_{\rm reio}$ + Commander (dotted lines) at low-$\ell$ changes the results on $\tau_{\rm reio}$ and $\mathrm{A_{\rm s}}$ and reduces their uncertainties, as expected. We observe small shifts on other parameters ($\omega_{\rm b}$, $\omega_{\rm cdm}$, $n_{\rm s}$), consistently for all three high-$\ell$ likelihoods, when fitting for $\mathrm{\Sigma m_{\nu}}$. They result from intrinsic correlations between ($\tau_{\rm reio}$, $\mathrm{A_{\rm s}}$) and the other cosmological parameters.

In the six-parameter $\mathrm{\Lambda{CDM}(3\nu)}$ case, hlpALL and hlpALLps give very similar results, but for a small difference on $n_{s}$. This is related to the more constraining point source model (we refer to the discussion in Couchot et al. 2017). The comparison, illustrated in Fig 12, shows the robustness of the cosmological parameters estimation with respect to the choice of the CMB (high-$\ell$ and low-$\ell$) likelihoods. The residual (small) differences between them illustrate the remaining systematic uncertainties. For example, the differences between Plik and HiLLiPOP  could be linked to the different choices made for masks, $\ell$ ranges and foreground templates used in both cases.

Finally, the values and uncertainties of the cosmological parameters fitted in the $\mathrm{\Lambda{CDM}(3\nu)}$ and $\mathrm{\Lambda{CDM}(1\nu)}$, with PlikALL, are very close to each other. This shows that the mass repartition has almost no effect on $\mathrm{\Lambda{CDM}}$ parameters when $\mathrm{\Sigma m_{\nu}}$ is fixed to 0.06  eV.

6 Conclusions

We have addressed the question of the propagation of foreground systematics on the determination of the sum of the neutrino masses through an extensive comparison of results derived from the combination of cosmological data including Planck CMB likelihoods with different foreground modelisations.

For this comparison we have worked within the $\mathrm{\nu\Lambda{CDM}(3\nu)}$ model assuming three mass-degenerate neutrinos, motivated by oscillations results. We have justified this approximation, showing that it leads to the same results as those obtained when considering normal or inverted hierarchy.

We have shown that the details of the foreground residuals modelling play a non-negligible role in the $\mathrm{\Sigma m_{\nu}}$ determination, and affect the results in two different ways. Firstly, they are unveiled by different $\mathrm{A_{L}}$ values for the various likelihoods, up to $2\sigma$ away from $\mathrm{\Lambda{CDM}}$. This impacts the $\mathrm{\Sigma m_{\nu}}$ limit: The higher the $\mathrm{A_{L}}$ value favoured by the data, the lower the upper bound on $\mathrm{\Sigma m_{\nu}}$. For this reason we have added the CMB lensing in the combination of data, providing a way to derive a limit with an $\mathrm{A_{L}}$ value fully compatible with the $\mathrm{\Lambda{CDM}}$ model. Secondly, it introduces a spread of the profile likelihoods, resulting in various limits on $\mathrm{\Sigma m_{\nu}}$, from which a systematic uncertainty was derived. We have compared CMB temperature and polarisation results, as well as their combination, and showed that the results are very consistent between themselves.

We have also discussed the impact of the low-$\ell$ likelihoods. We have shown, through the use of an auxiliary constraint on $\tau_{\rm reio}$ (derived from the latest Planck reionisation results) combined with Commander, that a better determination of the uncertainty on $\tau_{\rm reio}$ led to a reduction of the upper limit on $\mathrm{\Sigma m_{\nu}}$, of the order of a few $10^{-2}$  eV with respect to the lowTEB case.

We have also addressed the question of the neutrino hierarchy. We have shown that the profile likelihoods are identical in the normal and inverted hierarchies, proving that the current data are not sensitive to the details of the mass repartition. Still, cosmological data could rule out the inverted hierarchy if they lead to a low-enough $\mathrm{\Sigma m_{\nu}}$ limit. However, today, the $\mathrm{\Sigma m_{\nu}}$ upper bound is still too high to get to this conclusion.

Combining the latest results from CMB anisotropies with Planck (both in temperature and polarisation, and including the last measurement of $\tau_{\rm reio}$), with BAO, SNIa, and the CMB lensing, we end up with:

[TABLE]

The values of the $\chi^{2}$ of the profile likelihoods are also given for further use in neutrinos global fits. For the first time, all the following effects have been taken into account:

• •

Systematic variations related to foreground modelling error,

• •

a value of $\mathrm{A_{L}}$ compatible with expectations,

• •

a lower value for $\tau_{\rm reio}$ compatible with the latest measurements from Planck,

• •

the new version of the BAO data (DR12),

making our final $\mathrm{\Sigma m_{\nu}}$ limit a robust result. For all these reasons, we think that this is the lowest upper limit we can obtain today using cosmological data.

As far as cosmology is concerned, the uncertainty on the neutrino mass will be improved in the future: It could be reduced by a factor $\simeq 5$ if one refers, for instance, to the forecasts on the combination of next-generation ‘Stage 4’ B-mode CMB experiments with BAO and clustering measurements from DESI (Audren et al. 2013; Font-Ribera et al. 2014; Allison et al. 2015; Abazajian et al. 2016; Archidiacono et al. 2017). Nevertheless, the proper propagation of systematics, in particular coming from the modelling of foregrounds, is a more important topic than ever in today’s cosmology.

Acknowledgements.

We gratefully acknowledge the IN2P3 Computer Center (http://cc.in2p3.fr) for providing the computing resources and services needed to this work.