© ESO 2021

1. Introduction

Next-generation surveys (stage IV) of the cosmic large-scale structure will greatly improve both the amount and quality of data for cosmological investigations. For instance, in the coming decade the surveys carried out by Euclid¹, the Vera C. Rubin Observatory², and the Nancy Grace Roman Space Telescope³ will probe scales and redshifts that were previously inaccessible. The correct interpretation of such a large amount of high-quality data, however, poses a challenge for our theoretical modelling.

In this paper we explore the current modelling limitations for one of the most promising probes: cosmic shear, the measurement of the apparent distortions of galaxy shapes caused by the weak lensing (WL) effect of intervening matter between us and distant sources (see Kilbinger 2015, for a recent review). It provides a direct way to trace the distribution of matter, and as a result it can be used to infer the total matter power spectrum, P_δδ(k, z). In contrast, the galaxy power spectrum, P_gg(k, z), estimated from the clustering of galaxies, depends on the galaxy bias and how galaxies occupy non-linear structures.

A complication is that much of the constraining power of the cosmic shear signal relies on the ability to interpret scales far into the non-linear regime, corresponding to wavenumbers k ≈ 7 h Mpc⁻¹ (e.g., Huterer & Takada 2005; Semboloni et al. 2011; Taylor et al. 2018b). On those scales, perturbations of the matter density field are no longer small, and linear theory cannot be used to predict the evolution of large-scale structures.

There are theoretical approaches for predicting clustering beyond the limit of linear theory, such as: standard perturbation theory (Blas et al. 2014; see also Bernardeau et al. 2002 for a detailed review) that includes higher-order terms; renormalised perturbation theory (Crocce & Scoccimarro 2006; Crocce et al. 2012; Blas et al. 2016); response functions (Nishimichi et al. 2016); effective field theory (Baumann et al. 2012); and the reaction method (Cataneo et al. 2019). These methods are able to achieve accuracies on power spectrum predictions of ≈1% with respect to numerical simulations, up to scales k ≲ 0.3 h Mpc⁻¹ (see e.g., Foreman & Senatore 2016; Beutler et al. 2017; D’Amico et al. 2020). This is sufficient for modelling the mildly non-linear regime, where the baryonic acoustic oscillation peak is located, but these techniques cannot be used to predict the signal in the highly non-linear regime that WL analyses will probe. The reaction method, however, can in principle achieve an accuracy of ≈1% up to scales k ≲ 10 h Mpc⁻¹ when using the results of an emulator as a baseline (Giblin et al. 2019; Cataneo et al. 2020).

The common approach of modelling the non-linear part of the power spectrum instead relies on fitting formulae determined from comparisons to N-body simulations of cold dark matter particles (e.g., Halofit; Smith et al. 2003). While very economical in terms of computational time, these fitting formulae have a limited range of applicability because the simulations they are based on assume a specific model – usually the cosmological constant cold dark matter (ΛCDM) model or minimal extensions with a constant dark energy equation of state parameter allowed to be different from −1. Thus, applying these corrections to models outside the range constrained by the simulations, for example a more general dark energy fluid, may lead to biased results (see e.g., Casarini et al. 2011a; Seo et al. 2012).

Recently, so-called emulators – based on a large suite of simulations, such as the Coyote Universe (Heitmann et al. 2010, 2009, 2014; Lawrence et al. 2010), the Mira Universe (Heitmann et al. 2016; Lawrence et al. 2017), the Euclid Emulator Project (Knabenhans et al. 2019), and the BACCO simulation project (Angulo et al. 2020)– have been proposed as an alternative to fitting formulae. Emulators interpolate between high-resolution simulation runs at key points (nodes) in the cosmological parameter space. The main advantage of an emulator with respect to a fitting formula is that it does not degrade the accuracy of the corrections over the parameter space sampled by the simulations, such as the range in redshift and scales.

Another way to predict the matter power spectrum on small scales is provided by HMCode⁴ (Mead et al. 2015), which is based on the analytical halo model (Peacock & Smith 2000; Seljak 2000; Cooray & Sheth 2002) and tuned to match the Coyote Extended Emulator (Heitmann et al. 2014) results. It has subsequently been improved to include effects of neutrinos, chameleon and Vainshtein screened models, and dynamical dark energy (Mead et al. 2016) ⁵.

All these methods depend on the quality of the simulations on which they are based, in addition to the declared quality of the method itself (for example the accuracy of the fitting formula or the accuracy of the interpolation of the emulators). Restricting ourselves to just the non-linear clustering, the agreement between purely dark matter simulations is limited by the size of the simulated volume, the number of particles employed in the simulation, and the choice of initial conditions (see e.g., Casarini et al. 2015; Schneider & Teyssier 2015). Hence, part of the differences between the methods described above may be attributed to the simulation parameters on which they are based (e.g., when dimension and resolution are insufficient) rather than the methods themselves.

Moreover, our ability to extract cosmological information from WL measurements on small scales is limited further by baryonic feedback processes (Semboloni et al. 2011) because gas cooling, star formation, galactic winds, supernova explosions, and feedback from active galactic nuclei (AGNs) modify the expected distribution of matter on small scales (Jing et al. 2006; Rudd et al. 2008; Casarini et al. 2011b, 2012; van Daalen et al. 2011; Castro et al. 2018; Debackere et al. 2020). Accurate predictions of the matter power spectrum on those scales require hydrodynamical simulations that not only need to reproduce the non-linear clustering of cold dark matter particles, but should also reliably describe the baryonic component. Such hydrodynamical simulations are much more demanding in terms of computational resources, and the impact of baryonic feedback extracted from hydrodynamical runs has to be modelled and incorporated in the reconstruction of the cosmic shear signal (see e.g., Schneider & Teyssier 2015; Schneider et al. 2016, for a method to include the impact of feedback in the data analysis pipeline). In order to reduce the computational power requested to include such effects, methods relying on ’baryonification’ algorithms that mimic the effects of astrophysical processes induced by baryons have been proposed (see e.g., Aricò et al. 2021, 2020).

Given the cost of simulating large volumes with high resolution for every necessary point in the parameter space for each specific cosmological model, we are particularly interested in techniques that can drastically reduce the number of simulations (see Linder & White 2005; Francis et al. 2007; Casarini et al. 2009). As an example, in this work we use the PKequal⁶ method, which allows us to determine the non-linear power spectrum of a dynamical dark energy model at a particular redshift with an ensemble of non-linear spectra of simpler constant w models (Casarini et al. 2009, 2016).

In this paper, we investigate the impact of different implementations of the non-linear corrections and baryonic feedback prescriptions on cosmological parameter estimation, adopting cosmic shear measurements from Euclid (Laureijs et al. 2011) as our baseline (for an analysis of the impact of non-linear corrections on other LSS observables see e.g., Safi & Farhang 2020). We note that we do not investigate here other effects that might affect the parameter estimation pipeline, such as the common assumption of a Gaussian likelihood, which is only an approximation. The paper is organised as follows. After describing the various non-linear recipes in Sect. 2, in Sect. 3 we summarise the Euclid specifications relevant for our analysis. As a first example of the impact of non-linear prescriptions, we assess in Sect. 4 the constraining power of the Euclid survey on the dark energy parameters, here assumed to be described by the so-called CPL parameterisation (Chevallier & Polarski 2001; Linder 2003). In Sect. 5 we quantify the shifts in cosmological parameters when a wrong correction pipeline is used, focusing on the combination of Planck and mock Euclid data. In Sect. 6 we examine the impact of baryonic feedback.

2. Available non-linear prescriptions

In this section we describe the techniques that we use in this study to compute the matter power spectrum in the deeply non-linear regime.

2.1. Halofit

One of the first widely used prescriptions to model the non-linear part of the power spectrum, called Halofit, was developed by Smith et al. (2003). The authors measured the non-linear evolution of the matter power spectra using a large library of cosmological N-body simulations with power-law initial spectra (Jenkins et al. 1998).

The Halofit approach is based on the halo model (Peacock & Smith 2000; Seljak 2000; Ma & Fry 2000), in which the density field is described in terms of the distribution of isolated dark matter haloes. The correlations in the field are assumed to arise from the clustering of haloes with respect to each other on large scales, and through the clustering of dark matter particles within the same halo on small scales. The total non-linear power spectrum, P_NL(k), can then be decomposed into

$$\begin{array}{c}\hfill P_{\mathrm{NL}}(k)=P_{\mathrm{Q}}(k)+P_{\mathrm{H}}(k),\end{array}$$(1)

where P_Q(k) is the quasi-linear term related to the large-scale contribution to the power spectrum, and P_H(k) describes the contribution from the self-correlation of haloes. These terms are also known as the two-halo and the one-halo term, respectively, and we discuss them in this order below.

Seljak (2000), Ma & Fry (2000), Scoccimarro et al. (2001) proposed to use linear theory filtered by the effective window that corresponds to the distribution of haloes as a function of mass, n(M), convolved with their density profiles, $\stackrel{\sim }{\rho }(k,M)$, and a prescription for their bias with respect to the underlying mass field, b_H(M). The quasi-linear term can then be expressed as

