© The Authors 2024

1. Introduction

Over the last two decades, the continuous improvement of the precision and accuracy of cosmological observations has opened a window to constraining neutrino properties via cosmological data. In particular, the sum of the three neutrino masses can be well constrained due to its effect of suppressing the clustering of cold dark matter after the neutrino non-relativistic transition.

Neutrino oscillation experiments demonstrated that at least two neutrinos are massive by measuring two squared mass differences Δm_ij² = m_i² − m_j² as

$$\begin{array}{cc}\hfill \mathrm{\Delta }{m}_{21}^2& =(7.{42}_{-0.20}^{+0.21})\times {10}^{-5}{\mathrm{eV}}^2,\hfill \\ \hfill \mathrm{\Delta }{m}_{3i}^2& =\{\begin{array}{cc}(2.{517}_{-0.028}^{+0.026})\times {10}^{-3}{\mathrm{eV}}^2\hfill & (NO),\hfill \\ (-2.{498}_{-0.028}^{+0.028})\times {10}^{-3}{\mathrm{eV}}^2\hfill & (IO),\hfill \end{array}\hfill \end{array}$$(1)

where the sign of Δm_3i² depends on the mass ordering, which can be either normal (NO: m₁ < m₂ < m₃, Δm_3i² = Δm₃₁² > 0) or inverted (IO: m₃ < m₁ < m₂, Δm_3i² = Δm₃₂² < 0), see Esteban et al. (2020), Gonzalez-Garcia & Yokoyama (2022). According to these results, the minimum value of the neutrino mass sum is either 0.06 eV¹ in NO or 0.10 eV in IO. However, neutrino oscillation experiments cannot constrain the absolute neutrino mass sum. On the other hand, β-decay experiments are sensitive to the effective electron anti-neutrino mass. Recently, the Karlsruhe Tritium Neutrino Experiment (KATRIN) constrained m_ν < 0.45 eV (90% CL, see Aker et al. 2024). Future β-decay experiments, such as Project 8 (Project 8 Collaboration 2022), could potentially reach a sensitivity of 40 meV.

Despite the great progress made in the precision of β-decay experiments², cosmology provides the most stringent (albeit model-dependent) constraints to date on the absolute neutrino mass scale. The Planck Collaboration reported an upper bound of ∑m_ν < 0.24 eV (95% CL) by combining the temperature,

polarisation, and lensing of the cosmic microwave background (CMB). Adding external baryon acoustic oscillation (BAO) data further improves the bound to ∑m_ν < 0.12 eV (95% CL, see Planck Collaboration VI 2020). Including Pantheon type Ia supernovae, baryon acoustic oscillations and redshift-space distortions from the Sloan Digital Sky Survey and weak-lensing measurements from the Dark Energy Survey leads to an upper bound of 0.11 eV (Alam et al. 2021). Recently, the DESI Collaboration reported an upper limit of ∑m_ν < 0.072 eV at the 95% CL from the combination of DESI BAO and CMB data. This bound mildly favours the NO scenario and shows that there might be a possibility of indirectly constraining the neutrino mass hierarchy using cosmological data. Current data on the full-shape galaxy power spectrum do not yet improve the constraint on ∑m_ν compared to the combination of CMB and the BAO geometric information (Ivanov et al. 2020a), but future galaxy surveys are expected to improve it. Finally, it is important to reiterate that changing the underlying cosmological model could potentially relax the constraints obtained from many of the cosmological probes (Lambiase et al. 2019; di Valentino et al. 2022).

The ability of cosmology to constrain the neutrino mass paved the way to investigate, via cosmological data, additional neutrino properties and particularly the number density of neutrinos. The standard model of particles physics predicts three active neutrino species, as confirmed by accelerator experiments (Mele 2015). In cosmology, this would correspond to an effective neutrino number N_eff = 3 in the instantaneous decoupling limit, or (when sticking to the same definition of this effective number) to N_eff = 3.044 when neutrino decoupling is modelled accurately (Froustey et al. 2020; Akita & Yamaguchi 2020; Bennett et al. 2021). Any deviation from this value is a hint at non-standard physics in the neutrino sector (e.g. low-temperature reheating, non-thermal corrections to the neutrino distribution generated by new physics, non-standard neutrino interactions or sterile neutrinos, Dvorkin et al. 2022; Archidiacono & Gariazzo 2022; Dasgupta & Kopp 2021), or exotic radiation components (e.g. axions or other forms of dark radiation, Marsh 2016; Kawasaki & Nakayama 2013; Di Luzio et al. 2020).

The continuous improvement of cosmological constraints on neutrino properties has opened up the question of whether upcoming cosmological surveys will be able to deliver a detection of a non-zero neutrino mass, to distinguish between the IO and NO scenarios due to a tight upper limit, and to investigate additional neutrino properties. One of the primary science goals of Euclid is to improve the cosmological constraints on the neutrino mass, possibly delivering evidence for a non-zero value. The aim of this work is to assess the sensitivity of Euclid to the neutrino mass, as well as its robustness against variations of the number of neutrinos or the modelling of dark energy.

Euclid (Euclid Collaboration 2024d) is a medium-class mission of the European Space Agency, which will map the Local Universe to improve our understanding of the expansion history and of the growth of structures. The satellite will observe roughly 14 000 deg² of the sky through two instruments, a visible imager (VIS, Euclid Collaboration 2024b) and Near-Infrared Spectrometer and Photometer (NISP, Euclid Collaboration 2024c), delivering the images of more than one billion galaxies and the spectra of between 20 and 30 million galaxies out to a redshift of about 2. The combination of spectroscopy and photometry will allow us to measure galaxy clustering and weak gravitational lensing, aiming at a 1% level of accuracy in the corresponding power spectra.

The paper is organised as follows. In Sect. 2, we review the main effects of massive neutrinos in cosmology and discuss their impact on the theoretical predictions for Euclid observables in Sect. 3. In Sect. 4 we introduce additional probes like cluster number counts from Euclid itself and external CMB data, which are potentially very synergistic with Euclid primary (weak lensing and galaxy clustering) probes. The methodology used to perform the forecast is described in Sect. 5, the code is validated in Sect. 6, and the results are presented in Sect. 7. Finally, we draw our conclusions on the sensitivity of Euclid to neutrino physics in Sect. 8. Along the main text, we point the reader to several Appendices A, B, C, D, E, F, G, and H, providing technical details on the analysis, as well as additional figures and consistency checks.

2. Neutrino cosmology

In the following, we briefly review the main effects of the neutrino masses and number density on structure formation. For a thorough description, we refer the readers to Bashinsky & Seljak (2004), Hannestad (2006), Lesgourgues & Pastor (2006, 2012), Lesgourgues et al. (2013), Lattanzi & Gerbino (2018), Archidiacono et al. (2017, 2020), and Lesgourgues & Verde (2022).

2.1. Massive neutrinos

Massive neutrinos play an important role in the distribution of large-scale structures. They decouple from the primordial plasma around T ∼ 1 MeV, while still being relativistic. The redshift z_nr at which they enter the non-relativistic regime depends on the individual mass, m_ν, of each neutrino as (1 + z_nr)∼2 × 10³(m_ν/1 eV). After decoupling, neutrinos travel with a thermal velocity, v_th, that defines the neutrino free-streaming wavenumber (Lesgourgues & Pastor 2006) as

$$\begin{array}{c}\hfill k_{\mathrm{FS}}={\left(\frac{4\pi G_{\mathrm{N}}\overline{\rho }{a}^2}{v_{\mathrm{th}}^2}\right)}^{1/2},\end{array}$$(2)

where $\overline{\rho }$ is the total background density, a is the scale factor, and G_ontN is the Newton constant. After the non-relativistic transition the thermal velocity decays as v_th = ⟨p⟩ m_ν⁻¹ ∝ a⁻¹, therefore, the free-streaming wavenumber goes through a minimum value of $k_{\mathrm{nr}}=k_{\mathrm{FS}}(z=z_{\mathrm{nr}})\sim \mathcal{O}({10}^{-2})h{\mathrm{Mpc}}^{-1}$ at the time of the non-relativistic transition; here h = H₀/(100 km s⁻¹ Mpc⁻¹) is the reduced Hubble parameter. Neutrinos cannot cluster in regions smaller than their free-streaming length (λ_FS = 2πa/k_FS). On scales larger than the maximum free-streaming length, corresponding to k < k_nr, massive neutrinos always behave as cold dark matter and the power spectrum of matter density fluctuations, P_mm(k), does not change with respect to the power spectrum of a cosmology that assumes the neutrinos are massless, but their non-relativistic density is added to that of cold dark matter. On the other hand, on scales smaller than the current value of the free-streaming length (k > k_FS⁰), massive neutrinos induce a suppression of the linear matter spectrum, P_mm(k), at a redshift of z = 0 by approximately, ΔP_mm(k)/P_mm(k)≃ − 8f_ν, where f_ν = Ω_ν, 0/Ω_m, 0 is the fraction of neutrino with respect to matter density (Hu et al. 1998; Lesgourgues et al. 2013). The reason is twofold: on scales where neutrinos are free-streaming, they do not contribute to gravitational clustering; and they slow down the growth of cold dark matter perturbations, which is given by a^{1 − 3f_ν/5} during matter domination (Bond et al. 1980). The second effect is responsible for the majority of the overall suppression and the scale dependence of the linear growth factor induced by massive neutrinos must be taken into account in the modelling of large-scale structure observables (see Sect. 3).

An important consideration is that the suppression of the power spectrum (as well as the effect of neutrino masses on the CMB spectrum) is almost the same independently of which neutrino ordering (IO or NO) is present in nature. This has been shown, for example, in Lesgourgues et al. (2004), Lesgourgues & Pastor (2006) and Archidiacono et al. (2020). Indeed, the cosmological observables are mainly sensitive to f_ν, that is, to Ω_ν, 0 h² ≈ ∑m_ν/(93.12 eV), and thus, to the summed neutrino mass, ∑m_ν (see Mangano et al. 2005)³.

The step-like suppression induced by neutrino free streaming is best seen when varying ∑m_ν while fixing the density of total matter, the density of baryonic matter, and the fractional density of the cosmological constant, namely, the parameters {Ω_m, 0h², Ω_b, 0h², Ω_Λ}. This is illustrated in the top panel of Fig. 1. However, such a transformation also changes characteristic times and scales that are strongly constrained by CMB experiments, such as the redshift of equality between radiation and matter or the angular diameter distance to the last-scattering surface. To understand what is left for such experiments as Euclid to measure, given our knowledge of the CMB spectrum, it makes more sense to look at the variation in the matter power spectrum when floating ∑m_ν, while fixing the redshift of equality, baryonic matter density, and angular scale of the sound horizon at decoupling (Archidiacono et al. 2017; Lesgourgues & Verde 2022). In particular, this means fixing the parameters {Ω_c, 0h², Ω_b, 0h², θ_s}, where Ω_c, 0 accounts for the fractional density of cold dark matter only. We note that neutrinos with a realistic mass were still relativistic at the time of radiation-to-matter equality, so that for a fixed CMB temperature and effective neutrino number, the redshift of equality is fixed by (Ω_c, 0 + Ω_b, 0)h² – and not by Ω_m, 0h² = (Ω_b, 0 + Ω_c, 0 + Ω_ν, 0)h². In that case, the effect of the neutrino mass is displayed in the bottom panel of Fig. 1. The increase in ∑m_ν with a fixed θ_s implies a decrease of h that suppresses the large-scale power spectrum roughly by the same amount as neutrino free-streaming suppresses the small-scale power spectrum. As a result, in a combined fit to Euclid data and CMB data from (for instance) Planck, the neutrino mass would manifest itself mainly as an overall decrease of the amplitude of the matter power spectrum relative to that measured at the time of decoupling by CMB data. This decrease is slightly redshift dependent: the smaller the redshift, the more the power spectrum gets suppressed. Besides that, variations in ∑m_ν induce a small shift in BAO scales, responsible for the oscillations shown in the bottom panel of Fig. 1.

Fig. 1. Effects of the neutrino mass on the linear total matter (solid) or cold-plus-baryonic matter (dashed) power spectrum, presented as ratios with respect to the ΛCDM spectrum with massless neutrinos for four different values of the summed neutrino mass. Top: Singling out the effect of neutrino free-streaming, we keep the parameters {Ωm, 0h2, Ωb, 0h2, ΩΛ} fixed. Bottom: Quantities best constrained by CMB data, that is, {zeq, Ωb, 0h2, θs}, are fixed to show what is left for Euclid to measure. The mass is always equally split between the three neutrino species.

The Euclid weak-lensing probe traces total matter, whereas (to a very good approximation) the galaxy clustering probe traces only the fluctuations of cold and baryonic dark matter (Castorina et al. 2014). The two panels in Fig. 1 also show the impact of varying ∑m_ν on the cold-plus-baryonic matter power spectrum, P_cc (dashed line). The effect of the neutrino mass on the amplitude of the linear small-scale spectrum, P_cc, is slightly reduced compared to the total matter case, with a suppression given roughly by ΔP_cc(k)/P_cc(k)≃ − 6f_ν.

2.2. Number of relativistic species

The contribution of ultra-relativistic species to the background density of radiation, ρ_r (including that of neutrinos before their non-relativistic transition) can be parameterised in terms of an effective number of neutrinos N_eff as

$$\begin{array}{c}\hfill {\rho }_{\mathrm{r}}=[1+\frac{7}{8}{\left(\frac{4}{11}\right)}^{4/3}{N}_{\mathrm{eff}}]{\rho }_{\gamma },\end{array}$$(3)

where ρ_γ stands for the photon background density. In the above expression, the factor 7/8 accounts for the Fermi–Dirac statistics of neutrinos, and (4/11)^4/3 accounts for the neutrino-to-photon temperature ratio in the instantaneous decoupling approximation (Lesgourgues & Pastor 2006). We would expect N_eff = 3 for three standard model neutrinos that thermalised in the early Universe and decoupled well before electron-positron annihilation. However, theoretical predictions set N_eff = 3.044 in the standard cosmological model (Mangano et al. 2002; Froustey et al. 2020; Bennett et al. 2021; Drewes et al. 2024) because neutrinos decouple gradually with residual scatterings during this time⁴. This value would change in the presence of non-standard neutrino features (see for example Lesgourgues et al. 2013, for the case of a large leptonic asymmetry, non-thermal distortions from low-temperature reheating, non-standard interactions, etc.), or additional relativistic relics contributing to the energy budget (see for example Dvorkin et al. 2022).

At the level of background cosmology all these effects and models are captured by a single parameter N_eff. However, at the level of perturbations some model-dependent features may arise, for instance if the additional relics carry small masses or feature (self-)interactions. In this work, for simplicity, we stick to the class of models where additional non-relativistic relics are decoupled, massless and free-streaming. However, we must simultaneously take the effect of neutrino masses into account – in most of our forecasts, {∑m_ν, N_eff} are two free model parameters. In this work, we will consider two ways to distribute the total mass ∑m_ν over the different species.

• (a)

  By default, we considered three active neutrino species degenerate in mass (m_ν = ∑m_ν/3), thermally distributed and with a temperature chosen in such a way that each species contributes to N_eff by 3.044/3. Then, the additional free-streaming and massless dark radiation mentioned before enhances the total effective neutrino number as N_eff = 3.044 + ΔN_eff with ΔN_eff ≥ 0. In particular, this assumption is used throughout the results presented in Sect. 7. Considering three degenerate massive neutrinos offers the advantage of predicting a matter power spectrum nearly indistinguishable from that of the realistic NO and IO scenarios with the same total mass, while models with only one or two massive species provide poorer approximations (Lesgourgues et al. 2004; Lesgourgues & Pastor 2006; Archidiacono et al. 2020).

• (b)

  For the specific validation presented in Sect. 6, we adopted the same model as in most previous forecasts, both for the sake of comparison with earlier work and because the calculation of Fisher matrices is easier when parameters can vary symmetrically around their fiducial value. Then, like in the baseline analysis of Planck Collaboration VI (2020), we stuck to a single massive neutrino (m_ν = ∑m_ν), with the same thermal distribution and temperature as in the previous model. We considered additional free-streaming and massless species, contributing to the effective neutrino number by N_eff − 3.044/3, such that N_eff can be either bigger or smaller than 3.044.

We note that some of the models mentioned before (that feature non-standard physics in the neutrino sector) would require other model-dependent schemes for the splitting of the total mass across species and for the phase-space distribution of each massive species. However, the sensitivity of Euclid to the parameters ∑m_ν and N_eff is expected to depend only weakly on each particular scheme, such that our sensitivity forecast assuming case (a) is representative of other cases.

The effect of varying N_eff on the linear matter power spectrum depends on which other parameters are kept fixed. For instance, the top panel of Fig. 2 shows the impact of increasing N_eff with fixed Ω_m, 0, Ω_b, 0, and h. The leading effect is then a shift in the redshift of equality, inducing a suppression of the power spectrum. However, like for neutrino masses, in order to understand what is left for Euclid to measure, it is interesting to fix the cosmological parameters that are best constrained by CMB data, namely the redshift of radiation-to-matter equality z_eq, the baryon matter density, and the sound horizon angular scale (Hou et al. 2013; Lesgourgues et al. 2013; Lesgourgues & Verde 2022). With such a choice, the effect of varying N_eff is displayed in the bottom panel of Fig. 2. Since z_eq is given by the matter-to-radiation density ratio Ω_m, 0/Ω_r, 0, increasing the radiation density with a fixed baryon matter density implies an increase in the cold dark matter density, and thus, a decrease of the baryon-to-cold matter density ratio Ω_b, 0/Ω_c, 0. This has two well-known consequences: an enhancement of the amplitude of the matter power spectrum on scales k > k_eq, where $k_{\mathrm{eq}}\sim \mathcal{O}({10}^{-2})h{\mathrm{Mpc}}^{-1}$ is the wavenumber crossing the Hubble radius at equality; and a decrease in the amplitude of BAO oscillations. The scale of BAO peaks is also slightly shifted due to an enhanced neutrino drag effect (Bashinsky & Seljak 2004; Lesgourgues et al. 2013; Baumann et al. 2019). In the bottom panel of Fig. 2, one can clearly identify the power enhancement at k > k_eq, as well as the oscillations produced by the shift in BAO phase and amplitude. The small power suppression on large scales comes from the fact that the matter fractional density, Ω_m, 0, which controls the overall amplitude of the matter power spectrum on those scales, decreases slightly when we increase Ω_m, 0h², while fixing θ_s.

Fig. 2. Effect of the effective number of relativistic degrees of freedom Neff = 3.044 + ΔNeff on the linear total matter power spectrum, presented as ratios with respect to the ΛCDM spectrum with Neff = 3.044. Top: Density and Hubble parameters {Ωm, 0, Ωb, 0, h} fixed only. Bottom: Quantities best constrained by CMB data, that is, {zeq, Ωb, 0h2, θs}, are fixed to show what is left for Euclid to measure. Here, we assume only massless neutrinos.

3. Theoretical predictions for Euclid observables

The modelling of Euclid observables follows the recipes presented in Euclid Collaboration (2020), EC20 hereafter, and subsequently updated in Euclid Collaboration (2024a). However, the equations must be modified to account for the presence of massive neutrinos. Below we will briefly review the main equations, highlighting the relevant differences needed in massive neutrino cosmologies.

3.1. Photometric survey

The information from the Euclid imaging and photometric surveys can be embedded in three primary observables: the weak gravitational lensing (hereafter WL); the photometric reconstruction of galaxy clustering (hereafter GC_ph); and their cross-correlation. For the purpose of our forecast, we decompose the observables into spherical harmonics leading to 2D angular power spectra (C_ℓ). Assuming the Limber (see also Kilbinger et al. 2017) and the flat-sky approximations, the equation for the C_ℓ is expressed as

$$\begin{array}{cc}\hfill C_{\mathit{ij}}^{\mathit{XY}}(\ell )& =c{\displaystyle {\int }_{z_{\min }}^{z_{\max }}}\mathrm{d}z\frac{W_i^X(k_{\ell },z)W_j^Y(k_{\ell },z)}{H(z)r^2(z)}{P}_{\delta \delta }^{\mathit{XY}}(k_{\ell },z)\hfill \\ \hfill & +N_{\mathit{ij}}^{\mathit{XY}}(\ell ),\hfill \end{array}$$(4)

where X and Y each stand either for L (referring to WL) or G (referring to GC_ph), i and j denote the redshift bins⁵, and H(z) is the Hubble rate. The nonlinear power spectrum of density fluctuations is denoted by $P_{\delta \delta }^{\mathit{XY}}(k,z)$ and is evaluated at k_ℓ = (ℓ+1/2)/r(z) (Kilbinger et al. 2017), with r(z) being the comoving distance. Finally, W_i^X(z) is the window function of the i-th redshift bin at z for the X observable⁶, which can be written as

$$\begin{array}{cc}\hfill W_i^{\mathrm{L}}(k,z)& =\frac{3}{2}{\mathrm{\Omega }}_{\mathrm{m},0}\frac{H_0^2}{c^2}(1+z)r(z){\displaystyle {\int }_z^{z_{\max }}}\mathrm{d}{z}^{\prime }\frac{n_i(z^{\prime })}{{\overline{n}}_i}\frac{r(z^{\prime }-z)}{r(z^{\prime })}\hfill \\ \hfill & +W_i^{\mathrm{IA}}(k,z),\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill W_i^{\mathrm{G}}(k,z)& =\frac{H(z)}{c}{b}_i(k,z)\frac{n_i(z)}{{\overline{n}}_i},\hfill \end{array}$$(6)

for WL and for GC_ph, respectively. Here, $n_i(z)/{\overline{n}}_i$ is the normalised galaxy distribution and b_i(k, z) is the galaxy bias in the i-th redshift bin. The contribution of the intrinsic alignment is embedded in W_i^IA(k, z) and we adopt the modelling through nonlinear alignment method, the so-called extended nonlinear alignment (eNLA) model, used in EC20⁷. The intrinsic alignment window function is then expressed as

$$\begin{array}{c}\hfill W_i^{\mathrm{IA}}(k,z)=-\frac{{\mathcal{A}}_{\mathrm{IA}}{\mathcal{C}}_{\mathrm{IA}}{\mathrm{\Omega }}_{\mathrm{m},0}{\mathcal{F}}_{\mathrm{IA}}(z)}{D(z,k)}\frac{n_i(z)}{{\overline{n}}_i(z)}\frac{H(z)}{c},\end{array}$$(7)

where D(z, k) is the scale-dependent linear growth factor. The function ℱ_IA(z) depends on the luminosity function and is given by

$$\begin{array}{c}\hfill {\mathcal{F}}_{\mathrm{IA}}(z)={(1+z)}^{{\eta }_{\mathrm{IA}}}{\left[\frac{⟨L⟩(z)}{L_{\star }(z)}\right]}^{{\beta }_{\mathrm{IA}}},\end{array}$$(8)

where ⟨L⟩(z) and L_⋆(z) are the redshift-dependent mean and characteristic luminosities of source galaxies. The parameters 𝒜_IA and η_IA are nuisance parameters, and are allowed to vary around the fiducial values {𝒜_IA, η_IA}={1.72, −0.41}; the parameters 𝒞_IA = 0.0134 and β_IA = 2.17 are kept fixed.

Finally, $N_{\mathit{ij}}^{\mathit{XY}}$ is the shot-noise term, which is zero for the cross-correlation between observables [$N_{\mathit{ij}}^{\mathrm{GL}}(\ell )=0$], while for the auto-correlation it can be written as

$$\begin{array}{c}\hfill N_{\mathit{ij}}^{\mathrm{GG}}(\ell )=\frac{1}{{\overline{n}}_i}{\delta }_{\mathit{ij}},N_{\mathit{ij}}^{\mathrm{LL}}(\ell )=\frac{{\sigma }_{ϵ}^2}{{\overline{n}}_i}{\delta }_{\mathit{ij}},\end{array}$$

where ${\overline{n}}_i$ is the average number of galaxies per redshift bin expressed in steradians and obtained by dividing the expected total number of observed galaxies, $n_{\mathrm{gal}}=30{\mathrm{arcmin}}^{-2}$, by the number of redshift bins, N_bins = 10, and σ_ϵ = 0.30 is the variance of the observed ellipticities. Following EC20, we neglect any subdominant contributions from redshift-space distortions and lensing magnification.

The primary impact of the neutrinos can be summarised through their scale-dependent growth and their impact on the expansion history. Obviously, the neutrino mass enters as an additional component in the computation of H(z). Besides this trivial difference, the first relevant difference in the modelling of Euclid observables in massive neutrino cosmologies with respect to ΛCDM concerns $P_{\delta \delta }^{\mathit{XY}}(k_{\ell },z)$. Indeed, it has been shown (Castorina et al. 2014, 2015) that in the presence of massive neutrinos the tracer of galaxy clustering is given by the clustering of cold dark matter and baryons, neglecting the contribution of neutrinos. Defining the power spectrum in terms of cold dark matter only [$P_{\delta \delta }^{\mathrm{GG}}(k_{\ell },z)=P_{\mathrm{cc}}(k_{\ell },z)$ with c = CDM + baryons], makes the bias less scale dependent at linear and at mildly nonlinear scales. Therefore, we can approximate b_i(k, z)≃b_i(z) in Eq. (6) and we can model the bias with only one nuisance parameter for each redshift bin. We still take into account the scale-dependent growth in all relevant terms, as, for example, in Eq. (7). An erroneous definition of the bias in terms of total matter power spectrum in massive neutrino cosmologies leads to an overestimate of the sensitivity of future galaxy surveys to ∑m_ν (Vagnozzi et al. 2018).

On the other hand, massive neutrinos do contribute to the gravitational potential, which is the source of the weak lensing effects. Therefore, we set $P_{\delta \delta }^{\mathrm{LL}}(k_{\ell },z)=P_{\mathrm{mm}}(k_{\ell },z)$. For the cross-correlation of these two probes, we assume $P_{\delta \delta }^{\mathrm{GL}}(k,z)=\sqrt{P_{\delta \delta }^{\mathrm{GG}}(k,z)P_{\delta \delta }^{\mathrm{LL}}(k,z)}$. Finally, we stress that the recipe used here to model the photometric observables implies that all the power spectra are nonlinear; thus, the nonlinear corrections described in Sect. 3.3 were applied to both P_mm(k_ℓ, z) and P_cc(k_ℓ, z).

The final log-likelihood can be modelled as a Gaussian with respect to the angular power spectrum and written as

$$\begin{array}{c}\hfill {\chi }^2=f_{\mathrm{sky}}{\displaystyle \sum _{\ell }}(2\ell +1)(\frac{d_{\ell }^{\mathrm{mix}}}{d_{\ell }^{\mathrm{th}}}+ln\frac{d_{\ell }^{\mathrm{th}}}{d_{\ell }^{\mathrm{fid}}}-N),\end{array}$$(9)

with

$$\begin{array}{cc}\hfill d_{\ell }^{\mathrm{th}}& =det\left[C_{\mathit{ij}}(\ell )\right],\hfill \\ \hfill d_{\ell }^{\mathrm{fid}}& =det\left[C_{\mathit{ij}}^{\mathrm{fid}}(\ell )\right],\hfill \\ \hfill d_{\ell }^{\mathrm{mix}}& ={\displaystyle \sum _k}det\left[\{\begin{array}{cc}{C}_{\mathit{ij}}(\ell )\hfill & \mathrm{if}j\ne k\hfill \\ C_{\mathit{ij}}^{\mathrm{fid}}(\ell )\hfill & \mathrm{if}j=k\hfill \end{array}\right].\hfill \end{array}$$(10)

Here ‘fid’ denotes the values of C_ij(ℓ) computed for the fiducial model parameters, and f_sky = 0.3636 is the sky fraction covered by the wide survey⁸. Following EC20 and Euclid Collaboration (2024a), we assume two different sets of specifications for photometric probes, a pessimistic and an optimistic one, as listed in Table 1. We will use both settings in our validation tests, but in Sect. 7 all our final forecast results rely solely on the pessimistic settings, for which multipoles are included in the likelihood only up to ℓ_max = 1500. We note that at higher multipoles, it becomes important to include baryonic feedback effects in the modelling of the observable power spectrum. These effects could be potentially degenerate with those of the neutrino mass and number density. However, as hinted in previous works as Spurio Mancini & Bose (2023) and as explicitly shown in Appendix H, this is not the case when the data is cut at ℓ_max = 1500. Thus, we neglect baryonic feedback in the rest of this work.

3.2. Spectroscopic survey

The observable extracted from the Euclid spectroscopic survey is a 3D galaxy power spectrum, which can be written as:

$$\begin{array}{cc}\hfill P_{\mathrm{obs}}(k_{\mathrm{fid}},{\mu }_{\mathrm{fid}};z)& =\frac{1}{q_{\perp }^2(z)q_{\Vert }(z)}\hfill \\ \hfill & \times \left\{\frac{{[b{\sigma }_8(z)+f(k,z){\sigma }_8(z){\mu }^2]}^2}{1+f{(k,z)}^2{k}^2{\mu }^2{\sigma }_{\mathrm{p}}^2(z)}\right\}\hfill \\ \hfill & \times \frac{P_{\mathrm{dw}}(k,\mu ;z)}{{\sigma }_8^2(z)}{F}_z(k,\mu ;z)+P_{\mathrm{s}}(z),\hfill \end{array}$$(11)

where μ = k · r/(kr) is the cosine of the line-of-sight angle with respect to the wavenumber k (with absolute value k = |k|), and the subscript “fid” denotes the quantities computed with the fiducial cosmology (which is used to convert angles and redshifts to physical distances). In the following, we will briefly review the modelling of the main observational effects that are taken into account in Eq. (11) to convert the cold dark matter and baryon power spectrum P_cc(k_ℓ, z) into the observed galaxy power spectrum P_obs(k_fid, μ_θ, fid; z). As already explained in Sect. 3.1, since the tracers of this observable are galaxies, we always refer to P_cc(k_ℓ, z), thereby removing the contribution of neutrinos, rather than refer to the total matter power spectrum. We briefly summarise the impact of the neutrinos on the various terms below.

The first term on the right-hand side is the Alcock–Paczyński effect (Alcock & Paczynski 1979), arising from the assumption about the underlying cosmology applied in the conversion of redshifts and angles into parallel and perpendicular distances as

$$\begin{array}{c}\hfill q_{\perp }(z)=\frac{D_{\mathrm{A}}(z)}{D_{\mathrm{A},\mathrm{fid}}(z)},q_{\Vert }(z)=\frac{H_{\mathrm{fid}}(z)}{H(z)},\end{array}$$(12)

where D_A(z) is the angular diameter distance, and H(z) is the Hubble rate. Both quantities are obviously affected by the impact of neutrinos on the expansion rate. We also note that as usual μ = μ_fidq_⊥/q_∥/G as well as k = k_fid G/q_⊥, where G² = 1 + μ_fid²(q_⊥²/q_∥² − 1). Here ‘fid’ denotes that the values are set to their fiducials, following the recipe of EC20 and Euclid Collaboration (2024a).

The term in the numerator in curly brackets, [b(z)σ₈(z)+f(k,z)σ₈(z)μ²], accounts for redshift-space distortions, which are anisotropic perturbations appearing in redshift-space, due to the Doppler effect being an additional source of redshift beyond the cosmological one. The effect is modelled according to the Kaiser formula (Kaiser 1987). Here b(z) is the bias (see Sect. 3.1 for a justification of dropping the scale dependence) and f(k, z) is the scale-dependent growth factor of CDM+baryons, explicitly excluding neutrinos since we are interested in clustered objects (galaxies). The growth factor is computed with respect to the CDM+baryon component as $f(z,k)=\frac{1}{2}\frac{\text{d}ln\left[P_{\mathit{cc}}(z,k)\right]}{\text{d}lna}$⁹. Finally, σ₈(z) is defined also using CDM+baryons only. As such, the entire numerator is impacted by the mass of the neutrinos through the scale dependence of the growth as well as the reduction of the amplitude, σ₈(z).

The term in the denominator in curly brackets, [1+f(k,z)²k²μ²σ_p²(z)], represents the Fingers of God effect, due to the additional redshift coming from galaxy peculiar velocities. The effect was modelled as a Lorentzian factor, where σ_p(z) is the distance dispersion corresponding to the velocity dispersion σ, or explicitly σ_p(z) = σ/[H(z)a(z)]. Given the uncertainty in the modelling and in the redshift dependence of σ_p(z), we treated it as four additional nuisance parameters: one for each redshift bin.

P_dw(k, μ_θ; z) is the partially de-wiggled power spectrum, computed starting from the linear $P_{\delta \delta }^{\mathrm{GG}}$¹⁰, which is P_cc, and accounting for the smearing of the baryon acoustic oscillations due to nonlinear effects; it can be written as (Wang et al. 2013):

$$\begin{array}{cc}\hfill P_{\mathrm{dw}}(k,\mu ;z)& =P_{\delta \delta }^{\mathrm{GG}}(k;z){\mathrm{e}}^{-g_{\mu }{k}^2}\hfill \\ \hfill & +P_{\mathrm{nw}}(k;z)(1-{\mathrm{e}}^{-g_{\mu }{k}^2}),\hfill \end{array}$$(13)

where P_nw is the no-wiggle power spectrum obtained by removing the BAO from $P_{\delta \delta }^{\mathrm{GG}}$ (Boyle & Komatsu 2018), and

$$\begin{array}{c}\hfill g_{\mu }(k,\mu ,z)={[{\sigma }_{\mathrm{v}}(z)]}^2\{1-{\mu }^2+{\mu }^2{[1+f^{\mathrm{fid}}(z,k)]}^2\},\end{array}$$(14)

where σ_v(z), which has dimensions of length, reflects the galaxy velocity dispersion and is being treated as four additional nuisance parameters, one for each redshift bin, as for σ_p(z).

The term F_z(k, μ; z) = exp[−k²μ²σ_r²(z)] represents the damping due to the spectroscopic redshift errors along the line of sight, where ${\sigma }_{\mathrm{r}}(z)=\frac{c}{H(z)}{\sigma }_z$ with a redshift-independent error σ_z = 0.002.

The last term in Eq. (11), P_s(z), is the shot noise caused by the Poissonian distribution of measured galaxies on the smallest scales due to the finite total number of observed galaxies. We modelled it in the same way as in EC20 and Euclid Collaboration (2024a), with one contribution given by the fiducial inverse number density in each bin, and a second contribution accounting for residual shot noise that we treat as a nuisance parameter in each bin.

The final likelihood is modelled as a Gaussian, comparing the observed power spectrum to a fiducial one. The χ² is given by

$$\begin{array}{cc}\hfill {\chi }^2& ={\displaystyle \sum _i}{\displaystyle {\int }_{-1}^1{\int }_{{10}^{-3}}^{k_{\max }}}\mathrm{d}{k}_{\mathrm{fid}}\mathrm{d}{\mu }_{\mathrm{fid}}{k}_{\mathrm{fid}}^2\frac{V_i^{\mathrm{fid}}}{8{\pi }^2}\hfill \\ \hfill & \times {\left\{\frac{P_{\mathrm{obs}}[k(k_{\mathrm{fid}},{\mu }_{\mathrm{fid}},z_i),\mu ({\mu }_{\mathrm{fid}},z_i),z_i]-P_{\mathrm{obs}}^{\mathrm{fid}}(k_{\mathrm{fid}},{\mu }_{\mathrm{fid}},z_i)}{P_{\mathrm{obs}}[k(k_{\mathrm{fid}},{\mu }_{\mathrm{fid}},z_i),\mu ({\mu }_{\mathrm{fid}},z_i),z_i]}\right\}}^2,\hfill \end{array}$$(15)

where i denotes the index of the redshift bin and V_i is the comoving volume of the spherical shell of the redshift bins probed by the experiment and k_max is in Table 2. The partial sky coverage of the survey is approximately taken into account through multiplication of the volume by a sky fraction f_sky = 0.3636.

Following EC20 and Euclid Collaboration (2024a) we assumed two different sets of specifications for spectroscopic galaxy clustering, as reported in Table 2.

Over the past few years, more accurate modelling of galaxy clustering including higher order perturbations using the effective field theory of large scale structures has been developed (Senatore & Zaldarriaga 2017; Schmittfull et al. 2021; Chudaykin et al. 2021) and applied to real data from existing surveys (Chen et al. 2022; D’Amico et al. 2020; Ivanov et al. 2020b). However, the result of a forecast based on synthetic data depends weakly on the specific modelling; additionally, our final results (Sect. 7) combine galaxy clustering with photometric probes, thus, further reducing the impact of the higher order terms on the prediction of the overall sensitivity of Euclid to neutrino parameters.

3.3. Nonlinear modelling

Modelling nonlinear clustering in massive neutrino cosmologies is essential to achieve robust constraints on their mass. A thorough comparison of the performance of several N-body codes and emulators is shown in Adamek et al. (2023). Their findings show that, in cosmologies where only the total neutrino mass is varied, the most up-to-date emulators (HMcode, EuclidEmulator2, BACCOemulator) agree with simulations at the 1 − 2% level for the matter power spectrum. However, none of the aforementioned emulators has been explicitly trained on or built for models with a varying number of neutrino-like species. We present below a comparison of these emulators to N-body simulations in order, firstly, to confirm that they accurately capture the effect of neutrino mass, and secondly, to check whether they can also account for the impact of varying N_eff.

In Fig. 3 (left column) we show the accuracy in the prediction of the total matter power spectrum of HALOFIT (Takahashi et al. 2012) with the neutrino corrections of Bird et al. (2012), HMcode (Mead et al. 2021), EuclidEmulator2 (Euclid Collaboration 2021a)¹¹, and BACCOemulator (Angulo et al. 2021). The accuracy is computed with respect to the Dark Energy and Massive Neutrino Universe (DEMNUni) N-body simulations (Carbone et al. 2016; Parimbelli et al. 2022) at redshifts z = 0, 1, 2. We find that HALOFIT is the least accurate in reproducing the results from the simulations, both at intermediate (BAO) scales and at smaller scales for any redshift. Instead, HMcode and the two emulators show a very similar performance, apart from a few exceptions. At z = 0 HMcode underestimates the power at small scales with respect to the two emulators, which describe the neutrino mass suppression with 1 − 2% accuracy. At z = 1 HMcode slightly overestimates power at BAO scales, while deeply in the nonlinear regime it seems to reproduce the results of the simulations more accurately. Finally, at z = 2 HMcode performs slightly better than EuclidEmulator2, while BACCOemulator is not yet trained up to this redshift. Overall, both HMcode and the emulators are within 2% accuracy at any redshift and at any scale where the simulations can be trusted (k < 1 h Mpc⁻¹).

Fig. 3. Relative percentage difference of the nonlinear power spectrum computed with various recipes with respect to the DEMNUni simulations for the case ∑mν = 0.16 eV. We note that here, contrary to the forecast analysis, and to be consistent with the DEMNUni simulations, we assume the total neutrino mass to be equally shared among the three neutrino species m1 = m2 = m3 = ∑mν/3. The left column shows the total matter power spectrum, while the right column shows the cold dark matter power spectrum. The spectra are evaluated at z = 0, 1, 2 (first, second and third rows, respectively). The theoretical predictions are provided by HALOFIT (dashed orange line), HMcode (solid blue line), EuclidEmulator2 (dashed red line), and BACCOemulator (dot-dashed green line). The purple shaded area represents the shot noise of the simulations, and the vertical dotted line the maximum wavenumber. In the third row the predictions of BACCOemulator are missing because the emulator is trained only up to z = 1.5. The predictions for the cold dark matter power spectrum (neglecting the contribution of neutrino perturbations) is computed according to the approximate formula Eq. (16) (see text for details).

In our usage case, it is not only important to accurately model the clustering of the total power spectrum but also of the CDM+baryon power spectrum (see Sect. 3). In the right column of Fig. 3 we show the accuracy of the same nonlinear recipes in predicting the cold dark matter and baryons power spectrum P_cc(k, z) extracted from the DEMNUni simulations. In order to convert the nonlinear P_mm(k) into the nonlinear P_cc(k) we use the formula

$$\begin{array}{c}\hfill P_{\mathrm{mm}}(k)=f_{\mathrm{c}}^2{P}_{\mathrm{cc}}(k)+2f_{\nu }{f}_{\mathrm{c}}{P}_{\mathrm{c}\nu }(k)+f_{\nu }^2{P}_{\nu \nu }(k),\end{array}$$(16)

where P_mm(k) is the total matter auto-correlation power spectrum, P_νν(k) the neutrino auto-correlation power spectrum, P_cν(k) the cross-correlation cold dark matter–neutrino power spectrum, f_c the cold fraction of dark matter, and f_ν = (1 − f_c) the hot one. The formula is correct as long as all the power spectra appearing both in the left-hand side and in the right-hand side are either linear or nonlinear. An approximation arises when mixing linear and nonlinear power spectra. Here we consider nonlinear P_mm(k) and P_cc(k), while P_cν(k) and P_νν(k) are assumed to be linear. This assumption is accurate for P_νν(k), because for ∑m_ν ≲ 0.6 eV the neutrino free-streaming length is larger than the nonlinear scale today (and even more so at higher redshifts). On the other hand, while we expect the cold dark matter–neutrino cross-power spectrum P_cν(k) to exhibit some nonlinearities, these are found to be negligible.

The right panel of Fig. 3 shows precisely the accuracy of this approximation: we compare results from the DEMNUni simulations against the nonlinear P_cc(k) of the various emulators obtained by inverting Eq. (16), where only P_mm(k) is assumed to be fully nonlinear. Concerning the accuracy on P_cc(k) of the different prescriptions adopted to compute the theoretical P_mm(k) [from which we derive P_cc(k)] the same considerations drawn for the total matter power spectrum hold true here. While HALOFIT fails to accurately reproduce the results of the simulations, the accuracy of the emulators is better than the one of HMcode, especially at low redshift. However, both HMcode and the emulators remain within 2% accuracy with respect to the DEMNUni simulations at any redshift and scale.

In Fig. 4, we show the massive neutrino (∑m_ν = 0.16 eV and ∑m_ν = 0.32 eV) induced suppression in the total matter power spectrum (left plot) and in the CDM+baryon power spectrum (right plot), with respect to a ΛCDM cosmology for the DEMNUni simulations and for the nonlinear predictions described above. For larger neutrino masses, HMcode outperforms not only HALOFIT but also the emulators in the precision of modelling this suppression. We note that a neutrino mass of about 0.3 eV, although it is excluded in the minimal ΛCDM + ∑m_ν scenario by the most stringent cosmological constraints to date (e.g. Alam et al. 2021; Palanque-Delabrouille et al. 2020), remains within reach in extended models (Lambiase et al. 2019).

Fig. 4. Theoretical predictions of the suppression of the total matter power spectrum (left) and of the cold dark matter power spectrum (right) in massive neutrino cosmologies (∑mν = 0.16 eV and ∑mν = 0.32 eV) with respect to the pure ΛCDM case with massless neutrinos at z = 0. The colour coding of the theoretical predictions is the same as in Fig. 3. In the case ∑mν = 0.32 eV, EuclidEmulator2 is not shown because the value of the neutrino mass is too far away from the range of validity of the emulator. The DEMNUni simulations are depicted as black pluses.

Finally, we are also going to vary the effective number of neutrino-like species (N_eff) in our main analysis. To understand the nonlinear clustering in cosmologies in this case, we performed simulations with 1024³ particles in a box with size L = 1024 h⁻¹ Mpc varying N_eff between 0.2 and 0.4. In Fig. 5 we checked the accuracy of HALOFIT, HMcode, EuclidEmulator2, and BACCOemulator with respect to these N_eff simulations. The two emulators fail in reproducing the variations of the number of neutrino species (green and red lines, BACCOemulator and EuclidEmulator2, respectively); this was expected since they are not trained on cosmologies with non-standard N_eff. In order to overcome the parameter extension with the emulators, we tried to remap the variation in N_eff on large scales to a variation in the total matter and cold dark matter density parameters {Ω_m, 0, Ω_c, 0} and of the reduced Hubble constant h using the relationship derived in Rossi et al. (2015). However, tweaking the parameters does not improve the accuracy of the reconstruction of the N_eff variations by means of emulators to a level competitive with HMcode. On the other hand, HMcode provides a good fit of the relative difference in the nonlinear clustering induced by variations of the number of neutrino-like particles, both for the phase-shift in the BAO scale, and for the overall suppression at small scales. Indeed, even though the HMcode parameters are not fit to simulations with varying N_eff, the model is based on the convolution with the linear power spectrum, which naturally embeds the effect of N_eff.

Fig. 5. Theoretical predictions of the suppression of the total matter power spectrum in non-standard Neff cosmologies with respect to the standard one, at z = 0, and assuming fixed values of {Ωm, 0, Ωb, 0, h}, like in the top panel of Fig. 2. The colour coding of the theoretical predictions is the same as in Fig. 3. We add here two more prescriptions, obtained by rescaling Ωc, 0, Ωb, 0 and h to mimic changes in Neff. We refer to it as “ΛCDM equivalent”: the dot-dashed pink line represents results for EuclidEmulator2; and the dotted cyan line represents the BACCOemulator prediction. Black crosses represent results from N-body simulations.

To summarise, the emulators reach 1% accuracy in reproducing the CDM+baryon power spectrum in massive neutrino cosmologies out to k = 1 h Mpc⁻¹, thus performing slightly better than HMcode, as already noted in Adamek et al. (2023). However, given their range of validity in terms of redshift (BACCOemulator) and neutrino mass (EuclidEmulator2) and especially given that they are not trained in varying N_eff cosmologies, we opted for HMcode as our recipe to compute nonlinear corrections.

4. Additional probes

Through the history of the Universe neutrinos evolve from behaving like radiation to contributing to the total matter density. In order to capture neutrino effects at different epochs and remove parameter degeneracies, it is crucial to combine primary probes with additional secondary probes on larger scales from the Euclid survey, such as cluster number counts, and with external data probing the early Universe, such as from CMB anisotropies.

4.1. Cluster number counts from the Euclid survey

Clusters of galaxies are potentially a strong cosmological probe (Oukbir & Blanchard 1997; White et al. 1993; Bahcall & Fan 1998; Reiprich & Boehringer 2002; Mantz et al. 2008; Vikhlinin et al. 2009). Measurements of their abundance and evolution allow us to place constraints on both the geometry of our Universe and the growth of its density perturbations (Mohr 2005; Vikhlinin et al. 2009; Allen et al. 2011; Kravtsov & Borgani 2012). This could be achieved in particular from galaxy cluster number-counts experiments where, assuming that a halo number density mass function is predicted in a given cosmological model, one can confront the number of observed clusters computed in a given survey volume with its theoretical prediction and, from it, constrain cosmological parameters (Bocquet et al. 2019; Kirby et al. 2019; Costanzi et al. 2021; Sakr et al. 2022; Lesci et al. 2022). Euclid will study galaxy clusters abundance, an independent and complementary probe to the two primary ones, allowing the optical detection of clusters previously unattainable in terms of the depth and area covered. Under optimistic assumptions about the calibration of the mass-observable relation, Carbone et al. (2012), Basse et al. (2014), Cerbolini et al. (2013), Sartoris et al. (2016) find that the promising increase in the number of detected clusters will provide tight constraints on models of dark energy or non-minimal massive neutrinos.

To forecast constraints from cluster number counts, we have followed a similar approach to Sartoris et al. (2016) on the modelling of the mass-observable scaling relation and selection function. The cluster’s mass is derived from a scaling relation with an observable quantity, such as the cluster richness, luminosity, velocity dispersion or shear from gravitational lensing. We call M_obs this mass. The estimated number counts of these clusters for a redshift bin l and mass bin m, corresponding to z_l and M_obs, m, can then be expressed as

$$\begin{array}{cc}\hfill N_{l,m}& ={\displaystyle {\int }_{{\mathrm{\Omega }}_{\mathrm{tot}}}}\mathrm{d}\mathrm{\Omega }{\displaystyle {\int }_{z_l}^{z_{l+1}}}\mathrm{d}z\frac{\mathrm{d}V}{\mathrm{d}z\mathrm{d}\mathrm{\Omega }}\hfill \\ \hfill & \times {\displaystyle {\int }_0^{+\infty }\mathrm{d}M\frac{\mathrm{d}n(M,z)}{\mathrm{d}M}}\frac{1}{2}[\mathrm{erfc}(x_{m,l})-\mathrm{erfc}(x_{m+1,l})],\hfill \end{array}$$(17)

under the assumption of a log-normal observed mass distribution. Here dn(M, z)/dM is the cluster mass function of Euclid Collaboration (2023) defined below, dΩ is the solid angle element in steradians, dV/(dz dΩ) is the derivative of the comoving volume with respect to the redshift and solid angle element,

$$\begin{array}{c}\hfill \frac{\mathrm{d}V}{\mathrm{d}z\mathrm{d}\mathrm{\Omega }}=\frac{c{(1+z)}^2{D}_{\mathrm{A}}^2(z)}{H(z)},\end{array}$$(18)

with D_A(z) being the angular diameter distance and c the speed of light, while erfc(x) is the complementary error function, with the argument x being the (biased) logarithm of the mass (see below). We can explicitly write this as x_m, l = x (M_{obs, m, l}) defined in each mass bin m and redshift bin l as

$$\begin{array}{c}\hfill x_{m,l}=x(M_{\mathrm{obs},m,l})=\frac{ln(M_{\mathrm{obs},m,l}/M)-{\mathcal{B}}_{\mathrm{bias},l}}{\sqrt{2{\sigma }_{M_{\mathrm{obs},m,l}}^2}}·\end{array}$$(19)

We define the bias and the variance to be

$$\begin{array}{c}\hfill {\mathcal{B}}_{\mathrm{bias},l}=B_{M,0}+b_{\mathrm{E}}ln(\frac{M}{M_{\mathrm{pivot}}})+{\alpha }_{\mathrm{E}}ln(1+z_l),\end{array}$$(20)

and

$$\begin{array}{c}\hfill {\sigma }_{M_{\mathrm{obs},m,l}}^2={\sigma }_{M_{\mathrm{obs}},0}^2-1+{(1+z_l)}^{2{\beta }_{\mathrm{E}}},\end{array}$$

where B_M, 0, b_E, and α_E are free parameters to quantify, respectively, the change in calibration, slope, and redshift dependence in the mass-observable-biased log-normal mass distribution, with M_pivot = 3 × 10¹⁴ h⁻¹ M_⊙, while σ_Mobs, 0² and β_E serve also to respectively calibrate and account for the redshift dependence in the variance. The predicted number density of halos of mass M at redshift z (the mass function) is given by

$$\begin{array}{c}\hfill \frac{\mathrm{d}n(M,z)}{\mathrm{d}M}=\frac{{\rho }_{\mathrm{bc},0}}{M}\mathcal{F}[\sigma (R,z)]\frac{\mathrm{d}ln([\sigma (R,z){]}^{-1})}{\mathrm{d}M},\end{array}$$(21)

where σ(R, z) is the variance of the density field within a sphere of radius R at redshift z, and it is computed following Costanzi et al. (2013) neglecting the massive neutrino component, from the linear CDM+baryon power spectrum P_cc(k, z) as

$$\begin{array}{c}\hfill {\sigma }^2(R,z)=\frac{1}{2{\pi }^2}{\displaystyle {\int }_0^{\infty }}\mathrm{d}k{k}^2{P}_{\mathrm{cc}}(k,z)W^2(kR),\end{array}$$(22)

where W(kR) is the top-hat filter in k-space,

$$\begin{array}{c}\hfill W(kR)=3\frac{[sin(kR)-kRcos(kR)]}{{(kR)}^3},\end{array}$$(23)

R = R(M) = (3M/4πρ_bc, 0)^1/3 is the radius of a sphere enclosing a mass M, and ρ_bc, 0 = ρ_crit, 0 Ω_bc, 0 is the mean CDM+baryon energy density at z = 0. Here we have used the critical density ρ_crit, 0 and the CDM+baryon density fraction Ω_bc, 0 = Ω_m, 0 − Ω_ν, 0. Finally, the multiplicity function reads, according to the modified Press–Schechter formalism of Euclid Collaboration (2023):

$$\begin{array}{c}\hfill \mathcal{F}(\nu ,z)=\nu A(p,q)\sqrt{\frac{2a}{\pi }}{\mathrm{e}}^{-a{\nu }^2/2}(1+\frac{1}{{(a{\nu }^2)}^p}){(\nu \sqrt{a})}^{q-1},\end{array}$$(24)

where ν(M, z) = δ_c(z)/σ(R, z), and δ_c(z) is the linear density contrast for spherical collapse. This can be computed following the prescription of Weinberg & Kamionkowski (2003):

$$\begin{array}{c}\hfill {\delta }_{c,\mathrm{lin}}(z)\simeq \frac{3}{20}{(12\pi )}^{2/3}[1+0.0123{\log }_{10}{\mathrm{\Omega }}_{\mathrm{m}}(z)].\end{array}$$(25)

Note that this formula does not take into account any massive neutrino, except for Ω_m. Here the matter density fraction Ω_m(z) at redshift z is computed from the present day total matter density Ω_m, 0 as

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{m}}(z)={\mathrm{\Omega }}_{\mathrm{m},0}{(1+z)}^3\frac{H_0^2}{H^2(z)}·\end{array}$$(26)

The redshift and scale-dependence of the parameters of Eq. (24) can be written as

$$\begin{array}{cc}\hfill A(p,q)& ={\left\{\frac{2^{-0.5-p+q/2}}{\sqrt{\pi }}[2^p\mathrm{\Gamma }(q/2)+\mathrm{\Gamma }(-p+q/2)]\right\}}^{-1},\hfill \end{array}$$(27)

$$\begin{array}{cc}\hfill q& =q_{\mathrm{R}}{[{\mathrm{\Omega }}_{\mathrm{bc}}(z)]}^{q_z},\hfill \end{array}$$(28)

$$\begin{array}{cc}\hfill q_{\mathrm{R}}& =q_1+q_2[\frac{\mathrm{d}ln\sigma (R,z)}{\mathrm{d}lnR}+0.5],\hfill \end{array}$$(29)

$$\begin{array}{cc}\hfill p& =p_1+p_2[\frac{\mathrm{d}ln\sigma (R,z)}{\mathrm{d}lnR}+0.5],\hfill \end{array}$$(30)

$$\begin{array}{cc}\hfill a& =a_{\mathrm{R}}{[{\mathrm{\Omega }}_{\mathrm{bc}}(z)]}^{a_z},\hfill \end{array}$$(31)

$$\begin{array}{cc}\hfill a_{\mathrm{R}}& =a_1+a_2{[\frac{\mathrm{d}ln\sigma (R,z)}{\mathrm{d}lnR}+0.6125]}^2,\hfill \end{array}$$(32)

where the parameters a_i, p_i, and q_i are calibrated to the simulation. The adopted values of the seven fitted parameters are listed in Table 3.

The corresponding likelihood function is based on Poisson statistics (Cash 1979; Holder et al. 2001; Bonamente 2020):

$$\begin{array}{cc}\hfill ln{\mathcal{L}}_{\mathrm{CC}}& =ln\mathcal{P}(n_{l,m}|N_{l,m})\hfill \\ \hfill & ={\displaystyle \sum _{l=1}^{N_l}\sum _{m=1}^{N_m}}[n_{l,m}ln(N_{l,m})-N_{l,m}-ln(n_{l,m}!)],\hfill \end{array}$$(33)

where 𝒫(n_l, m|N_l, m) is the Poisson distribution probability of finding n_l, m clusters given an expected number of N_l, m in each bin in redshift and mass. The N_l, m are computed from Eq. (17), whereas the n_l, m are the fiducial values of the N_l, m. To estimate the number counts, we have considered equally spaced redshift bins in z ∈ [0.2, 1.8], with a width of Δz = 0.1. As for the limiting mass, we have followed the selection function used in Sartoris et al. (2016) in the pessimistic case where the lower mass for clusters is defined as the one corresponding to the significance of detection threshold, or the ratio between the cluster galaxy number count and the field RMS, N_c/σ_field is above 5. In our analysis, we also conservatively vary the nuisance parameters. In order to compute the fiducial mocks the latter are set to (B_M,0, b_E, α_E, σ_ln,M,0, β_E) = (0.0, 1.0, 0.0, 0.2, 0.125).

4.2. Cosmic microwave background

In Sect. 7 we forecast the sensitivity to neutrino parameters of the Euclid probes combined with CMB data from the Planck satellite. In the context of a forecast, it is easier to describe the Planck constraining power not through the actual data and likelihood, but through mock data and a synthetic likelihood mimicking the sensitivity of Planck. This allows us to assume the same underlying cosmology for our Euclid and Planck mock data, and perform our forecast in ideal conditions – of course, one should keep in mind that this is a very different exercise from a real data analysis, which can bring surprises, especially if different data sets are in tension due to statistical flukes, incomplete systematic modelling, or physical effects that have not been accounted for.

As in many previous works, we use for this purpose the fake_planck_realistic likelihood of the public MontePython package, which accounts for the measurement of CMB temperature, polarisation, and lensing by Planck, with a sky coverage of 57% and noise spectra very close to the actual ones. This Gaussian likelihood accounts for three correlated observables: the CMB temperature map; the E-mode polarisation map, and the reconstructed CMB lensing potential map. In principle, these observables are slightly correlated with Euclid observables because the same clusters can shear high-redshift Euclid galaxies and distort patterns on the last scattering surface. This correlation will be taken into account in the analysis of real Euclid data, but, for simplicity, we neglect it in the present forecasts.

Next, we estimate the potential of Euclid data in combination with future CMB data from the LiteBIRD satellite, optimised for large angular scales, and from the CMB Stage-IV survey (CMB-S4), optimised for smaller angular scales. As in Brinckmann et al. (2019) we model this combination with two mock likelihoods in the public MontePython package: the litebird_lowl likelihood accounts for LiteBIRD temperature and polarisation data in the multipole range where the survey is the most constraining, 2 ≤ ℓ ≤ 50; and the cmb_s4_highl likelihood for temperature, polarisation, and lensing data from CMB-S4 at ℓ > 50. Details on the assumed sensitivity can be found directly in the numerical package or in Brinckmann et al. (2019). In that case, we follow the same methodology as for Planck: the CMB mock data account for the same fiducial model as Euclid and we neglect correlations between the CMB and large-scale structure data, although this correlation will be more important for a survey very sensitive to CMB lensing like CMB-S4.

For CMB experiments, the predicted error on N_eff depends mainly on the sensitivity of data to small angular scales. This means that in our modelling of LiteBIRD+CMB-S4 the assumptions regarding the CMB-S4 sensitivity at large ℓ are crucial. As in the rest of this forecast paper, we stick to conservative assumptions. Similarly to Brinckmann et al. (2019), we assume that foregrounds can be removed from CMB temperature and polarisation data only up to ℓ_max = 3000, that the Galactic cut imposes a sky coverage f_sky = 0.40, and that the CMB-S4 instrument can be modelled with a beam width of θ = 3.0 arcmin and a sensitivity to temperature of σ_T = 1.0 μK arcmin. With such assumptions, the combination of LiteBIRD+CMB-S4 predicts σ(N_eff) = 0.038 in the ΛCDM+∑m_ν+N_eff case – see Table 5 of Brinckmann et al. (2019). Since in this case the posterior is nearly Gaussian, this implies a 95% CL upper bound ΔN_eff < 0.076. This is consistent with the forecasts of the CMB-S4 white paper Abazajian et al. (2016), but that paper explores many other assumptions on CMB-S4 data, with a maximum multipole for the temperature data in the range ℓ_max^TT ∈ [3000, 5000], a sky fraction f_sky ∈ [0.1, 0.8], a beam width of θ ∈ [1.0, 3.0] arcmin, and a sensitivity of σ_T ∈ [1.0, 3.0] μK arcmin – see Figs. 22–24 of Abazajian et al. (2016). We should stress that in that regard the summary plot of Fig. 23 (right panel) – which can be compared with our Fig. 18 – CMB-S4 forecasts are performed under extremely optimistic assumptions: ℓ_max^TT = 5000, f_sky = 0.5 (respectively 0.7), θ = 1.0 arcmin, and σ_T = 1.0 μK arcmin, which gives a 95% CL upper bound ΔN_eff < 0.054 (respectively 0.046). While these numbers provide an example of what CMB Stage-IV may ideally achieve, we want to give more conservative estimates of what will be possible with the combined surveys.

5. Forecast method

Having described our modelling of observables and likelihood functions in Sects. 3 and 4, we now explain how such likelihoods can be used for the purpose of Euclid forecasts. Our forecast methods are based on Euclid Collaboration (2024a), adjusted to account for the presence of massive neutrinos. The validation of the four different pipelines described in this section against each other is detailed in Sect. 6 and one of them (i.e. the MCMC pipeline) was used in deriving our main results, presented in Sect. 7.

If the parameters are well constrained, we can approximate the posteriors with a multivariate Gaussian. Since such a distribution can be described by a mean vector and a covariance matrix, with such an approximation all the information concerning the sensitivity of the experiment to the parameters of the model is contained in the Hessian matrix F_αβ of the log-likelihood at the best fit, also called the Fisher Information (FI) matrix:

$$\begin{array}{c}\hfill F_{\alpha \beta }\equiv -⟨{\partial }_{\alpha }{\partial }_{\beta }{ln\mathcal{L}|}_{\mathrm{b}}estfit⟩,\end{array}$$(34)

where the indices {α, β} run over model (cosmological or nuisance) parameters. Our version of the MontePython package implements mock likelihood functions ℒ(θ_α) for each of the Euclid probes described in Sect. 3. The pipeline referred to later as MP/Fisher (the MontePython package used in Fisher mode) directly evaluates the FI matrix from Eq. (34). First, the code calls the Einstein–Boltzmann solver (EBS) CLASS to evaluate the cosmological observables required by the likelihood calculation (power spectra, growth rate, etc.). Then, it performs a numerical calculation of two-sided second-order derivatives, based on the evaluation of the likelihood at the best fit and in 2N² neighbouring points, where N is the number of model parameters¹².

Due to the computational cost and complexity of evaluating second-order derivatives, one can use the fact that the likelihoods have a relatively simple Gaussian dependence on the power spectra (but not on the parameters) to express the Fisher matrix in terms of first derivatives only. For the photometric probe, the FI matrix can be expressed as

$$\begin{array}{c}\hfill F_{\alpha \beta }=\frac{1}{2}{f}_{\mathrm{sky}}{\displaystyle \sum _{\ell }}(2\ell +1)\mathrm{Tr}\{{\left[{\mathsf{C}}^{\mathrm{f}}id(\ell )\right]}^{-1}\left[{{\partial }_{\alpha }{\mathsf{C}}^{\mathrm{t}}h(\ell )|}_{\mathrm{fid}}\right]\\ \hfill {\left[{\mathsf{C}}^{\mathrm{f}}\mathsf{id}(\ell )\right]}^{-1}\left[{{\partial }_{\beta }{\mathsf{C}}^{\mathrm{t}}h(\ell )|}_{\mathrm{fid}}\right]\},\end{array}$$(35)

where for each multipole ℓ the matrices C^fid and C^th are built out of the angular power spectra $C_{\mathit{ij}}^{\mathit{XY}}$ defined in Eq. (4), evaluated, respectively, at fiducial parameter values in the case of C^fid, or at two points (per dimension) in parameter space to compute numerically the first-order derivative ${\partial }_{\alpha }{C}^{\mathrm{th}}{|}_{\mathrm{fid}}$. Since with this method only the first-order derivatives are evaluated, the total number of points needed reduces to 2N (using a first-order double-sided finite difference). We refer the reader to Euclid Collaboration (2024a) for more details on the definition of the matrices C^fid and C^th. For the spectroscopic probe, the FI matrix becomes instead

$$\begin{array}{cc}\hfill F_{\alpha \beta }& ={\displaystyle \sum _i}\frac{1}{8{\pi }^2}{\displaystyle {\int }_{-1}^1}\mathrm{d}\mu {\displaystyle {\int }_{{10}^{-3}}^{k_{\max }}}\mathrm{d}k{k}^2[{\partial }_{\alpha }\ln (P_{\mathrm{obs}}){|}_{\mathrm{fid}}\hfill \\ \hfill & \times {\partial }_{\beta }\ln (P_{\mathrm{obs}}){|}_{\mathrm{fid}}{V}_i^{\mathrm{fid}}],\hfill \end{array}$$(36)

where P_obs is the observable power spectrum defined in Eq. (11), k_max is given in Table 2 and the index i runs over four redshift bins with effective volume V_i^fid. Details on the assumed survey characteristics and binning strategy are specified in Euclid Collaboration (2024a).

The pipeline that we refer to as CF, that is, the CosmicFish code, uses Eqs. (35) and (36) to compute the FI matrix. It calculates the matrices C(ℓ) and the spectra P_obs at the fiducial point and at a given number of neighbouring points in parameter space, and numerically evaluates the first derivatives of these observables with respect to model parameters. By default, CosmicFish uses two-sided first-order derivatives (two neighbouring points). To check the numerical stability of this method, we have checked that we get similar result when using a 4-point forward stencil derivative (see Appendix A). The version of CosmicFish used in this work internally calls one of the EBSs (CAMB or CLASS ) to evaluate the cosmological observables, providing us with two forecast pipelines CF/CAMB and CF/CLASS¹³.

Our fourth pipeline called MP/MCMC is the MontePython package used in its default mode, which is running Markov chain Monte Carlo (MCMCs) with a Metropolis-Hasting algorithm to explore the parameter space and infer the posterior from the likelihood according to Bayes theorem. We always assume top-hat priors on model parameters, with prior edges chosen far in the exponential tails of the posteriors, such that the confidence intervals inferred from the MCMCs at the 95% CL are unaffected by the choice of prior edge. The only exception is the summed neutrino mass ∑m_ν, for which we impose a theoretical prior ∑m_ν > 0. The two main advantages of MCMC forecasts is that they rely neither on a Gaussian approximation to the likelihood, nor on the calculation of numerical derivatives. We note that the MP/Fisher and MP/MCMC pipelines call the very same Euclid mock likelihoods, computing the same functions ℒ(θ_α), which are the only parts of the pipeline depending on the cosmology; thus, it is impossible by construction to introduce a mistake in the modelling of cosmology (or of the Euclid survey) when passing from one pipeline to the other. This means that the validation of one pipeline automatically implies the validation of the other pipeline (as long as the FI matrix is well conditioned and the MCMC has converged).

This statement could be questioned in the case in which one of the two approaches (FI matrix or MCMC) requires a more accurate calculation of the cosmological observables or of the likelihood functions. In the case at hand, we know that the highest precision is required by the FI matrix approach, because a robust evaluation of numerical derivatives with small step sizes requires high accuracy and low numerical noise. However, our approach consists of validating the MP/Fisher pipeline first, which is more demanding in terms of accuracy; this in turn validates the MP/MCMC pipeline, which is less demanding. We also explicitly compare the results from the two pipelines in Sect. 6 to make sure that they agree.

6. Code validation

In order to validate the four forecast pipelines described in Sect. 5 against each other, we need to choose a framework (fiducial model and set of parameters) in which the FI formalism can be applied in a robust way. In the results section, we will be interested in extensions of the minimal five-dimensional (5D) flat ΛCDM model with up to four additional free parameters {∑m_ν, N_eff, w₀, w_a}¹⁴. However, floating these nine parameters simultaneously opens up strong degeneracies in the parameter space and leads to very non-Gaussian posteriors. In that case, the Gaussian approximation breaks down. Then, the FI matrix calculation becomes unstable, since it depends heavily on the choice of numerical derivative step sizes. Besides, the FI matrix no longer provides a good estimate of the experimental sensitivity.

To overcome this issue, we decided to validate our forecasting pipelines by looking separately at three 7-parameter models for which the posterior remains close to Gaussian (as confirmed by contour plots from MCMC runs displayed in Sect. 6.2). This strategy does not limit the range of validity of our final results, since at the end, after validating the pipelines against each other, we will use the MCMC pipeline for our main results. This pipeline does not rely on any Gaussian approximation and can be used even for full 9-parameter forecasts.

The three seven-parameter models used for validation are defined as follows. Firstly, we stick to a w₀w_aCDM cosmology, with a time-varying dark energy equation of state parameter w(a) = w_a(1 − a)+w₀ (Chevallier & Polarski 2001; Linder 2003) but fixed {∑m_ν, N_eff}. This model was already used in the forecasts of Euclid Collaboration (2024a) and allows us to cross-check that the previous validation still holds. Secondly, we investigate cosmologies where we vary the neutrino mass and the effective number of relativistic degrees of freedom, ΛCDM+∑m_ν+N_eff, but stick to a cosmological constant. This framework allows us to cross-check our implementation on neutrino effects on the observables. Thirdly, we consider models with varying neutrino mass and dark energy equation of state, w₀CDM+∑m_ν, with fixed N_eff = 3.044 and w_a = 0. This last set up will prove that our modelling of nonlinear corrections is consistent across the different pipelines even when neutrino and DE parameters are varied simultaneously.

The fiducial values of cosmological parameters used in Sects. 6 and 7 are summarised in Table 4. Despite identical fiducial values, the models used for validation in this section and for forecasts in the next section feature two subtle differences with respect to each other. First, as already explained in Sect. 2.2 (case b), we assume here a single massive neutrino species with mass m_ν = ∑m_ν, which allows for easier comparison with earlier work and for simpler FI matrix calculations (given that the N_eff posterior remains symmetric around the fiducial value). Second, we fix here the nonlinear modelling parameters σ_p and σ_v to their fiducial values in order to keep the contours of the spectroscopic probe more Gaussian. In Sect. 7, we instead consider three degenerate massive neutrino species, each with mass m_ν = ∑m_ν/3 (see Sect. 2.2, case a) and we will marginalise over σ_p and σ_v.

We performed validation tests for each of the three seven-parameter models, each probe (photometric or spectroscopic), and each case (pessimistic or optimistic). Thus, the number of validation tests amounted to 12. In each of these 12 situations, we proceeded as follows:

1. Validate the predictions and accuracy settings of the EBSs by comparing the CF/CAMB and CF/CLASS pipeline results;

2. Validate our MontePython likelihood and our FI matrix implementation by comparing the results from the CF/CLASS and MP/Fisher pipelines;

3. Check the validity and stability of our FI forecasts (choice of step size in numerical derivatives and validity of Gaussian approximation) by comparing the previous two pipelines with MP/MCMC runs.

In practice, we performed steps 1 and 2 simultaneously, by comparing all three derived FI matrices. For steps 1 and 2, our test consisted of comparing the marginalised and unmarginalised errors on each cosmological and nuisance parameter¹⁵. Validation was achieved when all errors were within 10% of the median. This subjective criterion¹⁶ was already adopted in previous Euclid validation papers like EC20 and Euclid Collaboration (2024a). The results of these tests are shown in Sect. 6.1. For step 3 we compared 1D posteriors and 2D confidence contours in triangle plots in Sect. 6.2.

6.1. Validation of Fisher pipelines

We show the marginalised and unmarginalised errors predicted by each pipeline (CF/CAMB, CF/CLASS, MP/Fisher) used to compute the FI matrix, each probe, and each case, respectively, in Fig. 6 for the w₀w_aCDM model, in Fig. 7 for the ΛCDM+∑m_ν+N_eff model, and in Fig. 8 for the w₀CDM+∑m_ν model. In all cases, the errors are within 10% of the median, leading to the validation of the pipelines.

Fig. 6. For the w0waCDM model, comparison of each of the Fisher marginalised (light grey) and unmarginalised (dark grey) errors on each cosmological and nuisance parameter for: photometric optimistic probe (top left); photometric pessimistic probe (top right); spectroscopic optimistic probe (bottom left); and spectroscopic pessimistic probe (bottom right). For the spectroscopic probe the labels lbsi are short for the parameters ln(b σ8)i.

Fig. 7. Same as Fig. 6, but for the ΛCDM+∑mν+Neff model.

Fig. 8. Same as Fig. 6, but for the w0CDM+∑mν model.

Uncertainties from CF/CAMB and CF/CLASS are usually within 2% of each other, which shows the excellent agreement between the two EBSs. Such a level of agreement was achieved after a careful setting of the parameters describing underlying assumptions on the cosmological models (especially those controlling the modelling of the neutrino sector) and of accuracy settings. Our choice of input and accuracy settings for CAMB and CLASS is described in Appendix B. The errors from CF/CAMB and CF/CLASS differ by more than 2% only in the photometric case, in which the σ₈ and ∑m_ν errors are 3–5% away from each other. We believe this is due to the tiny but irreducible difference in the neutrino treatment between CLASS and CAMB. As shown in Fig. B.1, the suppression of the power spectrum on small scales due to neutrinos seems to start earlier in CLASS than CAMB, inducing differences of the order of 0.1% in the linear matter power spectrum. Still, the final error bars achieve a level of agreement that is remarkable and sufficient for Euclid purposes.

Errors from CF/CLASS and MP/Fisher are within 10% of each other. The precise level of agreement is sensitive to the choice of step sizes, especially in the MP/Fisher pipeline. For this choice, we adopted the following strategy

• For CF/CAMB and CF/CLASS, we adopt by default relative steps of 1% in each parameter (or 0.01 for w_a). However, we observed that this choice is inappropriate for ∑m_ν, since this parameter has a smaller impact on the matter power spectrum. The effect of a 1% variation in the neutrino mass on both P_mm and P_cc is so small that numerical derivatives would be dominated by numerical noise from the EBSs. Thus, for ∑m_ν, we adopt a step size of 10%, that is, Δ(∑m_ν) = 6 meV. We checked that the results are relatively stable against small variations in the step sizes (see Appendix C) and against a change in the numerical derivative method (see Appendix A).

• In MP/Fisher, we stick to the recommendation of Euclid Collaboration (2024a) and choose the step size in relation to the marginalised 1D parameter uncertainty for each individual parameter for a given model. To get the uncertainties, we could either run an MCMC or use the error obtained from CosmicFish. In this case we use the uncertainties obtained from CosmicFish. In the case of the photometric probe, we choose a step size of 10% of the posterior error, while for the spectroscopic probe, we used a step size of 5% of the error.

These results are very sensitive to the modelling of nonlinear corrections. In Sect. 3.3, we argue that the most recent version of HMcode (Mead et al. 2021) reproduces simulations featuring massive neutrinos significantly better than HALOFIT (Takahashi et al. 2012; Bird et al. 2012). We checked that using HALOFIT instead of HMcode does affect the marginalised forecasted errors by a large factor of 1.5 (see Appendix D). The assumption that galaxies trace cold-plus-baryonic matter instead of total matter is also important for Euclid, but to a lesser extent. When replacing P_cc by P_mm everywhere in the spectroscopic likelihood, the uncertainties vary by at most 5% on ∑m_ν (for our fiducial value ∑m_ν = 60 meV) (as shown in Appendix E).

6.2. Comparison of Fisher ellipses and MCMC contours

For each of the three validation models and each of the two probes, we show the 1D marginalised posteriors and the 2D marginalised confidence contours of our MP/MCMC runs in the triangle plots of Figs. 9 to 11. For this test, we only consider the optimistic case, which is the most sensitive to details in the modelling of the observables, and the most likely to reveal possible discrepancies in the numerical implementation of the likelihoods. This case is also the most constraining, and thus closer to the Gaussian approximation and the most suited for comparison with FI matrix predictions. In Figs. 9 to 11, we actually plot the Gaussian posteriors and elliptical contours of the CF/CAMB and MP/Fish pipelines on top of the MP/MCMC results.

Fig. 9. Comparison of the 1D marginalised posteriors and 2D (95% C.L. and 68% C.L.) marginalised contours obtained with the MP/MCMC pipeline to the Gaussian posteriors and elliptical contours from the FI pipelines CF/CAMB and MP/Fisher for the w0waCDM model. Left: Spectroscopic optimistic probe. Right: Photometric optimistic probe.

First and foremost, we can see in Fig. 9 that the previous validation of photometric probes presented in Euclid Collaboration (2024a) for the w₀w_aCDM model is reproducible, despite the changes in the likelihood code and nonlinear recipe triggered by the assumption of massive neutrinos. Concerning the spectroscopic probes, here we have corrected for an inconsistency in the unit conversion from Mpc⁻¹ to h Mpc⁻¹ that affected previous forecasts; thus, the sensitivity to h degrades with respect to Euclid Collaboration (2024a), while the sensitivity to the other cosmological parameters is only mildly altered due to their correlations with h. Nevertheless, the conclusions drawn in Euclid Collaboration (2024a) about the validation of MP/MCMC for this case still hold: the posteriors of the photometric and spectroscopic probes are very close to multivariate Gaussians. We get another confirmation of this, since the MCMC and FI posteriors or contours are hardly distinguishable by eye and agree up to typical MCMC convergence errors.

In the case of the ΛCDM+∑m_ν+N_eff model, we see in Fig. 10 that non-Gaussian effects come into play. This is particularly obvious when looking at the marginalised neutrino posteriors, which are truncated at ∑m_ν = 0 by the physical prior, and whose tails fall quicker than a Gaussian distribution for large masses. This non-Gaussianity propagates to the parameters that are most correlated with ∑m_ν. Indeed, the orientation of the 2D contours shows that h, Ω_m, 0, n_s, and σ₈ are significantly correlated with ∑m_ν, while Ω_b, 0 and N_eff are not. The former parameters are also the ones whose posteriors deviate the most from a Gaussian shape. The fact that the MCMC contours are overall tighter than the FI ellipses is also a consequence of the non-Gaussianity with respect to ∑m_ν > 0, which cuts off a region of the parameter space (with either negative or too large neutrino mass) that would otherwise be reachable by FI contours. Deviations from Gaussianity are found to be even stronger within the photometric probe, both because the neutrino mass posterior is more suppressed for large masses and due to the very strong correlation between ∑m_ν, Ω_m, 0, and σ₈.

Fig. 10. Same as Fig. 9, but for the ΛCDM+∑mν+Neff model.

Finally, in the case of the w₀CDM+∑m_ν model, we see some qualitatively similar trends in Fig. 11. For the spectroscopic probe, we observe a ‘banana’-shaped contour in the n_s-∑m_ν plane and slightly non-elliptic contours in other planes. Once again, the non-Gaussianity is even more pronounced with the photometric probe, due to the cut-off of negative or too large neutrino masses.

Fig. 11. Same as Fig. 9, but for the w0CDM model+∑mν model.

As in Sect. 5, we stress that the MP/Fisher and MP/MCMC pipelines call the very same Euclid mock likelihoods functions ℒ(θ_α). Given that these functions have been validated by the tests presented in the previous subsection, and that the MCMC runs presented here have reached a high convergence level, the difference between the MCMC and FI contours can only be attributed to the intrinsic non-Gaussianity of the problem. This shows that for extended cosmological models, particularly in the case of a free neutrino mass parameter, robust Euclid forecasts (and even more so with respect to analyses of future real data) require a full MCMC analysis. Having validated the MontePython likelihoods against CosmicFish forecasts, we have proved that our MP/MCMC pipeline offers a robust and accurate implementation of the physical assumptions performed in this work. Thus, we can use it in the next section in order to derive our main forecast results.

7. Results

Here, we present the results of our forecast analysis. We use the previously validated MP/MCMC pipeline based on the MontePython sampler, interfaced with the Boltzmann solver CLASS. The underlying cosmology is either flat ΛCDM with the standard cosmological constant or w₀w_aCDM, where dark energy is a fluid with a time-varying equation-of-state parameter according to Chevallier & Polarski (2001) and Linder (2003).

As already explained in Sect. 2.2 (case a), from the many possible ways of distributing the total mass ∑m_ν among the different species contributing to the effective number of neutrinos N_eff (which could lead to potentially different results), we made a choice that is representative of a large class of well-motivated scenarios. We assume that the three active neutrino species are degenerate in mass with a fiducial value ∑m_ν = 60 meV close to the minimum allowed by normal ordering. We further assume that massive neutrinos contribute to N_eff through the standard value 3.044 and that additional free-streaming massless particles account for ΔN_eff = N_eff − 3.044, with ΔN_eff > 0 (see Appendix B.1 for details on the specific settings in CLASS and CAMB). Sticking to these assumptions, we consider four extended cosmological scenario with an increasing number of free parameters:

• ΛCDM + ∑m_ν (baseline);

• ΛCDM + ∑m_ν + ΔN_eff;

• w₀w_aCDM + ∑m_ν;

• w₀w_aCDM + ∑m_ν + ΔN_eff.

In Table 4 we list the fiducial values for the parameters of our models.

Here (in contrast to the validation of the MCMC pipeline) we use pessimistic specifications, additionally also varying the nonlinear nuisance parameters σ_v and σ_p of the GC_sp recipe. The reasoning behind this choice is that additional theoretical errors and systematic effects might have been overlooked in the modelling of the observables adopted in EC20 and Euclid Collaboration (2024a). Therefore, applying optimistic specifications might overestimate the sensitivity to the neutrino mass, while our aim is to provide a conservative forecast of the Euclid discovery potential. Moreover, the cut-off at ℓ_max = 1500, implied by pessimistic specifications, prevents our results from being biased due to the uncertainties on the modelling of baryonic feedback, as we discuss in Appendix H.

We fit both Euclid and CMB synthetic data (for the latter, see Sect. 4.2). For Euclid, the observables are the 3 × 2pt angular power spectra from photometric data (photometric galaxy clustering, weak lensing, and their cross-correlation, see Sect. 3.1), the 3D galaxy power spectrum from spectroscopic data (Sect. 3.2), and cluster number counts (CC, Sect. 4.1).

The cosmological parameter fiducial values were already specified in Table 4. The parameters w₀ and w_a are kept fixed to the fiducial value when the underlying cosmological model is ΛCDM. The parameters of the baseline model are always varied, while the parameters of the extended models (such as the number of relativistic degrees of freedom and the dark energy equation of state parameters) are kept fixed, unless specified. The optical depth to reionisation is fixed to the Planck best fit value τ = 0.0543 when only Euclid observables are included in the analysis, while we allow it to vary in the joint analysis with CMB probes.

In Table 5, we show the 1 σ uncertainty (or the 95% CL upper limit) on each cosmological parameter inferred from our MCMC forecasts, assuming different cosmological models and varying the combination of synthetic data¹⁷. In Fig. 12, we show the same results expressed in terms of a dimensionless sensitivity, which is the 1 σ uncertainty relative to the fiducial value¹⁸. In each panel we compare the two underlying cosmologies, which are ΛCDM + ∑m_ν (solid colour bars) and w₀w_aCDM + ∑m_ν (hatched colour bars). The effective number of neutrino species is additionally varied in the bottom panel. The different colours show the sensitivity of three different data combinations: WL+GC_ph+XC_ph+GC_sp (blue); Euclid +Planck (orange); and Euclid +CMB-S4+LiteBIRD (green).

Fig. 12. Sensitivity (relative errors with respect to the fiducial values) to the cosmological parameters, comparing ΛCDM cosmologies (solid colour bars) and w0waCDM cosmologies (hatched colour bars), for models with (bottom panel) and without (top panel) varying Neff. The different colours correspond to the sensitivity of different data combinations: WL+GCph+XCph+GCsp (blue); Euclid +Planck (orange); and Euclid +CMB-S4+LiteBIRD (green). We note that for wa we show the absolute errors, while for ΔNeff we show the 95% upper limits. For the neutrino mass we always show the relative 1 σ uncertainty as computed from the posterior variance in order to facilitate the comparison between different data sets.

Comparison with previous forecast. When comparing the results of the present forecast to those of EC20, one has to consider that we assume the same specifications labelled as pessimistic there. However, our baseline case includes the variation in the neutrino mass, which was not done in EC20. In light of these considerations, we observe that the sensitivity of the photometric probes (WL+GC_ph+XC_ph) to w₀ and w_a is consistent with the previous Euclid forecast from Table 11 of EC20 (note that they cite relative sensitivities, while we cite absolute sensitivities). On the other hand, the sensitivity to the other cosmological parameters degrades because of the correlation with the neutrino parameters; for example, the sensitivity to σ₈ degrades by a factor of 1.6 in w₀w_aCDM+∑m_ν+ΔN_eff due to the free neutrino mass, while varying ΔN_eff does not have any impact. A direct comparison of the results that include GC_sp with those of EC20 is not feasible because, as already mentioned in Sect. 6.2, here we have corrected an inconsistency in the unit conversion from h Mpc⁻¹ to Mpc⁻¹ that was present in the pipeline of EC20. This issue produced a fake signal on H₀, and, thus, a spurious enhancement of the sensitivity to h, which affected also the sensitivity to other cosmological parameters that are correlated with h. Note that we checked that reintroducing the bug would make our results consistent with those of EC20. Our results are also consistent with the older Euclid -like sensitivity forecast of Audren et al. (2013), which predicted a 1 σ sensitivity to ∑m_ν of 32 meV for Euclid WL  +  Planck or 25 meV for Euclid GC_sp + Planck, while we obtain 23 meV for the combined WL + GC_ph + XC_ph + GC_ph + Planck data set. Interestingly, to account for the imperfect modelling of nonlinear scales, that forecast included a marginalisation over a theoretical error smoothly growing towards smaller scales, instead of cutting the data abruptly at some k_max or ℓ_max like in the present work. Audren et al. (2013) adopted a theoretical error amplitude suggested by the covariance of the nonlinear power spectra predicted by different N-body codes for the same cosmology. The fact that our results are comparable suggests that their theoretical error and our pessimistic settings are equally conservative. A slightly different modelling of the theoretical error in Sprenger et al. (2019) led to a sensitivity to ∑m_ν of 24 meV for Euclid GC_sp + WL + Planck with pessimistic specifications (there labelled as conservative), consistent with our results. Other works in the literature present Euclid -like forecasts (Yankelevich & Porciani 2019; Chudaykin & Ivanov 2019; Sailer et al. 2022) including only galaxy clustering with a one-loop modelling. However, the differences in assumptions, methodologies, and survey specifics make direct comparisons challenging.

Complementarity of Euclid probes. In Fig. 13, we show the constraints on the extended w₀w_aCDM+∑m_ν+ΔN_eff model expected from Euclid-only probes. The same figure for the baseline model ΛCDM + ∑m_ν can be found in Appendix F. We note that in the extended w₀w_aCDM+∑m_ν+ΔN_eff model the constraints on h and ΔN_eff are dominated by the spectroscopic probes and those on σ₈, ∑m_ν, w₀, and w_a by the photometric probes. The complementarity of the photometric and spectroscopic probes is impressive: the GC_sp contours exhibit parameter degeneracies in different directions than the WL+GC_ph+XC_ph ones in several planes of the parameter space. In particular, the correlations between the dimensionless Hubble constant h and some of the other cosmological parameters that are present in GC_sp are broken by the photometric information. This leads to a remarkable sensitivity increase in the joint analysis of the two primary Euclid probes.

Fig. 13. Marginalised 1D distribution (diagonal) and 2D (off-diagonal) 68% and 95% contours for a selection of parameters of the w0waCDM+∑mν+Neff model obtained by fitting GCsp (grey), WL+GCph+XCph (red), and their combination (blue). The horizontal and vertical dashed lines mark the fiducial values.

As far as the total neutrino mass is concerned, photometric probes are potentially more powerful than spectroscopic ones in constraining ∑m_ν. As long as we are not including CMB information, the observable effect of the total neutrino mass is best described by the top panel of Fig. 1, that is, a steplike suppression of the matter and cold-plus-baryonic power spectra. The longer lever arm of WL data thus provides the best constraint. However, the Euclid probes alone cannot provide significant evidence for a non-zero neutrino mass, neither in the w₀w_aCDM cosmology (∑m_ν ≲ 0.22 eV, 95% CL, see Table 5) nor in ΛCDM [σ(∑m_ν) = 56 meV, see Table 5].

When using CMB data, there is a well-known correlation between ∑m_ν and h, due to their opposite effects on the angular diameter distance to last scattering. This strong correlation does not show up when fitting Euclid probes alone, mainly thanks to the sensitivity of galaxy clustering data to the unique scale dependence of the growth factor induced by massive neutrinos. Note that this conclusion does not depend on the specific underlying cosmological model, whether it is ΛCDM or w₀w_aCDM, and whether ΔN_eff is varying or fixed. We anticipate that the different directions of degeneracy between ∑m_ν and h in Euclid and CMB probes will induce a significant improvement on the ∑m_ν constraints in the Euclid +CMB joint analysis. Finally, note that we do not find any specific correlation between the neutrino mass and the bias parameters, neither in GC_ph nor in GC_sp, as we show in Appendix G.

Concerning dark radiation, as already mentioned, a variation in ΔN_eff impacts the matter and cold-plus-baryonic matter power spectra both through a change in the overall shape (caused by a shift of the time of equality) and in the phase and amplitude of BAO peaks (due to neutrino drag), as can be checked in the top panel of Fig. 2. The first effect induces a correlation between ΔN_eff and h, since both parameters can be increased while maintaining a fixed redshift of equality. This correlation, combined with the prior volume effect induced by the physical prior ΔN_eff > 0 (Hadzhiyska et al. 2023), leads to a small bias in the mean value of h and of n_s reconstructed from photometric probes. We have explicitly checked that this bias disappears in the models where ΔN_eff is fixed. The Euclid constraints on ΔN_eff are dominated by the spectroscopic probe, because GC_sp data partially break the correlation between ΔN_eff and h thanks to a better measurement of BAO peaks. Nevertheless, the photometric probes contribute to the ΔN_eff constraints by breaking the degeneracy in the ∑m_ν versus ΔN_eff plane. Normally, one would expect a negative correlation between these parameters, due to their opposite effects on the BAO peak scale¹⁹. This negative correlation is clearly visible on the GC_sp-alone contour of Fig. 13. Thanks to the large lever arm of WL data, the photometric probe allows us to distinguish the difference in overall shape and growth rate induced by variations of ∑m_ν or ΔN_eff in the matter and cold-plus-baryonic power spectra, and to lift the degeneracy. As a consequence, a combination of the probes improves the sensitivity to both ΔN_eff and ∑m_ν.

On the other hand, under our assumptions concerning the modelling of the Euclid cluster count data and likelihood, we find no significant improvement on cosmological parameter sensitivity when adding clusters. Table 5 shows the result for WL+GC_ph+XC_ph+GC_sp+CC, where CC stands for cluster count, in the ΛCDM+∑m_ν and w₀w_aCDM+∑m_ν+N_eff cases. The difference between the results with and without CC lies within MCMC convergence errors. Additionally, in the w₀w_aCDM+∑m_ν+N_eff case, the upper bound on ΔN_eff gets a tiny enhancement from a Bayesian prior volume effect related to small parameter degeneracies between ΔN_eff and the nuisance parameters of the CC likelihood. We recall that we performed a very conservative modelling of CC data with wide priors on several nuisance parameters. With a more optimistic modelling, the same data may lead to tighter constraints.

Combination with CMB. In Figs. 14 and 15, we show the constraints placed on the parameters of the w₀w_aCDM+∑m_ν+ΔN_eff and ΛCDM+∑m_ν models for the combination of Euclid probes with CMB experiments, Euclid +Planck (orange contours), and Euclid +CMB-S4+LiteBIRD (green contours). The same figure for the other cosmological models can be found in Appendix F. These figures show that the Euclid probes will play a crucial role in improving constraints on the cosmological parameters considered in this work. As a matter of fact, we can identify several cases in which future CMB data sets allow for parameter degeneracies that are broken by Euclid probes.

Fig. 14. Marginalised 1D distributions (diagonal) and 2D (off-diagonal) 68% and 95% contours for a selection of parameters of the w0waCDM+∑mν+Neff model obtained by fitting Euclid +Planck (orange), and Euclid +CMB-S4+LiteBIRD (green). The horizontal and vertical dashed lines mark the fiducial values.

This is particularly true for the degeneracy between ∑m_ν and h. We already discussed above why these parameters are correlated in a CMB-only analysis, and why the Euclid probes can break this degeneracy, as can be checked again by comparing the {∑m_ν, h} panel of Fig. 13 with the same panel in Fig. 14 (see also Fig. 15). Additionally, when CMB data are taken into account, the observable impact of the total neutrino mass is best described by the bottom panel of Fig. 1, and amounts mainly to an overall (redshift-dependent) suppression of the matter and cold-plus-baryonic power spectra. Then, the role of CMB data is to fix the amplitude of the primordial spectrum of fluctuations at high redshift (modulo some uncertainty on the optical depth to reionisation that depends on CMB polarisation measurements on large angular scales), while the Euclid WL probe fixes the amplitude of the matter and cold-plus-baryonic power spectra at low redshift. Therefore, the combination of Euclid and CMB complementary data yields much stronger bounds on ∑m_ν than either of the two data sets taken separately, allowing eventually for a neutrino mass detection. Even in the minimal hierarchy scenario, which we assume as our fiducial model, the joint Euclid plus Planck forecast returns a 2.6 σ evidence for a non-zero mass with σ(∑m_ν = 60 meV) = 23 meV (Euclid +Planck). Finally, replacing Planck with a combination of future CMB surveys leads to the detection of a non-zero neutrino mass at almost 4 σ with σ(∑m_ν = 60 meV) = 16 meV (Euclid +CMB-S4+LiteBIRD). One of the reasons behind this improvement is the tightening of the constraints on the reionisation optical depth from LiteBIRD. Indeed, Fig. 16 shows how the joint fit of LSS and CMB data leads to a strong degeneracy between ∑m_ν and τ (Allison et al. 2015; Archidiacono et al. 2017). Therefore, the sensitivity to the neutrino mass can be further improved by means of an independent measurement of the reionisation optical depth. Such measurement can arise from future 21-cm experiments (Sailer et al. 2022)²⁰, from probing the kinetic Sunyaev-Zeldovich effect in the CMB (Park et al. 2013), and from CMB polarization (Allys et al. 2023).

Fig. 15. Same as Fig. 14, but for the ΛCDM+∑mν model.

Fig. 16. Marginalised 2D contours of ∑mν and τreio for the ΛCDM+∑mν model.

Given the fantastic ability of Euclid to constrain the neutrino mass sum in combination with CMB data, it is worth investigating whether Euclid data will help to reconstruct the neutrino mass hierarchy. Figure 17 shows the neutrino mass sum required by the two neutrino mass ordering schemes as a function of the mass of the lightest neutrino. We see that, if the neutrino mass sum has indeed the minimum value allowed by neutrino oscillations in normal ordering (NO), which coincides with the fiducial value of our forecast (∑m_ν = 60 meV), then the Euclid +CMB-S4+LiteBIRD sensitivity to the neutrino mass sum will also imply a 2.5 σ indication in favour of NO. Of course, the statistical significance of the evidence in favour of NO will decrease if the best-fit neutrino mass turns out to be larger than the minimum value of NO. It is important to stress that this measurement is not a direct probe of the neutrino mass ordering, but rather a consequence of the strong sensitivity to ∑m_ν.

Fig. 17. Sum of neutrino masses ∑mν in the inverted ordering (red solid line) scenarios and in the normal ordering (blue solid line), including the uncertainties on the mass squared differences from oscillation experiments, as a function of the lightest neutrino mass. The grey dot-dashed lines show the 95% CL upper bounds from Planck only (dark grey), and in combination with external data (Planck +ext., light grey). For the latter we use the value quoted in Alam et al. (2021) for Planck +Pantheon-SNe+SDSS(BAO+RSD)+DES 3 × 2pt, which includes Planck, Pantheon type Ia supernovae, baryon acoustic oscillations and redshift-space distortions from the Sloan Digital Sky Survey, and weak lensing measurements from the Dark Energy Survey. The region above the blue dashed line is excluded at 95% CL by Euclid only primary probes. The orange shading shows the 68% and 95% CL limits on ∑mν from Euclid +Planck in the baseline ΛCDM+∑mν model, assuming a fiducial value ∑mν = 60 meV. The 68% and 95% CL constraints from Euclid +CMB-S4+LiteBIRD are shown in green. We note that we split the two constraints into two different ranges of the lightest neutrino mass only to avoid an overlap of the two shadings. Assuming that the true value of the neutrino mass sum is indeed the minimum allowed by neutrino oscillation experiments in normal ordering, the combination Euclid +CMB-S4+LiteBIRD will rule out the inverted ordering at 2.5 σ.

Contrary to the case of ∑m_ν, the constraints on ΔN_eff are dominated by CMB probes. Nevertheless the Euclid probes play an important role in breaking some of the degeneracies with other cosmological parameters, and in particular between ΔN_eff and h. Euclid improves the Planck sensitivity to the Hubble constant by a factor of 2.5 (σ(h) = 0.0021, Euclid +Planck, in ΛCDM+∑m_ν, compared with Table 2 of Planck Collaboration VI 2020); as already discussed and shown in Fig. 13, this is mainly due to GC_sp data. This significantly increases the sensitivity of combined CMB and Euclid data to ΔN_eff. Additionally, as shown in the bottom panel of Fig. 2, once we include CMB information (thus fixing the angular size of the sound horizon at recombination), an increase in N_eff leads to an enhancement of the matter and cold-plus-baryonic power spectra at small scales. This effect provides additional sensitivity to N_eff. Figure 18 shows that adding Euclid improves the current limit (Planck Collaboration VI 2020) on the presence of new light particles by more than a factor of 2 (ΔN_eff < 0.144, 95% CL, Euclid +Planck, in ΛCDM+∑m_ν+ΔN_eff). On the longer term, the Euclid probes even have the potential to improve the sensitivity of future CMB data sets like LiteBIRD+CMB-S4 to N_eff – once more, not because of a better sensitivity to the effect of N_eff itself, but thanks to a significantly more accurate determination of the correlated parameter h. Under the conservative assumptions described in Sect. 4.2, LiteBIRD+CMB-S4 alone would achieve ΔN_eff < 0.076 (95% CL, ΛCDM+∑m_ν+N_eff). The combination Euclid +LiteBIRD+CMB-S4 would increase the sensitivity to ΔN_eff < 0.063 (95% CL, ΛCDM+∑m_ν+N_eff)²¹. Such bounds approach the asymptotic limits for particles of spin 1 and 1/2 at temperatures above the electroweak phase transition, and thus can almost exclude their existence, independently of their decoupling temperature (see Fig. 18).

Fig. 18. Contribution of new light particles beyond the Standard Model to ΔNeff as a function of their decoupling temperature. Assuming natural units, we report the temperatures in GeV. As a reference we show the contribution of a Goldstone boson (solid brown line), a Weyl fermion (solid red line), and a massless Gauge boson (solid purple line). The grey vertical bands denote the range of temperatures of the QCD phase transition (around 200 MeV) and neutrino decoupling (around 1 MeV), while the thin grey vertical lines mark the annihilation temperature of Standard Model particles. The grey horizontal dot-dashed line represents the current 95% CL limit from Planck TT,TE,EE+lensing+BAO (Planck Collaboration VI 2020). The horizontal dashed lines mark the 95% upper bounds on ΔNeff from Euclid only (in blue) and in combination with current (Planck) and future (CMB-S4+LiteBIRD) CMB surveys (in orange and green, respectively). As explained in Sect. 4.2, we choose to adopt more conservative assumptions regarding the CMB-S4 sensitivity and foreground removal than in a similar figure in the CMB-S4 white paper, Fig. 23 of Abazajian et al. (2016).

The sensitivity of Euclid to the Hubble parameter will bring strong indications concerning the puzzle raised by the Hubble tension in the near future. We note that the sensitivity of the combination of Euclid spectroscopic and photometric probes is such that, in this case, future CMB observations will not lead to a significant improvement. For instance, in the ΛCDM+∑m_ν model, the sensitivity will marginally increase from σ(h) = 0.0021 with Euclid +Planck to σ(h) = 0.0016 with Euclid +CMB-S4+LiteBIRD. In the w₀w_aCDM+∑m_ν model, Euclid data are essential to break the correlation between the dark energy parameters and the other cosmological parameters, including h, therefore, there is no improvement at all in the sensitivity to h due to future CMB experiments (from σ(h) = 0.0036 with Euclid +Planck to σ(h) = 0.0035 with Euclid +CMB-S4+LiteBIRD).

Moreover, Euclid will provide some key information regarding the σ₈ tension, having a 1 σ error on S₈ = σ₈(Ω_m/0.3)^0.5 of 0.0049 from Euclid only in the ΛCDM+∑m_ν case. This is an improvement of nearly 3 (2.8) with respect to the DES Y3 3 × 2pt result (Amon et al. 2022; Secco et al. 2022).

Finally, the sensitivity of Euclid to dynamical dark energy parameters, which was already illustrated in previous forecasts assuming a fixed neutrino mass and no extra relativistic species (EC20,Casas et al. 2023), is confirmed in the framework of extended cosmologies with varying ∑m_ν and ΔN_eff. As a matter of fact, the contours derived from Euclid probes do not show significant correlations between the dark energy and neutrino parameters (see Fig. 13). On the other hand, the combination of Euclid with CMB brings back the correlation between ∑m_ν and {w₀, w_a}, as we shall comment in the next paragraph on the neutrino mass sensitivity in extended cosmologies.

Neutrino mass sensitivity in extended cosmologies. The constraints on the total neutrino mass from the Euclid primary probes alone are nearly independent of the assumed underlying model, as can be checked in Table 5 and in Fig. 12. However, once combined with external data sets, the constraints become more model dependent.

Floating the number of relativistic relics ΔN_eff has a limited impact on ∑m_ν bounds, but allowing for time-dependent dark energy degrades the mass sensitivity by a significant amount. For instance, the Euclid +Planck sensitivity decreases from σ(∑m_ν = 60 meV) = 23 meV in the ΛCDM case to 38 meV in the w₀w_aCDM case. Similar trends are observed with the Euclid +CMB-S4+LiteBIRD data set. Previous studies have shown that a correlation between ∑m_ν and {w₀, w_a} is already present in the analysis of CMB data alone (Brinckmann et al. 2019). As a matter of fact, the effect of the total neutrino mass and of the dark energy parameters on CMB lensing are difficult to disentangle. The effects of ∑m_ν (with fixed sound horizon scale, like in the bottom panel of Fig. 1) and of {w₀, w_a} on the overall amplitude of the matter power spectrum and on the scale of CMB and BAO peaks are also relatively similar, and do not allow us to remove this correlation. Such degeneracies can be removed by including in the analysis additional statistics, such as the bispectrum, as explained in Hahn et al. (2020) (see also Ajani et al. 2020; Chudaykin & Ivanov 2019; Uhlemann et al. 2020). It is nevertheless interesting to see that the constraint on ∑m_ν never degrades by more than a factor of 2 when the cosmological constant is replaced by dynamical dark energy.

When the dark energy and radiation density parameters {w₀, w_a, ΔN_eff} are varied simultaneously, the sensitivity remains close to that in the w₀w_aCDM+∑m_ν model. In the w₀w_aCDM+∑m_ν+ΔN_eff model, it is still possible to get about 2 σ evidence for a non-zero mass from the Euclid +CMB-S4+LiteBIRD combination.

8. Conclusions

Building on the work of Euclid Collaboration (2024a), we have implemented and validated a pipeline to consistently forecast the constraints on the neutrino parameters (such as mass and abundance) that will be measurable from the upcoming Euclid mission of the European Space Agency. We have studied how to model the observables of the photometric survey (allowing for cosmic shear and galaxy clustering measurements) and those of the spectroscopic survey (allowing for a reconstruction of the galaxy clustering power spectrum) for cosmologies including massive neutrinos, additional relativistic degrees of freedom, and variations in the dark energy equation of state.

We find that the Fisher matrices computed with the CosmicFish and MontePython pipelines agree to a 10% level even after marginalising over all corresponding nuisance parameters, independently of which EBS (CAMB or CLASS) is used. For the purposes of this implementation, we have also compared various nonlinear treatments. We conclude that for our purpose the predictions of HMcode (Mead et al. 2021) are the closest to numerical simulations for the cosmological models considered here, even outperforming some of the emulators that we have considered.

The validated forecast pipeline predicts a great deal of complementarity between the various probes measured by Euclid, particularly between the photometric and spectroscopic surveys. The combination of these probes will help to break the degeneracies between the parameters {H₀, σ₈, ΔN_eff, ∑m_ν, w₀, w_a}. There is a strong complementarity with additional probes as well, which can be used to lift the residual degeneracy among the small-scale fluctuation amplitude, Hubble constant, and the neutrino parameters. Once combined with data from CMB anisotropy experiments, the Euclid data will lead to extremely tight constraints on {σ₈, H₀, ∑m_ν, N_eff}.

Indeed, we find that while Euclid data alone will already allow us to constrain σ(∑m_ν = 60 meV) = 53 meV in the standard ΛCDM+∑m_ν model, its combination with CMB anisotropy data from Planck will lead to a sensitivity of σ(∑m_ν = 60 meV) = 23 meV. Thus, we can expect Euclid +Planck data to raise evidence for a non-zero neutrino mass at least at the 2.6 σ level. If the combined data turns out to be best fit with ∑m_ν ≃ 60 meV, it would indicate a preference for normal ordering over inverted ordering roughly at the 2 σ level. In combination with future CMB data from LiteBIRD and CMB Stage-IV, Euclid will reach a sensitivity of around σ(∑m_ν = 60 meV) = 16 meV, which will allow for a neutrino mass detection at almost 4 σ. Under the same assumption, we would obtain a clear disambiguation of the neutrino mass hierarchy at 2.5 σ. These constraints are only mildly sensitive to the addition of ultra-relativistic degrees of freedom (ΔN_eff > 0) and end up being degrading only by a factor of roughly 2 when considering a flexible model of time-varying dark energy.

Euclid is also expected to measure the Hubble parameter with exquisite precision in these cosmologies, with an uncertainty expected to reach 0.1–0.3 km s⁻¹ Mpc⁻¹ in combination with CMB data. This will bring very strong clues regarding the puzzle raised by the current Hubble tension.

The determination of the effective neutrino number N_eff by CMB experiments is limited by a parameter degeneracy between N_eff and H₀. The Euclid probes alone will be sensitive to N_eff and, even more importantly, when combined with CMB data, they will be able to partially lift this degeneracy. Assuming a fiducial model with no additional relativistic degrees of freedom beyond standard neutrinos (N_eff = 3.044), the combination of Planck data and Euclid primary probes will provide a 95% CL upper bound ΔN_eff < 0.144 in the ΛCDM+∑m_ν+N_eff case. The bound will be further tightened to ΔN_eff < 0.063 with future CMB data from LiteBIRD and CMB Stage-IV. These predictions remain stable when marginalising over dynamical dark energy parameters. Such bounds will constrain many well-motivated models of additional relativistic particles such as Goldstone bosons, Weyl fermions, or massless gauge bosons. Not only will decoupling times after the QCD transition be potentially ruled out by Euclid alone, but in combination with CMB experiments, decoupling times even slightly before the QCD transition can be constrained. Indeed, in combination with future CMB surveys, Euclid will come close to excluding the presence of these particles at any decoupling temperature (even far before the QCD transition, see Fig. 18).

In summary, Euclid will put tight constraints on the parameters of the ΛCDM model as well as some of its most commonly considered extensions, while bringing strong indications concerning the Hubble constant tension. In combination with either current or future CMB data, Euclid will constrain neutrino physics, potentially providing the first detection of a non-zero neutrino mass, along with some strong hints about the mass ordering. Moreover, the combined data will have an enhanced sensitivity to the effective number of neutrinos and a significant potential to rule out many candidates for light relics beyond the Standard Model. As such, Euclid will consolidate the ability of precision cosmology to probe neutrino physics.

Acknowledgments

The Euclid Consortium acknowledges the European Space Agency and a number of agencies and institutes that have supported the development of Euclid, in particular the Agenzia Spaziale Italiana, the Austrian Forschungsförderungsgesellschaft funded through BMK, the Belgian Science Policy, the Canadian Euclid Consortium, the Deutsches Zentrum für Luft-und Raumfahrt, the DTU Space and the Niels Bohr Institute in Denmark, the French Centre National d’Etudes Spatiales, the Fundação para a Ciência e a Tecnologia, the Hungarian Academy of Sciences, the Ministerio de Ciencia, Innovación y Universidades, the National Aeronautics and Space Administration, the National Astronomical Observatory of Japan, the Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the Research Council of Finland, the Romanian Space Agency, the State Secretariat for Education, Research, and Innovation (SERI) at the Swiss Space Office (SSO), and the United Kingdom Space Agency. A complete and detailed list is available on the Euclid web site (https://www.euclid-ec.org). We warmly thank Julien Bel for providing matter power spectra measurements from the DEMNUni simulations. N.F. is supported by the Italian Ministry of University and Research (MUR) through Rita Levi Montalcini project “Tests of gravity at cosmological scales” with reference PGR19ILFGP. N.F. and F.P. also acknowledge the FCT project with ref. number PTDC/FIS-AST/0054/2021. F.P. acknowledges partial support from the INFN grant InDark and the Departments of Excellence grant L.232/2016 of the Italian Ministry of University and Research (MUR). S.P.’s simulations were performed with computing resources granted by RWTH Aachen University under project thes1329. Z.S. acknowledges funding from DFG project 456622116 and support from the IRAP and IN2P3 Lyon computing centers. N.S. acknowledges support from the Maria de Maetzu fellowship grant: CEX2019-000918-M, financiado por MCIN/AEI/10.13039/501100011033. The DEMNUni simulations were carried out in the framework of “The Dark Energy and Massive-Neutrino Universe” project, using the Tier-0 IBM BG/Q Fermi machine and the Tier-0 Intel OmniPath Cluster Marconi-A1 of the Centro Interuniversitario del Nord-Est per il Calcolo Elettronico (CINECA). We acknowledge a generous CPU and storage allocation by the Italian Super-Computing Resource Allocation (ISCRA) as well as from the coordination of the “Accordo Quadro MoU per lo svolgimento di attività congiunta di ricerca Nuove frontiere in Astrofisica: HPC e Data Exploration di nuova generazione”, together with storage from INFN-CNAF and INAF-IA2.

References