$$\begin{array}{c}\hfill P_{\mathrm{Q}}(k)=P_{\mathrm{L}}(k){\left[\frac{1}{\overline{\rho }}{\displaystyle \int }\mathrm{d}M{b}_{\mathrm{H}}(M)n(M)\stackrel{\sim }{\rho }(k,M)\right]}^2,\end{array}$$(2)

with $\overline{\rho }$ the homogeneous background matter density, and P_L(k) the linear power spectrum.

A simpler approach was proposed by Peacock & Smith (2000), who assumed that the quasi-linear term corresponds to pure linear theory, P_Q(k) = P_L(k). However, quasi-linear effects must modify the relative correlations of haloes away from linear theory, irrespective of the allowance made for the finite size of the haloes (see Smith et al. 2003, and references therein). Halofit takes then an empirical approach, allowing the quasi-linear term to depend on n(M), and truncating its effects at small scales. If we define the dimensionless power spectrum as

$$\begin{array}{c}\hfill {\mathrm{\Delta }}^2(k)\equiv \frac{k^3}{2{\pi }^2}P(k),\end{array}$$(3)

the quasi-linear term in Halofit is given by

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\mathrm{Q}}^2(k)={\mathrm{\Delta }}_{\mathrm{L}}^2(k)\frac{{[1+{\mathrm{\Delta }}_{\mathrm{L}}^2(k)]}^{{\beta }_n}}{1+{\alpha }_n{\mathrm{\Delta }}_{\mathrm{L}}^2(k)}{\mathrm{e}}^{-f(y)},\end{array}$$(4)

where y ≡ k/k_σ, k_σ is a non-linear wavenumber related to the spherical collapse model (Press & Schechter 1974; Sheth & Tormen 1999; Jenkins et al. 2001), α_n and β_n are coefficients sensitive to the input linear spectrum, and f(y) = y/4 + y²/8 governs the decay rate at small scales.

To describe the clustering of matter on small scales, we need a description for P_H, which is given by

$$\begin{array}{c}\hfill P_{\mathrm{H}}(k)=\frac{1}{{\overline{\rho }}^2{(2\pi )}^3}{\displaystyle \int }\mathrm{d}Mn(M){|\stackrel{\sim }{\rho }(k,M)|}^2.\end{array}$$(5)

In order to model this term, we can use an expression that looks like a shot-noise spectrum on large scales, but progressively vanishes on small scales by the filtering effects of halo profiles and the mass function. A good candidate is

$$\begin{array}{c}\hfill {{\mathrm{\Delta }}_{\mathrm{H}}^2}^{\prime }(k)=\frac{a_n{y}^3}{1+b_{n}y+c_n{y}^{3-{\gamma }_n}},\end{array}$$(6)

where a_n,  b_n,  c_n, and γ_n are dimensionless numbers that depend on the input spectrum.

Cooray & Sheth (2002), however, showed that with this expression the halo model disagrees with low-order perturbation theory in some cases. To solve this, Halofit modifies Eq. (6) to obtain a spectrum steeper than Poisson on the largest scales

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\mathrm{H}}^2(k)=\frac{{{\mathrm{\Delta }}_{\mathrm{H}}^2}^{\prime }(k)}{1+{\mathrm{\mu }}_n{y}^{-1}+{\nu }_n{y}^{-2}},\end{array}$$(7)

where μ_n and ν_n are, again, coefficients that depend on the input spectrum. Smith et al. (2003) showed that Halofit is able to reproduce the measurements from simulations more accurately, and down to smaller scales, than the halo model.

2.2. Halofit with Bird and Takahashi corrections

All the coefficients used in the Halofit recipe were determined by Smith et al. (2003) from a fit to cold dark matter simulations in boxes of lengths of 256 h⁻¹ Mpc containing 256³ particles (Jenkins et al. 1998). As a consequence of the relatively large particle mass, Halofit may not be suitable if we want to test cosmologies with massive neutrinos, or go down to very small scales, where the impact of baryonic interactions is non-negligible. Moreover, the limited simulation volume results in large sample variance (see Casarini et al. 2015; Schneider & Teyssier 2015), and thus may lead to inaccurate results even for a ΛCDM cosmology (White & Vale 2004; Casarini et al. 2009, 2012; Hilbert et al. 2009; Heitmann et al. 2010). Finally, as the simulations were performed for the standard cosmological model, using Halofit with a dark energy equation of state w ≠ −1 may yield an incorrect estimate of the power spectrum (Casarini et al. 2011a; Seo et al. 2012).

To address these limitations, Bird et al. (2012) investigated the impact of massive neutrinos, and performed several N-body simulations of the matter power spectrum incorporating massive neutrinos with masses between 0.15 and 0.6 eV. They focussed on non-linear scales below 10 h Mpc⁻¹ at z < 3, and extended the Halofit approximation to account for massive neutrinos. They found that in the strongly non-linear regime Halofit over-predicts the suppression due to the free-streaming of the neutrinos. In particular, the asymptotic behaviour of the non-linear term in Halofit is given by ${\mathrm{\Delta }}_{\text{H}}^2\sim y^{{\gamma }_n}$, and therefore Bird et al. (2012) adjusted γ_n to their ΛCDM simulations with massive neutrinos. Moreover, they modified the non-linear power spectrum with the ansatz

$$\begin{array}{c}\hfill {({\mathrm{\Delta }}_{\mathrm{NL}}^{\nu })}^2={\mathrm{\Delta }}_{\mathrm{NL}}^2(1+Q_{\nu }),\end{array}$$(8)

with

$$\begin{array}{c}\hfill Q_{\nu }=\frac{l{f}_{\nu }}{1+k^{3}m},\end{array}$$(9)

where f_ν = Ω_ν/Ω_m is the ratio between the neutrino and total matter energy densities, and l and m are fitted to the simulations. They also modified Eq. (4) to

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\mathrm{Q}}^2(k)={\mathrm{\Delta }}_{\mathrm{L}}^2(k)\frac{{[1+{\stackrel{\sim }{\mathrm{\Delta }}}_{\mathrm{L}}^2(k)]}^{{\stackrel{\sim }{\beta }}_n}}{1+{\alpha }_n{\stackrel{\sim }{\mathrm{\Delta }}}_{\mathrm{L}}^2(k)}{\mathrm{e}}^{-f(y)},\end{array}$$(10)

with

$$\begin{array}{cc}& {\stackrel{\sim }{\mathrm{\Delta }}}_{\mathrm{L}}^2={\mathrm{\Delta }}_{\mathrm{L}}^2(1+\frac{p{f}_{\nu }{k}^2}{1+1.5k^2}),\hfill \end{array}$$(11)

$$\begin{array}{cc}& {\stackrel{\sim }{\beta }}_n={\beta }_n+f_{\nu }(r+n^{2}s),\hfill \end{array}$$(12)

where p, r, and s are fitted to the simulations.

Another important improvement to the original Halofit was introduced in Takahashi et al. (2012), who updated the fitting parameters using high-resolution N-body simulations with box sizes of L = 300−2000 h⁻¹ Mpc with n_p = 1024³ particles each, for 16 cosmological models around the best fitted cosmological parameters from WMAP data (Hinshaw et al. 2013), including dark energy models with a constant equation of state. This revised version of Halofit provides an accurate prediction of the non-linear matter power spectrum down to k ∼ 30 h Mpc⁻¹ and up to z ≥ 3 with an accuracy ∼5–10%. In the remainder of this paper we refer to the non-linear prescription that includes the improvements from Takahashi et al. (2012) and Bird et al. (2012) as Halofit, for simplicity.

2.3. Halofit with PKequal

One of the limitations of the standard Halofit approach, even after the corrections by Bird et al. (2012) and Takahashi et al. (2012) are considered, is that it is based on a fit to N-body simulations with a constant value for the dark energy equation of state. However, given the precision of stage IV surveys, we are particularly interested in determining whether the data prefer an evolving dark energy equation of state. To avoid biases in our non-linear predictions, we could run N-body simulations that include a time dependence for w, such as the CPL parameterisation (Chevallier & Polarski 2001; Linder 2003) given by

$$\begin{array}{c}\hfill w(a)=w_0+w_a(1-a).\end{array}$$(13)

This approach implies the need for significant computational resources. Another option, however, is to map the non-linear power spectra of dark energy models with a constant equation of state to those with a time varying one. In this work we consider the PKequal code (Casarini et al. 2016), which implements the spectral equivalence from Casarini et al. (2009). Francis et al. (2007) showed how predictions for constant w models at z = 0 can be related to the power spectra of cosmologies with an evolving equation of state w(a) given by Eq. (13) with an accuracy ∼0.5% up to k ≃ 1 h Mpc⁻¹. The PKequal technique achieves this precision also at z > 0 for a general equation of state w = w(a) by imposing the equivalence of the distance to the last scattering surface and requiring that the amplitudes of the density fluctuations at the redshift of interest are the same. For a given set of values of w₀ and w_a, these two conditions yield at each z a unique value of w_eq and σ_8, eq for the constant w model.

The performance of this method has been tested for several dark energy models (Casarini et al. 2009; Casarini 2010), and also in the presence of gas cooling, star formation and SN feedback (Casarini et al. 2011b). These studies find differences in power spectra between the mapped dynamical dark energy models and the ensemble of equivalent constant w models that are within 1 % up to k ≃ 2–3 h Mpc⁻¹. With this method it is possible to extend both emulators (see Casarini et al. 2016) and fitting formulae (as in this work) to dynamical dark energy models if they are valid for constant w models.

2.4. HMCode

An alternative approach to predict the non-linear matter power spectrum, called HMCode, was proposed by Mead et al. (2015). It introduces physically motivated free parameters into the halo model formalism, instead of using empirical fitting functions. Mead et al. (2015) fit these to N-body simulations with box sizes L = 90–1300 h⁻¹ Mpc and n_p = 512³–1024³ particles for a variety of ΛCDM and wCDM models (Heitmann et al. 2010). HMCode also accounts for the effects of baryonic feedback on the power spectrum by fitting the halo model to hydrodynamical simulations that include gas cooling, star formation, as well as supernova and AGN feedback (Schaye et al. 2010; van Daalen et al. 2011).

In Mead et al. (2016) HMCode was updated to account for massive neutrinos, chameleon and Vainshtein screening mechanisms, as well as evolving dark energy equations of state described by the CPL parameterisation. Throughout the rest of the paper we will consider this latest version when referring to HMCode. We note, however, that the PKequal approach can also be used with the original HMCode by applying it to the equivalent {w_eq, 0} constant w model at any redshift z, and imposing the same σ₈(z). This yields the prediction for a given {w₀, w_a} CPL model, with an accuracy similar to that of the fit and the simulations on which HMCode itself is based.

2.5. Comparison

Before we study how the different methods affect the inference of cosmological parameters and related confidence regions, we show by how much the power spectra differ as a function of the wavenumber k for the prescriptions described above. We note that theoretical approaches based on perturbation theory or the effective field theory of large-scale structure do not yet provide accurate power spectra for the small scales we consider here.

In Fig. 1 we show the relative difference between the various predictions and the one from Halofit. The power spectra were computed using CAMB (Lewis et al. 2000) at redshift 0, where discrepancies are most pronounced. Hence, these represent a worst case scenario as WL probes mostly intermediate redshifts. The upper-left panel shows the relative differences between the linear power spectrum (orange), Halofit with PKequal (blue), and HMCode (red) for the ΛCDM model using the Planck TTTEEE+lowE+lensing+BAO 2018 mean values for the various cosmological parameters (Planck Collaboration VI 2020). The linear power spectrum diverges from Halofit for k > 0.1 h Mpc⁻¹, while the disagreement between Halofit and HMCode is below 7.5% for k < 10 h Mpc⁻¹. While in general non-linear corrections boost the matter power spectrum, Fig. 1 shows how the linear power spectrum takes larger values than the corrected one in the range 0.005 ≲ k ≲ 0.1. This feature, which is slightly counter-intuitive though well known, is due to the effect of the quasi-linear term of the halo model, relevant at these scales, that can lead to a suppression of power with respect to the linear case (Smith et al. 2003).

Fig. 1. Relative difference between the various predicted matter power spectra and the Halofit prediction at z = 0 for different cosmologies. Top-left panel: using the mean values of the ΛCDM parameters given by Planck 2018 (Planck Collaboration VI 2020). Top-right panel: differences when the ΛCDM parameters are drawn from a five-dimensional Gaussian distribution with the diagonal given by five times the 1σ uncertainty quoted by Planck TTTEEE+lowE+lensing+BAO 2018. Bottom-left panel: differences when also the constant dark energy equation of state is drawn from a Gaussian distribution with a dispersion of 0.3. Bottom-right panel: difference for evolving dark energy models where, alongside w0, also wa is Gaussian distributed with a dispersion of 1.0. In the bottom panelsw0 and wa (if present) are always chosen in such a way that w(z) <  − 1/3. The feature visible at k ≈ 0.005 h Mpc−1 corresponds the scale at which the non-linear corrections are turned on in CAMB.

The upper-right panel shows the variation on small scales when we consider parameters around the mean values of ΛCDM. For this comparison we generated 100 power spectra with cosmological parameters drawn from a five-dimensional Gaussian distribution with the diagonal given by five times the 1σ uncertainty quoted by Planck TTTEEE+lowE+lensing+BAO 2018⁷. These results highlight that the discrepancy between Halofit and HMCode can be larger than 10% for k > 3 h Mpc⁻¹.

The bottom-left panel considers models beyond ΛCDM, allowing for a constant equation of state parameter, w, that can differ from −1. The different lines correspond to power spectra when we draw parameters from a six-dimensional Gaussian distribution where we adopt a dispersion of 0.3 for w (but we do require w < −1/3). The bands showing the discrepancy between the different non-linear corrections increase slightly (a 10% discrepancy is reached at scales of k = 1–2 h Mpc⁻¹), but the overall shape remains the same.

Finally, in the bottom-right panel we consider dynamical dark energy models with a dark energy equation of state given by Eq. (13). In this case we add w_a (so that we draw parameters from a seven-dimensional Gaussian) with a dispersion of 1.0 and we require the dark energy equation of state to be always smaller than −1/3, that is w₀ + w_a < −1/3. The overall shape for the bands is the same for HMCode and the linear spectra, although the discrepancies are significantly larger (already 10% at scales of k = 0.6 h Mpc⁻¹). We also note that, since we allow w_a to vary, there is a difference between Halofit and Halofit+PKequal.

Lawrence et al. (2017) presented an updated version of COSMIC EMU that includes massive neutrinos. Their high-resolution simulations are interpolated with an accuracy of ∼4%. Lawrence et al. (2017) compared their COSMIC EMU predictions with Halofit, and HMCode. When massive neutrinos are considered, the different approaches show differences of ∼20% and ∼15%, respectively, in the power spectra for scales above k = 0.1 Mpc⁻¹, indicating the need for further improvements.

3. The Euclid Cosmic Shear survey

Our aim is to investigate how the expected constraints on cosmological parameters from Euclid data depend on the recipe that is used to predict the matter power spectrum on non-linear scales, although we note that our finding are also relevant for other stage IV experiments. Euclid is an M-class mission of the European Space Agency (Laureijs et al. 2011) that will carry out a spectroscopic and a photometric survey of galaxies over an area of 15 000 deg². The cosmic shear measurements use high-quality imaging at optical wavelengths, supported by multi-band optical ground-based photometry and near-infrared observations by Euclid. The telescope is designed so that (residual) instrumental sources of bias in the observed cosmic shear signal are subdominant compared to the statistical uncertainties (e.g., Cropper et al. 2013; Euclid Collaboration 2020b). However, to achieve its objectives, it is essential that the signal can be accurately predicted in the non-linear regime. Although this is also relevant to fully exploit the data from the clustering of galaxies and the cross-correlations with the lensing signal, we focus on the cosmic shear case in this paper and defer a more comprehensive study to future work.

We adopt the baseline specifications for the Euclid data, which are described in Euclid Collaboration (2020a, hereafter EC19). The redshift distribution of the sources is given by

$$\begin{array}{c}\hfill n(z)\propto {\left(\frac{z}{z_0}\right)}^{2}exp[-{\left(\frac{z}{z_0}\right)}^{3/2}],\end{array}$$(14)

with $z_0=0.9/\sqrt{2}$, resulting in a mean redshift of ⟨z⟩ = 0.96. The sample is divided into 10 equi-populated redshift bins n_i(z) (with the index i indicating the tomographic bin). We assume an average number density of galaxies with precise shape measurements of ${\overline{n}}_{\mathit{g}}=30{\mathrm{arcmin}}^{-2}$. To capture the noise arising from the intrinsic galaxy ellipticities we adopt a dispersion of σ_ϵ = 0.21 for each of the two ellipticity components. We also use the same approach as EC19 to compute the WL power spectrum, defined as

$$\begin{array}{c}\hfill C_{\mathit{ij}}^{ϵϵ}(\ell )=c{\displaystyle \int }\mathrm{d}z\frac{W_i^{ϵ}(z)W_j^{ϵ}(z)}{H(z)r^2(z)}{P}_{\delta \delta }[\frac{\ell +1/2}{r(z)},z],\end{array}$$(15)

where P_δδ is the non-linear matter power spectrum, r(z) is the comoving distance to redshift z, and the window function $W_i^{ϵ}$ is defined as

$$\begin{array}{c}\hfill W_i^{ϵ}(z)=W_i^{\gamma }(z)-\frac{{\mathcal{A}}_{\mathrm{IA}}{\mathcal{C}}_{\mathrm{IA}}{\mathrm{\Omega }}_{\mathrm{m},0}{\mathcal{F}}_{\mathrm{IA}}(z)H(z)n_i(z)}{D(z)c},\end{array}$$(16)

with n_i(z) normalised such that ∫dz n_i(z) = 1. In Eq. (16), the first term corresponds to the usual lensing kernel,

$$\begin{array}{cc}& W_i^{\gamma }(z)=\frac{3}{2}{\mathrm{\Omega }}_{\mathrm{m},0}{\left(\frac{H_0}{c}\right)}^2(1+z)r(z)\hfill \\ \hfill & \times {\displaystyle {\int }_z^{z_{\max }}\mathrm{d}{z}^{\prime }{n}_i(z^{\prime })[1-\frac{r(z)}{r(z^{\prime })}],}\hfill \end{array}$$(17)

whilst the second term models the effect of intrinsic alignments (IA). Here D(z) is the growth factor, 𝒞_IA = 0.0134 is a constant so that the normalisation of the model 𝒜_IA can be compared to the literature. To describe the dependence of the IA signal as a function of scale, redshift and galaxy luminosity, we adopt the extended non-linear alignment model (Joachimi et al. 2015), and the function ℱ_IA(z) is given by

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

where ⟨L⟩(z) and L_⋆(z) are the redshift-dependent mean and a characteristic luminosity of source galaxies as computed from the luminosity function, respectively. The parameters η_IA, β_IA, and 𝒜_IA are free parameters that can be determined observationally. We use {𝒜_IA, η_IA, β_IA}={1.72, −0.41, 2.17} as fiducial values, as was done in EC19, while they are allowed to vary in the parameter estimation.

Our approach for the modelling of IA is phenomenological, but we note that more physically motivated models have been proposed. For instance, Fortuna et al. (2021) used a halo model approach to link direct observations of IA to implications for cosmic shear. Finally, we note that the modelling of the IA signal is linked to non-linear structure formation, but exploring this further is beyond the scope of this paper.

We adopt the same approach presented in EC19 to model the covariance of the WL observable. As usual, it can be split into Gaussian and non-Gaussian contributions (Takada & Hu 2013; Cooray & Hu 2001; Hamilton et al. 2006; Hu & Kravtsov 2003; Kayo et al. 2012). The latter involve the convergence trispectrum, but there are significant uncertainties on how to model this quantity. Relying on the halo-model formalism (see e.g., Cooray & Hu 2001), it was shown in EC19 that the signal-to-noise ratio of the shear power spectrum (with the specifications used in this work) decreases by 30% at ℓ_max = 5000 when we add the non-Gaussian contribution. This loss of information corresponds to an effective cut at ℓ_max = 1420 when using only a Gaussian covariance. Following EC19, we decide to avoid uncertainties related to the non-Gaussian modelling and use a Gaussian covariance with a cut at ℓ_max = 1500. This corresponds to the pessimistic scenario. We also consider the impact of changing this cut in the optimistic settings, where we set it at ℓ_max = 5000. It is important to note that an ℓ-cut preserving the signal-to-noise ratio is not directly linked to constraints on the parameters (see e.g., Copeland et al. 2018), but the amplitude of the impact of the non-Gaussian modelling used here is small compared to other assumptions.

4. Impact of non-linear corrections on forecast constraints

Given our desire to use measurements on small scales to estimate cosmological parameters, it is essential to assess how the different methods to model the power spectrum on highly non-linear scales affect the results. In the following, we limit ourselves to a simple extension of the ΛCDM model: We assume that dark energy is dynamical with its equation of state parameterised according to Eq. (13), and that density perturbations for this component are well described by the parameterised post-Friedmann approach, which assumes that the dark energy field remains smooth with respect to matter at the scales of interest (Hu & Sawicki 2007; Hu 2008; Fang et al. 2008).

To evaluate the impact of the different recipes for the non-linear power spectrum on the final parameter estimation from Euclid, we use the Fisher matrix approach (see EC19, for an extensive review of the methodology). For the different methods we compute the figure of merit (FoM) for the parameters w₀ and w_a, where the FoM is defined as

$$\begin{array}{c}\hfill \mathrm{FoM}=\sqrt{\det {\stackrel{\sim }{\mathsf{F}}}_{w_0{w}_a}},\end{array}$$(19)

and ${\stackrel{\sim }{\mathsf{F}}}_{{\mathit{w}}_0{\mathit{w}}_a}$ is the Fisher matrix marginalised over all the cosmological parameters except for w₀ and w_a. This allows us to examine whether degeneracies between cosmological parameters differ when switching from one method to the other. We can also quantify the biases in dark energy parameters and changes in the FoM, which captures the performance of Euclid.

The free parameters for this analysis are: the matter and baryon density parameters Ω_m, 0 and Ω_b, 0; the dark energy parameters, w₀ and w_a; the spectral index of primordial perturbations, n_s; the dimensionless Hubble parameter, h; and σ₈, that is the root mean square of present-day linearly evolved density fluctuations in spheres of 8 h⁻¹ Mpc radius. The fiducial values of these parameters are listed in Table 1.

We used the specification of Euclid for WL observations that were detailed in Sect. 3 and computed the FoM from the WL power spectra $C_{\mathcal{l}}^{ϵϵ}$ for different values of the maximum multipole, ℓ_max. We considered two different values for ℓ_max, namely 1500 and 5000, with the latter probing deeper into the non-linear regime. We note that here the same ℓ_max applies to all redshift bins, thus leading to a different cut in scales (k_max) for each bin. We investigated the difference between this analysis and one considering a k_max rather than a multipole cut in Appendix A. In both cases the minimum multipole was fixed to ℓ_min = 10, following EC19. We note that the cuts we used here should be considered an approximation; in general, there is no direct mapping between ℓ_max and k_max, and a more refined approach would be needed to convert a cut in multipoles into a cut in wavenumber (Taylor et al. 2018b).

Table 2 lists the resulting FoM values for our baseline w₀w_aCDM cosmology for the three different non-linear recipes. The variation in the predicted FoMs is substantial, which is not that surprising given the differences we see in Fig. 1; at large k the differences exceed 10% when the parameters are allowed to vary with respect to the fiducial model.

The large variation is caused by two separate effects. Firstly, the fiducial $C_{\ell }^{ϵϵ}$ is obtained by integrating up to k ≃ 30 h Mpc⁻¹ and then inverted when computing the covariance matrix that enters the Fisher matrix forecast. Small differences can become large in the inversion process. Secondly, what is important is not so much the fiducial model itself, but rather the derivatives of the power spectrum with respect to the model parameters. Different non-linear corrections predict different derivatives, thus leading to different Fisher matrix elements. This is supported by the fact that the FoM for HMCode and Halofit+PKequal are very similar. These two prescriptions explicitly take deviations of w₀ and w_a from the fiducial ΛCDM values into account, while this is not the case for the standard Halofit, which is designed to describe ΛCDM cosmologies and extended to cases where dark energy is described by a constant equation of state. As a consequence, the amplitude of the derivatives of the matter power spectrum with respect to {w₀,  w_a} are similar between HMCode and Halofit+PKequal, but different for Halofit. This explains why the FoM values in the third and fourth columns are so similar. We also tested the impact of the differences in the covariance matrix and in the derivatives with respect to w_a, obtaining the FoM through the described procedure, but keeping fixed to the a reference non-linear prescription either the first or the latter.

It is also worth noting that the FoM is highest for Halofit because in that case w_a impacts the power spectra at non-linear scales through its effect on the linear matter power spectrum, while in the other two methods the non-linear corrections are affected by this parameter as well. Overall, the non-linear corrections dampen the derivatives with respect to w_a, leading to weaker constraints on this parameter and to a lower FoM. These results demonstrate that the FoM depends critically on the non-linear model that is used, highlighting the need for (more) accurate prescriptions.

In addition to this, Table 2 shows also that the difference in the FoM between HMCode and HaloFit decreases in the optimistic case with respect to the pessimistic one. While this might seem counter intuitive, this effect is due to the fact that the derivatives with respect to w_a, which are the main responsible for the variation in the FoM value, tend to vanish at high multipoles. For such a reason, including smaller scales in the analysis reduces the impact of such derivatives, and the main difference between the recipe is now in the covariance matrices. However, they are computed at the parameter fiducial value, where the differences between recipes are less significant. Instead, if we compare HaloFit with HaloFit+PKequal, the change in moving from the pessimistic to the optimistic case is negligible, due to the different construction of the latter recipe, which is a correction of HaloFit itself, rather than being based on a different approach.

5. Bias on cosmological parameter estimates

A major concern is that inaccuracies in the theoretical predictions on non-linear scales translate into shifts in the inferred cosmological parameters. To quantify the impact of such biases, we create a mock dataset of Euclid WL observations, using the specifications listed in Sect. 3. The fiducial cosmology that we use to generate these mock data is summarised in Table 3, where A_s is the amplitude of the primordial power spectrum, ω_b, 0 = Ω_b, 0h², ω_c, 0 = Ω_c, 0h² where Ω_c, 0 is the present time density parameter for cold dark matter, while for the intrinsic alignment nuisance parameters, also free in this analysis, we assume the fiducial values {𝒜_IA, η_IA, β_IA}={1.72, −0.41, 2.17}. Again, we follow EC19 in neglecting the uncertainties on the galaxy distributions and shear bias, for simplicity. However, we refer the reader to Tutusaus et al. (2020) and Euclid Collaboration (2020b) for analyses focussed on these systematic uncertainties. In contrast to what was done in Sect. 4, we assume here that the expansion history is provided by a dynamical dark energy, assuming w₀ = −0.9 and w_a = 0.1. This allows us to assess the impact of the different non-linear recipes when a non-standard model is the true underlying cosmology. This is motivated by the fact that the recipes in Sect. 2 differ in how they account for such extensions of the standard cosmological model.

We adopt the non-linear correction provided by Halofit+PKequal as the reference to which we compare the parameter estimates for the other recipes. We stress that we do not advocate the use of a particular prescription, but rather wish to quantify the shifts in the estimated cosmological parameters that may arise from using a different prediction.

We analyse this mock dataset with a Markov chain Monte Carlo method (MCMC), using the MontePython⁸ suite (Audren et al. 2013; Brinckmann & Lesgourgues 2019), with a Euclid lensing mock likelihood as presented in Sprenger et al. (2019)⁹. With this setting we sample the cosmological parameters listed in tab:fiducialmcmc using either Halofit (switching off PKequal) or HMCode. When using HMCode, the baryonic feedback parameters are fixed to the values fitting the COSMIC EMU dark-matter-only simulations (Heitmann et al. 2014).

As done in Sect. 4, we performed the analysis for ℓ_max = 1500 and 5000 in order to explore the impact of the non-linear corrections on the results when including high multipoles. Also in this case, we explore the difference with a scale cut in Appendix A. We complemented the WL mock dataset with TT, TE, EE, and lensing data from a mock Planck likelihood, thus reproducing the sensitivity of the full mission. We note that we did not use real Planck data to avoid a mismatch between the fiducial Euclid cosmology and the actual best fitted Planck values.

This approach enables us to determine the ‘bias’ (B) on cosmological parameters with respect to our fiducial cosmology, that is the offset of the mean of the estimated posterior distribution from the true value, which in turn allows us to quantify the impact of using a different prescription for the non-linear evolution of density perturbations. To assess the significance of a particular bias, it is useful to compare it to the expected statistical uncertainty. We therefore consider the (relative) bias on all the cosmological parameters given by

$$\begin{array}{c}\hfill B(\theta )=\frac{|{\theta }^{\star }-{\theta }^{\mathrm{fid}}|}{\sigma },\end{array}$$(20)

where θ^⋆ is the mean value of the parameter found with the MCMC analysis, σ is the 68% uncertainty estimated from the chains, and θ^fid is its fiducial value. We note that this estimate implicitly assumes that the posterior distribution can be approximated by a Gaussian, which is indeed valid for constraints based on the combination of Planck and Euclid data.

The statistical uncertainties of Euclid are determined by the survey design (see Sect. 3). To draw reliable conclusions from these data, it is essential that systematic uncertainties are sufficiently small. Ideally biases should vanish, but generally it is too costly to achieve this. A reasonable compromise, however, is to adopt B_thr ≲ 0.1 (see e.g., Sect. 4.1 in Massey et al. 2013), which is what we do here¹⁰. We note that the intrinsic variance of the mean estimated by the MCMC also contributes to our estimate of B(θ). We quantified this contribution by computing the scatter of the mean value by bootstrapping the chains. We find that it is always within 1% of the final error on cosmological parameters.

The results are reported in Table 4 and the posteriors are presented in Fig. 2. Although the optical depth τ is a free parameter in the model, it is effectively constrained by the Planck measurements alone. We therefore do not report its value here. Oppositely, as a comparison, we can see how constraints on other parameters tighten thanks to the inclusion of Euclid data. For instance, Planck Collaboration VI (2020) report for σ₈ a 68% confidence level (C.L.) interval equal to 1.5% or 1% of the parameter value in a w₀w_aCDM cosmology, when Planck is respectively combined with RSD and WL data or BAO and SN data. Instead, with Planck-TT and Euclid-WL only we obtain here 0.7%.

Fig. 2. One-dimensional posterior distributions, and 68% and 95% confidence level marginalised contours for the dark energy parameters (w0 and wa) and the parameters h and σ8. The left panel refers to the ℓmax = 1500 case, while the right panel uses measurements deeper into the non-linear regime, with ℓmax = 5000. The mock data of Euclid cosmic shear assume Halofit+PKequal non-linear corrections as reference, while the parameter estimation is performed with either HMCode (blue), or Halofit (orange). Black dashed lines mark the fiducial model.

As expected, the biases are larger for HMCode compared to Halofit without PKequal. We also find that for Halofit, increasing the range from ℓ_max = 1500 to ℓ_max = 5000 does not increase the bias significantly, whereas the bias strongly depends on the ℓ-range for HMCode, both in amplitude and in sign. This can be explained by looking at Fig. 1: at scales larger than a few h Mpc⁻¹, HMCode systematically over-predicts the power with respect to Halofit.

We find that in the Halofit case, the biases in the cosmological parameters approximately satisfy B(θ) ≲ B_thr when ℓ_max = 1500, while for HMCode the biases for almost all the parameters exceed this threshold. When setting ℓ_max = 5000, the parameters estimated in the HMCode case are all biased significantly more than the acceptable threshold (except for ω_b, 0), and now also the Halofit case exhibits biases larger than B_thr for h, Ω_m, 0 and σ₈.

In order to correct for the significant mismatch in the non-linear prescriptions, HMCode increases ω_c, 0, h, ln(10¹⁰A_s), and w_a, while at the same time the values for n_s and w₀ are decreased. This tweaking of parameters increases the amplitude of the linear matter power spectrum at scales 0.2 h Mpc⁻¹ ≲ k ≲ 2 h Mpc⁻¹, where HMCode has a lack of power with respect to Halofit+PKequal (see Fig. 1). As the scales around 0.2 h Mpc⁻¹ are those that mainly contribute to the estimate of σ₈, this explains the large bias observed for this (derived) parameter, B(σ₈) ∼ 5.

Overall, the Δχ² ≲ 1 indicates that replacing Halofit + PKequal with Halofit-only does not have a strong impact on the results as it is well within the range of the statistical uncertainties¹¹. On the other hand, using HMCode leads to a signficantly higher Δχ², highlighting how the difference between the two non-linear prescriptions cannot be fully compensated by modifying the background quantities and the linear growth.

It is worth to noting that for both Halofit and HMCode the parameters that are most significantly biased are H₀ and σ₈. These are the parameters that currently show tension between high- and low-redshift measurements (e.g., Riess et al. 2019; Hildebrandt 2020; Spurio Mancini et al. 2019). Our results imply that the Euclid cosmic shear measurements have the statistical power to resolve this, but only if we can accurately model the non-linear scales.

6. Impact of baryons

Up to this point, we have limited our study to the impact of changing the recipe that is used to compute the non-linear evolution of cold dark matter perturbations. On the small scales of interest, however, baryons collapse into the dark matter haloes to form stars, or are heated up, or even expelled into the intergalactic medium. These processes modify the matter distribution, and it is therefore important to account for baryonic physics when computing the matter power spectrum P_δδ(k, z) (e.g., van Daalen et al. 2011; Casarini et al. 2012; Castro et al. 2018; Debackere et al. 2020). This can be done by multiplying the non-linear power spectrum, P_δδ(k, z) – computed using one of the non-linear prescriptions discussed above – with ℬ(k,z), a ‘baryon correction model’ (BCM) that captures the baryonic effects (e.g., Semboloni et al. 2011), so that

$$\begin{array}{c}\hfill P_{\mathrm{c}+\mathrm{b}}(k,z)=P_{\delta \delta }(k,z)\mathcal{B}(k,z),\end{array}$$(21)

where P_c + b(k, z) is the corrected power spectrum. The function ℬ(k,z) can be estimated by fitting Eq. (21) to power spectra obtained from hydrodynamical simulations that include baryons.

The challenge is that baryonic effects cannot (yet) be incorporated into cosmological simulations from first principles. The different implementations that have been used, not surprisingly, lead to a variety of possible BCM prescriptions. Here we consider three recent proposals.

The first one presented in Harnois-Déraps et al. (2015, HD15 hereafter), is based on three scenarios of the OverWhelmingly Large hydrodynamical simulations (Schaye et al. 2010). These were used to calibrate the power spectra for z < 1.5. It is able to reproduce the simulated results with an accuracy better than 2% for scales k < 1 h Mpc⁻¹. The functional form of ℬ(k,z) is given by

$$\begin{array}{cc}& \mathcal{B}(k,z)=1-A_{\mathrm{HD}15}(z)exp\left\{{[B_{\mathrm{HD}15}(z)x(k)-C_{\mathrm{HD}15}(z)]}^3\right\}\hfill \\ \hfill & +D_{\mathrm{HD}15}(z)x(k)exp\left[E_{\mathrm{HD}15}(z)x(k)\right],\hfill \end{array}$$(22)

with x(k) ≡ log₁₀(k/[h Mpc⁻¹]) and X_HD15(z) are polynomial functions of redshift given in Harnois-Déraps et al. (2015).

As a second model, we consider the results obtained by Schneider & Teyssier (2015, ST15 hereafter), who accounted for the effects of baryons following a different approach. They start from a suite of DM only N-body simulations and modify the density field in such a way that it mimics the effects of a particular feedback recipe. They achieve this by explicitly modelling the main constituents of the haloes, which are dark matter, hot gas in hydrostatic equilibrium, ejected gas and stars. The model parameters are set to resemble SZ and X - ray observations. The resulting modifications to the power spectrum are shown to be well reproduced by defining ℬ(k,z) as

$$\begin{array}{c}\hfill \mathcal{B}(k,z)=\frac{1+{(k/k_{\mathrm{s}})}^2}{{[1+k/k_{\mathrm{g}}(z)]}^3}\mathcal{G}(z)+[1+{(k/k_{\mathrm{s}})}^2][1-\mathcal{G}(z)],\end{array}$$(23)

with k_g(z) and 𝒢(z) auxiliary functions provided in Schneider & Teyssier (2015). We set the model parameters to the following fiducial values

$$\begin{array}{c}\hfill \{k_{\mathrm{s}},log{M}_{\mathrm{c}},z_{\mathrm{b}},{\eta }_{\mathrm{b}}\}=\{67h{\mathrm{Mpc}}^{-1},13.8,2.3,0.17\}.\end{array}$$(24)

We note that these are different from those of Schneider & Teyssier (2015) since we have updated them to the best fitting values obtained using the more recent Horizon-AGN simulations (as done in Chisari et al. 2018).

Chisari et al. (2018, hereafter Ch18) found that the ST15 model performs well at low redshift, but its accuracy degrades for larger z. Fitting the Horizon-AGN simulations, they therefore proposed the third model we will consider here, with

$$\begin{array}{c}\hfill \mathcal{B}(k,z)=\frac{{[1+k/k_{\mathrm{s}}(z)]}^2}{{[1+k/k_{\mathrm{s}}(z)]}^{1.39}},\end{array}$$(25)

where k_s is no longer a constant, but a function of z instead. The detailed form is given in Chisari et al. (2018).

We can now quantify the impact of the choice of BCM on the FoM by comparing it to the results for the dark-matter-only forecasts. To this end, we use Halofit as our benchmark model to compute P_δδ(k, z), which is consistent with what is done in the quoted papers. Our results are presented in Table 5.

The more recent Ch18 model yields the smallest change in the FoM, but the differences are never larger than ∼15%, even in the scenario with the largest change, that is the HD15 model with ℓ_max = 5000. This is the consequence of two opposite effects that partially cancel each other out.

On the one hand, at k ∼ 3–13 h Mpc⁻¹, gas ejection due to AGN feedback suppresses the power spectrum, while for larger k the effect of stars is to increase it. These very small scales (k > 15 h Mpc⁻¹) are, however, weighted down by the lensing kernel so that the overall effect is to reduce the signal in the C^ϵϵ(ℓ), which tends to reduce the FoM. However, reducing C^ϵϵ(ℓ) also decreases the Gaussian covariance that is used to estimate the uncertainty in the WL signal as the baryonic effects also change the power spectra used to compute the covariance matrix. As a result, the inverse covariance boosts the FoM. By moving the FoM in opposite directions, these two effects roughly compensate each other, so that the choice of the BCM model impacts the constraining power of the WL signal on dark energy parameters only marginally. Qualitatively similar results are obtained when we compare the marginalised uncertainties for individual parameters. In general, we find that the parameter bounds are less affected than the FoM, in accordance with the argument given above.

Although the detailed implementation of the BCM does not affect the constraining power much, it is nonetheless necessary to investigate whether an incorrect choice biases the cosmological parameter estimates. To this end, we adopt a procedure similar to that of Sect. 5: We generate three simulated Euclid WL datasets assuming that the BCM describing the true effect of baryons is CH18, ST15, or HD15, respectively. We then analyse these datasets, together with simulated Planck CMB data, with theoretical predictions that neglect baryonic effects.

As done in the previous sections, we perform the MCMC analysis for both ℓ_max values, 1500 (pessimitic) and 5000 (optimistic), which allows us to assess the relevance of the BCM for high multipoles. Our results are reported in Table 6 and shown in Fig. 3. For ℓ_max = 1500 we notice that the bias exceeds the threshold B_thr for all parameters except for ω_b, 0 and ω_c, 0 in the CH18, while in ST15 also the bias on ln(10¹⁰A_s) is lower than B_thr. HD15 instead only has this latter parameter within our acceptable threshold, a result confirmed by the highest Δχ² among the three cases. The most significantly biased parameters in all three analysis are h, w₀ and w_a. Overall, however, all three cases produce very similar results, as can be seen in the left panel of Fig. 3.

Fig. 3. One-dimensional posterior distributions, and 68% and 95% marginalised joint two-parameter contours for w0 and wa, and the parameters h, ns and σ8 from the MCMC analysis. The results are obtained by neglecting baryon effects when fitting mock datasets created without baryonic effects (green) and with baryonic effects (orange for CH18, blue for ST15 and purple for Hd15). The left panel refers to the ℓmax = 1500 case, while the right panel goes deeper into the non-linear regime, with ℓmax = 5000.

For ℓ_max = 5000, we find that the biases are very large when BCM effects are neglected; B(θ) > B_thr for all parameters, with B ≳ 5 for the dark energy parameters w₀ and w_a, B ≳ 3 for h, and B ≳ 4.5 for n_s. As expected, the biases in the power spectrum amplitude, the baryon density, and the cold dark matter density are the less significant because these are all well constrained by the Planck measurements.

Ignoring BCM effects could lead to a false detection of a time-varying dark energy equation of state. Moreover, with the current tension between H₀ measurements between CMB and late-time probes, an unbiased measurement of H₀ will also be crucial. Our results confirm earlier work (e.g., Semboloni et al. 2011) that correctly modelling the impact of baryonic feedback on the power spectrum is essential for the analysis of Euclid data.

7. Conclusions

Forthcoming surveys of the large-scale structure will deliver datasets of exquisite quality that will allow us to pursue what is usually referred to as precision cosmology. To interpret correctly these measurements, significant improvements in the underlying theoretical predictions are called for: We need to ensure that errors in the modelling of density fluctuations on non-linear scales do not introduce biases in the inferred cosmological parameters that are larger than the expected statistical uncertainties. This is a particular concern for cosmic shear tomography, given that one has to integrate the matter power spectrum P_δδ(k, z) deep into the non-linear regime, where baryonic physics complicates matters even further. Motivated by these concerns, we have investigated how different popular prescriptions that account for these effects influence both the accuracy and the precision with which Euclid can infer cosmological parameters using cosmic shear alone.

We used Fisher matrix forecasts to quantify the impact of three different non-linear recipes on the dark energy FoM. The recipes that we considered are the revised implementation of Halofit, Halofit+PKequal, and the HMCode prescription. These differ significantly from one another when the cosmological parameters are left free to take values other than the fiducial ones. As a consequence, the derivatives of P_δδ(k, z) that enter the determination of the Fisher matrix are changed, leading to quite discrepant FoM values. In particular, we find that the Halofit case provides the higher FoM because of the different role played by w_a in the non-linear corrections. Although we have explicitly considered the case of Euclid, this result is generic for cosmic shear tomography analyses, although the size of this effect will depend on the details of the survey of interest. Hence our findings highlight the importance of choosing the most reliable non-linear model, in order to compute realistic estimates of the expected performance of a particular WL survey.

While it is important to quantify the precision, that is how tight the constraints will be, it is perhaps even more important to establish the accuracy of the results: We need to be confident that an incorrect choice of theoretical ingredients does not introduce an unacceptably large bias, that is a deviation from the (unknown) true value. Whether or not the bias is too large also depends on the precision with which that parameter can be measured. We therefore define B = |θ − θ^fid|/σ, and adopt a theshold of B < B_thr ≈ 0.1.

To study the accuracy with which cosmological parameters can be determined, we created mock data with a given prescription for non-linearities and/or baryon physics, and fitting these with a different model. This allowed us to address the issues of the choice of the non-linear recipe and the baryon correction model separately. To examine the impact of the recipe used to compute the power spectra on non-linear scales, we created a mock dataset using Halofit+PKequal that comprises Euclid cosmic shear and Planck CMB data. We emphasise that the choice of Halofit+PKequal is arbitrary as we do not know which of the non-linear corrections better describes the true small-scale evolution. We fitted these with Halofit or HMCode.

We find that B ≲ B_thr if Halofit is used for ℓ_max = 1500, while for ℓ_max = 5000, some cosmological parameters, namely h, Ω_m, 0 and σ₈, exceed the threshold. This is not surprising given the similarities of the two models when one only looks at P_δδ rather than at its derivatives. In contrast, comparing the theoretical spectra obtained using HMCode based corrections with the assumed Halofit+PKequal fiducial leads to strong biases, with almost all parameters already biased by more than B_thr if we restrict the analysis to ℓ_max = 1500. Including very non-linear scales (ℓ_max = 5000), we find that B > 1 for all parameters, except for ω_b, 0, which is actually tightly constrained by Planck. In particular, the estimate of w₀ shifts towards its ΛCDM value even if the mock data were created using {w₀,  w_a}={−0.9,  0.1}. What is even more interesting is that the most biased parameters are h (B = 3.02) and σ₈ (B = 4.62); the values of both of these are currently debated. The sensitivity of these parameters to the adopted prescription for the non-linear power spectrum highlights the need for further improvements, which may already be needed to correctly interpret current data.

It is also essential that the changes to the power spectrum caused by baryon physics have to be taken into account. For a fixed non-linear recipe and baryon correction model prescription, severe biases are found when fitting the mock data with the right non-linear correction, but not accounting for the presence of baryons, in line with earlier work (e.g., Semboloni et al. 2011). In the most constraining setting (ℓ_max = 5000), for all the three cases we considered, we find significant biases for all the cosmological parameters, except for {ω_b, 0, ln(10¹⁰A_s)}, which are actually constrained by the Planck data rather than by cosmic shear. While this work was near to completion, Schneider et al. (2020) presented a similar analysis, but their method to account for baryons is very different from the one we have adopted here. They use instead a model for the baryonification of dark matter only simulations to determine the matter power spectrum (Schneider et al. 2019). Notwithstanding these differences, which make a straightforward comparison impossible, their conclusions are in agreement with what we found here.

As a final remark, we remind the reader that our results refer to the case where cosmic shear is used as the only probe. This is, however, only part of the information that future surveys will provide. Indeed, the same data used to do cosmic shear tomography (WL) can and will be used to compute the photometric galaxy clustering (GCph) and to cross correlate the shear and density fields (XC). As was shown in EC19 and further investigated in Tutusaus et al. (2020), it is the joint use of WL+GCph+XC that is needed to achieve ∼{1, 10}% errors on {w₀, w_a}, rather than any single probe by itself. It is therefore worth considering whether, and to which extent, the results we have obtained here change if all the three probes are considered.

For instance, EC19 limited the GCph and XC to smaller multipoles compared to what was used for WL, reducing their sensitivity to the small-scale corrections. This diminishes the impact of errors in the prediction of P_δδ(k, z) in the large k regime, but at the expense of larger statistical uncertainties. This motivates extending our work to quantify the impact of modelling the small-scale power spectra for the joint probes. Such a study, which is beyond the scope of our initial exploration, would provide guidance on how to proceed in order to exploit the high-quality data that stage IV surveys will provide.

Acknowledgments

D. F. M. thanks the Research Council of Norway for their support, and the resources provided by UNINETT Sigma2 – the National Infrastructure for High Performance Computing and Data Storage in Norway. M. M. has received the support of a fellowship from “la Caixa” Foundation (ID 100010434), with fellowship code LCF/BQ/PI19/11690015, and the support of the Spanish Agencia Estatal de Investigacion through the grant “IFT Centro de Excelencia Severo Ochoa SEV-2016-0597”. This paper is based upon work from the COST action CA15117 (CANTATA), supported by COST (European Cooperation in Science and Technology). M. A. acknowledges the computing support of INFN-CNAF. Stefano Camera acknowledges support from the Italian Ministry of Education, University and Research (MIUR) through Rita Levi Montalcini project ‘PROMETHEUS – Probing and Relating Observables with Multi-wavelength Experiments To Help Enlightening the Universe’s Structure’, and the ‘Departments of Excellence 2018–2022’ Grant awarded by MIUR (L. 232/2016). The work of SC is supported by the Belgian Fund for Research FNRS-F.R.S. I. T. acknowledges support from the Spanish Ministry of Science, Innovation and Universities through grant ESP2017-89838-C3-1-R. I. T. and T. K. acknowledge funding from the H2020 programme of the European Commission through grant 776247. AP is a UK Research and Innovation Future Leaders Fellow, grant MR/S016066/1. 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 Academy of Finland, the Agenzia Spaziale Italiana, the Belgian Science Policy, the Canadian Euclid Consortium, the Centre National d’Etudes Spatiales, the Deutsches Zentrum für Luft- und Raumfahrt, the Danish Space Research Institute, the Fundação para a Ciência e a Tecnologia, the Ministerio de Economia y Competitividad, the National Aeronautics and Space Administration, the Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, 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 (http://www.euclid-ec.org).

References

Appendix A: Comparison between scales and multipole cuts

In the analyses performed in this paper, we considered two cases, with different cuts in multipoles used; the optimistic ℓ_max = 5000 cut represents the situation in which data for all scales are included in the analysis, while the ℓ_max = 1500 mimics the removal of non-linear scales that could be performed in the data analysis in order to reduce the impact of small-scale modelling.

However, our approach can be seen as a first approximation as a constant multipole cut corresponds to different scales for different redshifts. In this subsection we investigate a less simplistic approach, implementing a k_max scale cut, rather than a multipole one, with the purpose of removing part of the non-linear scales.

We assume k_max = 0.25 h Mpc⁻¹ and, using the Limber approximation considered to express the C^ϵϵ(ℓ) in Eq. (15), we convert this cut in a maximum multipole for each redshift bin:

$$\begin{array}{c}\hfill {\ell }_{\max }(z_i)=k_{\max }r(z_i)-\frac{1}{2},\end{array}$$(A.1)

where z_i is the mean redshift of each redshift bin.

We apply this approach to our MCMC analysis, and compare the results we obtain with those of the ℓ_max = 1500 case. In Fig. A.1 we show the comparison in the results for the analysis performed using Halofit (left panel) and HMCode (right panel). In both cases we find that the k_max analysis leads to broader constraints with respect to the ℓ_max case; this is due to the fact that the chosen k_max translates into a much more agressive cut in multipoles, especially at low redshift. Therefore, a significative amount of information is lost with respect to the ℓ_max = 1500 case.

Fig. A.1. One-dimensional posterior distributions, and 68% and 95% confidence level marginalised contours for the dark energy parameters (w0 and wa) and the parameters h and σ8. Here we compare results obtained with a multipole cut at ℓmax = 1500 (green contours) with those related to the analysis with the scales cut kmax = 0.25 h Mpc−1 (red contours). The mock data of Euclid cosmic shear assume Halofit+PKequal non-linear corrections as reference, while the parameter estimation is performed with either HMCode (right panel), or Halofit (left panel). Black dashed lines mark the fiducial model.

Concerning the bias found on cosmological parameters, B(θ) does not change significantly in the Halofit analysis, while for HMCode the biases are slightly enhanced with respect to the ℓ_max = 1500 case, with w₀,  w_a,  h and σ₈ now reaching B ≈ 1. This apparent increase however is mostly due to the loss of constraining power when the scale cut is implemented; as it can be seen in Fig. A.1, the marginalised posterior distributions for the cosmological parameters now exhibit non-Gaussian features while Eq. (20) implicitly assumes Gaussian distributions. In Fig. A.1 it is shown how the scale cut case produces contours closer to the expected fiducial values, a result supported also by the change in Δχ² moving from the ℓ_max to the scale cut, which changes from 32 to 14 in the HMCode case.

We performed the same analysis at the Fisher matrix level finding a still significant dependence of the results on the adopted non-linear recipe. We indeed get FoM = (2.93, 3.59, 1.59) for Halofit, HMCode, and Halofit + PKequal, respectively. While the severe decrease in the FoM is expected given that we are removing a large part of the data, it is somewhat surprising to still find such a variety of values. This is, however, a consequence of the integrated nature of the WL $C_{\mathit{ij}}^{ϵϵ}(\ell )$. To understand what is going on, we focus on the case i = j = 5 giving ℓ_max ≃ 500. Because of the photo - z broadening of the lensing kernel, the integral giving C₅₅(ℓ_max) gets contributions from the redshift range (0.1, 2.6). Over this range, the argument k_ℓ(z) of the matter power spectrum P_δδ(k, z) feeding the integral is larger than k_max for z < 0.84 so that which non-linear recipe is adopted still matters. Such an argument can be repeated for all the bins combinations and the multipoles thus explaining why the FoM is still dependent on the non-linear recipe even with these very conservative scale cuts, a result in agreement with what was discussed in Taylor et al. (2018a). Therefore, in order to remove completely the dependence on the non-linear description from the analysis, different approaches are needed, for example using band powers rather than a C^ϵϵ(ℓ) analysis (Joachimi et al. 2021).

Appendix B: MCMC results validation

The results of this paper have been obtained using both Fisher matrix and MCMC codes. The Fisher analysis relies on one of the codes used in EC19, which has passed through a careful validation procedure that included intensive comparisons between different Fisher matrix codes. Our MCMC analysis relies on a public MontePython likelihood for Euclid WL (Brinckmann & Lesgourgues 2019), adapted to the Euclid specifications of EC19. It has been further modified to include different models of baryonic feedback effects. This likelihood code was first used in Sprenger et al. (2019). In contrast to the Fisher code, it has not been validated against other codes. In this Appendix, we therefore present a comparison between validated Fisher forecasts and the MCMC ones.

In the Euclid-only case, our analysis reveals some deviations that are attributed to the intrinsic limitation of any Fisher analysis, due to the non-Gaussianity of the posterior distribution, as well as some important parameter degeneracies. Nevertheless, the impact on the forecasts becomes negligible when Planck constraints are included in the analysis. In this case, we find that the forecasts on cosmological parameters obtained with the Fisher and MCMC methods agree very well. This therefore validates our MCMC approach against the Fisher codes used in EC19.

The Euclid forecasts for WL in EC19 are accompanied by a series of public Fisher matrices, corresponding to different setups and cosmologies. In order to validate our MontePython likelihood and MCMC analysis, we compared forecasts obtained for both the pessimistic and optimistic setups, and using the same cosmological parameters (in particular, Ω_b, 0 and Ω_m, 0 instead of ω_b and ω_c). We performed these comparisons for Euclid only, and in combination with Planck. In the latter case, we used the mock Planck likelihood available in Montepython that accurately reproduces the Planck limits on cosmological parameters. We construct a covariance matrix from the MCMC chains in the Planck-only case. Its inverse provides a Fisher matrix that can be added to the validated Euclid Fisher matrices. We checked that the Planck-only case constrains the standard cosmological parameters well, with close-to-Gaussian two-dimensional posterior distributions. This is a good indication that one can safely use the Planck covariance matrix for the Fisher analysis. In contrast to the standard cosmological parameters, most of the constraining power for the dark energy parameters w₀ and w_a comes from Euclid, not Planck. As a consequence, the corresponding entries in the Planck Fisher matrix are not relevant and do not significantly impact the Euclid + Planck forecasts for these parameters.

For Euclid only, the marginalised two-dimensional posterior distributions and the Fisher contours obtained in the case ℓ_max = 1500 for five varying cosmological parameters (Ω_m, 0, Ω_b, 0, h, n_s and σ₈), with DE and IA parameters fixed to their fiducial values, are shown in Fig. B.1. We find that the directions and widths of all the Fisher ellipses are well recovered by the MCMC approach. However, we also find two significant differences. First, the MCMC contours in the plane (Ω_b, 0, h) display a banana shape, leading to a more stringent constraint on Ω_b, 0 compared to the Fisher method. In turn, this affects two-dimensional contours between Ω_b, 0 and other parameters. Second, the marginalised posterior distribution for h has a significant non-Gaussian shape, falling more rapidly at lower values than at higher values. Consequently, some of the contours are less extended on one side compared to the Fisher ellipses. This also arises from the degeneracy between h and Ω_b, 0, with these parameters being poorly constrained with Euclid WL only. These features cannot be recovered by the Fisher analysis that relies on the assumption that the posterior distribution is Gaussian, which is not the case for strong variations of h and Ω_b, 0. Nonetheless, these results show that both the MCMC and Fisher methods provide consistent results, even if the agreement between them is limited by the intrinsic limitation of the Fisher analysis. Adding dark energy parameters in the analysis further degrades the agreement between the Fisher and MCMC forecasts, for similar reasons. Here, we have only shown the case in which w₀ and w_a remain fixed, for a better illustration of the effect of parameter degeneracies on the two method comparison.

Fig. B.1. Fisher matrix vs MCMC forecasts, for Euclid WL only with ℓmax = 1500. In red, one-dimensional marginalised posterior distributions, with two-dimensional 68% and 95% confidence level marginalised contours, for varying cosmological parameters and fixed DE and IA parameters. Black curves represent the corresponding Fisher ellipses, from the validated Fisher matrix of EC19.

For Planck + Euclid WL, we find a very good agreement between Fisher and MCMC forecasts, even when the DE and IA parameters are varying. This is shown in Fig. B.2 representing the MCMC marginalised two-dimensional posterior distributions, which fit well to all the Fisher ellipses. We find that the differences in the forecasts are less than 10% for all the cosmological parameters. This therefore validates the approach we have used throughout the paper. It also illustrates the importance of adding CMB data to break parameter degeneracies. Compared to the Euclid WL only case, the differences between the Fisher and MCMC methods are suppressed. The few remaining differences can be either due to a limited convergence of the MCMC chains, to uncertainties induced by the binning method, or to a slightly non-Gaussian likelihood function for Planck-only that induces small differences between MCMC posteriors and the corresponding Planck Fisher ellipses extracted from the covariance matrix of the MCMC chains. Similar results have been obtained for the optimistic setup with ℓ_max = 5000.

Fig. B.2. Fisher matrix vs MCMC forecasts, for Planck+Euclid WL, with ℓmax = 1500. In red, one-dimensional marginalised posterior distributions, with two-dimensional 68% and 95% confidence level marginalised contours, for varying cosmological parameters including w0 and wa and IA parameters. Black curves represent the corresponding Fisher ellipses, from the validated Fisher matrix of EC19 combined with a Planck-only Fisher matrix.

All Tables

All Figures

Fig. 1. Relative difference between the various predicted matter power spectra and the Halofit prediction at z = 0 for different cosmologies. Top-left panel: using the mean values of the ΛCDM parameters given by Planck 2018 (Planck Collaboration VI 2020). Top-right panel: differences when the ΛCDM parameters are drawn from a five-dimensional Gaussian distribution with the diagonal given by five times the 1σ uncertainty quoted by Planck TTTEEE+lowE+lensing+BAO 2018. Bottom-left panel: differences when also the constant dark energy equation of state is drawn from a Gaussian distribution with a dispersion of 0.3. Bottom-right panel: difference for evolving dark energy models where, alongside w0, also wa is Gaussian distributed with a dispersion of 1.0. In the bottom panelsw0 and wa (if present) are always chosen in such a way that w(z) <  − 1/3. The feature visible at k ≈ 0.005 h Mpc−1 corresponds the scale at which the non-linear corrections are turned on in CAMB.

In the text

Fig. 2. One-dimensional posterior distributions, and 68% and 95% confidence level marginalised contours for the dark energy parameters (w0 and wa) and the parameters h and σ8. The left panel refers to the ℓmax = 1500 case, while the right panel uses measurements deeper into the non-linear regime, with ℓmax = 5000. The mock data of Euclid cosmic shear assume Halofit+PKequal non-linear corrections as reference, while the parameter estimation is performed with either HMCode (blue), or Halofit (orange). Black dashed lines mark the fiducial model.

In the text

Fig. 3. One-dimensional posterior distributions, and 68% and 95% marginalised joint two-parameter contours for w0 and wa, and the parameters h, ns and σ8 from the MCMC analysis. The results are obtained by neglecting baryon effects when fitting mock datasets created without baryonic effects (green) and with baryonic effects (orange for CH18, blue for ST15 and purple for Hd15). The left panel refers to the ℓmax = 1500 case, while the right panel goes deeper into the non-linear regime, with ℓmax = 5000.

In the text

Fig. A.1. One-dimensional posterior distributions, and 68% and 95% confidence level marginalised contours for the dark energy parameters (w0 and wa) and the parameters h and σ8. Here we compare results obtained with a multipole cut at ℓmax = 1500 (green contours) with those related to the analysis with the scales cut kmax = 0.25 h Mpc−1 (red contours). The mock data of Euclid cosmic shear assume Halofit+PKequal non-linear corrections as reference, while the parameter estimation is performed with either HMCode (right panel), or Halofit (left panel). Black dashed lines mark the fiducial model.

In the text

	Fig. B.1. Fisher matrix vs MCMC forecasts, for Euclid WL only with ℓmax = 1500. In red, one-dimensional marginalised posterior distributions, with two-dimensional 68% and 95% confidence level marginalised contours, for varying cosmological parameters and fixed DE and IA parameters. Black curves represent the corresponding Fisher ellipses, from the validated Fisher matrix of EC19.
In the text

Fig. B.2. Fisher matrix vs MCMC forecasts, for Planck+Euclid WL, with ℓmax = 1500. In red, one-dimensional marginalised posterior distributions, with two-dimensional 68% and 95% confidence level marginalised contours, for varying cosmological parameters including w0 and wa and IA parameters. Black curves represent the corresponding Fisher ellipses, from the validated Fisher matrix of EC19 combined with a Planck-only Fisher matrix.

In the text