© Euclid Collaboration 2020

1. Introduction

Euclid¹ will explore the expansion history of the Universe and the evolution of large-scale cosmic structures by measuring shapes and redshifts of galaxies, covering 15 000 deg² of the sky, up to redshifts of about z = 2. It will be able to measure up to 30 million (Pozzetti et al. 2016) spectroscopic redshifts, which can be used for galaxy clustering measurements and 2 billion photometric galaxy images, which can be used for weak lensing observations (for more details, see Amendola et al. 2018 and Laureijs et al. 2011). The Euclid telescope is a 1.2 m three-mirror anistigmat and it has two instruments that observe in optical and near-infrared wavelengths. In the optical, high-resolution images will be observed through a broadband filter with a wavelength range from 500 to 800 nm (the VIS band) and a pixel resolution of 0.1 arcsec that will result in a galaxy catalogue complete down to a magnitude of 24.5 AB in the VIS band. In the near-infrared, there will be imaging in the Y, J, and H bands with a pixel resolution of 0.3 arcsec and a grism (slitless) spectrograph with one “blue” grism (920−1250 nm) and three “red” grisms (1250−1850 nm; in three different orientations). This will result in a galaxy catalogue complete to magnitude 24 in the Y, J, and H bands. The optical and near-infrared instruments share a common field-of-view of 0.53 deg².

This paper is motivated by the challenging need to have reliable cosmological forecasts for Euclid. This is required for the verification of Euclid’s performance before launch in order to assess the impact of the real design decisions on the final results (for an updated technical description of the mission design, see Racca et al. 2016). It is also necessary that the different cosmological probes in Euclid use a consistent and well-defined framework for forecasting and that the tools used for such forecasts are rigorously validated and verified. These forecasts will then provide the reference for the performance of Euclid to the scientific community. Forecasting in this context refers to the question of how well Euclid will perform when using the survey data to distinguish the standard cosmological model from simple alternative dark energy scenarios, given its current, up-to-date specifications. This paper represents the outcome of an intense activity of comparison among different and independent forecasting codes, from various science working groups in the Euclid Collaboration. This work was conducted by a group known as the inter-science taskforce for forecasting (IST:F) and it presents the results of the code comparison for forecasts based on Fisher-matrix analyses.

Euclid forecasts were previously made in the Definition Study Report (hereafter “Red Book”, Laureijs et al. 2011), using Fisher matrix predictions, similar to those presented in this paper. The description of the Red Book forecasts was not detailed enough due to the length constraints of that document. In contrast, in this paper we specify in detail the methodology and the procedures followed in order to validate the results and verify their robustness. The Red Book did specify the Euclid instrument, telescope, and survey specifications in detail, and these have been used by the community to create a suite of predictions for the expected uncertainties on cosmological measurements for a range of models. Many forecasts have been made for several scenarios beyond the parameters of the standard ΛCDM model (see in particular Amendola et al. 2018): some, for example Hamann et al. (2012), have demonstrated the expected constraints on neutrino masses; others have examined the prospect of constraining the properties of different dark matter candidates; and most commonly the different parameterisation of dark energy and modified gravity properties have been considered (e.g. di Porto et al. 2012; Majerotto et al. 2012; Wang 2012; de Putter et al. 2013; Casas et al. 2017, among many others).

In this paper we focus on the primary cosmological probes in Euclid: weak lensing (WL); photometric galaxy clustering (GC_ph); spectroscopic galaxy clustering (GC_s); their combination; and the addition of the data vector of cross-correlation (XC) between GC_ph and WL. For each of the probes we define a recipe, that details the methodology used to make forecasts, which should serve as a useful documentation to enable the scientific community to be able to implement and reproduce the same methodology. The comparison presented here has involved several different codes. Therefore, rather than providing a single open source code, in a specific programming language, we grant a validation stamp for all the codes that took part in the code comparison and satisfied the requirements. These are listed in Sect. 4.1.

Since the Red Book the Euclid baseline survey and instrument specifications have been updated according to the progress in our knowledge of astrophysics, and in the evolution of the design of the telescope and instruments. In this paper we also update the Red Book forecasts by including some of these changes, which directly impact the statistical precision with which Euclid can constrain cosmological parameters. Such effects include changes to the survey specifications, changes to the instrument that impact the populations of galaxies probed, and changes in the understanding of astrophysical systematic effects. Some of these changes, although not all, have an impact on Euclid’s ability to constrain cosmological parameters (Majerotto et al. 2012; Wang et al. 2013).

For the GC observables a significant change since the Red Book has been an increased understanding of our target galaxy sample (Pozzetti et al. 2016), which has influenced the instrument design and the consequent changes in the optimisation of the survey strategy (Markovič et al. 2017). One of the difficulties in predicting the performance of Euclid’s spectroscopic survey accuracy has always been the poor knowledge of the number density and evolution of the galaxy population that will be detected as Hα emitters, and whose redshifts will be measured by the survey. More recent observations imply lower number densities at z >  1 than originally assumed based on information available at the time of the Red Book. This has resulted in a re-assessment of performance and the resulting cosmological parameter constraints. The main consequence of a reduced number density n(z) of observed galaxies is an increase in the shot noise and a consequent increase of the parameter uncertainties. The science reach of an experiment depends heavily on the “effective” volume covered, which in the case of GC studies is a function combining the cosmological volume with the product of the mean density of the target population and the amplitude of the clustering power spectrum; this is denoted nP. A further limitation of the baseline forecasts in the Red Book was that the observing strategy consisted of just two rotations, that is to say observing position angles (PA – also called dithers) for each of the two red and blue channels. This had the disadvantage that when a single PA observation is lost, the spectrum corresponding to the affected grism (covering a specific blue or red redshift range) could not be reliably cleaned of possible contaminating signals from adjacent spectra. In this paper we update the GC forecasts to reflect the changes in the expected PAs and the n(z).

For WL one of the main difficulties in making cosmological forecasts has been the ability to accurately model the intrinsic alignments of galaxies – a local orientation of galaxies that acts to mimic the cosmological lensing signal. In the Red Book a non-parametric model was chosen in redshift and in scale, and prior information on the model’s nuisance parameters included the expected performance of spectroscopic galaxy-galaxy lensing results at the time of Euclid’s launch. Since the time of Red Book there has been an increased attention in this area and several physical models have been proposed (Joachimi et al. 2015; Kirk et al. 2015; Kiessling et al. 2015) that model intrinsic alignments in a more realistic manner. A second area of attention has been in the appreciation of the impact of the small-scale (high k-mode) matter power spectrum on WL two-point statistics (e.g. Taylor et al. 2018a; Copeland et al. 2018). This has led to improved models of the impact of baryonic feedback, neutrino mass, and non-linear clustering on small scales. In this paper we update the WL forecasts to reflect these improved models.

As well as Euclid, there are several other large-scale cosmological experiments and survey projects that have produced forecasts since the Red Book, often including Euclid forecasts for comparison. For example the Square Kilometer Array² (SKA) is expected to deliver a wealth of radio wavelength data. Its Phase 1 Mid array (SKA1-MID) is going to be commissioned on similar timescales to Euclid and is offering a unique opportunity for multi-wavelength synergies. As discussed in Kitching et al. (2015), a multitude of cross-correlation statistics will be provided by comparing Euclid and SKA clustering and weak lensing data, with the additional advantage of the cross-correlations being less affected by systematic effects that are relevant for one type of survey but not the other (for more recent reviews, see also Weltman et al. 2020; Bacon 2020). Cross-correlation forecasts for an SKA1-MID 21 cm intensity mapping survey and Euclid galaxy clustering have been performed in Fonseca et al. (2015) and Pourtsidou et al. (2017). Forecasts for the cross-correlations of shear maps between Euclid and Phase 2 of SKA-MID assuming 15 000 deg² of sky overlap have been performed in Harrison et al. (2016). Three large optical surveys that will join Euclid are the Dark Energy Spectroscopic Instrument³ (DESI), the Large Synoptic Survey Telescope⁴ (LSST), and the Wide Field Infrared Survey Telescope⁵ (WFIRST). Joint analyses of their data, combined with data coming from cosmic microwave background (CMB) missions can give new insights to a broad spectrum of science cases, ranging from galaxy formation and evolution to dark energy and the neutrino mass. Some of the possibilities of combining Euclid, LSST, and WFIRST have been described in Jain et al. (2015) and Rhodes et al. (2017); LSST and WFIRST have also made their own forecasts e.g. Chisari et al. (2019). Concentrating on the merits for large-scale structure measurements, we note that the complementarity of Euclid, LSST and WFIRST results in significant improvement in cosmological parameter constraints. Similarly, joint ventures with γ-ray experiments such as the Fermi satellite or the proposed e-ASTROGAM mission may lead to an improvement in our understanding of the particle nature of dark matter (Camera et al. 2013; De Angelis et al. 2018). In this paper we do not make forecasts for other surveys, but the results presented here are an update for the Euclid forecasts presented in previous inter-survey comparisons.

This paper is arranged as follows. In Sect. 2 we describe the cosmological models considered in our forecasts and the parameters that characterise them. In Sect. 3 we introduce the Fisher-matrix formalism to estimate errors on cosmological parameters and describe our methodology to apply it to the different cosmological probes in Euclid. In Sect. 4 we present the forecasting codes included in our analysis and describe in detail the code-comparison procedure that we implemented, as well as the level of agreement found for the different cases and parameter spaces considered. In Sect. 5 we present the final cosmological parameter forecasts for the different probes in Euclid, considered separately and in combination. Finally, we present our main conclusions in Sect. 6. In this paper we focus on cosmological forecasts that explicitly do no include systematic effects relating to instrument design and performance of data reduction algorithms (although we do include astrophysical systematic effects); the assessment of the impact of these on cosmological parameter estimation is subject to an exercise known as science performance verification (SPV) in the Euclid Consortium, and will be presented in a series of separate papers.

2. The cosmological context

In the following subsections we describe different cosmological models and their associated parameters. We start by introducing and discussing the most commonly used background cosmological quantities.

2.1. Background quantities

Applying Einstein’s field equations of general relativity to the metric line element⁶ of a homogeneous and isotropic universe allows us to derive the Friedmann equations, which describe the time evolution of the scale factor a(t) as a function of the curvature parameter K and the energy content of the universe, characterised by the total energy density ρ and pressure p. They read

(1)

(2)

Here the overdot denotes differentiation with respect to the cosmic time t, c is the speed of light, G the gravitational constant and K can assume negative, zero, or positive values for open, flat or closed spatial geometry respectively. A more convenient parameter than the time t to describe the evolution of the scale factor is the redshift, z = a₀/a − 1, where a₀ corresponds to the present-day value of the scale factor, normalised as a₀ = 1; we convert the t variable into a or z. The Hubble expansion rate, H(z), can be expressed as

(3)

with H₀ ≡ H(z = 0) being the Hubble parameter today, which is commonly written as

(4)

where h is the dimensionless Hubble parameter. The function E(z) will be later specified for each cosmological model we use. For a given value of H(z), there is a value of ρ that results in a spatially flat geometry (K = 0). That is the critical density

(5)

For a generic component labelled i, we define the density parameter Ω_i(z) ≡ ρ_i(z)/ρ_crit(z). Based on Eq. (1), we can also introduce an effective curvature density parameter Ω_K(z) = − Kc²/[a²(z)H²(z)]. With these definitions, Eq. (1) takes the form

(6)

where the sum is over all N species considered in the model. We use present-day values of the density parameters, which, unless specified otherwise, are indicated with a subscript 0.

Equations (1) and (2) can be combined into an energy conservation equation that specifies the relation between ρ (and p) and the scale factor. A solution of this equation requires to specify the properties of each energy component in the form of an equation of state, p = p(ρ); we specify the latter in terms of the equation of state parameter w ≡ p/ρc², which can be redshift-dependent. For the case in which the equation of state parameter is constant in time the energy conservation equation implies

(7)

Once the relations ρ_i(a) are known, these can be used in Eq. (1) to find a solution for a(t). We also note that −3[1 + w_i] always gives ∂ln ρ_i(a)/∂ln a.

The matter energy density at late times is mainly in the form of baryons and cold dark matter (CDM) particles, which are described by w_b = w_c = 0. The photon radiation density is characterised by w_γ = 1/3. A contribution from massive-neutrinos can be described by a varying equation of state parameter w_ν, which matches w_γ at early times and w_c when they become non-relativistic. For the purpose of galaxy clustering and weak lensing measurements, we can consider radiation density to be negligible, effectively setting Ω_γ, 0 = 0, and treat the massive-neutrinos as part of the total matter contribution, with Ω_m, 0 = Ω_c, 0 + Ω_b, 0 + Ω_ν, 0, since they are non-relativistic at the low redshifts relevant for these probes.

In the context of general relativity, cosmic acceleration requires a fluid, dubbed “dark energy” (DE), with an equation of state w_DE <  −1/3. The standard model of cosmology, commonly referred to as the ΛCDM model, assumes that this phenomenon is due to the presence of a cosmological constant, referred to as Λ, described by a constant equation of state w_Λ = −1, which, according to Eq. (7), corresponds to a time-independent energy density ρ_Λ. The ΛCDM model currently fits observations very well (see Planck Collaboration VI 2020, and references therein) but it suffers from several fundamental problems: the observed value of the cosmological constant is many orders of magnitude smaller than the theoretical predictions (the cosmological constant problem); the fact that the Λ and CDM densities are similar today while they have evolved very differently through time marks our epoch as a special time in the evolution of the Universe (the coincidence problem); and in general the model and its parameters cannot be predicted from physical principles.

A more general scenario for the component responsible for cosmic acceleration postulates a dynamical DE, with a redshift-dependent equation of state parameter w_DE(z). A commonly used and well-tested parameterisation of the time dependence is

(8)

where w₀ is the present (z = 0) value of the equation of state and w_a is a measure of its time variation. In this case, the evolution of the DE density obeys

(9)

Using Eqs. (7) and (9) in Eq. (1), the function E(z) defined in Eq. (3) becomes

(10)

with the current DE density, Ω_DE, 0, satisfying the relation Ω_DE, 0 = 1 − Ω_m, 0 − Ω_K, 0. The ΛCDM model can be recovered by setting w₀ = −1 and w_a = 0, in which case the function E(z) takes the form

(11)

2.2. Distance measurements

The comoving distance to an object at redshift z can be computed as

(12)

Although this quantity is not a direct observable, it is closely related to other distance definitions that are directly linked with cosmological observations. A distance that is relevant for our forecasts is the angular diameter distance, whose definition is based on the relation between the apparent angular size of an object and its true physical size in Euclidean space, and is related to the comoving distance by

(13)

Also relevant for our forecasts is the comoving volume of a region covering a solid angle Ω between two redshifts z_i and z_f, which is given by

(14)

for a spatially flat universe (K = 0), this becomes

(15)

These expressions allow us to compute the volume probed by Euclid within a given redshift interval.

2.3. Linear perturbations

The structure we see today on large scales grew from minute density fluctuations generated by a random process in the primordial Universe. The evolution of these fluctuations, for non-relativistic matter on sub-horizon scales, can be described by ideal fluid equations (Peebles 1980). Density fluctuations for a given component i are characterised by the density contrast

(16)

which quantifies the deviations of the density field ρ_i(x, z) around the mean spatial density over space, where x is a three-dimensional comoving coordinate at a redshift z. To describe these fluctuations statistically, it is convenient to work in Fourier space by decomposing δ into plane waves,

(17)

The power spectrum, P_i(k, z), for the generic component i is defined implicitly as

(18)

where δ_D is the Dirac delta function. Under the assumptions of statistical homogeneity and isotropy, the power spectrum can only depend on k = |k| and z.

The dimensionless primordial power spectrum of the curvature perturbation ζ generated by inflation is parameterised as a power law

(19)

where A_s is the amplitude of the primordial scalar perturbation, the scalar spectral index n_s measures the deviation (tilt) from scale invariance (n_s = 1), and k₀ is a pivot scale. The corresponding power spectrum defined in Eq. (18) in terms of density perturbations is related to the primordial one via the transfer function 𝒯_i through

(20)

for a generic component i. In the early Universe during radiation domination, curvature perturbations with comoving scales smaller than the horizon are suppressed, whereas super-horizon fluctuations remain unaffected, until they enter the horizon. In the matter dominated era, curvature perturbations on all scales remain constant. This implies that a characteristic scale corresponding to the epoch of matter-radiation equality is imprinted on the shape of the transfer function, and hence on the matter power spectrum.

The growth of density fluctuations obeys a second-order differential equation. At early enough times, when those fluctuations are still small, the fluid equations can be linearised. During matter domination, considering matter as a pressureless ideal fluid, the equation for the evolution of the density contrast becomes

(21)

In the ΛCDM scenario with no massive-neutrinos, this equation can be written in terms of the redshift as

(22)

where we neglect dark energy perturbations, the prime refers to the derivative with respect to z, and Ω_m(z) is given by

(23)

The solutions δ_m(k, z) of Eq. (22), at late times, are scale-independent, which motivates the introduction of the growth factor D(z) through

(24)

where z_i is an arbitrary reference redshift in the matter-dominated era. A useful quantity is the growth rate parameter, defined as

(25)

Using these definitions⁷, the growth rate satisfies a first-order differential equation

(26)

with initial condition f(z = z_i) = 1. Using Eq. (25) the solution for D(z) can be expressed in terms of f(z) by the integral

(27)

In ΛCDM the late-time matter growth is scale-independent, so that the transfer function 𝒯_m(k, z) can be split into a scale-dependent part T_m(k) (normalised so that T_m → 1 for k → 0) and the scale-independent growth factor D(z) introduced above. A convenient way to express the power spectrum defined in Eq. (20) for matter is

(28)

with the normalisation constant

(29)

where W_TH(x) = 3(sin x − x cos x)/x³ is the Fourier transform of the top-hat filter, and R₈ = 8 h⁻¹ Mpc. The pivot scale k₀ and the amplitude of the scalar mode A_s are absorbed into the normalisation, which is designed to give a desired value of σ₈, the root mean square (rms) of present-day linearly evolved density fluctuations in spheres of 8 h⁻¹ Mpc, which is given by

(30)

Hence, the generic power spectrum described in Eq. (20) relates to the matter power spectrum of Eq. (28) when the transfer function of species i is

(31)

Since we consider models beyond ΛCDM, we need to specify the behaviour of the dark energy perturbations. We use a minimally-coupled scalar field, dubbed quintessence (Wetterich 1988; Ratra & Peebles 1988), which is often considered as the standard “dynamical dark energy”. At the level of linear perturbations, this choice corresponds to a fluid with a sound speed equal to the speed of light (), and no anisotropic stress (σ_DE = 0). This implies that it is smooth deep inside the horizon, and does not develop significant fluctuations, in which case the same form of the growth and power spectrum equations given above are still valid. In addition, we want to allow for the possibility that w_DE(z) crosses w = −1 (which is not allowed in quintessence, but could happen in multi-field scenarios, Peirone et al. 2017), and to achieve this we use the Parameterised Post-Friedmann (PPF) prescription (Hu & Sawicki 2007; Fang et al. 2008; Planck Collaboration XIV 2016). In the following, we drop the subscript DE in w and , but keep it, for clarity, in σ_DE.

2.4. Parameterising the growth of structure: choice of γ

Outside of the ΛCDM framework, both the background and the perturbations can be modified (see Planck Collaboration XIV 2016, for a collection of scenarios and constraints). A simple way that was extensively used in the past to model the modified growth of perturbations is based on the observation that in ΛCDM the growth rate of Eq. (25) is well approximated by

(32)

with a constant growth index parameter γ ≈ 0.55 (Lahav et al. 1991; Linder 2005). A scenario with modified growth then corresponds to a different value of γ. However, for most realistic modified gravity models, a constant growth index γ is too restrictive if scale dependence exists. More importantly, the prescription of Eq. (32) is incomplete as a general description of the evolution of perturbations – which in general requires, under the assumption of adiabaticity, at least two degrees of freedom as a function of time and space (Amendola et al. 2008) – and so we prefer a more complete parameterisation, as discussed below. Nonetheless, we include γ as a parameter in our analysis in order to facilitate comparison with the Euclid Red Book (Laureijs et al. 2011) where this case was explored. We notice that the assumption done in this case of no modifications to the lensing potential, detailed further in the rest of this section, effectively reduces the degrees of freedom to a single one. Generalisations or alternative parameterisations are left for future work.

As mentioned, we need two functions of time and space to describe the evolution of the perturbations in general. To ensure that the recipes for different probes consistently implement the same assumptions for the linear evolution of perturbations, we need to relate γ explicitly to these two functions. In order to understand this relation, it is convenient to make these two free functions explicit in the perturbation equations. Since these functions are free, there are different choices (equally general) that can be made to define them. Here we express the temporal and spatial metric potentials, Ψ and Φ, in terms of the (μ_MG, Σ_MG)⁸ parameterisation (see e.g. Planck Collaboration XIV 2016). The first free function is μ_MG, which parameterises the growth of structure and is implicitly defined as a function of scale factor and scale via

(33)

where ρ_i is the energy density of the generic component i, p_i its pressure and δ_i its comoving density perturbation, while σ_i is the anisotropic stress, which is non-vanishing for relativistic species. The second free function Σ_MG describes the deflection of light and is defined through

(34)

An alternative option is the pair (μ_MG, η), with η defined as

(35)

which reduces to⁹

(36)

for negligible shear, as happens at the (low) redshifts of interest for Euclid. Equivalently, combining Eqs. (33) and (35), one has

(37)

These definitions differ (at high redshift) from the ones used in Planck Collaboration XIV (2016), since they make the anisotropic term explicit in the equations, separating it from contributions beyond ΛCDM that may be contained in μ_MG, Σ_MG, and η. The case for μ_MG = 1, Σ_MG = 1 (or equivalently η = 1) corresponds to ΛCDM, while any deviation – either due to modified gravity or to extra relativistic species – is encoded in at least one of the functions being different from 1. At low redshifts, the anisotropic stress terms are negligible, and the equations above are simplified, matching the definition used in Planck Collaboration XIV (2016).

Linear perturbations for a model are then fixed only once we make a choice for two of these free functions (such as (μ_MG, Σ_MG) or (μ_MG, η)). When assuming Eq. (32), the function μ_MG defined in Eq. (33), which in general is a function of time and space, can be converted into a function of the scale factor and of γ. In particular, (see Mueller et al. 2018, where μ_MG, Σ_MG are called G_M, G_L, respectively) μ_MG is related to γ via the following expression:

(38)

where the prime denotes differentiation with respect to ln a and Ω_m = Ω_m(a); for a constant γ, the last term vanishes and μ_MG only depends on the scale factor a through Ω_m and H. Fixing γ, however, does not fix the second free function defining linear perturbations (such as Σ_MG or η described above), for which there is still a choice to be made to define the model.

A possible choice for the second condition that fixes linear perturbations is Σ_MG = 1, since in this case light deflection (and therefore the lensing potential) is the same as in ΛCDM. This is also a reasonable choice from a theoretical point of view, as it is realised in standard scalar tensor theories, such as Brans–Dicke and f(R) (see Amendola et al. 2008; Pogosian & Silvestri 2016). In the limit Σ_MG → 1, the usual equations for the lensing power spectrum are valid (see Eq. (122)), which means the weak lensing description used for ΛCDM can still be applied. This choice is therefore also the one typically implicitly adopted in past analyses of γ, such as the one described in the Red Book. We restrict our analysis to the choice (γ, Σ_MG = 1) here as well to facilitate the comparison with previous analysis.

Finding a value of γ different from ΛCDM (i.e. from a value around 0.55) is then only related to having μ_MG differ from 1 (i.e. from the expected value of this function in ΛCDM), which in turn physically means that the matter spectrum P_δδ(k, z) is affected by a different growth rate. Details on how γ is specifically treated in different observational probes, and in the non-linear regime, is further discussed, separately, in Sects. 3.2.3–3.3.4.

2.5. Impact of neutrinos

In the presence of massive neutrinos, the definition of the linear growth rate in Eq. (25) needs to be modified to allow for a dependence on both redshift and scale, f(z, k), even at the linear level. To handle such dependencies in a semi-analytical way, approximations to the growth rate f(z, k) as a function of the neutrino fraction f_ν = Ω_ν, 0/Ω_m, 0 have been developed. One approach is the fitting formulae of Kiakotou et al. (2008), where

(39)

with

(40)

and the functions A(k), B(k), and C(k) have been obtained via fits of power spectra computed using the Boltzman code CAMB (Lewis et al. 2000).

Figure 1 shows the effect of massive neutrinos on f in the linear regime, as ratios between the massive and massless neutrino cases, that it is to say that we show the function μ_ν(k), which is independent of redshift. As can be seen, in the presence of massive-neutrinos the linear growth rate acquires a scale dependence, which decreases with decreasing neutrino mass, as the suppression of perturbations due to neutrino free streaming decreases, and as the free streaming scale becomes larger.

Given the current upper limits on the total neutrino mass (Palanque-Delabrouille et al. 2020), the growth suppression is of the order of 0.6% at maximum, and mainly affects scales k >  0.1 h Mpc⁻¹. For simplicity, we ignore this sub-percent effect in our forecasts, assuming f to be scale independent for values ∑_im_ν, i = 0.06, 0.15 eV tested, and computed in the massless limit. For ∑_im_ν, i = 0.06 eV solar oscillation experiments constrain neutrinos to be in a normal hierarchy, and for ∑_im_ν, i = 0.15 eV either an inverted or a normal hierarchy (Jimenez et al. 2010); we choose a normal hierarchy in this case.

Concerning the linear matter power spectrum P_δδ(k, z), in the presence of massive neutrinos, Eq. (28) in Sect. 2.3 can be modified by replacing the transfer function T(k) by a redshift-dependent one, 𝒯(k; z), while keeping the scale-independent linear growth factor D(z) as given in the absence of massive-neutrino free-streaming (Takada et al. 2006; Takada 2006; Eisenstein & Hu 1997; see Eq. 25 of Eisenstein & Hu 1997)¹⁰. The fiducial value of D(z) at each redshift can be computed through numerical integration of the differential equations governing the growth of linear perturbations in the presence of dark energy (Linder & Jenkins 2003). The linear transfer function 𝒯(k; z) depends on matter, baryon, and massive-neutrino densities (neglecting dark energy at early times), and is computed via Boltzmann solver codes in each redshift bin. The forecasts presented in this work will assume a fixed Ω_m, 0, that is when Ω_ν, 0 is varied, Ω_cb, 0 ≡ Ω_c, 0 + Ω_b, 0 varies as well, in order to keep Ω_m, 0 unchanged.

2.6. Magnification bias

We note that we have not included the effects of magnification bias in our forecasts. Magnification bias refers to the distortion of the observed galaxy number counts caused by gravitational lensing (see e.g. Duncan et al. 2016; Thiele et al. 2020; Tanidis et al. 2020), and it can be important for Stage IV galaxy surveys. Concerning galaxy clustering analyses with samples of photometrically detected galaxies, previous works have shown that neglecting the magnification effects can lead to significant shifts for several cosmological parameters (see e.g. Cardona et al. 2016; Montanari & Durrer 2015). Very recently it was shown that magnification bias is also relevant for galaxy clustering analyses of galaxies detected with spectroscopic techniques. More specifically neglecting it may lead to biased parameter estimation, like an incorrect inference of the growth rate of structures (Jelic-Cizmek et al. 2020). Furthermore, magnification can also be important for the cross-correlations between weak lensing and galaxy clustering of photometrically-detected galaxies (Ghosh et al. 2018).

In our forecasts we have only focused on the uncertainties on the different parameters. However, it is important to investigate the effect of magnification, in both the shift of the best-fit and the uncertainties on the parameters, in the context of model validation for Euclid’s final data analysis pipeline.

2.7. The standard ΛCDM model and its extensions

The spatially flat ΛCDM model is the baseline case considered in this paper, and corresponds to having a cosmological constant (w₀ = −1, w_a = 0). For a spatially flat cosmology Ω_K, 0 = 0 and the value of Ω_DE, 0 is a derived parameter, since Ω_DE, 0 = 1 − Ω_m, 0. The baseline is then described by a minimal set of 6 parameters:

• Ω_b, 0 and Ω_m, 0, the baryon and total matter energy densities at the present time;

• h, the dimensionless Hubble parameter, describing the homogeneous background evolution;

• σ₈, describing the amplitude of density fluctuations;

• n_s, the spectral index of the primordial density power spectrum;

• ∑m_ν: the sum of neutrino masses.

In this work, we further study the power of Euclid primary probes for constraining deviations from the ΛCDM model by also analyzing extensions of this baseline parameter space. In particular, we consider the following extensions:

• spatially non-flat models, by varying {Ω_DE, 0} (equivalently one could vary Ω_K, 0, however we choose the first option in the numerical implementation);

• dynamical DE models, by varying the background values of w₀ and w_a;

• modifications in the growth of structures, by varying the growth index γ.

The specific fiducial choice of all parameters will be discussed in Sect. 3.1.5.

3. Fisher matrix formalism for forecasting

In this paper we use a Fisher matrix formalism to estimate errors for cosmological parameter measurements. In this section we describe the general formalism and define some specific quantities that will be used throughout. We also present the detailed recipes used to implement the Fisher matrix formalism to compute forecasts of the different cosmological probes in Euclid.

3.1. General formalism

The aim of the analysis presented here is to obtain estimates on the uncertainties on the cosmological parameter measurements, i.e. the posterior distribution P(θ|x) of the vector of (model) parameters θ, given the data vector x. Using Bayes’ theorem, this can be obtained as

(41)

where P(θ) is the prior information on our parameters, P(x) is the Evidence, and L(x|θ) is the likelihood of the data vector given the parameters.

The Fisher matrix (Bunn 1995; Vogeley & Szalay 1996; Tegmark et al. 1997) is defined as the expectation value of the second derivatives of the logarithmic likelihood function (the Hessian) and can be written in the general form

(42)

where α and β label the parameters of interest, θ_α and θ_β, and the derivatives are evaluated in the point θ_ref of the parameter space, which coincides with the maximum of the likelihood distribution. The Fisher matrix thus corresponds to the curvature of the logarithmic likelihood, describing how fast the likelihood falls off around the maximum. For a Gaussian likelihood function this has an analytic expression that depends only on the expected mean and covariance of the data, i.e.

(43)

where μ is the mean of the data vector x and C = ⟨(x − μ)(x − μ)^⊺⟩ is the expected covariance of the data. The trace and sum over p or q here represent summations over the variables in the data vector. Often, for Gaussian distributed data x with mean ⟨x⟩=μ and covariance C, either the mean is zero, or the covariance is parameter-independent. In both cases, one of the two terms in Eq. (43) is non-vanishing. In Sects. 3.2–3.4 we specify which of these terms are used, and in the cases that either or both can be used we show both expressions.

Once the Fisher matrix is constructed, the full expected error covariance matrix of the cosmological parameters is the inverse of the Fisher matrix¹¹,

(44)

The diagonal elements of the error covariance matrix contain the marginalised errors on the parameters. For example, the expected marginalised, 1-σ error on parameter θ_α (i.e. having included all the degeneracies with respect to other parameters), is

(45)

The unmarginalised expected errors, or conditional errors, can be computed by , i.e. the square root of the reciprocal of the appropriate diagonal element of the Fisher matrix. We can also define the correlation coefficient between the errors on our cosmological parameters as ρ, which contributes to the off-diagonal elements in the parameter covariance matrix,

(46)

It is important to note that ρ_αβ = 0 if α and β are completely independent. In order to “marginalise out” a subset of the parameters and obtain a smaller Fisher matrix, one should simply remove the rows and columns in the full parameter covariance matrix that correspond to the parameters one would like to marginalise over. Re-inverting will then give the smaller, marginalised Fisher matrix . This is equivalent to taking the Schur complement (Haynsworth 1968; Zhang 2005; Kitching & Amara 2009) of the Fisher matrix for the smaller subset of parameters.

3.1.1. Visualising confidence regions

The Fisher matrix can be used to plot the marginalised joint posterior probability, or projected confidence region, of two parameters (θ_α, θ_β), assuming that the joint posterior probability is well approximated by a Gaussian. The log-likelihood function is then locally close to the function of a hyper-dimensional ellipsoid, defined by the Fisher matrix of Eq. (42), and the projection on any two-dimensional sub-space is simply an ellipse. The parametric form of the confidence level ellipses is defined by the semi-minor and semi-major axes of the ellipses, a and b, and the angle of the ellipse ϕ; a and b are expressed in terms of the larger and smaller eigenvalues of the covariance matrix, and ϕ is related to the ratio of the y-component to the x-component of the larger eigenvector (as given by Eq. (40.72) of the Particle Data Group’s Review of Particle Properties section on statistics¹²). They are:

(47)

where A is a constant factor defined as A² = 2.3,  6.18,  11.8 for two-parameter contours at 1-, 2-, and 3-σ confidence level (C.L.), respectively (see e.g. Press et al. 2007), where the C_αβ are defined in Eq. (44).

3.1.2. Figure of merit

When considering a particular experiment, its performance in constraining specific parameters, e.g. those related to the dark energy model, can be quantified through a figure of merit (FoM). We adopt the FoM as defined in Albrecht et al. (2006) which is inverse proportional to the area of the 2-σ contour in the marginalised parameter plane for two parameters θ_α and θ_β, under the usual Gaussianity assumption inherent in the Fisher forecast formalism. Therefore the FoM can be calculated simply from the marginalised Fisher submatrix for those two parameters. That is

(48)

The above formula was utilised in the original dark energy FoM definition from Albrecht et al. (2006), which considers the w₀, w_a parameterization defined in Sect. 2.1:

(49)

which is in fact equivalent (within a constant factor) to other definitions used in the literature (see, e.g., Rassat et al. 2008; Wang 2008; Amendola et al. 2012; Majerotto et al. 2012; Wang 2010). Different FoMs can also be defined for any arbitrary set of parameters by simply taking the determinant of the appropriate Schur complement. As an example, for a model where the linear growth rate f(z) of density perturbations is parameterised as in Eq. (32), the FoM would read , and similarly for any other cosmological parameter pair (e.g. Majerotto et al. 2012). Throughout this paper we refer to the “FoM” as being that defined in Eq. (49).

3.1.3. Correlation matrix and the figure of correlation

As shown in Casas et al. (2017) one can define a correlation matrix P for a d-dimensional vector p of random variables as

(50)

where C is the covariance matrix. By definition, P is equal to the unit matrix if all parameters are uncorrelated, while it differs from it if some correlation is present. We plot this matrix in Sect. 5 to visualise the impact on the attainable cosmological constraints of cross-correlations in the photometric survey. We further adopt the introduction of a “figure of correlation” (FoC), first defined in Casas et al. (2017), defined (here without the logarithm) as:

(51)

which is 1 if parameters are fully uncorrelated. Off-diagonal non-zero terms (indicating the presence of correlations among parameters) in P will correspond to FoC >  1. The FoC and the FoM are independent quantities, (see Casas et al. 2017, for a geometrical interpretation).

3.1.4. Projection into the new parameter space

The Fisher matrix is defined for a set of parameters θ, but a new Fisher matrix can be constructed for an alternative set of parameters p(θ). In this case the new Fisher matrix S is related to the original Fisher matrix F by a Jacobian transform

(52)

where the Jacobian matrices ∂θ_α/∂p_i relate the original to the new parameterization. Notice that if p(θ) is not a linear function then the likelihood in the new parameters is in general not Gaussian, even if it was in the original parameter space θ. Hence the Gaussian approximation inherent in the Fisher formalism may be valid for one choice of parameters but not for another.

3.1.5. Fiducial parameter values

In this section we detail the choice of fiducial model, about which the derivatives and cosmological quantities used in the Fisher matrices considered in this paper are computed (for the ΛCDM case and for its extensions). Since we are combining information from the Euclid galaxy clustering and weak lensing probes, our final Fisher matrix should have consistent rows and columns across all the probes. Within the assumption of a minimal cross-covariance between the probes, the combination of the Fisher matrices related to the two observables can then be achieved through a direct addition of the matrices. The final minimal set of cosmological parameters for both probes is defined as (see Sect. 2.7)

(53)

The (minimal) dark energy and modified gravity models discussed in Sects. 2.1 and 2.4 involve the following parameters in addition to the minimal ΛCDM set:

(54)

(55)

(56)

These appear in the expansion history and growth of structure respectively (see Eqs. (8) and (32)). The chosen parameter values of our fiducial cosmological model are shown in Table 1.

In addition to the fiducial cosmology we also need to define the survey specifications. These are described in more detail in the following sections, but we summarise here some of the main characteristics that are common to both GC and WL.

In Tables 2–4 we describe the specifications used to compute the forecasts for GC and WL observations within the Euclid mission. In this paper we adopt the specifications of the Euclid Red Book (Laureijs et al. 2011) as closely as possible, to allow for a direct comparison with those forecasts. We note that during the writing of this paper the Euclid Consortium underwent a process of science performance verification (SPV) – a review process that also produced cosmological parameter forecasts that included systematic effects induced by instrumental and data reduction procedures. The forecasts in the SPV used the validated IST:F codes, but slightly different areas and depth specifications were used that were motivated from the internal Euclid simulations. The SPV results are not publicly available, and since these simulation-derived quantities may change before publication of any SPV results, we choose to remain with the Red Book specifications; however in future work these will also be updated.

3.1.6. Required accuracy in Fisher matrix forecasts

Fisher matrices with a large number of correlated parameters can be difficult to compute and to invert accurately due to the numerical complexity involved. Therefore one requires accurate estimates of the parameter uncertainties, i.e. “the error on the errors”. One aspect is simply the accuracy of the Fisher matrix approximation itself (i.e. the fact that it assumes a Gaussian likelihood), but others are the precision with which the Fisher matrix is computed. Typically, problems arise in the evaluation of the derivatives, the approximation of a sum over redshift planes rather than an integral, and the precision of the inversion of the Fisher matrix to obtain the parameter covariance matrix giving the parameter uncertainties and their correlations. The assessment of the accuracy with which one needs to compute a Fisher matrix is simple in one dimension, with

(57)

where F is the Fisher matrix and C is the correlation matrix, as defined above. The computed fractional error on the parameter is then of the same order as the numerical fractional error on the Fisher matrix.

In the following we assume a 10% requirement on the uncertainty in the parameter error arising from numerical inaccuracies. However we note that this is a somewhat arbitrary requirement. We see that for one dimension this does not impose a tight requirement on numerical accuracy, since one would expect Fisher matrix numerical fractional errors < 10⁻⁴.

For more than one parameter, the covariances between them are important and this will be crucial for determining the sensitivity. We can obtain the main effects by considering the amplitude and orientation of the parameter confidence regions and ellipses for two parameters, θ_α and θ_β. In two dimensions,

(58)

(59)

(60)

and

(61)

where for a given α, β pair where α ≤ 2 and similarly for β¹³. Since

(62)

then we see that the first term in Eq. (61) is just what we had in the one parameter case and does not impose tight restrictions on the Fisher matrix precision. Therefore, it is the determinant term, involving covariances and the potential for degeneracies between parameters, that requires care. In two dimensions, for a given α and β pair

(63)

where the correlation coefficient ; we also note that from Eq. (46). From this, one has

(64)

Again, we see that the first term indicates that the precision of the Fisher matrix propagates to errors on the parameter uncertainty of the same order, and hence only the second term can cause a change. The final result, in two dimensions, is

(65)

where δ_max is the maximum fractional imprecision on an element of the Fisher matrix. We note that for two dimensions r_αβ is a single number and so the summation convention does not apply to the above equation (one could write σ_α | β, i.e. the error on α given a fixed β). Assuming r_αβ ≤ 0.999 for realistic (well chosen) parameters and a requirement that δσ_α/σ_α ≤ 0.1 we find a required numerical accuracy of δ_max ≤ 10⁻⁴, which is a reasonable aim for a good numerical implementation.

We note that the FoM accuracy is simply given by

(66)

for any parameter combination given in Eq. (48). Thus if we want the accuracy of FoM to better than 10%, this becomes the most restrictive constraint.

Another important element is the orientation of the joint parameter contour ellipse. This is especially relevant for combining Euclid’s constraints with other probes. The orientation ϕ is defined in Eq. (47), and we find its uncertainty is

(67)

For many parameter pairs, for example (w₀, w_a), one of the parameter uncertainties is much larger than the other. Assuming the error on θ_β to be much larger than the one on θ_α (σ_β ≫ σ_α) and using Eq. (46), we can rewrite

(68)

The dominant term overall is the first term in the square brackets, and is suppressed with respect to δC_ββ/C_ββ itself by σ_α/σ_β, hence it is smaller than the effect on the error on the amplitude of the parameter uncertainty. If we require the parameter uncertainty to be δσ_α/σ_α ≃ 0.1, then the orientation angle will be accurate to δϕ ≈ 2%.

The conclusion is that, for ≤10% precision on parameter errors and FoM and ≤2% precision on the orientation of the contours, the required fractional Fisher matrix precision (both in calculation of elements, and its inversion) is approximately 10⁻⁴. In other words, one certainly does not need to push near machine precision. We note here that for two-sided or three-point finite differences for evaluation of derivatives with a parameter stepsize ϵ, the errors go as ϵ², while five-point differences go as ϵ⁴.

3.2. Recipe for spectroscopic galaxy clustering

This section describes the forecasting procedure recommended for our galaxy clustering observable – the full, anisotropic, and redshift-dependent galaxy power spectrum, P(k, μ; z). We spend most of the section describing how to model our observable. We conclude by specifying how to use our observable in the Fisher matrix calculation, as described generally in Sect. 3.1.

The initial goal is to calculate a Fisher matrix for the following full set of cosmological parameters:

(69)

The first set of parameters determines the shape of the linear matter power spectrum (as defined in Sects. 2.1 and 2.3): the physical densities of baryons, the physical density of total matter, the dimensionless Hubble parameter, and the scalar spectral index. In particular, these parameters control the transfer function and the power law of the spectrum of primordial density perturbations (see Sect. 2.3 for details).

The second set carries the uncertainty in our theoretical knowledge about late-time non-linearities and redshift-space distortions (RSDs), which damp the baryon acoustic oscillation (BAO) signal in the galaxy power spectrum, and cause the so-called fingers-of-God effect. These are considered nuisance parameters, to be marginalised over. Such marginalisation necessarily increases the uncertainty on the other parameters (in the Fisher matrix). The non linear recipe adopted in this paper will be discussed in detail later (in Sect. 3.2.2).

Like the second set, the redshift-dependent set of parameters is also specific to galaxy clustering. Anisotropies in the power spectrum with respect to the line-of-sight direction due to the Alcock–Paczynski (AP) effect are parameterised in terms of the angular diameter distance D_A(z) and the Hubble parameter H(z), while the pattern of RSDs is described by the combination fσ₈(z) ≡ f(z)σ₈(z) defined in Sect. 2.3; note that here we generalise Eq. (30) to be redshift-dependent by changing the redshift argument to the power spectrum in that integral. The additional parameters bσ₈(z) ≡ b(z)σ₈(z) and P_s(z), characterise the bias of our target galaxy sample, and respectively the shot noise residual, as we describe below. Each parameter in this set is varied freely in each redshift bin of the survey. In this way, the constraints on D_A(z), H(z), and fσ₈(z) obtained at each redshift provide us with dark-energy-model-independent constraints on the expansion and structure growth histories of the Universe. The key cosmological information contained in these measurements can later be transformed into parameters specific to our chosen models. For practical ease, we calculate the natural logarithms of these parameters where possible.

These three sets of parameters capture all the main contributions of uncertainty in our model, some to be marginalised over, and some to be propagated into our final, general parameter set. After marginalising over our observational and theoretical “nuisance parameters” ln[bσ₈(z)], P_s(z), σ_p(z_i), and σ_v(z_i), we transform the constraints on the remaining set of redshift-dependent and “shape” parameters to those on the parameters of our chosen specific DE and modified gravity model by using a Jacobian matrix (see Eq. (52)), in terms of the final set of cosmological parameters θ = {Ω_b, 0,  h,  Ω_m, 0,  n_s,  σ₈,  ∑m_ν} with or without the extensions {Ω_DE, 0, w₀,  w_a,  γ}, as described in Sect. 3.1. Finally, we must determine the FoM (see Sect. 3.1.2) on the DE parameters of Eqs. (55) and (56), having marginalised over the rest.

3.2.1. The observable galaxy power spectrum

We work in Fourier space to describe the distribution of our target sample – Hα-emitting galaxies – given the statistics of the underlying dark matter field. We work in terms of P(k, μ; z), where k is the modulus of the wave mode in Fourier space, and is the cosine of the angle θ between the wave-vector k and the line-of-sight direction , and z is the redshift of the density field¹⁵. We give the final, full model for the observable galaxy power spectrum in Eq. (87) and describe all the individual effects in the preceding sections.

In total, five main observational effects must be modelled beyond simply calculating the matter power spectrum, the linear model of which is given in Eq. (28). In this brief section we show how we model each effect in the linear, as well as in the non-linear, regime of perturbation theory.

There are five main effects that we need to model:

1. the galaxy bias (in the case of Euclid spectroscopy the bias of Hα-emitting galaxies) with respect to the total matter distribution,

2. anisotropies due to RSD induced by the peculiar velocity component of the observed redshift causing discrepancies in the clustering strength measured at different angles with respect to the line-of-sight.

3. the residual shot noise that remains even after the known noise – due to Poisson sampling by target galaxies of a (theoretically) smooth underlying matter density field – has been subtracted,

4. the redshift uncertainty that suppresses the correlation between galaxy positions by smearing out the galaxy field along the line of sight,

5. distortions due to the AP effect, which introduces an anisotropic pattern in the power spectrum by rescaling wavemodes in the transverse and radial directions by different factors.

The combination of the above five effects are applied to the power spectrum model. We now describe how we model these effects, focusing first on the linear theory predictions. We describe (Sect. 3.2.2) how we extend this recipe to take into account the impact of non-linearities to derive a final model of the observed galaxy power spectrum.

Effective galaxy bias. We do not discuss the issue of galaxy bias in detail here, but refer the reader to the comprehensive review of Desjacques et al. (2018) and use the simple linear relation:

(70)

where P_δδ(k; z) (see Eq. (28)) is the linear matter power spectrum, which, when multiplied by the square of the redshift-dependent effective bias of the galaxy sample, b²(z), yields the galaxy linear power spectrum. In our case, the sample in question contains what will be detected by Euclid as Hα-line emitter galaxies.

Anisotropies due to RSD. The measured galaxy redshifts contain a non-cosmological contribution due to the line-of-sight component of the peculiar velocity

(71)

where v_∥ is the galaxy peculiar velocity along the line of sight. If one assumes that the observed redshifts are entirely cosmological when estimating the distances to each galaxy, this distorts the density field in a way that imprints a specific pattern of anisotropies onto the observed power spectrum; this is known as the RSD. As the peculiar velocity displacements are sourced by the true underlying density field, the pattern of the RSDs provides us with an additional source of cosmological information on the relation between the density and peculiar velocity fields, which in the linear regime depends on the growth rate parameter f(z).

In the linear regime, the relation between the real- and redshift-space galaxy power spectra (here labelled as “zs” for redshift-space), is given by Kaiser (1987)

(72)

where σ₈(z) is defined in Eq. (30), but generalised to be redshift-dependent. Written in this way, Eq. (72) illustrates that the quantities controlling the anisotropies in the two dimensional power spectrum are the combinations b(z)σ₈(z) and f(z)σ₈(z), which we treat as free parameters.

Redshift uncertainty and shot noise. The effect of redshift uncertainties results in a modification of the power spectrum in the form

(73)

where the factor F_z(k, μ; z) is given by

(74)

and accounts for the smearing of the galaxy density field along the line of sight direction k_∥ = kμ due to possible redshift errors. This error propagates into a comoving distance error (see e.g. Wang et al. 2013),

(75)

where σ_z(z) = σ_0, z(1 + z), and σ_0, z is a linear scaling of the redshift error. We note that the damping due to redshift errors does not vary with changes in the expansion history, since k_∥ ∝ H(z) and σ_r ∝ H⁻¹(z) (Wang et al. 2013).

As described before, P_s(z) is a scale-independent offset due to imperfect removal of shot noise, which is taken to be P_s(z) = 0 at all redshifts in our fiducial model, but is allowed to vary with redshift. It is standard in galaxy clustering analysis to treat the shot-noise term as a free parameter (see e.g. Gil-Marín et al. 2016; Beutler et al. 2017; Grieb et al. 2017). We therefore include it here to assess the impact that marginalising over this parameter might have on the expected parameter constraints from Euclid.

AP projection effects. The measurement of the galaxy power spectrum requires the assumption of a reference cosmology to transform the observed redshifts into distances. Assuming an incorrect cosmology leads to a rescaling of components of the wave-vector k in the directions parallel and perpendicular to the line-of-sight, k_∥ and k_⊥, as

(76)

where the coefficients q_⊥ and q_∥ are given by the ratio of the angular-diameter distance and the Hubble parameter in the true and reference cosmologies as

(77)

We label the redshift dependence here, but will suppress this in the following equations for clarity. Eqs. (76) and (77) can be used to convert from the known reference k_ref and μ_ref to the unknown, true k and μ as

(78)

These relations can be used to model the impact of the particular choice of the reference cosmology on the observed power spectrum by mapping the values of k_ref and μ_ref, to the true ones as in Ballinger et al. (1996)

(79)

where the overall scaling of the power spectrum by follows from the dilation of the volume element due to the AP effect.

The full linear model for the observed power spectrum of Hα emitters. We can now write down the observed, linear anisotropic galaxy power spectrum as

(80)

where all k = k(k_ref, μ_ref) and μ = μ(μ_ref).

3.2.2. Non-linear scales for galaxy clustering

At late times and small scales, the matter distribution is affected by the non linear evolution of density fluctuations, which distorts the shape of the power spectrum beyond the predictions of linear perturbation theory. Most of the Euclid spectroscopic sample will lie in the range 0.9 <  z <  1.8. We expect that we may gain constraining power by including mildly non-linear scales into the analysis. However, this requires modelling non-linear effects that become relevant at the scales of the main feature we are trying to measure, namely the BAO signature. We follow a phenomenological model for non-linear effects in galaxy clustering similar to other recipes considered in the literature (e.g. Seo & Eisenstein 2007; Wang et al. 2013). We separately model two main consequences of non-linearities, which are relevant to our calculation of the observed galaxy power spectrum. Our final model of the observed power spectrum deviates from the linear recipe described in Sect. 3.2.1 in two separable ways, which we parameterise using two “non linear parameters”, σ_p(z) and σ_v(z). Both these parameters have the same physical meaning, and represent the linear galaxy velocity dispersion , with u_∥(r, a) = v_∥(r, a)/(aHf), and can thus be computed from the linear matter power spectrum as

(81)

However, we keep them separate to allow for some theoretical uncertainty in our practical model of non-linearities.

The first effect of non-linearity is that the description of RSD requires an additional factor that we model as a Lorentzian, which accounts for the finger-of-God (FoG) effect under the assumption of an exponential galaxy velocity distribution function (Hamilton et al. 1998). Secondly, we must account for the damping of the BAO feature, which we calculate with another exponential damping factor, acting in both directions: radial and transverse. This is accounted for by using the so-called “de-wiggled” power spectrum, P_dw(k, μ; z), (see Wang et al. 2013). Therefore we model the true non-linear underlying redshift-space power spectrum as (cf. Eq. (72) above)

(82)

where the combination 2f(z)σ_p(z) would correspond to the pairwise peculiar velocity dispersion along the line of sight, and P_dw(k, μ; z) is given by

(83)

We can imagine this sum as a simple selective suppression of amplitude of the wiggles only, applied first by suppressing the amplitude of the whole linear power spectrum in the first term, and the restoration of the power in the broadband power spectrum (excluding the BAO wiggles) in the second term. Therefore, P_δδ(k; z) is still the linear power spectrum, whereas P_nw(k; z) is a “no-wiggle” power spectrum with the same broad band shape as P_δδ(k; z) but without BAO features, which we compute using the formulae of Eisenstein & Hu (1998). The function

(84)

is the non-linear damping factor of the BAO signal derived by Eisenstein et al. (2007). The parameter σ_v(z) controls the strength of the non-linear damping of the BAO signal in all directions in three-dimensions.

We use Eq. (81) to define the fiducial values of σ_p(z) and σ_v(z). In practice, in Eqs. (82) and (84) we write the σ_v(z) and σ_p(z) as

(85)

(86)

where z_mean is the central redshift of the survey, and D(z) is the growth factor from Eq. (27), treating the values of σ_v(z_mean) and σ_p(z_mean) as free parameters. We note that we assume in Eq. (81) that , so that the above equation is most likely not the correct scaling in detail; we are making the assumption that we understand the increase with time of the non-linear impact of the parameters as being related to the growth of structure. This is in line with the intuition that σ_v(z) and σ_p(z) are related to the velocity dispersion of galaxy populations

An additional effect of the non-linear evolution of density fluctuations is to modify the broadband shape of the power spectrum, increasing the power on small scales. In the range of scales and redshifts in which we apply our recipe (see Sect. 4.2), this effect can be approximately described by changing the value of the shot noise term P_s(z) that we include in our recipe (see e.g. Ho et al. 2012). A more detailed treatment of small scale non linearities is out of the scope of our analysis. We note that this is a conservative treatment of non-linear scales; the implementation of a more detailed modelling of non linearities in which we could accurately predict these broadband distortions as a function of the cosmological parameters would allow us to extract additional information from the shape of the observed power spectrum. By marginalizing over P_s(z) we are conservatively discarding this information from our forecasts.

The full non-linear model for the observed power spectrum of Hα emitters. We can now write down the full observed, non-linear anisotropic galaxy power spectrum as

(87)

where all k = k(k_ref, μ_ref) and μ = μ(μ_ref). All the quantities in this expression can be calculated using equations given in this section. A number of alternatives have been considered in literature, but so far, the phenomenological model of Eq. (87) has proved to be at least practical, all models being in any case only an approximation (see e.g. Scoccimarro & Frieman 1996; Bernardeau et al. 2012; Baldauf et al. 2015; Casas et al. 2017, for further discussion).

3.2.3. Modelling beyond ΛCDM

In order to calculate the DE FoM (see Sect. 3.1.2), we must transform our original parameter set (see Eq. (69)) into the subset of parameters whose FoM we need. We choose to freely vary the redshift-dependent parameters in Eq. (69) to avoid limiting ourselves to any particular DE or MG model throughout most of our galaxy clustering analysis. We then project our initial Fisher matrix using a Jacobian matrix of derivatives into the parameter space specific to our non-ΛCDM model, as described in Sect. 3.2.5.

In particular, the redshift-dependent parameters D_A(z) and H(z) will carry most of the information on the expansion history of our universe and therefore contribute most to the constraints on w₀ and w_a. On the other hand, it is the redshift-dependent parameter fσ₈(z) that most strongly depends on γ via Eq. (32). Therefore, this parameter (in all its redshift bins) will contribute most to the constraints we achieve on γ. However, when it comes to γ, this simple picture of dependence is not the whole story. Since σ₈(z) = σ₈D(z) depends on the growth factor as well, it directly depends on the γ parameter. The same is also true for the matter power spectrum itself. For this reason, we always use the ratio . In this way the dependence on γ is cancelled as long as the approximation is valid. Secondly, the parameter σ_p is calculated from an integral over P_δδ(k; z), and hence also depends on γ. This means that when σ_p is considered simply as a nuisance parameter and marginalised over, no constraints on γ come from it. On the other hand, if the value of this parameter is computed through Eq. (81), it contributes to the constraint on γ as well as the shape parameters¹⁶. Finally the function g_μ from Eq. (84) depends on the cosmological parameters (γ and shape) through P_δδ in σ_v (calculated in Eq. (81)) and f(z), which we also include in the modelling.

3.2.4. Covariance matrix for galaxy clustering

The statistical constraining power of a galaxy survey depends mainly on the abundance of its target galaxy sample and how the survey strategy determines the survey’s selection function. Modelled simply, this means that the statistical errors on the power spectrum can be estimated by knowing the observed number density of its galaxy targets (in our case Hα emitters) and the survey volume. The expected number of Hα emitters per unit area and redshift intervals that will be observed by Euclid is given in Table 3. This can be converted to the number density in each redshift bin,

(88)

where the width of the redshift bin is given in terms of the minimum and maximum redshift of the bin, Δz = z_max − z_min. The minimum and maximum redshifts, and the volume of the redshift bin, V_s, are given in Table 3, while the survey area A_survey is specified in Table 4.

Assuming that the observed power spectrum follows a Gaussian distribution, its covariance matrix can be approximated by

(89)

where δ_D is the Dirac delta function (reflecting the independence of different Fourier modes) and V_eff(k) is the effective volume of the survey given by

(90)

Here, P_obs(k, μ; z) is the full observed two-dimensional power spectrum, and all quantities are evaluated for the fiducial cosmology.

3.2.5. Applying the Fisher formalism to galaxy clustering

Following Tegmark (1997) and Seo & Eisenstein (2003) the Fisher matrix for the observed galaxy power spectrum in a redshift bin centred on z_i is:

(91)

where α and β run over the cosmological parameters we wish to vary, which in turn enter implicitly into P_obs(z; k, μ) (see Eq. (87)). The values of k and μ in Eq. (91) are those obtained in the reference cosmology, k_ref and μ_ref, but we drop the subscript “ref” from now on for clarity. The total Fisher matrix is then calculated by summing over the redshift bins:

(92)

We note here that for the redshift-dependent parameters, we assume that the redshift bins are independent, therefore the sum in Eq. (92) reduces to a single term. For example, if we take α  =  H(z_j)  =  H_j, we have

(93)

since

(94)

The sum in Eq. (92) then reduces to a single term with i = j:

(95)

This includes both the off-diagonal terms containing the covariance between different redshift bins of the same redshift-dependent parameter, and the covariance between different z-bins of different parameters, leaving only the non-zero covariances between redshift-dependent parameters within the same redshift bin. In other words,

(96)

if α(z_i) and β(z_j) are redshift-dependent parameters, and in the following we abbreviate these subscripts to ij.

Explicitly, the model-independent Fisher matrix is given by Eq. (91) for one single redshift bin, where the subscripts α and β run over the shape- and redshift- dependent parameters. In this way we have a matrix for one single redshift with 6 + (5 × 1) rows and columns in the full non linear case as reported in Eq. (69) (in the linear case, the number of z-independent parameters falls to 4). Schematically, let us assume that the Fisher matrix for one redshift bin is

(97)

where each term corresponds to a block matrix, with each element calculated from derivatives with respect to the parameters of the row and column. The first 6 × 6 block is composed of derivatives with respect to the shape parameters only, the 5 × 5 block by those with respect to the redshift dependent parameters only, and the 5 × 6 block (which is the transpose of the 6 × 5) contains the mixed derivatives. The total Fisher matrix, including all redshift bins, is then a matrix of dimension 6 + 5 × n_bin. Making use of the formalism in Eq. (97) we then have

(98)

where we must sum over the contributions to the shape parameter elements from each redshift, but treat the redshift bins as creating separate redshift-dependent parameters.

Derivatives. In order to convert the data covariance into the desired parameter covariance in the Fisher matrix calculation in Eq. (91), we need the derivatives of the observable power spectrum P_obs with respect to the parameters of the chosen parameter space. From Eq. (87) one may calculate all derivatives numerically using finite differencing methods. In practice, a derivative step should be chosen such that the convergence level towards the true value (generally the analytical one) is reached. However given that the analytical values are inaccessible for most of the derivatives computed here, the best possible test consists in checking the step range for which the derivatives are stable. Choosing large steps raises the truncation error due to the discretization of the parameter values, while using too small steps increases roundoff errors. We note that for the residual shot noise, the following analytical equation may be used instead:

(99)

(see Sect. 4.2 for practical guidelines on how to calculate the derivatives).

Projection to a new parameter space. As discussed in Sect. 3.2, we first consider as independent parameters the values of D_A(z_i), H(z_i), fσ₈(z_i), bσ₈(z_i) and P_s(z_i), in each of the N_z redshift bins. We use two additional nuisance parameters, σ_v(z_mean) and σ_p(z_mean), to account for the uncertainty in our modelling of non-linear effects. The full list of parameters, {θ_i}, of the initial Fisher matrix then includes redshift-independent parameters, {θ_zi}, and redshift-dependent ones {θ_zd}: {θ_i}={θ_zd}+{θ_zi}. These parameters are listed in Eq. (69) and recalled below for convenience:

• 5 × N_z redshift-dependent parameters, {θ_zd}={ln D_A(z_i), ln H(z_i), ln[fσ₈(z_i)],  ln[bσ₈(z_i)],  P_s(z_i)}

• 4 cosmological redshift-independent parameters and 2 non-linear model parameters, {θ_zi}={ω_b,0, h, ω_m,0, n_s, σ_v(z_mean),σ_p(z_mean)}.

We must first decide which of these parameters to treat as nuisance parameters, and which we assume contribute to our knowledge of our final parameter set. We marginalise out the following nuisance parameters:

• 4 nuisance parameters to be marginalised over, {θ_nuisance}={ln[bσ₈(z_i)],  P_s(z_i),σ_v(z_mean),σ_p(z_mean)}.

We finally wish to forecast the errors on the following final cosmological parameter set in order to be able to combine the GC-only results with the constraints expected from the weak lensing measurements:

• 9 cosmological parameters as chosen in Sect. 3.1.5, {θ_final}={Ω_b, 0,  h,  Ω_m, 0,  n_s,  Ω_DE, 0 w₀,  w_a,  σ₈,  γ}.

As described in Sect. 3.1.4, in order to transform our Fisher matrix from one parameter set to another, we must understand the mutual dependence of the remaining initial and final parameters, which is encapsulated in the Jacobian matrix, J of derivatives of the initial parameters, α, with respect to the desired final parameters, κ, as

(100)

and so

(101)

which defines the final Fisher matrix that we use.

FoM. The specific Fisher submatrices needed to compute the DE FoM (see Sect. 3.1.2) are calculated through a combination of marginalising out the DE-independent parameters and projecting others into the DE parameter space of choice. The exact choices we make, in terms of scales included and cosmology considered, can impact the final result. We have tested the GC recipe for a variety of different choices and in Sect. 5 we identify a pessimistic and optimistic case. Different marginalization approaches are described in the literature, (e.g. Wang et al. 2010). We choose to follow what could be referred to as the “full P(k) method, with growth information included”, where, after marginalizing out the observational nuisance parameters, b(z_i)σ₈(z_i) and P_s(z_i), the remaining Fisher matrix is projected into the final parameter space of {θ_final}, and the FoM is obtained after again marginalizing out everything except w₀ and w_a. We therefore first marginalise out the observational nuisance parameters, bσ₈(z_i) and P_s(z_i) by removing their corresponding rows and columns from the parameter covariance matrix, (as usual, see Sect. 3.1 and Coe 2009). We are then left with 3 × N_z redshift-dependent parameters and 4 shape parameters. We project these 3 × N_z + 4 parameters into the required final parameter space using Eq. (101) above.

3.3. Recipe for weak lensing

This section describes the forecasting procedure recommended for weak lensing, using the observable tomographic cosmic shear power spectrum. We first describe the calculation of this observable, and then how it is applied in the Fisher matrix context; that is described in general in Sect. 3.1. We note that we refer to weak lensing as the physical phenomenon, and cosmic shear as the observable and summary statistic that uses weak lensing to extract cosmological parameters.

The large-scale cosmic structure deflects the path of photons from distant galaxies, which induces distortions in the images of these galaxies. Locally such distortions can be decomposed at the linear level into convergence and a (complex) shear distortion, which are respectively related to the magnification and shape distortion of the image. More specifically, these are the trace and trace-free parts of the (inverse) amplification matrix that describes the linearised mapping from lensed (image) coordinates to unlensed (source) coordinates (see e.g. Kilbinger 2015, for a general review). Convergence is a change in the observed size of a galaxy, and shear is a change in the observed third flattening or third eccentricity (known as “ellipticity”, or polarisation). In this paper we only consider the cosmological signal in the shear field (see e.g. Alsing et al. 2015, for a discussion of the use of the convergence in cosmological parameter estimation). The shear field caused by large-scale structure has a zero mean because we assume isotropy and homogeneity of the Universe, but its two-point correlation function and its power spectrum (as well as higher-order correlations) contain cosmological information that probes both the background evolution of the Universe and the growth of cosmic structure; these two-point statistics are known as “cosmic shear”.

There are several ways in which the cosmic shear can be exploited for cosmological parameter inference. The real/configuration space measurement of the two-point statistic, as a function of angular seperation on the celestial sphere, is known as the cosmic shear correlation function. On the other hand, the angular spherical-harmonic measurement of the two-point statistic is known as the cosmic shear power spectrum. Both of these statistics can be computed in a series of redshift slices/bins that capture the geometry of the three-dimensional shear field. This redshift binning is known as “tomography” and is required in order to achieve high-precision dark energy measurements (Hu 1999; Castro et al. 2005; Casas et al. 2017; Spurio Mancini et al. 2018). Depending on which observable one is interested in, different approximations may be used to make the equations more tractable, but a thorough description of all methods is beyond the scope of this paper. In the following, we shall focus on the cosmic shear (tomographic) power spectrum.

In this section we describe this observable in detail. In addition we quantify the major astrophysical contamination of cosmic shear, which is the intrinsic (local) alignment (IA) of galaxies. The IA signal can be modelled as a projected power spectrum, similar to the cosmic shear observable itself.

As described in Sect. 3, θ_final = {Ω_b, 0,  h,  Ω_m, 0,  n_s,  Ω_Λ,  0,  w₀,  w_a,  σ₈,  γ} is the full set of cosmological parameters of interest. In contrast to the GC Fisher matrix, derivatives of the weak lensing observable (cosmic shear) can be computed directly in this parameter space, and therefore no pre-computation, or projection from a different parameter set is required.

3.3.1. The observable tomographic cosmic shear power spectrum

In this section we define the primary weak lensing observable. While the cosmological signal with which we are concerned is the additional ellipticity caused by the lensing of the large-scale structure, that we summarise in the cosmic shear power spectrum, we also wish to model the primary astrophysical systematics. We consider five main quantities that must be modelled in order to recover the observable cosmic shear power spectrum: Firstly, the (theoretical) cosmic shear power spectrum, namely the primary cosmological power spectrum; secondly, the intrinsic alignment power spectrum, modelling the local alignment of galaxies, representing the main astrophysical systematics; as a third quantity, we consider the small-scale part of the matter power spectrum, including a halo model describing the clustering of dark matter on small scales beyond linear theory (k <  7 h Mpc⁻¹) to which the cosmic shear power spectrum is particularly sensitive (Taylor et al. 2018a), which includes the impact of baryonic feedback (e.g. Semboloni et al. 2011; Copeland et al. 2018); the fourth item includes photometric redshifts and number density, which model the inferred uncertainty in galaxy positions due to the broad-band estimates of the redshifts, where the redshift distribution affects the signal part of the cosmic shear power spectrum but the total number of galaxies does not; as fifth, and last point, the shot noise due to Poisson sampling by galaxy positions of the shear field, which is affected by the total number of galaxies observed. We now describe each of these effects and show how the observed cosmic shear power spectrum is constructed. The discussion of the non-linear modelling of the matter power spectrum is postponed to Sect. 3.3.2.

The cosmic shear power spectrum. As defined above, shear is the change in the ellipticity of the image of a background galaxy, caused by the lensing effect of large-scale structure along the line-of-sight. For an individual galaxy we express this, to linear order, as

(102)

where γ is the cosmological shear and ϵ^I is the intrinsic (unlensed) ellipticity. This ellipticity is a spin-2 quantity with zero mean over a large survey area and a non-zero two-point correlation function or power spectrum, encoding information on both the expansion history of the Universe and the matter power spectrum. The spherical harmonic transform of the two-point correlation function is the angular power spectrum. While the full computation of the cosmic shear power spectrum is relatively laborious (see e.g. Taylor et al. 2018b) due to spherical Bessel functions and several nested integrals, the Limber approximation (expected to be suitable for angular scales of ℓ ≳ 100; see Kitching et al. 2017; Kilbinger et al. 2017; Lemos et al. 2017; Kaiser 1992; Giannantonio et al. 2012) allows one to simplify its expression to

(103)

where i and j identify pairs of redshift bins, E(z) is the dimensionless Hubble parameter of Eq. (11), r(z) is the comoving distance, and P_δδ(k, z) is the matter power spectrum evaluated at k = k_ℓ(z) ≡ (ℓ + 1/2)/r(z) due to the Limber approximation; we define the weight function W^γ(z) as

(104)

(105)

with z_max the maximum redshift of the source redshift distribution. With respect to the standard formalism, we introduce the reduced window function (also colloquially known as the lensing efficiency), which will be important later. In the above we assume that the changes in the formalism due to non-flat spatial geometries can be captured via changes in the cosmological parameters, power spectrum and reduced window function rather than changes to the Limber approximation caused by modifying the spherical Bessel functions in the non-Limber approximated case to hyper-spherical Bessel functions; this assumption is well justified as show in Taylor et al. (2018b).

The intrinsic alignment power spectrum. Tidal interactions during the formation of galaxies, and other processes (e.g. spin correlations) may induce a preferred intrinsically correlated orientation of galaxy shapes, affecting the two-point shear statistics (Joachimi et al. 2015; Kiessling et al. 2015; Kirk et al. 2015). As a consequence, the observed shear power spectrum will include these additional terms, which cannot be neglected. This effect, referred to as intrinsic alignment (IA), can be considered to be an astrophysical systematic error, contributing to the observed signal and that cannot be removed efficiently by observational strategies. The contribution of IA can be seen by considering the two-point correlation function of Eq. (102), which results in four terms

(106)

where i <  j, represents the correlation between background shear and foreground intrinsic alignment, while represents correlations between foreground shear and background ellipticity. The latter is zero, because a foreground shear should not be correlated with a background ellipticity except if the galaxy redshifts are misassigned. The former, together with the intrinsic-intrinsic (II) alignment autocorrelation power spectrum, can be written in a very similar way to the cosmic shear expression Eq. (103), using a model known as the linear-alignment model:

(107)

The weight function for IA can be conveniently written as

(108)

where the term c/H(z) = dr/dz reflects the choice of integrating with respect to z instead of r.

Intrinsic alignment power spectra can be expressed as a function of the matter power spectrum. There are many models available in the literature for describing the effect of IA (e.g. Joachimi et al. 2015); here we follow a simplified yet observationally motivated approach with

(109)

(110)

where the function ℱ_IA(z) reads

(111)

We refer to this as the extended non-linear alignment (eNLA) model. For ℱ_IA(z) = 1, the model reduces to the so-called non-linear alignment (NLA) model widely used in the literature (Bridle & King 2007). In Eq. (111) ⟨L⟩(z) and L_⋆(z) are the redshift-dependent mean and the characteristic luminosity of source galaxies, as computed from the luminosity function. With respect to NLA, the extended model therefore includes luminosity dependence of the IA, as hinted at by low redshift studies and hydrodynamical simulations (Tenneti et al. 2016; Hilbert et al. 2017; Chisari et al. 2015). The parameters η_IA, β_IA and 𝒜_IA are free parameters of the model and should be determined by fitting the data and/or via carefully designed simulations. Here, we fix them to the following fiducial values: {𝒜_IA, η_IA, β_IA} = {1.72, − 0.41, 2.17}, while 𝒞_IA = 0.0134 remains fixed in the Fisher analysis, as it is degenerated with 𝒜_IA.

Photometric redshifts and number density. One of the key ingredients in Eqs. (105) and (108) is the number density distribution n_i(z) of the observed galaxies in the ith bin. For photometric redshift estimates, this can be written as

(112)

where are the edges of the ith redshift bin. The underlying true distribution n(z) appearing in this expression is chosen to be in agreement with the Euclid Red Book¹⁷ (Laureijs et al. 2011; Smail et al. 1994):

(113)

where with z_m being the median redshift. According to the Euclid Red Book, z_m = 0.9 and the surface density of galaxies is . With this choice, we can set the edges of the 10 equi-populated bins, which are

(114)

with and . We note that, since we use 10 equi-populated bins, the number density in each bin is simply . We note that the minimum redshift range is lower than in Laureijs et al. (2011) where a minimum of 0.2 was used as a conservative limit to avoid potential catastrophic redshift outliers. We use a more optimistic value of 0.001, but note that since the lensing kernel peaks around z = 0.3 the change in inferred parameter errors is small.

The number density n(z) is convolved in Eq. (112) with the probability distribution function p_ph(z_p|z), describing the probability that a galaxy with redshift z has a measured redshift z_p (Kitching et al. 2009). A convenient parameterisation for this quantity is given by

(115)

which allows us to include both a multiplicative and additive bias in the redshift determination of a fraction (1 − f_out) of sources with reasonably well measured redshifts, and a fraction f_out of catastrophic outliers (i.e., systems with severely incorrect estimate of the redshift). By modifying the parameters of this function, one can mimic different cases of interest. Our choice is summarised in Table 5.

We note that we have fixed these quantities in the Fisher matrix estimate and do not explore their impact on the forecast; in a cosmological parameter inference on data, such parameters should be varied and self-calibrated. These are nevertheless not entirely flexible and can be inferred independently, with some errors, for instance with clustering-based methods (Gatti et al. 2018). The subsample of galaxies for which Euclid will provide spectroscopic redshift measurements will also help to calibrate the photometric redshift error distribution. The maximal degradation of the FoMs, when redshift errors are considered as free parameters and are self-calibrated from data, has been recently investigated in Tutusaus et al. (2020).

The shot noise. The uncorrelated part of the intrinsic (unlensed) ellipticity field acts as a shot noise term in the observed power spectrum. This is non-zero for auto-correlation (intra-bin) power spectra, but is zero for cross-correlation (inter-bin) power spectra, because ellipticities of galaxies at different redshifts should be uncorrelated. This term can be written as

(116)

where is the galaxy surface density per bin, and has to be consistently expressed in inverse steradians; is the Kronecker delta symbol; and is the variance of the observed ellipticities whose square root value is given in Table 4 which is taken from Massey et al. (2004). We further discuss sources of errors in Sect. 3.3.4.

The observed cosmic shear power spectrum. Finally in the flat-sky and Limber approximations, we can write the total observed tomographic shear angular power spectrum as

(117)

which, in an expanded form reads

(118)

where the first three terms give the contributions from the cosmological (theoretical) shear, cosmological shear-IA, and IA–IA correlations, respectively. The ij label the tomographic redshift bin combinations into which the redshift distribution of sources have been divided. Each term in Eq. (118) has the same general structure: an integral over redshift of the product of weight functions (depending only on the cosmic evolution) and a power spectrum (probing the growth of structures). The integration formally extends up to the redshift of the horizon, but in practice there is no contribution for z ≥ z_max with z_max the maximum redshift of the source redshift distribution. Moreover, the lower limit must be different from zero, to avoid divergence caused by assumming the Limber approximation in Eq. (103). Therefore, the integration in Eq. (118) is in the range z_min ≤ z ≤ z_max with (z_min, z_max) = (0.001, 2.5) as a reasonable choice.

Compact notation. For convenience we rewrite the constituents of Eq. (118) as

(119)

having defined the following kernel functions

(120)

(121)

(122)

It is worth noticing that the IA kernel function depends on the product n_i(z)n_j(z) so that is expected to be non-zero only over the redshift range where n_i(z) and n_j(z) are both non-vanishing.

3.3.2. Non-linear scales for cosmic shear

The standard cosmic shear power spectrum is sensitive to k-modes in the range 0 <  k ≤ 7 h Mpc⁻¹ (Taylor et al. 2018a), which therefore necessitates a modelling of high-k behaviour of the matter power spectrum. Whilst it is possible to exclude some high-k information by performing a cut in ℓ-mode this procedure is not a “clean” way to remove modes and is inefficient, particularly at low redshift. As shown in Taylor et al. (2018c), it is possible to employ a weight function to exactly remove scales above a particular maximum k-mode, a method known as “k-cut cosmic shear”; however, in this work we wish to present a standard formalism so that results can be more readily compared to the literature, and Laureijs et al. (2011). On a related matter, Copeland et al. (2018) have recently found that baryonic effects on the non-linear power spectrum can decrease the FoM by up to 40%, depending on the choices of model and prior. We therefore need to model the high-k behaviour and we discuss this modelling here.

The linear power spectrum of Eq. (28) is able to match the one measured from numerical simulations only over a limited range in k. For cosmic shear, one needs to model the matter power spectrum deep into the non-linear, high-k, regime that corresponds to projected angular modes in up to multipoles ℓ ∼ 1000, and in fact already important at ℓ ∼ 100 (Taylor et al. 2018a). A common approach to model the non-linear matter power spectrum is to define it as a general function of the linear power spectrum, P_lin,δδ(k,z), that translates it into the full, non-linear power spectrum P_δδ(k,z). Various recipes are available in the literature, each of which has a different performance on different scales. Here, we implement the two currently most popular recipes, which we shall denote as “HF” (for halofit, see Smith et al. 2003; Takahashi et al. 2012) and “HM” (for halo model, see Cooray & Sheth 2002; Mead et al. 2016). These are both analogous formalisms, and both define the so-called 1- and 2-halo terms. The 1-halo term accounts for correlations between dark matter particles within the same dark matter halo, and dominates on small scales (smaller than the size of the halo); the 2-halo term describes correlations between distinct dark matter haloes, dominates on scales larger than the ones where the 1-halo term dominates and is proportional to the linear matter power spectrum multiplied by the effective bias of the species under consideration (e.g. dark matter haloes, galaxies etc.). The two HF and HM methods differ in how the 1- and 2-halo terms are computed. Additionally, both methods depend on parameters that are fitted to reproduce the matter power spectrum measured from state-of-the-art (at the time they were proposed) simulations.

The HF method dates back to Smith et al. (2003), but here we adopt the revised version from Takahashi et al. (2012), that extends the original recipe (designed for a CDM-dominated universe) to models with constant equation of state dark energy. The 1- and 2-halo terms are related to the linear matter power spectrum by a set of empirically motivated fitting functions, whose parameters have been set in order to fit the measured power spectrum from 16 different high-resolution N-body simulations. The refined HF method provides an accurate prediction of the non-linear matter power spectrum over a wide range of wavenumbers k ≤ 30 h Mpc⁻¹ at redshifts 0 ≤ z ≤ 10: the precision is 5% for k ≤ 1 h Mpc⁻¹ at 0 ≤ z ≤ 10 and 10% for 1 h Mpc⁻¹ ≤ k ≤ 10 h Mpc⁻¹ at 0 ≤ z ≤ 3. Since massive neutrinos are included in the analysis, we consider the standard HF method with the extension for massive-neutrinos described in Bird et al. (2012).

The HM method has recently been proposed by Mead et al. (2016) for CDM with a cosmological constant and was later extended to some non-standard scenarios including massive-neutrinos and some modified gravity theories. It takes into account the halo mass function and the halo density profile for the 1-halo term. The total power spectrum is then modified through the use of empirically motivated functions, whose parameters are determined by fits against the outcome of N-body simulations. The resulting non-linear matter power spectrum is accurate to ∼5% up to k ≤ 10h Mpc⁻¹ for 0 ≤ z ≤ 2. This is a smaller range compared to the HF technique, but it is worth stressing that the simulations used have higher resolution so that the overall performance is actually comparable if not better than for HF. A further advantage of the HM recipe is that it automatically includes nuisance parameters that can be used to correct for the impact of baryonic effects on the small scale power spectrum as verified by a comparison to hydrodynamical simulations.

We summarise here the non-linear recipes available at the time of this work, which also include massive-neutrino effects.

• Bird. Bird et al. (2012) modified the non-linear HF version from Smith et al. (2003) to account for massive neutrino effects in the non-linear evolution of the total matter power spectrum. This is a fitting to N-body simulations that include a massive neutrino component as collisionless particles. However, this version is not used in the present work, since it has been updated by the fitting below.

• Mead. Mead et al. (2016) provided a different fitting from the previous ones, to account for massive-neutrino effects at the non-linear level. Also in this case, N-body simulations, which include a massive-neutrino component as a collisionless species, are used, but the parameterisaton is based on a different approach than Bird et al. (2012). The HM method accounting for massive-neutrinos is not adopted in this work.

• TakaBird. As reported in the documentation of the CAMB webpage¹⁸, on March 2014 modified massive-neutrino parameters were implemented in the non-linear fitting of the total matter power spectrum, to improve the accuracy with the updated version of halofit from Takahashi et al. (2012). This is the current version adopted in this work when massive-neutrinos are included. The fitting parameters that account for non-linear corrections in the presence of massive-neutrinos are different from the ones implemented by Bird et al. (2012) in the original HF version from Smith et al. (2003), and described in the point above.

Of the three approaches listed above, we decided to use only the TakaBird implementation since we deem this to be the most reliable one. The Bird et al. (2012) correction to HF explicitly refers to a ΛCDM scenario thus preventing us from exploring deviations from it when constraining the DE (w₀, w_a) parameters. On the other hand, the Mead et al. (2016) recipe, although based on more recent results, has not been tested against simulations which include neutrinos with a mass as small as the one we are choosing as a fiducial value. Hereafter, therefore, only the Takahashi-Bird, TakaBird, power spectrum is used to correct for non-linearities in the matter power spectrum used as input for cosmic shear tomography.

3.3.3. Modelling beyond ΛCDM

Following the discussion outlined in Sect. 2.4, we can recast the growth factor as a function of the growth parameter γ according to

(123)

In the presence of a non-standard growth of structures, as is the case of modified gravity, we effectively deal with a growth factor D_MG(z; γ) that is different from the ΛCDM one, and depends on γ, free parameter of the theory (see e.g. Casas et al. 2016, 2017). However we find that small (percent) deviations from γ_fid ≃ 6/11, lead to changes of less than 0.1% in the growth factor over the redshift range of interest, and so one can safely set D_MG(z; γ) = D(z). Nonetheless, we rescale it as

(124)

where A here refers to one of the observables in the set {δδ, δI, II}. The shear power spectra can now be computed following the cosmic shear recipe detailed above, provided one replaces with its modified version .

In doing this, we are making two assumptions. First, Eq. (124) must be considered as a first order approximation. Indeed, while it is exact in the linear regime, its validity in the non linear one depends on the particular MG model being considered. Indeed, by using Eq. (124), we are implicitly stating that MG only rescales the growth of structures, while scale dependent modifications are actually expected. As a second issue, one should also note that, given the way that the IA spectra P_δI and P_II are related to the matter one P_δδ, we are here assuming that MG does not change this relation, which is an assumption that has not been tested.

3.3.4. Covariance matrix for cosmic shear

The error on the observed cosmic shear angular power spectrum can be expressed as

(125)

where f_sky is the fraction of surveyed sky, and Δℓ is the multipole bandwith. The last term in of Eq. (118) is a Poisson noise term, whereas the first three, that is all terms in Eq. (119), represent the contribution from cosmic variance. Lastly, the expression under the square root in Eq. (125) accounts for the limited number of available independent ℓ modes, with f_sky being the surveyed fraction of the sky. In the case that one measures the power spectrum directly from the data, a four point function should be defined to encapsulate the error on the measurement. In this case Eq. (125) the four point function obtained assuming Gaussianity and applying Wick’s theorem can be recast as

(126)

Both of these are used in the construction of the Fisher matrix.

A critical ingredient in the estimate of the Fisher matrix is the covariance matrix accounting for the errors in the measured observable. For the shear power spectrum, 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. 2013). The non-Gaussian terms involve the convergence trispectrum, and there are large theoretical uncertainties on how to model this quantity, both in the quasi-linear and in the non-linear regime. An approximate analytical way to carry out this calculation is to rely on the halo model formalism (see e.g. Cooray & Hu 2001). Following this approach with the shear specifications given in Sect. 3.3.1, we find that the information content of the shear power spectrum, defined as the total signal-to-noise ratio (Rimes & Hamilton 2005; Sato et al. 2009; Takada & Jain 2009; Kayo et al. 2013; Tegmark et al. 1997), decreases by 30% at ℓ_max = 5000 when we add the non-Gaussian contribution. This loss of information content corresponds to an effective cut at ℓ_max = 1420 in a forecast that only uses the Gaussian covariance given by Eq. (125). Motivated by this result and in order to avoid uncertainties due to non-Gaussian modelling, we decide not to implement full non-Gaussian terms in our analysis and instead investigate the impact of a cut at ℓ_cut = 1500 (defined as our pessimistic choice) in our Gaussian forecasts. The impact of a different choice in this scale cut (e.g. with respect to the alternative value of ℓ_cut = 5000, corresponding to our optimistic choice) is discussed in our numerical results in Sect. 5.

We note that our approach is practical, but has limited applicability. For example an equivalent ℓ-cut that preserves signal-to-noise is not directly linked to constraints on parameters beyond a global amplitude; parameter sensitivities are generally scale-dependent (Copeland et al. 2018) and therefore respond differently to small-ℓ cuts and the correlation introduced by non-Gaussian covariances. However, the amplitude of the impact of non-Gaussian modelling used here is small compared to other assumptions made.

3.3.5. Applying the Fisher matrix formalism to cosmic shear

For cosmic shear the mean of the shear field is zero, therefore in Eq. (43) one can either use the first term on the right-hand side where the signal is in the covariance, or one can redefine the mean as being the power spectrum itself. In the former case only the covariance (power spectrum) appears in the expression, while in the latter case the covariance of the covariance (i.e. the four point function) appears. If one were measuring the cosmic shear signal from data the former would represent a measurement of the spherical harmonic coefficients, and a likelihood construction where the signal was in the covariance; the latter would then represent the case that the power spectrum was measured directly. We define both here, since two of the codes tested cover both these choices.

In the case that one assumes the signal is the covariance (i.e. a measurement of the spherical harmonic coefficients) the Fisher matrix reads

(127)

where C^ϵϵ(ℓ) is defined in Eq. (118). No cross-correlation between modes with different ℓs is included, which is true under the Gaussian covariance assumption. We shall discuss the quantities ℓ_min, ℓ_max, and Δℓ in Sect. 4.3. Here is the ij element of the inverse of the matrix defined in Eq. (125).

In the case that one assumes the signal is the mean power spectrum (i.e. a measurement of the power spectrum coming directly from data) the Fisher matrix reads

(128)

Both of these definitions are used in the code comparison, and we find agreement between the two approaches; although this may be expected to hold analytically in any case (Carron 2013).

Derivatives. What is needed to compute the Fisher matrix elements is then not only the tomographic matrix of the ellipticity signal, C^ϵϵ(ℓ), but also its derivatives with respect to the cosmological (and nuisance) parameters. Using the compact notation introduced in the previous section Eq. (119), we can now compute the shear power spectrum derivatives:

(129)

where {A, a}∈{{δδ, γγ},{δI, Iγ},{II, II}}. We note that the first term is only present when differentiating with respect to the background parameters {Ω_m, w₀, w_a, h}, whilst in the second term the derivatives of the power spectrum may be further decomposed into

(130)

where we remind the reader that k_ℓ(z) = (ℓ + 1/2)/r(z). Again, we note that the second contribution is only present for derivatives with respect to {Ω_m, 0, w₀, w_a, h}. Finally, we stress that, different from the galaxy clustering case, the lensing Fisher matrix is directly computed with respect to the cosmological parameters of interest, so that there is no need for any projection.

3.4. Recipe for the combination of galaxy clustering and weak lensing

While the use of complementary probes is well known to be essential to provide tight limits on cosmological parameters, the joint analysis of correlated probes has been relatively little explored until now, and thus represents a rather new field in cosmology. In this section we describe our approach for combining galaxy clustering (from the spectroscopic redshift survey and photometric survey) with weak lensing for Euclid.

For the specific combination of weak lensing and galaxy clustering coming from the photometric redshift survey (GC_ph), we can use the angular power spectrum formalism described in Sect. 3.3.1: this allows us to show the impact of including or not cross-correlation terms (for every redshift bin) between WL and the GC_ph.

3.4.1. Observables for the photometric survey and probe combination

The maximum set of observables that may be used for probe combinations are all individual Euclid probes, the cross-correlation (XC) between those probes (which can be used either as observable or as covariance), and correlations between external data and Euclid probes. In the latter case the external data can also be used as individual probes in the combined set.

For this work however, we focus on the combination of only cosmic shear and galaxy clustering, the latter using both the photometric and spectroscopic galaxy surveys. In addition to simply combining these probes, we also include in our analysis the cross-correlation of photometric galaxy clustering and cosmic shear. Both latter probes are quantified as 2D projected observables, which allows us to compute the cross-observables and cross-covariance matrix in a simple, unified scheme. The cross-correlation of cosmic shear with spectroscopic galaxy clustering is expected to be small given that the spectroscopic Euclid survey selects galaxies at high redshift (z >  0.9, see Sect. 3.1.5) and thus has a small overlap with the lensing kernel. However, when we combine the photometric and spectroscopic galaxy clustering probes, the two galaxy distributions overlap and, therefore, we consider the following two different cases:

• In the pessimistic case we only consider galaxies below z = 0.9 in the photometric sample, and neglect the scatter from high redshifts into the photometric sample due to photometric redshift errors: these probes are uncorrelated, at the price of discarding some parts of the observations.

• In the optimistic case we do not perform any cut in redshift space for the photometric sample and we combine the spectroscopic and photometric galaxy clustering probes neglecting any cross-correlation that might appear between them.

For the photometric survey we use the specifications of Table 4 and the galaxy distribution function n(z) of Eq. (113), convolved with the probability distribution function p_ph(z) of Eq. (115). For the spectroscopic redshift survey we use instead the specifications of Table 2 and we follow the approach described in Sect. 3.2 to model this observable.

We first denote the angular density contrast of galaxies in redshift bin i by , the shear field of source tomography bin j as γ^j, and the normalised galaxy number density in each tomographic bin i of the photometric survey as n_i(z), given by Eq. (112). Notice that we assume here the same galaxy distribution for the photometric shear and galaxy clustering surveys; this is not guaranteed to be the case for the real photometric survey, since the sampled galaxies for GC_ph might be a subset of the full galaxy sample. Therefore, more realistic galaxy distributions are still being investigated.

We can now define the radial weight function for galaxy clustering as

(131)

where H(z) is the Hubble parameter defined in Eq. (3), and b_i(k, z) is the galaxy bias in tomography bin i. For the latter, we assume a constant value in each redshift bin and no dependence on the scale k; the values of b_i in each bin are considered as nuisance parameters over which we marginalise when producing our results. Their fiducial values are given by

(132)

with the mean redshift of each bin. The main assumptions made on the galaxy bias are therefore that it appears linearly in Eq. (131) and that it does not depend on the scale k. Removing these assumptions would require the addition of extra nuisance parameters to control the uncertainty on the galaxy bias and would lead to less constraining power from the photometric galaxy clustering survey. Given these assumptions, in the following we drop the k dependence from the expression of .

We use the same radial weight function for cosmic shear given by Eq. (105), and we include the effect of intrinsic alignment described in Sect. 3.3.1; combining Eq. (118) with Eqs. (109) and (110), we can define a cosmic shear power spectrum as

(133)

where we have defined the overall weight function , which includes intrinsic alignment contributions

(134)

with defined in Eq. (108). This is exactly equivalent to Eq. (119), except that we do not include the noise term for clarity. We change notation here to make the relationship between the galaxy clustering power spectra and the cosmic shear power spectrum clearer. Finally, using the Limber approximation, the other observables are given by

(135)

(136)

and represent, respectively, the cosmic shear power spectrum, the shear-galaxy power spectrum and the galaxy-galaxy power spectrum (using photometric galaxy clustering).

3.4.2. Non-linear scales for probe combination

For the modelling of non-linear scales in the cosmic shear, photometric galaxy clustering, and cross-correlations observables, we proceed as discussed in Sect. 3.3.2. We apply it for cosmic shear up to multipole ℓ_max = 5000 for the optimistic case and ℓ_max = 1500 for the pessimistic case (see Sect. 4.4.2 for the detailed description of the optimistic and pessimistic cases presented in the results). For the photometric galaxy clustering and cross-correlations, we consider multipoles up to ℓ_max = 3000 for the optimistic case and ℓ_max = 750 for the pessimistic case.

3.4.3. Modelling beyond ΛCDM

In order to calculate the DE FoM we extend our cosmological model beyond ΛCDM by adding a dark energy fluid parameterised by two equation of state parameters, w₀ and w_a, as was done in Sects. 3.2.3 and 3.3.3. However, when we take cross correlations into account we do not consider cosmological models with deviations from GR, that is to say we do not consider the γ parameterisation of the growth of structure that is considered for single probes. There are two main reasons for this. The first one is that the codes used in this work to include cross correlations into the analyses are coupled to Boltzmann codes (contrary to the codes for single probes, which use external input spectra) and they obtain the growth of structure directly from them. Therefore, it is not possible for these codes to include the γ parameterisation by only modifying the growth of structure and rescaling the standard matter power spectrum, as is done for single probes. The second main reason is that the public version of the Boltzmann code allowing for the γ parameterisation, MGCAMB, does not allow for this parameterisation together with an evolving dark energy equation of state, which is the case considered here. As a consequence, we do not consider models beyond GR for cross-correlations in this work.

3.4.4. Covariance matrix for probe combination

The contribution of spectroscopic galaxy clustering is added as a separate Fisher matrix (both in the optimistic and the pessimistic cases) and hence additional covariance terms are not included.

In the case of cosmic shear, photometric galaxy clustering and their cross-correlation, a joint covariance matrix is required, that can be defined in two different ways (see also Sect. 3.3.4); we present both and their outcome is compared in Sect. 4.4. One can define a covariance as

(137)

where A, B run over the observables L and G, Δℓ is the width of the multipoles bins used when computing the angular power spectra, and i, j run over all tomographic bins. Otherwise, one can define the fourth-order covariance as

(138)

As above, A, B, A′,B′ = L, G and i, j, k, l run over all tomographic bins. In both approaches, the noise terms take the form

(139)

(140)

(141)

with and introduced in Sect. 3.3.1. Notice that we assume here that the Poisson errors on cosmic shear and galaxy clustering from the photometric surveys are uncorrelated, which yields .

As for the cosmic shear case (see Sect. 3.3.4), the covariance matrix of photometric galaxy clustering can be split into the sum of Gaussian and non Gaussian contributions. The non Gaussian terms involve the galaxy trispectrum, whose modelling beyond the linear regime suffers from significant theoretical uncertainties. As for cosmic shear, we rely on the halo model formalism (Lacasa & Rosenfeld 2016; Lacasa 2018) to estimate the impact of the non Gaussian contribution to the photometric galaxy clustering covariance matrix. It is important to note that the photometric galaxy clustering is more non Gaussian than shear due to sampling high density regions instead of total matter, being integrated over a smaller volume (bin width versus whole light cone), and because shot noise is subdominant for clustering while it has a high impact on the shear Gaussian errors. Following this approach with the galaxy specifications given in Sect. 3.4.1, we find that the signal-to-noise ratio saturates for ℓ ≳ 500 when non Gaussian contributions are included. We mimic this effect using a conservative cut on the multipoles used for GC_ph, removing the contributions of multipoles ℓ >  750 from the analysis (pessimistic setting). We still show the impact of varying this cut by testing an ideal (unrealistic) case in which non-Gaussian contributions are entirely neglected up to ℓ_max = 3000 (optimistic setting).

3.4.5. Applying the Fisher matrix formalism to probe combination

The final expression of the combined Fisher matrix for the angular power spectra, including the contributions of photometric galaxy clustering, cosmic shear, and their cross-correlation, is given in the case of the fourth order covariance by

(142)

and in the case of the second-order covariance as

(143)

The block descriptors A, B, C, D run over the combined probes L and G, thus including the three observables described in Eq. (135), that is to say cosmic shear auto-correlation, cross-correlation between galaxy clustering and cosmic shear galaxy clustering auto-correlation. The indices in the sum ij and mn run over all unique pairs of tomographic bins (i ≤ j, m ≤ n) for the cosmic shear auto-correlation and the galaxy clustering auto-correlation, while all pairs of tomographic bins are considered to take into account all the cross-correlations between galaxy clustering and cosmic shear.

Finally, in addition to this combination, we include also the contribution of the spectroscopic galaxy clustering. Since we are neglecting its correlation with the photometric survey observable (both in the pessimistic and optimistic cases), we simply add the spectroscopic Fisher matrix computed as described in Sect. 3.2. Therefore, our final matrix is

(144)

4. Code comparison

In this section we describe the code comparison procedure and results for Euclid’s main cosmological probes, namely spectroscopic galaxy clustering and weak gravitational lensing, as well as the cross-correlation of photometric galaxy clustering with weak lensing. We have compared and validated six different Fisher matrix codes for galaxy clustering, seven for weak lensing, and five of the codes also implement cross-correlations. All the codes implement the recommended recipes for GC_s, WL, GC_ph and XC discussed in Sect. 3. For clarity and in order to isolate and compare different cases of interest (e.g. linear modelling versus non-linear corrections, Fisher matrix projection to different parameters, including dark energy and modified gravity parameterisations) we have validated the codes across various intermediate cases of increasing complexity.

This section aims at describing the implementation and code comparison among the different codes and represents our main deliverable: we include detailed instructions in Appendix A for the reader who may want to validate another Fisher Matrix forecast code; we make publicly available the recipe used, the reference output Fisher matrices, and plotting routines to reproduce your own comparison plot, in the form of Fig. 3, to check that your own code is validated. In Sect. 4.1 we describe the codes used in our comparison, in the WL, GC and XC cases. As the comparison was done among codes which had been already partially developed before this work, and was motivated by several codes providing different results in both amplitude and orientation, code ownership belongs to the authors who developed them. We have however encouraged all code owners to make these codes publicly available: CosmoSIS and CosmoMathica (used to produce the input files for this work) are already publicly available, the authors of four codes (FisherMathica, SOAPFish, CCCPy, TotallySAF) have expressed intention to make them publicly available, five codes (BFF, fishMath, CosmicFish, CarFisher and STAFF) can be available upon request. Requests can be made directly to the authors, as listed below in Sect. 4.1.

In Sects. 4.2 and 4.3 we describe specific implementation steps regarding the galaxy clustering and weak lensing code comparisons, in particular the construction of the background cosmological quantities, the implementation of the survey geometry and galaxy sample properties, and the perturbed quantities, in order to construct the observed power spectrum. We give specific and detailed information on the parameters we vary, and the calculation of the derivatives entering the Fisher matrix. We also describe how we implement non-linear corrections and how we project our Fisher matrices to derive constraints on commonly used dark energy and gravity parameters, such as w₀, w_a and γ. We define the different cases we use to validate our codes and finally show the results of the code comparison, that is to say how well the codes match. In Sect. 4.4 we repeat this procedure for the probe combination case. Overall, we find very good (percent level) agreement between the different codes for all cases. In Sect. 4.5 we comment on critical lessons we learnt during the code comparison validation, and give some advices for writing a Fisher matrix code for cosmology from scratch.

4.1. The codes used in the comparison

Here we describe each of the codes used in the comparison, for the cosmic shear, galaxy clustering and cross-correlation observables. The comparison started from codes that were already available at the time of this paper, summarised in Table 6. There are a variety of languages and approaches used, some are privately written codes and others use publicly available packages. We further provide in this paper instructions and reference input/output files to proceed with the validation of any other external code the reader may have. We emphasise that these codes were validated in their current versions, as of February 2019, and with regard to their specific Fisher matrix implementation and output.

BFF (BONN FISHER FORECAST) is a code developed by Victoria Yankelevich and is written in Fortran 90. It was designed to comply with the IST specifications by following the guidelines. However, it is based upon the code developed in the higher-order statistics work package to make forecasts for the galaxy bispectrum (and its combination with the power spectrum; see Yankelevich & Porciani (2019). Due to technical requirements, these two codes have different structures but give consistent results. The main advantage of BFF is computational time, which is about a couple of minutes on a one core laptop.

CarFisher is a Fisher matrix code written by Carmelita Carbone in IDL, which computes cosmological parameter forecasts for a 3D galaxy power spectrum. It has been used, in the linear regime version, to provide results in Laureijs et al. (2011) and Carbone et al. (2011a, b); Carbone et al. (2012). Now it also includes the non-linear modelling described in this work. It interfaces with CAMB to produce matter power spectra files, and with a F90 routine for the computation of scale independent growth factor and growth rate. The projection onto the final parameter space is performed by a separate module, that enables one to choose different settings (e.g. k-, z- or μ-intervals) without the need to rerun the Boltzmann code for power spectra production.

FisherMathica is a Fisher matrix code written by Santiago Casas in Mathematica in a modular way, including packages for different applications. Based on the code used in Amendola et al. (2012), it has been used for forecasts of modified gravity and dark energy (Casas et al. 2016, 2017, 2018).

The particular advantage of this code is flexibility in terms of cosmological parameters. It accepts any basis and can convert between them, such as A_s versus σ₈ or ω_m, 0 versus Ω_m, 0. It can accept also new parameters, which can be useful for extended models of modified gravity. The CosmologyTools package takes care of computing Hubble functions, distances, coordinate and unit transformations, analytical growth rate equations and other cosmological observables. The publicly available CosmoMathica package¹⁹ deals with the interaction with common codes in the cosmological community, such as the Boltzmann codes CLASS and CAMB and its modified gravity version MGCAMB, Eisenstein & Hu (1998) transfer functions, the Coyote Cosmic Emulator (Heitmann et al. 2014) and other fitting formulas for non-linear power spectra. CosmoMathica has been used to produce the input files for this work. The FisherGC package takes care of the galaxy clustering Fisher matrix, by computing the derivatives of the power spectrum with respect to the shape parameters and with respect to the redshift dependent parameters. It reads the GC specifications for different surveys from given text files and computes all needed quantities for the Fisher matrix. The integration settings are very flexible and are the ones provided by the Mathematica NIntegrate function. The FisherWL package computes the kernels, the window functions and the C(ℓ) functions for weak lensing. It also includes several models of intrinsic alignment. The Fisher matrix is calculated in a parallelised way for each of the parameters and redshift bins. The UsefulTools package takes care of input/output, creation of tables, Fisher ellipses plots and post-processing. It contains many useful matrix and vector operations that are needed to analyse and interpret Fisher matrices results. Finally, the code can be called interactively using Mathematica notebooks or from the command terminal, using .m file scripts.

fishMath is a code that computes galaxy clustering forecasts, written in Mathematica by Elisabetta Majerotto. A previous version of it, which did not include forecasts on shape parameters, was used in various papers to produce forecasts (Laureijs et al. 2011; Majerotto et al. 2012; Bianchi et al. 2012), and more recent versions of it have been used to make various cross-checks for other papers (Sapone et al. 2013, 2014; Majerotto et al. 2016). The code allows one to switch from linear modelling to the two different parameterisations of non-linear modelling. It consists of different modules, so that it is possible, for example to start by producing the Fisher matrix or to read the Fisher matrix from a file and marginalise it or project it over the preferred parameters. It also computes errors on parameters and the FoM for dark energy.

SOAPFish is a code implemented using a combination of Python 3 and Mathematica 10.0.0.0, by Domenico Sapone. It is a modified and improved version of the code used in Sapone et al. (2009, 2010, 2013, 2014), Majerotto et al. (2016), Hollenstein et al. (2009), Bueno Belloso et al. (2011) and Sapone & Majerotto (2012), where the authors forecasted the sensitivity with which the next generation of surveys will measure the cosmological parameters with different cosmological models. The code is divided into cells and three main parts. The first one specifies the values of the cosmological parameters and the survey properties, such as survey area, number of bins and error on redshift. Once these input quantities have been called, the code evaluates the cosmological functions, such as the Hubble parameter and the angular diameter distance, and the survey specifications, such as the bin volumes and the galaxy number density in each bin. The third and last part consists of evaluating the derivatives of the observed galaxy power spectrum in a fully and optimised numerical manner from which the Fisher matrices are evaluated.

TotallySAF (Totally Super Accurate Forecast) is a code that has been implemented using Python (compatible with Python 2 and Python 3) and C++, developed by Safir Yahia-Cherif and Isaac Tutusaus. The code consists on 2 different packages: SpecSAF and XSAF. The former is a modified and improved version of the code used in Tutusaus et al. (2016), where forecasts for the spectroscopic galaxy clustering probe of Euclid have been computed for theoretical models with varying dark matter and dark energy equation of state parameters. In order to compute the integrals this code implements a quick and precise method that consists of making sums over four cells in the k and μ space (the one under consideration, the right nearest-neighbour, the down nearest-neighbour and the right-down nearest-neighbour); this whole part of the code is parallelised. Concerning XSAF, it allows for the computation of weak lensing and photometric galaxy clustering forecasts, taking into account their cross-correlations. Another feature of this code is that the user can choose higher point numerical derivatives than the 3 and 5 points stencil. The 7, 9, 11, 13 and 15 point stencils are available for users who require high precision calculation. For each parameter used it is possible to choose the number of points and the step of the derivatives. The code can also provide the input files and produce the visualization plots.

CCCP (Camera’s Code for Cosmological Perturbations) is a code developed by Stefano Camera and has been used as a standalone code, or in order to compute either angular power spectra only or Fisher matrices only, in various analyses (e.g. Camera & Nishizawa 2013; Camera et al. 2015a,b,c; Bonaldi et al. 2016). At the time of the present code validation, it is written in Mathematica 11.2.0.0 and it is kept updated from version to version. It consists of two main routines: (i) one calculating cosmological background and perturbation quantities for a given set of cosmological parameters; and (ii) another computing the Fisher matrix for the parameters of interest. The former routine employs, whenever is possible, analytical expressions to avoid interpolations and the propagation of numerical errors. Nonetheless, CCCP also allows one to load the growth, the transfer function, the power spectrum or all functions from an external file, and it can be easily interfaced with a Boltzmann solver to obtain these quantities on the fly. Derivatives of angular power spectra with respect to cosmological parameters are computed according to the method discussed here in Sect. 4.5.2, and first proposed by Camera et al. (2017). It is found that it is in good agreement with MCMC results. Regarding the latter routine, that responsible for Fisher matrix computations, it can be used either with input data provided by the former routine, or with external spectra in form of tables of multipoles and values as inputs. The stability of the Fisher matrix is kept under control via several tests, and it is enforced, when necessary, for example by considering the logarithm of the parameter instead of the parameter itself: this reduces the dynamic range between the maximum and minimum eigenvalue, thus stabilising matrix inversion (see Fonseca et al. 2015). Despite the reduction of the dynamic range, matrix inversions may still propagate numerical instabilities into the resulting Fisher matrix. Therefore, CCCP checks the stability of the inverse matrix F⁻¹ (computed by default via one-step row reduction) by comparing it to the Moore–Penrose pseudo-inverse and to the inverse obtained via eigen-decomposition (Albrecht et al. 2006; Camera et al. 2012, 2018).

CCCPy was born as a Python version of CCCP. In its state at the time of the publication of this paper, it is based on the Core Cosmology Library (CCL, Chisari et al. 2019), a standardised library of routines to calculate basic observables used in cosmology. Specifically, CCCPy employs CCL for the computation of the various angular power spectra (WL, GC_ph, and XC), adjusted so to implement the same recipes as discussed in the previous sections, like binning, galaxy bias, and IA. Moreover, it has utilities implementing the SteM derivatives, with and without accuracy checks. It can be run as a single jupyter notebook, from the computation of the observables, to the derivatives, to the final Fisher matrices and the corresponding parameter constraints.

CosmicFish is a Fisher matrix code written in the Fortran90 and Python languages by Marco Raveri and Matteo Martinelli. It is interfaced with CAMB_sources and its generalizations are based on common CAMB modifications, such as MGCAMB (Zhao et al. 2009; Hojjati et al. 2011) and EFTCAMB (Hu et al. 2014; Raveri et al. 2014). CosmicFish uses as primary parameters for ΛCDM cosmology the set {Ω_b, 0h²,  Ω_c, 0h²,  h,  n_s,  log10¹⁰A_s,  (w₀,  w_a,  γ)} with the parameters in parentheses added when considering extensions to ΛCDM. Numerical derivatives are computed with with respect to these parameters, with the options of either 3 or 5 points stencil. Finally it projects the results on the common parameter basis combining the original Fisher matrix with a Jacobian matrix. CosmicFish has been used in previous works to investigate the information gain given by past and future cosmological surveys (Raveri et al. 2016a), as well as to obtain forecasts for proposed CMB missions (Abazajian et al. 2016). Also, a preliminary validation of the code was done reproducing results available in the literature at the time of the code release (Raveri et al. 2016b). The version validated in this paper contains key features and calculations not available in the current public code, e.g. the effect of intrinsic alignments and the possibility of working in the common parameter basis. These new features will be included in the next update of the public version of this code.

CosmoSIS is a publicly available cosmological parameter estimation code²⁰ (Zuntz et al. 2015). It has a modular structure and this enables the use of existing codes for observable predictions, and different experimental likelihoods, written in different languages, like Fortran90, Python and C/C++. The public version of CosmoSIS provides a standard library with different modules for the Boltzmann solver, like CAMB (Lewis et al. 2000; Howlett et al. 2012), CLASS (Lesgourgues 2011a,b; Blas et al. 2011; Lesgourgues & Tram 2011), or MGCAMB (Zhao et al. 2009; Hojjati et al. 2011), the background quantities, or nuisance parameters (like the galaxy bias, and intrinsic alignments). More importantly for us, it also contains a sampler that computes Fisher matrices. CosmoSIS runs as a pipeline. Some modifications have been made to the standard implementation by Isaac Tutusaus and Martin Kilbinger for the purposes of this paper, and the pipeline used works as follows: the Boltzmann solver used (CAMB) uses A_s instead of σ₈, fixing A_s to a dummy value and adding a third module into the pipeline to rescale the matter power spectrum. Afterwards the fourth module is added to take into account the non-linear correction to the matter power spectrum, and the fifth to consider the intrinsic alignments. Two more modules are needed to load the distribution of galaxies, n(z), and apply the galaxy bias model to the power spectrum. The projection to the C_ij(ℓ) values is performed for each of the considered observables, adding the intrinsic alignment power spectra and finally the full covariance matrix for the observables is computed. Once this is done, the Fisher matrix sampler calls the pipeline iteratively for all the variations of the parameters required to compute the numerical derivatives, and finally the Fisher matrix. Two derivative methods are available (3-point and 5-point stencils) with a user-defined step size per parameter. CosmoSIS has been widely used in the literature to analyse real data (Abbott et al. 2018a,b; Baxter et al. 2019), and to compute forecasts (Harrison et al. 2016; Park 2016; Olivari et al. 2018).

STAFF (Shear Tomography Accuracy from Fisher Forecast) is a Mathematica code developed by Vincenzo F. Cardone for the code comparison challenge. It comes as a single notebook (staff) performing all the necessary steps to output the final Fisher matrix, or as three separate notebooks for intermediate computations, namely pscalc for the estimate of the matter power spectrum and its derivatives, cijdercalc to obtain the shear, photometric galaxy clustering, and cross-correlation power spectra and derivatives, and fmcalc for the Fisher matrix evaluation. The main virtues of STAFF are its simplicity and versatility which allows the user to quickly modify the code for any non standard cosmological model of interest. Unfortunately, being written in Mathematica, the code is slower than some others, and hence it is not ideal from this point of view. Nevertheless, the most time consuming step is when computing the matter power spectrum and its derivatives and this operation needs only to be performed once. This is because since the semi-analytical approach implemented for the derivatives of the different power spectra of interest make it possible to express all of them as operations on the matter power spectra.

Table 6 presents a summary the forecasting codes involved in the comparison, including:

• name: label or name of the code;

• contact: contact person(s) for the corresponding code;

• GC_s: whether the code has been validated to forecast galaxy clustering from the spectroscopic survey;

• GC_ph: whether the code has been validated to forecast galaxy clustering from the photometric survey;

• WL: whether the code has been validated to forecast cosmic shear weak lensing from the photometric survey;

• XC^{(GC_ph, WL)}: whether the code has been validated to forecast cross-correlations between galaxy clustering and weak lensing, both from the photometric survey;

• curvature: whether the code has been validated to forecast non-flat models;

• Boltzmann code: which Boltzmann code(s) the code can interface with (if wished);

• external input: whether the code can interface with externally provided fiducial observables and derivatives;

• language: programming language(s) used in the code;

• public: whether the code is public at the time of the release of this paper.

Since the objective of this paper is to present validated forecasts that are independent of the code used, we do not express a preference for any of the codes described here. Users can decide which code to use based on their specific needs, considering aspects such as availability, cosmological probes considered, programming language, or any other specific feature(s).

4.2. Galaxy clustering code comparison

In this section we describe the procedure to compare the cosmological forecasts for GC measurements from the Euclid spectroscopic galaxy sample that are predicted by the different codes described in Sect. 4.1. All codes employ the recipe described in Sect. 3.2. The Fisher matrix for GC is given by Eq. (91) and it depends on the model for the observed galaxy power spectrum given by Eq. (87), and on the characteristics of the Euclid survey encoded in the effective volume, given by Eq. (90). Here, we describe additional issues to take care of when implementing this recipe. If the reader wishes to compare and validate their own GC Fisher matrix code for Euclid GC, they should follow these instructions in order to do so. As part of the code-comparison process, all codes were adapted to use common inputs. If the reader wishes to follow suit, we make these inputs, alongside our outputs, available along with this paper, describing them and providing instructions for their use in Appendix A.

We first comment on specific issues within the implementation, focussing particularly on the derivatives needed for the Fisher matrix. We then prepare a list of comparison cases with increasing complexity.

4.2.1. Implementation of the galaxy clustering recipe

The first step in the implementation of the GC recipe described in Sect. 3.2 is to construct background unperturbed quantities that enter into the Fisher matrix calculation. The Hubble parameter, H(z), the angular diameter distance, D_A(z), and the comoving volume, V(z_i, z_f), are given by Eqs. (3), (13) and (14), respectively. These functions should be general enough to allow for deviations from the standard flat ΛCDM cosmology, including non-zero curvature, dynamical DE models parameterised as in Sect. 2.1, and deviations from GR as described in Sect. 2.4.

Once these background functions are defined, we can proceed with the estimation of quantities describing the Euclid survey that are required to compute the effective volume V_eff(k_ref, μ_ref; z), of each redshift bin, defined in Eq. (90). To this end, we first evaluate their comoving volumes, V_s(z), using Eq. (14). This requires the knowledge of the area of the survey, which is given in Table 2. As a second step we evaluate the galaxy number density n(z), computed using Eq. (88) and the fractional number densities listed in Table 3.

We must then calculate the quantities that define our observable, namely the two-dimensional galaxy power spectrum, P_obs(k_ref, μ_ref; z), in Eq. (87). This is achieved using the full non-linear model discussed in Sect. 3.2.2. In particular, the following quantities are required to compute the observed power spectrum (see Sect. 3.2 for more details):

• fσ₈(z) and bσ₈(z) enter through the modelling of RSDs. The former is evaluated using the growth rate, f(z), defined in Eq. (32), and the values of σ₈(z) from our input files (see Appendix A.2 for details). The latter can be computed using the galaxy bias factors listed in Table 3 and the same values of σ₈(z).

• The pure matter power spectrum normalised by , . Although these redshift-independent ratios are again given as input files (see Appendix A.2), they can be estimated directly from the output of any Boltzmann code (e.g. Lewis et al. 2000 or Lesgourgues 2011a).

• The “no-wiggles” power spectrum, which we obtain using the formulae of Eisenstein & Hu (1998). The necessary ratios of are treated as inputs to the codes (see Appendix A.2).

• The absolute error in the distances expressed in terms of the redshift errors σ₀ = 0.001, given by Eq. (75).

• The residual shot noise term is represented by P_shot(k, z), which is set to zero at the reference cosmology. That is, our fiducial case assumes a perfect Poisson shot noise subtraction.

Our initial set of cosmological parameters is given in Eq. (69). At the linear level, they are divided into four redshift-independent shape parameters and five redshift-dependent parameters. The corresponding Fisher matrix has a dimension n_tot = 4 + 5 × N_z, where N_z is the number of redshift bins considered in the survey. When considering the full non-linear model, we add two parameters to the set and hence two rows and columns to the Fisher matrix. We then project these initial “model-independent” matrices into the final, model-specific parameter spaces following Sect. 3.2.5.

A note about derivatives. Derivatives of log₁₀P_obs(k, μ, z) with respect to cosmological parameters translate the covariance matrix of the power spectrum into that of the model parameters. Uncertainties in the derivatives can have a significant impact on the final results (as they enter the Fisher matrix that has to inverted). This is particularly true for derivatives that depend on numerical estimates of the slope of the power spectrum around the BAO wiggles. It is therefore imperative to pay close attention to the implementation of the numerical derivatives and to ensure that the parameter steps used lie well within the convergence region. More specifically, the parameter steps used should be small enough to be able to accurately represent the dependence of the underlying function (the power spectrum) on the various parameters, but large enough to avoid machine precision errors. We have checked this by ensuring we got stable results when the step size is varied in the range 0.0001−0.01 and with the different codes.

As already discussed, most of the derivatives of the observed galaxy power spectrum are calculated numerically. We start by considering the derivatives of the matter power spectrum with respect to the shape parameters. To this end, we use input files corresponding to matter power spectra computed for different parameter values. We evaluate the matter power spectrum at the fiducial cosmology and at individually incremented parameter values, keeping all other parameters fixed. The relative increments with respect to the reference value of a given shape parameter, θ_ref, are computed as θ_ref(1 + ϵ) and θ_ref(1 − ϵ). These matter power spectra are computed for all redshifts bins. In total, we then have N_z matter power spectra for the reference cosmology and 2 × N_z matter power spectra corresponding to the relative increments.

We evaluate derivatives with respect to the shape parameters using three-point rules:

(145)

where θ_shape is any of our four shape parameters, and for simplicity we wrote . We remind the reader that the only quantity depending on the shape parameters is the power spectrum with the wiggles removed, which means that the μ dependence cancels out.

Similarly for the two non-linear parameters, we also use 3-point derivatives to calculate the following:

(146)

We use different approaches for the derivatives with respect to redshift-dependent parameters. The derivative of the observed galaxy power spectrum with respect to the residual shot noise can be computed analytically as the inverse of the observed galaxy power spectrum, see Eq. (99). For the rest we use numerical rules at either 3 or 5 points in the parameter space.

For derivatives with respect to ln fσ₈(z) and ln bσ₈(z) we use the 3-point rule:

(147)

(148)

For ln H and ln D_A we choose 5-point derivative rules in order to better capture the oscillatory behaviour of the galaxy power spectrum near the BAO, which is modified by the AP distortions via these two parameters²¹. The derivatives of the galaxy power spectrum with respect to the Hubble parameter and the angular diameter distance are:

(149)

(150)

Here we omitted the direct dependence on k and μ for simplicity. It is important to remember that also the wavenumber k and the direction μ are modified through the cosmology (due to projection effects). This modification is taken into account using Eq. (78), where the subscript “ref” refers to the reference cosmology. It is clear from Eq. (78) that k and μ will differ from the corresponding reference values when H_ref(z)/H(z) or D_A(z)/D_A,ref deviate from unity, that it is when we take derivatives of the observed galaxy power spectrum with respect to ln H(z) and ln D_A(z). We would also like to notice that the independent variables we use are in fact logarithms of H(z) and D_A(z), therefore the input power spectra correspond to parameter increments given by, for example

(151)

Therefore the relative increment on the variable H(z) appearing explicitly in the observed galaxy power spectrum is ; the same discussion applies to D_A(z).

It is clear from the above expressions that the derivatives of the observed galaxy power spectrum need to be evaluated at each redshift bin.

4.2.2. Settings definition

In order to compare the outputs of the different forecasting codes in a controlled way, we consider a series of settings with increasing complexity, which we describe in this section. We compare our codes at each step and provide the resulting Fisher matrix files in Appendix A for the reader who wishes to perform the comparison with their own code. We recommend taking these steps as they have been optimised for easy de-bugging. We now describe each of these settings in detail. A summary of all settings is presented in Table 7.

Linear setting. In this case we only consider the observed power spectrum for scales k <  k_max = 0.25 h Mpc⁻¹ for all redshift bins and use the linear-theory prediction for P_obs(k, μ). This is achieved by fixing the non-linear parameters σ_p = σ_v = 0 in the non-linear recipe of Eq. (87), which is then reduced to its linear version of Eq. (80). We must compute derivatives of this model with respect to the remaining cosmological parameters. We note that, when we estimate derivatives with respect to ω_m,0 by varying the matter density, we must also vary the energy density of the remaining component, Ω_Λ, 0 in this case, in order to keep the cosmology flat. This does not mean that we cannot later project into a new parameter set allowing for curvature, since non-flat models should still be allowed, which we ensure by keeping D_A(z) and H(z) free. We use as input files the ratio , (calculated using CAMB), as well as for the growth rate f(z). This is described in detail in Sect. 4.2.1.

The Fisher matrix that we obtain has dimensions 4 + 5 × N_z, composed of the derivatives of the four shape parameters given by Eq. (145), and the derivatives of the redshift-dependent parameters given in Eqs. (147)–(150), as well as the analytical derivative with respect to the shot noise, given in Eq. (99).

Pessimistic setting. As a second case we implement the full non-linear model for the observed power spectrum of Eq. (87) and consider scales k <  k_max = 0.25 h Mpc⁻¹. As this is the same range of scales of above linear setting, the comparison of the results obtained in these two cases illustrates the degradation of the constraining power of P_obs(k, μ) due to non-linearities²². This configuration requires P_nw(k;  z), which is also given as a common input to all codes (see Sect. 4.2.1 for details).

In this step we also need to vary the non-linear parameters, σ_p(z_mean) and σ_v(z_mean), where z_mean is the mean redshift of the survey. In this case, we now allow for variations of the growth rate in the FoG term. To this end we rewrite the term containing the growth rate in Eq. (82) as

(152)

Here, we have assumed that σ_p can be written as σ_p(z_mean)D(z)/D(z_mean). Furthermore σ₈(z) can also be written in terms of a σ₈(z_mean), as it was the case for σ_p(z) and σ_v(z), and Eq. (152) then reads:

(153)

where we have now defined a new parameter . Hence the two independent parameters become σ_v(z_mean) and . The FoG term in Eq. (82) is written as

(154)

whereas the damping term remains as

(155)

We calculate the derivatives with respect to the two new non-linear parameters using the same three-point numerical method as those used for the shape parameters (see Eq. (146)). The Fisher matrix now has its largest dimensions: 4 + 2 + 5 × N_z rows and columns.

Intermediate setting. This case corresponds to the pessimistic setting but the maximum wavenumber considered in the analysis is extended up to k_max = 0.30 h Mpc⁻¹. Therefore, the dimension of the Fisher matrix is still 4 + 2 + 5 × N_z.

Optimistic setting. This case corresponds to intermediate setting where the maximum wavenumber considered in the analysis is k_max = 0.30 h Mpc⁻¹; however, we keep σ_p = σ_v fixed to their reference non-zero values, computed using Eq. (81). This means that non-linear corrections are included into our analysis, but we assume that we have perfect knowledge of them. The Fisher matrix in this step has dimensions 4 + 5 × N_z. In practice, the Fisher matrix in this case is the same than in the intermediate setting where we throw the rows and columns corresponding to the parameters σ_p and σ_v.

Projecting to the final parameter space. Once the full Fisher matrices have been computed, they can be projected into the same final cosmological parameter spaces as the other probes described in this paper (see Sect. 3.2.5).

The choice of the parameters in θ_final defines the cosmological models to be considered. We explore extensions of the ΛCDM parameter space by allowing for curvature, dynamical dark energy, and deviations from general relativity.

• PR1

  As a first step, we restrict the analysis to the ΛCDM cosmological model, extended to allow also for non-flat models. In this case, the Jacobian only contains the rows and columns corresponding to the derivatives of the original parameter set with respect to the following :

  (156)

  corresponding to a (4 + 3 × N_z)×6 matrix.

• PR2

  In this step we extend the parameter space considered in PR1 by allowing also for dynamical dark energy models, in which case

  (157)

  Now the Jacobian is given by a (4 + 3 × N_z)×8 matrix. The function σ₈(z) can be expressed as

  (158)

  and consequently fσ₈(z) is given by

  (159)

  where f(z) is the numerical solution of the first-order differential equation in Eq. (26) and the growth factor D(z) is given in terms of f(z) in Eq. (27).

• PR3

  Finally, we consider our most general cosmological models by extending the parameter space of PR2 by including the growth rate exponent parameter, γ:

  (160)

  In practice we substitute the expression of the growth rate given by Eq. (25) into Eq. (159). This case requires one to take into account the derivatives of the observed power spectrum with respect to γ through Eq. (25), which means that the dimensions of the Jacobian matrix are (4 + 3 × N_z)×9.

PRn-flat When we consider flat cosmologies, which correspond to the same parameter spaces as in cases PR1, PR2, and PR3, we eliminate Ω_DE,  0 by assuming into the equations that Ω_DE,  0 = 1 − Ω_m,  0.

As mentioned in Sect. 3.2.5, before projecting the Fisher matrix into a new parameter space, we need to marginalise over some parameters, in particular these are σ_p, σ_v, ln bσ₈(z), and P_shot. Hence, the dependence of the new parameters θ_final enter in the four shape parameters and in the redshift dependent parameters ln D_A(z), ln H(z), and ln fσ₈(z).

In order to project the marginalised Fisher matrix into a new parameter space, we need to create a Jacobian matrix; here we report the most general transformation (which correspond to the ΛCDM non-flat w₀ + w_a + γ cosmology), which is equal for all the four cases listed above²³. The Jacobian is explictly

(161)

The first four rows of the matrix correspond to the derivatives of the shape parameters with respect to the new set of parameters θ_final, whereas all remaining rows correspond to the derivatives or the redshift-dependent parameters with respect to θ_final, which are evaluated at the mean redshift of each bin.

We note that, since background functions do not depend on perturbed quantities, the derivatives of ln D_A and ln H with respect to n_s and γ are zero, as well as those with respect to σ₈. Moreover, since we are using ω_b, 0 and ω_m, 0 as independent parameters, H(z) and D_A(z) only depend on the total matter density, represented by ω_m, 0. Furthermore, the derivatives of the growth rate f(z) with respect to n_s and h are zero.

4.2.3. Results of the code comparison

In this section we compare the output of the six codes capable of estimating cosmological forecasts for GC measurements described in Sect. 4.1. These are BFF, CarFisher, fishMath, FisherMathica, SOAPFish, and TotallySAF. We have compared the results of these codes for each of the cases summarised in Table 7. The main conclusion of this exercise is that all codes are able to provide consistent cosmological forecasts for Euclid.

In Fig. 2, top-left panel, we show the relative difference between the full marginalised errors of the four shape parameters produced by the six codes considered in the current comparison. To produce this plot – and all the similar plots henceforth – we compare, for the constraint on each cosmological parameter, the performance of each code to the median of the constraints obtained by the set of all the available codes. Then, we show the relative difference between a given code and such median in percentage. We note that, by definition of median, in the case an odd number of codes is available, there will always be a code that actually provides the median constraints, meaning its corresponding marker will be placed at 0%. The agreement between the outcomes is better than a few percent for the Intermediate case. In Fig. 2 we also show the relative difference between the full marginalised errors for the redshift dependent parameters [ln D_A(z), ln H(z), ln fσ₈(z)] produced by the six codes considered in the current comparison. Also for these set of parameters the difference among all the outcomes is of the order of 4%.

Finally in Fig. 3 we also show the relative difference of the full marginalised projected parameters among six GC codes, for the non-flat ΛCDM + w₀ + w_a + γ case, using the non-linear recipe up to k_max = 0.30 h Mpc⁻¹, corresponding to the optimistic setting. The agreement between the outcomes are at most of the order of a few percent.

4.3. Weak-lensing code comparison

We now describe the code comparison and validation of the cosmic-shear weak-lensing Fisher-matrix codes available, from the implementation of the recipe to the final results, showing how well the output from the different codes now match. As for the GC case, we began the comparison process from codes whose results differed in both amplitude and orientation of the expected probability contours, with errors as large as 100%, and reduced them to be at least within 10%, as required and discussed in Sect. 3.1.6.

4.3.1. Implementation of the weak-lensing recipe

We implement the Fisher matrix definition in Eq. (127), with the Euclid specifications described in Table 4.

Binning. The sum over ℓ modes in the implementation of Eq. (127) is performed within the range ℓ_min = 10 and ℓ_max = 1500(5000) for the pessimistic (optimistic) case, as discussed in Sect. 3.3.4. We use N_ℓ = 100 logarithmically equi-spaced bins. All the quantities entering Eq. (127) are meant to be evaluated at the linear centre of the bin. We compute the cosmic shear power spectrum and its derivatives in each multipole bin k and define λ = log ℓ, such that

(162)

with λ_k = λ_min + (k − 1)Δλ, and Δλ = (λ_max − λ_min)/N_ℓ. The width of the bin is

(163)

which is therefore different depending on the ℓ value.

Derivatives. A fundamental aspect of each Fisher forecasting code is the computation of the derivatives, which has been discussed and introduced in Eqs. (129) and (130). A typical procedure to compute derivatives suggests that one should first compute the integral of the cosmic shear power spectrum, and then calculate the derivatives by numerical differentiation using, e.g. a finite differencing technique (i.e. a numerical implementation). However, since the integration limits do not depend on the parameters with respect to which one is deriving, the derivative operator can be taken inside the integral, as done in Eq. (129); in this case we can skip numerical differentiation to a large extent, since most of the derivatives can be performed analytically (what we call a “semi-analytical” implementation). Mathematically, the two approaches are equivalent, but in practice this is not always the case because numerical differentiation can introduce a certain degree of noise and some unphysical features. These effects can be reduced to an acceptable level by carefully setting the specifics of any particular implementation. The codes used here adopt different implementation choices, which are then compared.

We detail below the derivatives needed to implement the semianalytical approach. Equation (129) shows that there are two sets of derivatives one must compute to implement the semi-analytical recipe, namely derivatives of the kernel and of the matter and IA power spectra. We address them separately in the following discussion.

Semi-analytical approach: derivatives with respect to background quantities. Derivatives of the kernel and of the matter power spectrum will include derivatives of the Hubble rate and the comoving distance.

With regard to the Hubble rate, for the non-flat case, the latter can be written in terms of E(z) as in Eq. (3), whose expression as a function of the cosmological parameters is given in Eq. (10); we then obtain²⁴

(164)

(165)

(166)

(167)

where ℰ(z) ≡ ln(1 + z) − z/(1 + z). In the flat case we enforce the constraint Ω_DE = 1 − Ω_m, at all redshifts, so that Ω_K = 0. Hence, the (1 + z)² term in Eq. (10) drops, while the derivative with respect to Ω_m, 0 has a further contribution coming from expressing Ω_DE, 0 = 1 − Ω_m, 0. As a consequence, the derivatives ∂ln E(z)/∂Ω_DE, 0 do not appear and the derivative with respect to Ω_m, 0 becomes

(168)

With regard to the comoving distance, we give below the logarithmic derivatives of interest. For the general non-flat case, these are

(169)

(170)

(171)

(172)

(173)

where . Again, when the flatness constraint is enforced, the derivative with respect to Ω_m, 0 includes an extra term and becomes

(174)

while all the others remain the same, provided Ω_DE, 0 = 1 − Ω_m, 0.

Semi-analytical approach: kernel functions. The derivatives of the kernel functions with respect to the parameters p_μ, needed for Eqs. (129) and (130) are conveniently expressed as derivatives of the logarithm of the kernel functions:

(175)

where we recall that the index a ∈ {γγ, Iγ, II}. Let us consider first the lensing kernel (a = γγ). We then obtain

(176)

(177)

(178)

(179)

(180)

with

(181)

All the relevant derivatives of ln E(z) needed for such kernel derivative have been given above.

In order to compute the derivatives of the shear-IA kernel (i.e. a = Iγ), it is convenient to first split it as

(182)

with

(183)

(184)

The non-null derivatives can then be written as

(185)

(186)

(187)

(188)

(189)

Finally, the derivatives of the IA–IA kernel (a = II) are

(190)

where the first row holds for derivatives with respect to parameters Ω_m, 0, Ω_DE, 0, w₀, and w_a, and ∂_μ is a short-hand notation for ∂/∂p_μ.

Semi-analytical approach: matter and IA power spectrum derivatives. As from Eq. (130), the derivatives of the matter power spectrum ∂P_A(k, z)/∂p_μ (with A ∈ {δδ, δI, II}) contain two terms, the latter being present only for background parameters {Ω_m, 0, Ω_DE, 0, w₀, w_a, h}. The second factor in the second term is computed for fixed values of the model parameters, while the first factor in the second term uses the expressions derived in Eqs. (164)–(168) for the derivative of the comoving distance. On the other hand, the first term can only be evaluated numerically, and to do this we rely on a procedure that proves to be more stable than the usual finite differences. Therefore, we adopt the technique proposed by Camera et al. (2017) (see their Appendix A). For each k-z pair, we evaluate P_δδ(k, z) for different values of the parameter of interest p_μ. We choose values that span the ±10% range around the fiducial value of the parameter. If the neighbourhood is small enough, we can approximate the dependence of P_δδ on p_μ as linear (following the definition of derivative). Then we fitted a line through these values and take its slope as our final estimate of the ∂P_δδ(k, z)/∂p_μ derivative.

Once the matter power spectrum derivatives have been computed, the IA related ones can be analytically related to them. Again, it is convenient to move to logarithmic derivatives so that we obtain

(191)

(192)

and, from Eqs. (109) and (110), we find the following expression for the derivative with respect to the cosmological parameters {Ω_m, 0, Ω_DE, 0, Ω_b, 0, w₀, w_a, h, n_s, σ₈}:

(193)

(194)

We note that the second term in Eqs. (193) and (194) is only present for Ω_m, 0, Ω_DE, 0, w₀, w_a, and γ, since the growth factor D(z) depends on these parameters only.

The derivatives with respect to the IA parameters read

(195)

for the δI term, and

(196)

for derivatives of the IA–IA power spectrum.

4.3.2. Settings definition

For the purpose of validating different codes, we have compared their output in many different scenarios. A selection of the cases tested is in Table 8. We consider ΛCDM and the extensions described in Sect. 2 with neutrino mass m_ν = 0.06 eV. We implement the “TakaBird” (Bird et al. 2012; Takahashi et al. 2012) linear to non-linear mapping, as discussed in Sect. 3.3.2. We proceeded by steps of increasing complexity, starting with the simplest assumptions and then adding features after validation. We first compared spatially flat models by setting Ω_DE, 0 = 1 − Ω_m, 0 (hence Ω_K, 0 = 0), neglecting contribution from IA in Eq. (119) and keeping only standard cosmological parameters (plus neutrino mass). Then we tested extensions to non-flat models, as well as including modelling of IA systematics (eNLA model described in Sect. 3.3.1), and implemention of non-standard parameters such as w, w_a, γ discussed in Eq. (32).

The large number of cases offered the possibility to investigate how theoretical uncertainties impact the cosmic shear forecasts and how they degrade as the number of nuisance parameters increases.

4.3.3. Results of the code comparison

The codes used for the code comparison of the cosmic shear weak lensing are described in Sect. 4.1, and include, in particular, seven codes: CCCP, CCCPy, CosmicFish, CosmoSIS, FisherMatica, STAFF, TotallySAF. One of them is publicly available, three of them are planned to be made available, and the features of all codes are summarised in Table 6. Computing the Fisher matrix for the different cases in Table 8 is conceptually simple using all the relevant formulae in Sects. 2 and 3, but there many steps to be implemented, and also many potentials errors that can hamper the implementation. Moreover, the critical role played by the numerical derivatives led us to compare the results from codes implementing different recipes. Among the codes described in Sect. 4.1, for example, it is worth stressing that CCCP, CCCPy, and STAFF adopt the semi-analytical approach, while CosmicFish, CosmoSIS, and TotallySAF rely on fully numerical derivatives; FisherMatica follows a mixed approach (the derivatives of the shear power spectrum are computed by first differentiating and then integrating – as for the semi-analytical case – but the differentiation of the argument of the integrals is performed numerically). A further important difference among the seven codes concerns the way in which derivatives of the matter power spectrum are implemented. While CosmicFish and CosmoSIS compute them internally, the other five codes take them as external input (using the same input). As a consequence, we expect better agreement among these five codes. The external input for CCCP, CCCPY, FisherMatica, STAFF, and TotallySAF has been produced using CAMB with the same settings for CosmicFish and CosmoSIS, in order to be consistent and validated as described in Sect. 4.4.2. Any discrepancy among them should therefore not originate from the inputs, but rather from how the derivatives are computed.

In Fig. 4 we show how the output of the seven codes used for cosmic-shear weak lensing compare among each other, for two representative cases from Table 8. In particular, we plot the percentage difference on the error on each parameter, as found by each code, when marginalising (light grey band) or not marginalising (dark grey contours) on the other parameters. Dots correspond to the different codes used in the comparison. What matters is how close they are to each other, for each parameter, i.e. the hight of the grey bands for each parameter. In all cases (including the ones tested but not shown in the plots for brevity) the percentage difference on errors is within 10%, as required. We note that we plot the relative difference with respect to the median, as explained in Sect. 4.2.3.

4.4. Probe combination code comparison

In this section we compare the codes used to obtain results for the combination of Euclid’s probes, including galaxy clustering coming from the spectroscopic survey (GC_s), weak lensing (WL) and photometric galaxy clustering (GC_ph). For the two latter probes, we also include cross-correlation terms (XC). When cross terms are neglected the Fisher matrices are added. However in the following we discuss the specific implementation of the recipe for cross-correlation terms discussed in Sect. 3.4 and the validation of the four codes that took part in the comparison.

4.4.1. Implementation of the probe combination recipe

The codes used in the comparison, CosmicFish, CosmoSIS, STAFF, and TotallySAF, consider different sets of cosmological primary parameters. CosmoSIS, STAFF, and TotallySAF are able to work directly with the parameter basis of Eq. (53), while CosmicFish uses CAMB’s physical parameters ϑ = {Ω_bh², Ω_ch², h, n_s, A_s, w₀, w_a} as the primary set, and then projects on the chosen set of derived parameters. For CosmicFish, CosmoSIS, and TotallySAF we use the three-point stencil method to compute the numerical derivatives; in order to have the same step size for the different parameters in all codes, CosmicFish uses the steps for the ϑ parameters that are equivalent to the steps used by CosmoSIS and TotallySAF in the set of parameters of Eq. (53). A more detailed explanation of the problems that may arise from the change of the parameter basis is presented in Sect. 4.5. The STAFF code instead uses the semi-analytical approach described in Sect. 4.3, where the numerical part of the derivatives is obtained following the method proposed by Camera et al. (2017) (see their Appendix A).

It is also important to note that the four codes are completely independent of each other and use very different approaches. CosmicFish, STAFF, and TotallySAF implement a second order covariance given in Eq. (137) and build the Fisher matrix as in Eq. (143), while CosmoSIS considers a fourth order covariance given in Eq. (138) and builds the Fisher matrix as in Eq. (142).

Moreover, the codes differ also in how they obtain the cosmological quantities needed to compute the observables used to construct the Fisher matrix. STAFF and TotallySAF obtain these from external input, produced as described in Sect. 4.3. CosmicFish and CosmoSIS are instead directly connected to Boltzmann codes, CAMB_sources for the former and CAMB for the latter; while CosmoSIS fetches from CAMB the cosmological quantities and computes the spectra as described in Sect. 4.5.3, CosmicFish obtains the C(ℓ) directly from the internal computation of CAMB_sources and can therefore in principle avoid using the Limber approximation, which is instead implicit in the recipe of Sect. 4.5.3. The fiducial C(ℓ) obtained with these methods have been compared, finding an agreement to within approximately 1%.

4.4.2. Settings definition

For the purpose of the code validation, as for GC and WL single probes, we have investigated several different choices of parameters and settings. We discuss here some of them, as presented below. We use for the non-linear power spectrum the “TakaBird” modified HaloFit (HF) prescription of Takahashi et al. (2012), with additional corrections for massive-neutrinos following Bird et al. (2012), as discussed in Sects. 3.4.2 and 3.3.2. We use the fiducial discussed in Sect. 3.1.5 and, if not otherwise specified, a baseline neutrino mass ∑m_ν = 0.06  eV (with one massive-neutrino and two massless neutrinos) and we include intrinsic alignment in the weak-lensing calculations.

As discussed in Sects. 3.3.4 and 3.4.4, and as will be further clarified in Sect. 5, in the following we identify two configurations, corresponding to a stronger or weaker cut of the maximum 2D Fourier mode, ℓ_max, for the C_ℓ of weak lensing, photometric galaxy clustering, and the cross-correlation of the two probes; we further consider optimistic and pessimistic cuts when adding spectroscopic galaxy clustering as follows.

(197)

(198)

We note that when combining with the GC_s in the pessimistic case we cut GC_ph at redshift z <  0.9 in order to avoid overlap of the two surveys. This is done under the conservative assumption that we are not able to model the correlation of the two. In the optimistic case instead, we perform no cut in redshift and use the full GC_ph survey, assuming that the correlation of the two clustering surveys can be neglected given the relatively small common redshift range. Finally, in the optimistic case the GC_s Fisher matrix is obtained fixing the non-linear parameters σ_p(z_mean) and σ_v(z_mean) (see Sect. 4.2 for details), rather than marginalising over them as is done in the pessimistic case.

4.4.3. Results of the code comparison

In Fig. 5 we show the comparison of the four codes for the pessimistic and optimistic configurations (left and right panels respectively) for the combination of observables GC_ph+WL+XC and the full combination GC_s+GC_ph+WL+XC described in Sect. 4.4.2. In all cases investigated here we use the non-linear corrections given by the “TakaBird” formalism and marginalise over intrinsic alignment and galaxy bias nuisance parameters, included in the analysis as described in Sect. 3.4. We find a good agreement between the four codes, with the relative difference on the errors on cosmological parameters on each pair of codes always within 10%.

In addition to the results shown here, the comparison was performed using also a different neutrino mass (∑m_ν = 0.15  eV), in a ΛCDM cosmology, combining subgroups of observables, for nonflat models, and without intrinsic alignment nuisance parameters, obtaining similar results in all cases. We also compared the combination of GC_ph and WL without the inclusion of the cross correlation, again finding compatible results between the four codes.

4.5. Lessons learned

The process of comparison and validation among different codes involved a series of tests, choices, and lessons learned through the process. Here we describe some of these in detail, as a list of suggestions meant to facilitate comparison with any other Fisher matrix code. As usual, we highlight the main points, separating them among different probes.

4.5.1. Galaxy clustering

We list here some critical points that can be responsible for the different results one can find when comparing different galaxy clustering Fisher matrix codes, together with general tips to keep in mind when writing a new forecasting code.

Input. Input quantities are needed in order to evaluate the observed galaxy power spectrum and these are the matter power spectra produced with CAMB or CLASS. If input files are produced independently we suggest to switch off the default spline in CAMB. This is to avoid, where possible, unnecessary interpolations that may introduce instabilities. In the case of CLASS, one would have to modify the source code in order to do the same.

Interpolations. In order to compute reliable derivatives of P_obs(k, μ; z), needed in Eq. (87), an interpolation scheme is necessary: we can choose to interpolate in the normal basis {k,P_δδ(k,z)}, or in the logarithmic basis {log₁₀k,log₁₀P_δδ(k,z)}. We found that the latter method produces more reliable results, as it is the best way to follow quantities that vary over orders of magnitude. Moreover, we found that it is best to interpolate such logarithmic tables with cubic splines, and to avoid any other unnecessary interpolation.

Derivatives. The observed galaxy power spectrum in Eq. (87) is governed by different parameters with respect to which we need to take the derivatives. For some of the parameters an analytical derivative can be constructed. However, we suggest considering fully numerical derivatives for all the parameters. To explain the reason for this let us take the particular example of the derivatives with respect to the Hubble parameter and the angular diameter distance. In the literature it is found that it is possible to take the derivative of the matter power spectrum with respect to ln H(z) and ln D_A(z) by using the chain rule, i.e.

(199)

we note that k depends on H(z) via Eq. (78). This is a semi-analytical derivative because the derivative of P_δδ(k, z) with respect to k has to be computed numerically, whereas the dk/dln H(z) can be calculated analytically. During the code comparison we found that the semi-analytical derivatives of P_obs(k, μ; z) with respect to ln H(z) and ln D_A(z) do not necessarily agree with the numerical derivatives for k >  0.1 h Mpc⁻¹. Moreover, we found that the semi-analytical derivative is not stable when varying the k-binning, while the numerical derivative is quite stable. For the latter reason it is safer to produce numerical derivatives with respect to ln H(z) and ln D_A(z). This also has the advantage of facilitating the procedure of passing from linear modelling to non-linear corrections to P_obs(k, μ; z) in Eq. (87).

For a generic function f(x)²⁵, the three-point stencil derivative of f with respect to x at x₀ is given by

(200)

with ε ≪ 1 an increment, usually an arbitrary choice. The three-point stencil approximates the derivative up to order ε² errors. To be more precise, another choice is represented by the five-point stencil

(201)

When dealing with numerical derivatives, the choice of the step ε is crucial. On the one hand, if it is chosen too small, the subtraction yields a large rounding error. Actually, all the finite difference formulae are ill-conditioned, and due to cancellations they produce a value of zero if ε is small enough. On the other hand, if ε is too large, the calculation of the slope of the secant line is more accurately calculated, but the estimate of the slope of the tangent by using the secant could be worse. Hence, we recommend checking the convergence of the derivatives when choosing the steps, in order to obtain realistic errors and FoM values, and reach agreement among different implementations of the same fisher recipe. Moreover, in this respect, we have to warn the reader that single-probe results can be step-dependent, but they become step-independent when Euclid probes are combined together (see Sect. 4.5.3). This is physically due to the difficulty for the single probe to simultaneously constrain many parameters, for the GC case, in particular, a total of 24 different parameters, i.e. 4 shape parameters and 20 redshift-dependent ones.

Integrals. During the code comparison we did not find any issue concerning the integration method. However, in order to avoid unnecessary interpolations, we suggest to consider the same k binning as in the input files.

k-binning in input files. In order to validate the code we suggest to check stability of the code when using input files having different step sizes in log₁₀k space. We have tested the stability against several possible step sizes, namely 20 steps per decade in k, 50 steps, and 200 steps; as our default we have chosen the case of 50 steps per decade in k.

Stability of the results. As a final test, we performed several stability checks in order to ensure that the codes used in the code comparison give reliable results. One of the main challenges in the Fisher matrix approach is to consider an appropriate step for the derivatives of the observed matter power spectrum with respect to the parameters. In our analysis, we performed a systematic analysis considering, firstly the step in the derivatives with respect to the shape parameters ω_m, 0,  h, ω_b, 0, n_s, and, separately, we let the step for the derivatives vary with respect to ln D_A(z) and ln H(z) only²⁶.

In Fig. 6 left panel, we show the FoMs when varying the step for the derivatives of the observed galaxy power spectrum with respect to the shape parameters for the w₀w_aflat projection case; we neglected to show all the other projection cases as they manifest a similar behaviour. The analysis has been performed by varying one parameter at a time. In Fig. 6 (right panel), we show the FoMs when varying the step for the derivatives of the observed galaxy power spectrum with respect to the ln D_A for all the projection cases.

4.5.2. Weak lensing

Computing the weak lensing Fisher matrix involves several different steps, which must be dealt with using great care to avoid biasing the final result. The code comparison challenge has allowed us to identify several critical issues, and we briefly comment upon them in the following.

Fiducial shear power spectra. The tomographic angular power spectrum of cosmic shear, , is given in Eq. (118). Neither of the two integrands, the various kernels or the power spectrum, are analytic functions. Hence, each code in the code comparison has to compute them numerically, and then interpolate using an intermediate look-up table. We found that the lensing kernel of Eq. (120), , is a well-behaved function of redshift, and the results (in terms of spectra and Fisher matrices) are stable provided that the number of points in the table N_z ≥ 300.

The matter power spectrum, on the contrary, is a much more involved function to compute. This is because it depends on two variables, namely the scale k and the redshift z, and it presents an oscillating behaviour because of the baryon acoustic oscillation wiggles. Moreover, we evaluate the various P_A(k, z) (A = {δδ, δI, II}) at k_ℓ = (ℓ+1/2)/r(z), which means spanning a large range in k. On the other hand, for the various a = {γγ, Iγ, II}, the product is suppressed at both small and high z because the lensing kernel is only significantly non-zero over a limited redshift range. To equally well sample the support of the integrand, and to model accurately the baryon acoustic oscillations, we find that a good strategy is to create a table with {log₁₀k, z, log₁₀P_A(k, z)} over the range 0.0005 ≤ k/Mpc⁻¹ ≤ 35 and 0 ≤ z ≤ 2.5 with N_k × N_z entries. The results are stable with the same N_z as before and N_k = 800. We stress that we do not use h units, i.e. we take k in Mpc⁻¹ and P_A(k, z) in Mpc³ in order to avoid unnecessary complications when taking derivatives with respect to the dimensionless Hubble constant h.

Numerical versus semi-analytic derivatives. The derivatives of the shear power spectra are a critical ingredient of the WL Fisher matrix. Two different approaches are possible, as shown by the identity

(202)

(203)

According to this relation, one can either: (i) first compute the integral and then the derivatives, for instance through finite-difference methods; or (ii) first compute the derivatives of the integrand and then integrate the result. The former choice fully relies on numerical techniques, while the latter is partly analytic, since the derivatives of the kernel may be transformed into integrals of analytic quantities (at least for the models we have considered in this paper). Both approaches have their pros and cons. The fully numerical approach is quite generic and does not require any intervention from the user if a different cosmological model is adopted (apart from changes in the Hubble rate expression), but it is prone to the instability of numerical derivatives. On the contrary, the “semi-analytic” approach transforms almost all the derivatives into integrals (the only exception being those of the matter power spectrum). However, the user has to adjust the code to update the expressions for the derivatives of the Hubble parameter and the radial comoving distance each time a new cosmological model is considered. Of the seven codes participating to the challenge, three (CCCP, CCCPy, and STAFF) adopt the semi-analytic approach, three (CosmicFish, CosmoSIS, and TotallySAF) employ a fully numerical method, and a single code (FisherMathica) uses a mixed technique. Provided care is taken when computing the numerical derivatives, we find all the approaches provide consistent results. Choosing among the two schemes is therefore the users choice and does not affect results.

Numerical differentiation methods. Regardless of the choice one decides to adopt for derivative computation, the method implemented to compute the numerical derivatives of the quantities of interest plays a critical role in the estimate of the Fisher matrix. Finite differences methods are the typical choice.

As already mentioned in point 3. of Sect. 4.5.1, particular care should be taken in the choice of the step ε of the derivatives. In fact, when dealing for example with derivatives of the tomographic cosmic-shear angular power spectrum , the situation is much more complicated. This is due to the fact that even the most careful choice of ε is not, in general, hold true over the whole range of ℓ modes, and for any ij bin pair combination. As an example, consider the baryon fraction Ω_b, 0: its main impact on the matter power spectrum, and, consequently, on the shear signal, is at BAO scales, corresponding today to k_BAO, 0 ≃ 2π/150 Mpc⁻¹. However, this physical scale shows up at different angular multipoles, depending on the redshift bin considered. Therefore, a step ε optimised to capture at best the response of to Ω_b, 0 at i = j = 1 and ℓ = k_BAO, 0r(z = 0), may itself (and indeed typically will) lead to potentially significant numerical errors for a different combination of ij and ℓ mode.

As a possible solution to this issue, Camera et al. (2017) proposed an alternative method, first implemented in the CCCP code (see Sect. 4.1). Let us consider a parameter x: for each combination of ij and ℓ – in our example, the generic function f(x) – we sample the x-line in parameter space in, say, 15 x_i points around the fiducial value x₀ (including x₀). In the hypothesis that the neighbourhood is small enough, we interpolate the x_i values with a straight line. Then, by the definition of a derivative, the slope of the linear interpolation is the derivative of the function at the fiducial value, i.e. df/dx|_x0. We can test our ansatz by checking whether the spread between the linearly fitted f_fit(x_i) and the true values f(x_i) is less than some required accuracy, e.g. 1%. If this requirement is not met, we zoom in on the sampled x-range by cutting out a few values on the edges, until we reach the requested accuracy.

Of the seven codes participating in the comparison challenge, four of them implemented the aforementioned technique to compute derivatives at different levels. In particular, CCCP, CCCPy, and STAFF use the technique to compute the derivatives of the matter power spectrum, while the rest of derivatives are computed semi-analytically. FisherMathica uses this method for the matter power spectrum derivative, while the derivatives are computed with standard three- and five-point stencil methods. On the other hand, CosmicFish, CosmoSIS, and TotallySAF rely on three-point stencil.

The choice of derivative method impacts the ability to make statements on the stability of the results. For example, CCCP, CCCPy, and STAFF do not rely on the classical finite difference method, so there is no need to check their stability against the step ε used in the numerical differentiation; while some concern might be related to the range used in the fitting procedure, we have tested this and found to have no impact on the results. However, ε plays a key role in CosmicFish, CosmoSIS, and TotallySAF so that a detailed analysis of how the results change with ε has been carried out. We find that for 4%≤ε ≤ 10% both the marginalised errors on the cosmological parameters and the FoM are stable (i.e. unchanged over this range). A similar test has also been done for FisherMathica where stability is achieved for ε ≃ 3% no matter whether three- or five-point stencil derivatives are used; we note that the weaker requirement on ε is also related to the fact that CosmicFish, CosmoSIS, and TotallySAF are also used for the probe combination, which is more demanding in terms of stability than cosmic shear.

We find that the FoM and marginalised errors on {Ω_m, 0, w₀, w_a, σ₈} are stable with respect to the methodology and the step size, whereas this is not the case for {Ω_b, 0, h, n_s}. For the unstable parameters this is not surprising, since weak lensing is less sensitive to these parameters. However, what is really difficult to control is the relative orientation of the contours in the 2D parameter projections. Indeed, we find that these orientations – governed by the off-diagonal elements of the Fisher matrix – turn out to be quite sensitive to the choice of ε and to whether a three- or five-point stencil is used. However, provided the above quoted requirements on the step are met, stability of the results is achieved.

Number of multipoles and spacing. The WL Fisher matrix is obtained by summing up contributions from different multipoles ℓ spanning the wide range 10 ≤ ℓ ≤ ℓ_max. In determining the number and values of the ℓ modes to use in the calculation two options are possible: one can use N_ℓ linearly equispaced bins in ℓ; or one can use N_λ logarithmically equispaced bins in λ = log ℓ. We have explored both possibilities when evaluating the quantities of interest (i.e. and its derivatives)²⁷. We found that N_ℓ = N_λ = 100 gives consistent results among the two approaches for both ℓ_max = 1500 and ℓ_max = 5000. Furthermore, we found that reducing the number of logarithmic bins up to N_λ = 60 changes the FoM and the marginalised constraints by less than 1%, so that speeding up the computation is possible with a small degradation in accuracy. We nevertheless recommend using N_λ = 100 logarithmic bins, since they more finely sample the shear power spectra.

4.5.3. Probe combination

In this section we list some of the critical points that can lead to different results when comparing outputs from the probe combination codes. These involve both the construction of the Fisher matrix and the operations needed to obtain cosmological constraints from it, as well as possible issues arising when obtaining the relevant quantities from the Boltzmann solvers.

Matrix inversion. The final sensitivity expressed through the correlation matrix obtained by inverting the Fisher matrix is at the heart of the methodology. This has a drawback: matrix inversion can be tricky when the condition number is large i.e. the matrix is then ill-conditionned²⁸. In the case of weak lensing and galaxy clustering the condition number is of the order of 2 × 10⁶. In such a case, in principle, a precision better than 0.5 × 10⁻⁷ is needed for a secure inversion (i.e. no numerical instability) of the matrix with an accuracy better than 10%. However, in practice the required accuracy depends on the matrix. This means that each Fisher matrix needs a specific treatment to quantify the accuracy needed to avoid numerical instability. The ability to achieve this accuracy was tested by perturbing the Fisher matrix with a Gaussian-distributed random percentage of all the elements of the WL Fisher matrix (using the TakaBird non-linear prescription) according to the approach proposed by Yahia-Cherif et al. (2020). The 10% agreement level on the Fisher matrix is reached at > 99% degree of confidence when using a (relative) perturbation of 4 × 10⁻⁵. The 65% confidence level is reached at 10⁻⁴. More investigations were done in the spectroscopic GC linear case: all elements of the Fisher were perturbed one by one. Some of the elements can be perturbed by 100% without changing the FoM by more than 10%, whereas the most sensitive elements to the perturbations can give differences in the FoM higher than 10% for a perturbation of the order of 5 × 10⁻⁴.

The above results have implications on the accuracy required for the derivatives entering the Fisher matrix. Standard methods of derivation are the N-point stencils with a fixed step. A single step size for all parameters may not be appropriate. Two sources of errors exist: the truncation errors that increase with the step size; and the numerical noise that increases when the step size diminishes. The outputs of the Boltzmann code may be subject to a significant numerical noise for small steps: in the calculation of the power spectrum P_δδ(k) this problem was noted at a step-size below 10⁻⁵ using CAMB and 10⁻⁴ with CLASS. The fact that the outputs have a different sensitivity depending on the parameter does not ease the determination of an optimal step. An adaptive scheme can then be used but at the price of a significant increase of the computational time.

This issue is extremely relevant for parameters that are not well constrained by the probe in use; this is because poorly constrained parameters lead to elements of the Fisher matrix that have small numerical values, and therefore lead to complications when inverting the matrix. We therefore find that this problem is most significant when using the observational probes separately, while the combination of these probes has a weaker dependence on the choice of the step size.

Scale sampling. The computation of C(ℓ) for the considered observables is performed by obtaining the matter power spectrum from the Boltzmann solver, sampled in the redshift and scale ranges considered. The comparison performed highlighted how results can depend on this sampling. In particular, we noticed that non-linear corrections are particularly sensitive to the number of scales at which the matter power spectrum is computed; this is due to the fact that we compute the power spectrum up to k_max = 70 Mpc⁻¹ and the highest scales are sampled more sparsely due to the logarithmic step in k used in the Boltzmann solver. Since the codes interpolate between these k values when performing the integral needed to compute the C(ℓ), a step that is too large in the sampling leads to artificial oscillations in the matter power spectrum that can significantly impact the results.

We recommend therefore to check carefully that the chosen step is small enough so that no artificial oscillations are introduced.

Change of parameter basis. The codes investigated for probe combination (CosmicFish, CosmoSIS, STAFF, and TotallySAF) deal with the choice of the primary parameters in different ways; while the first uses CAMB parameters ϑ = {ω_b, 0,  ω_c, 0,  h,  n_s,  A_s} as the primary set and can then project on the chosen set of derived parameters, CosmoSIS is able to work both with this set and the target set of this paper, i.e. θ_final = {Ω_m, 0,  Ω_b, 0,  h,  n_s,  σ₈}. During our comparison effort we found at first that while the agreement of the two codes was optimal when both codes used the first basis, the same agreement was not reached when using the latter for CosmoSIS and using the projection of CosmicFish. This comes from the fact that varying one parameter in the θ_final basis, while keeping the others fixed, implies a variation of more than one parameter in the ϑ set; as an example, a variation of Ω_m, 0 would obviously imply a variation of ω_c, 0 in the ϑ set, with ω_b, 0 fixed in order not to vary Ω_b, 0, but it would also require a variation of A_s, needed to keep σ₈ at its fiducial value. In order to fix this issue, CosmicFish was modified in order to be able to effectively work in the θ_final basis, i.e. using as input the full set of ϑ parameters corresponding to each variation of the θ_final set. We therefore recommend to pay attention to this issue when the parameters varied in the Fisher code do not coincide with those varied in the Boltzmann code that is used.

5. Results

In this section we show the results of the validated Fisher-matrix codes for the cosmological parameters of interest (baseline and extensions), whose fiducial parameters have been illustrated in Table 1:

(204)

We do this for the following scenarios:

• ΛCDM within a spatially flat (and non-flat) geometry;

• w₀, w_a within a spatially flat (and non-flat) geometry;

• w₀, w_a, γ within a spatially flat (and non-flat) geometry.

In the ΛCDM case Ω_DE,  0 ≡ Ω_Λ, the equation-of-state parameter and its derivative are fixed to w₀ = −1 and w_a = 0, and the modified gravity parameter γ is fixed to γ = 6/11²⁹. For the flat case Ω_DE,  0 ≡ 1 − Ω_m, 0. All results refer to a fixed value of the neutrino mass of ∑m_ν = 0.06 eV and to the fixed value of τ reported below Table 1. The latter does not have any influence on the cosmological observables considered, and is simply mentioned here so that anyone interested in validating their own code can perform a consistent comparison.

As illustrated in Sect. 3.4, we compute cross-correlation terms only between the photometric galaxy clustering and the weak-lensing survey. We recall here for convenience that in Sect. 4.4.2 we have identified two configurations, a pessimistic (optimistic) setting, corresponding to: (a) a stronger (weaker) cut of the maximum angular mode ℓ_max, for the C(ℓ) of WL, GC_ph, and the cross-correlation of the two probes; and (b) a stronger (weaker) cut when adding GC_s. For the pessimistic case, the WL non-linear regime is cut at ℓ_max = 1500. We further limit the maximum multipole used for the analysis of GC_ph to ℓ_max = 750, with a further cut in redshift of z <  0.9 when this probe is combined with GC_s to safely neglect cross terms between photometric and spectroscopic data. Consequently, in the pessimistic case the cross-correlation analysis is limited to the lowest of these multipoles, with . For GC_s we reach a value of k_max = 0.25 h Mpc⁻¹. The optimistic case extends the regime of GC_s to k_max = 0.30 h Mpc⁻¹, fixing the parameters σ_p and σ_v; the range used for GC_ph is also extended up to with no cut in redshift. As discussed in Sect. 3.4.4, this choice of GC_ph is quite ideal, since it entirely neglects non-Gaussian terms, which are nevertheless important to consider for this probe and become important earlier than for WL; it does, however, show the potential of such a probe, if one were able to also include such high multipoles. The optimistic WL setting extends the non-linear regime up to ℓ_max = 5000. We can summarise this as follows :

(205)

(206)

For GC_s, we have further tested two more settings, as they can provide a valuable reference for comparison:

(207)

(208)

We have also verified that when considering a higher value in neutrino mass of ∑m_ν = 0.15 eV, the results stay substantially the same, with variations of less than 1% on the relative error for each parameter. We do not consider here a linear or intermediate setting for weak lensing, whose power is very much related to the inclusion of non-linear scales.

5.1. ΛCDM

We start by presenting results for a cosmological constant ΛCDM scenario. The marginalised 1σ errors on cosmological parameters, relative to the corresponding fiducial values, are shown in Fig. 7 for the pessimistic (left panels) and optimistic (right panels) settings. Top (bottom) panels refer to a spatially flat (non-flat) cosmology. Relative marginalised errors are also detailed in Table 9.

The combination of GC_s and WL is very powerful in breaking degeneracies among cosmological parameters and reduces the error on all parameters down to at least 0.3% in the flat case and at least the percent level in the non-flat case, both in the pessimistic and in the optimistic settings. When the information of GC_ph is added (green bars in Fig. 7) the errors are further reduced for all cases and parameters; on the other hand, the gain in adding cross-correlation terms (XC) among GC_ph and WL (yellow bars) is very small in the ΛCDM scenario, and could even degrade the results relative to an optimistic case when added in a pessimistic setting. However where results get worse when adding more information, we attribute this to numerical instability of the Fisher matrix formalism itself (even despite the high accuracy we achieve in this paper).

This situation changes for more complex models, where we see that XC has a much larger (positive) impact. When comparing the pessimistic and optimistic errors (Table 9), GC_s is the probe that gains most from the inclusion of higher modes. Finally, relative errors for a non-flat ΛCDM cosmology, shown in Fig. 7 (lower panels) and in Table 10 are, as expected, larger, since we are adding an extra parameter.

5.2. Changing the background: w₀, w_a

We now move to the case in which the equation of state of a dark-energy component is not w = −1 at all times (as in a cosmological constant case), but can vary in time. We consider the case in which such variation in time is parameterised by w₀, w_a, as defined in Eq. (8). Results on the marginalised errors forecasted for the Euclid survey are shown in Fig. 8. As before, left (right) panels in Fig. 8 refer to the pessimistic (optimistic) settings, while top (bottom) panels refer to a flat (non-flat) cosmology. Fisher forecasts assume a ΛCDM fiducial model, around which contours are drawn. Marginalised errors for different combinations of the observables and different settings are shown in more detail in Tables 11 and 12. The uncertainty on all cosmological parameters is reduced when GC_s and WL are combined. Furthermore, unlike what happens in the ΛCDM case, cross-correlations have a significant impact on the estimation of the errors, and help to further tighten uncertainties, in particular on w₀.

For this set of cosmologies, we can also estimate the FoM for w₀, w_a, as defined in Eq. (49), which is an estimate of the performance of the experiment in constraining a specific set of parameters. Results for the FoM for w₀, w_a are shown in Table 13 for both the single probes (GC_s, GC_ph, WL separately) and their different combinations, within the pessimistic and optimistic settings described in Sect. 5.

The combination of all three probes (without cross-correlations) reaches a FoM of 122 (376) in the pessimistic (optimistic) settings, for a spatially flat cosmology. Including cross-correlations has a substantial (positive) impact on the FoM, which reaches 441 in the pessimistic setting (i.e. enhancing it by a factor of 3.5) and improves up to 1245 in an optimistic setting (i.e. a factor of 3.3 larger than in absence of cross-correlations). This demonstrates the importance, for a w₀, w_a scenario, of including cross-correlations in order to fully exploit future Euclid data. This is confirmed when looking at the correlation matrix, plotted for this model (in the optimistic setting) in Fig. 9 (this plot is obtained following Fig. 4 in Casas et al. 2017). The left panel shows the photometric survey combination, without cross-correlation. The matrix is not diagonal, indicating the presence of correlations among cosmological parameters. When XC terms are taken into account (right panel), the correlation matrix non-diagonal terms are reduced, i.e. correlations among the corresponding parameters are removed; this in turn allows for better constraints on those parameters.

The impact of GC_ph is substantially more significant in the optimistic setting, i.e. ideally, there is much gain in being able to include as many multipoles as possible in the analysis, in order to retain information from the GC_ph. It is important to recall, however, as discussed in Sect. 3.4.4, that the optimistic setting for this probe entirely neglects non-Gaussian terms in the analysis and still includes all multipoles up to ℓ_max = 3000 for this specific probe; this might be regarded as too optimistic, as non-Gaussian terms are estimated to become important earlier than that, and earlier than for cosmic shear. We still include results for this ideal setting here to show the potential power that GC_ph could have, provided we could extend the analysis to such high multipoles.

A final remark on WL: not including IA would give FoM of (27, 6) (flat and non flat), slightly higher than the one shown in the table, in the pessimistic setting for the flat and non-flat cosmologies, respectively, and values (54, 15) in the optimistic setting. However the absence of IA makes the result more unreliable and unrealistic and therefore we do not include these values in the table, but only report them here in the text for comparison with other analyses.

For the w₀, w_a scenario we also show the contour plots in Fig. 10 (and an expanded in Fig. 11), for the optimistic setting and flat cosmology. As expected, parameters characterising dark energy, such as w₀, w_a, are degenerate with the spectral index n_s, since they both tilt the matter power spectrum at a certain scale. The combination of different probes, as in GC_s+WL, helps to break this degeneracy and to tighten constraints on both parameters.

Furthermore, in Fig. 12 we show the photometric survey alone (blue) including GC_ph and WL, and the impact of their cross-correlations (red), in the absence of any spectroscopic information; the ellipses overlap with the full (yellow) combination, highlighting the impact of adding GC_s to the photometric survey.

A comparison of yellow and red contours shows that the inclusion of cross-correlations in the photometric survey causes the impact of GC_s to remain quite significant on Ω_b, 0, h, n_s; however GC_s has less impact on w₀, w_a, Ω_m, 0, and σ₈, where the yellow contours almost overlap with the red ones. As discussed above, this is an ideal setting, but again shows the potential of the photometric probe if one could extend the analysis to high multipoles.

5.3. Changing the growth: w₀, w_a, γ

We now consider the case in which both the background and the growth of structure are changing, as described in Sect. 2.4. In the w₀, w_a, γ scenario the relative marginal errors for GC_s, WL, and their combination, are shown in Fig. 13 for flat and non-flat spatial cosmologies, for pessimistic (left panels) and optimistic (right panels) settings. The corresponding errors are also shown in Tables 14 and 15. Including the extra parameter γ into the analysis degrades the errors, as expected, since we are adding extra degeneracies. For a pessimistic setting, GC_s and WL are comparable, with the exception of Ω_m, 0 (better determined by WL) and h (better determined by GC_s). The combination of the two probes is able to push the accuracy down to the few percent level even in the pessimistic setting, allowing us to also determine w_a, which is nearly unconstrained by the single probes. For the optimistic case we have an improvement of approximately a factor of 3 for GC_s and of approximately 2 for WL, with respect to the pessimistic case.

5.4. Overview results for galaxy clustering and weak lensing

We have so far shown results for different probe combinations, for a given model. We finally collect here results for a given probe and different models, to have at hand how they perform when comparing different models and settings. In Table 16 we show the four settings tested for GC_s only, for all models and parameters. In Table 17 we show the overview for WL only, for which we recall that we only tested the non-linear settings, since this is essential to grasp any information from this probe.

Finally, the combination of all probes, including cross-correlations, is summarised in Table 18 for convenience, for the available models. If we increase the maximum wavenumber k, i.e. comparing the pessimistic with the intermediate settings, we clearly gain more information because more scales are taken into account. We observe a decrease on the errors of about 10% for the set of parameters in both ΛCDM flat and non-flat cases. For the other set of parameters, i.e. w₀w_a and w₀w_aγ, for both flat and non-flat cases, the errors decrease by about 25−30%, with the exception of γ for which we gain 15% information only.

Finally, comparing the intermediate with the optimistic setting, i.e. keeping k_max = 0.3 h Mpc⁻¹ but assuming perfect knowledge on the non-linear parameters σ_p and σ_v, we gain information over all the cosmological parameters. For example, considering the projection case w₀, w_a flat, the errors on the dark energy parameters decrease by about 100% and 150%, for w₀ and w_a, respectively. If curvature is also considered in the analysis, the errors on w₀ and w_a decrease by 28% only.

6. Conclusions

In this paper we have presented the validation of Fisher-matrix forecasts for the Euclid satellite, as a result of activity carried out by the inter-science task force for forecasting, joining expertise from different science working groups (SWGs) of the Euclid Collaboration. In particular, we have included in this analysis consideration of Euclid primary probes, spectroscopic galaxy clustering and weak lensing cosmic shear, as well as the cross-correlation between WL and photometric galaxy clustering. Our target was to provide reliable numerical codes for each probe, which can be then used for a variety of forecasts and specifications – including applications beyond the scope of this paper. We also provide full documentation and tools that will allow any new code owner to perform their own validation

At the start of this work, some numerical codes were already available, although they were providing different results, both in amplitude and in orientation. A few of them were producing similar results due to a similar validation carried out, separately, within the theory working group (TWG). However the recipe used for the different probes had not yet been developed or validated. Our task was therefore to join expertise from different SWGs and validate a common recipe for each probe, including GC_s, WL, GC_ph and the cross correlation across each redshift bin for the two photometric probes. The first result of this analysis is therefore the documentation of each of these recipes in Sect. 3, which are now defined for the Euclid Collaboration and for anyone else who would like to implement and/or validate their own numerical code. The development of the recipe allowed an improvement in all pre-existing codes, making sure that they are all aiming at describing the same features. In Sect. 3.1.6 we have also elaborated on the accuracy that is required in Fisher matrix forecasts, i.e. on the parameter uncertainties, on the accuracy of the Fisher matrix approximation itself, and on how it is computed. We have concluded that to reach ≤10% precision on parameter errors and FoM and ≤2% precision on the orientation of the contours, the corresponding required fractional Fisher matrix precision (both in calculation of elements, and its inversion) is approximately 10⁻⁴, which is a reasonable aim for numerical implementations.

As a second step, we validated the way in which the implementation was performed, for each recipe and numerical code. In Sect. 4 we provided a detailed documentation of each step, describing what works best, and how different codes compare to each other. If readers wish to compare and validate their own Fisher matrix code(s) for Euclid, they should follow these instructions to do so. As part of the code-comparison process, all codes were adapted to use common inputs. We make these inputs, alongside our outputs, available with this paper, describing them and providing instructions for their use in Appendix A. One main product of this work is then the validation of several codes, some of which will be made publicly available so that they can be reliably used, in their present version, for forecasts. A description of all validated codes is provided in Sect. 4 where we highlight their differences (in terms of which probe they consider, which programming language they use, which features they can deal with – or example external input), as summarised in Table 6. All validated codes match with an error on the parameters uncertainty of less than 10% which, as initially stated, was our requirements. We have further developed an appropriate visualisation that allows us to check how close the outputs of the codes are with respect to each other, for each probe: these are shown in Figs. 2, 3 for GC_s, Fig. 4 for WL, and Fig. 5 for probe combination. In Appendix A we also provide the plotting scripts, to allow readers to compare their own results or in order to overlap Euclid Fisher matrices with those of other experiments.

Many tests have been carried out along the course of this work for all probes, involving different numerical strategies (e.g. how to calculate the derivatives), as well as different recipes (e.g. how to model the non-linear contribution, and how to take different cuts depending on the emerging presence of non-Gaussian contributions or baryons). A summary of the lessons we learned while performing such validation has been presented in Sect. 4.5 for each probe. Further investigation will certainly be necessary to optimise the modelling at non-linear scales and will be an object of study for a new taskforce. In this paper we have identified pessimistic and optimistic settings, summarised in Sect. 5 and Eq. (206), for which we have shown final results.

We have focused on the cosmological context described in Sect. 2 which includes: ΛCDM; a (w₀, w_a) scenario in which the equation of state of dark energy is changing with time; the case (w₀, w_a, γ) in which the growth is parameterised by γ and allowed to change compare to its fiducial value in general relativity. All cosmologies were tested within a spatially flat and non-flat geometry. For all models we present a series of histogram bar plots (Fig. 7 for ΛCDM and Fig. 8 for w₀, w_a; Fig. 13 for w₀, w_a, γ) that allow marginalised 1σ errors on cosmological parameters to be visualised for different probe combinations.

In the (w₀, w_a) case we can define a FoM Eq. (49) for which we show our final results in Table 13 for different probe combinations. This table represents one of the main results of this paper. The FoM for the full combination of probes reaches 377 (1257) for a flat cosmology within a pessimistic (optimistic) setting. Values decrease as expected if a non-flat cosmology is allowed.

We have demonstrated that the impact of cross-correlations in redshift among GC_ph and WL is particularly relevant in models beyond ΛCDM: XC improves the FoM for a (w₀, w_a) scenario by a factor of ∼3.6 (∼3.3) in a pessimistic (optimistic) setting with respect to the case in which primary probes and GC_ph are combined as if they were independent probes. We have visually shown how including XC removes degeneracies among parameters (Fig. 9). Within the ΛCDM scenario, instead, XC can have a negligible impact.

All Euclid specifications used in this paper are detailed in Table 2 for the spectroscopic survey and in Table 4 for the photometric survey. In particular, the expected number density of observed Hα emitters for the spectroscopic survey shown in Table 3 has been updated since the Red Book (Laureijs et al. 2011) to match new observations of number densities and new instrument and survey specifications. Further modifications, such as the different areas and depth specifications motivated from the Euclid flagship simulations used within the Science Performance Verification (SPV), may still be possible in the future and will now be straightforward to implement, given a validated forecast pipeline. Since these simulation-derived quantities will change before publication of any SPV results, we chose to remain as close as possible to the Red Book specifications with regard to the photometric survey; in future work photometric values will also be updated.

Most importantly, we have provided a validated framework and reliable numerical codes to perform forecasts. Although we restricted the analysis to primary probes and cross correlations, reliable codes are now available as a baseline for the analysis of other probes (such as the CMB) and for further extensions beyond the cosmologies investigated in this paper. Our results and numerical codes can be implemented for all future updates in the specific settings chosen for Euclid, and are currently used, for example, for the process of SPV.

---

Acknowledgments

Authors’ contributions: All authors have significantly contributed to this publication and/or made multi-year essential contributions to the Euclid project which made this work possible. The paper-specific main contributions of the lead authors are as follows: Coordinators and corresponding contacts: Thomas Kitching, Valeria Pettorino, Ariel Sánchez. Galaxy Clustering Task: Ariel G. Sánchez, Domenico Sapone, Carmelita Carbone, Santiago Casas, Katarina Markovic, Alkistis Pourtsidou, Safir Yahia-Cherif, Elisabetta Majerotto, Victoria Yankelevich. Weak Lensing Task: Vincenzo F. Cardone, Santiago Casas, Stefano Camera, Matteo Martinelli, Isaac Tutusaus, Thomas Kitching, Martin Kilbinger. Cross-correlations Task: Matteo Martinelli, Isaac Tutusaus, Alain Blanchard, Martin Kunz, Fabien Lacasa, Vincenzo F. Cardone, Sebastien Clesse, Stéphane Ilić, Martin Kilbinger, Marco Raveri, Ziad Sakr. Cosmological context: Valeria Pettorino, Eric Linder, Martin Kunz, Matteo Martinelli, Martin Kilbinger, Alkistis Pourtsidou, Ariel G. Sánchez, Thomas Kitching, Carmelita Carbone, Ziad Sakr. Paper editing: Valeria Pettorino, Ariel G. Sánchez, Thomas Kitching, Domenico Sapone, Matteo Martinelli, Santiago Casas, Katarina Markovic, Isaac Tutusaus, Stefano Camera, Alkistis Pourtsidou, Vincenzo Cardone, Ziad Sakr. We acknowledge support of a number of agencies and institutes that have supported the development of Euclid. A detailed complete list is available on the Euclid web site (http://www.euclid-ec.org). 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- and Raumfahrt, the Danish Space Research Institute, the Fundação para a Ciênca e a Tecnologia, the Ministerio de Economia y Competitividad, the National Aeronautics and Space Administration, the Netherlandse Onderzoekschool Voor Astronomie, the Norvegian Space Center, 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. Stefano Camera is supported by 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 by the “Departments of Excellence 2018–2022” Grant awarded by MIUR (L. 232/2016). Carmelita Carbone acknowledges financial support from the European Research Council through the Darklight Advanced Research Grant (n. 291521), and from the MIUR PRIN 2015 “Cosmology and Fundamental Physics: illuminating the Dark Universe with Euclid”. Santiago Casas and Safir Yahia-Cherif acknowledge support from french space agency CNES. Sebastien Clesse aknowledges support from the Belgian Fund for Research F.R.S-FNRS. Martin Kunz and Fabien Lacasa acknowledge financial support from the Swiss National Science Foundation. Eric Linder acknowledges support by NASA ROSES grant 12-EUCLID12-0004. Katarina Markovic carried out some of the work at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004), and acknowledges support from the UK Science & Technology Facilities Council through grant ST/N000668/1, and from the UK Space Agency through grant ST/N00180X/1. Matteo Martinelli acknowledges support from the D-ITP consortium, a program of the NWO that is funded by the OCW. Alkistis Pourtsidou is a UK Research and Innovation Future Leaders Fellow and also acknowledges support from the UK Science & Technology Facilities Council through grant ST/S000437/1. Domenico Sapone acknowledges financial support from the Fondecyt Regular project number 1200171. Isaac Tutusaus acknowledges support from the Spanish Ministry of Science, Innovation and Universities through grant ESP2017-89838-C3-1-R, and the H2020 programme of the European Commission through grant 776247. Victoria Yankelevich acknowledges financial support by the Deutsche Forschungsgemeinschaft through the Transregio 33 “The Dark Universe”, the International Max Planck Research School for Astronomy and Astrophysics at the Universities of Bonn and Cologne and the Bonn-Cologne Graduate School for Physics and Astronomy. Thomas D. Kitching acknowledges support from a Royal Society University Research Fellowship. Stéphane Ilić acknowledges financial support from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007–2013)/ERC Grant Agreement No. 617656 “Theories and Models of the Dark Sector: Dark Matter, Dark Energy and Gravity.” Elisabetta Majerotto acknowledges financial support from the Swiss National Science Foundation.

References

Appendix A: How you can perform your own code comparison

In this appendix we address the following two questions that the reader may have:

• I have another Fisher matrix code, how can I validate my own code and its output?

• How can I overplot Euclid marginalised posterior contours with other contours (e.g. from another experiment)?

In Sect. 3 we described the recipes required to calculate forecasts for Euclid’s main cosmological probes, and their combination; in Sect. 4 we showed how to validate different codes against each other and the presented plots that allow to visualise the comparison. We now provide the tools needed to build a Fisher matrix and compare it with representative results of this paper. To address the second question, we also provide plotting routines and reference Fisher matrices to reproduce a plot like the one in Fig. 10. These tools can be found in the public repository³⁰. We describe its content in detail below.

A.1. Content of the repository

The public Euclid IST:F repository contains representative Fisher matrices obtained following the recipes of Sect. 3. One Euclid IST:F reference Fisher matrix is provided for each probe and cosmology.

These matrices are in the All_results folder, divided in the optimistic and pessimistic subfolders, corresponding to the two sets of specifications used in this paper. These subfolders, further contain a flat and non-flat folder each. Within these, the files containing the Fisher matrices follow the naming convention For example, a Fisher matrix whose file name is corresponds to the combined Fisher matrix for photometric galaxy clustering (GCph), weak lensing (WL) and their cross correlation (XC), obtained in the flat case (flat) for the flat w₀, w_a cosmology (w0wa), using the optimistic set of specifications (optimistic). Each Fisher matrix contains a header listing the cosmological parameter of the matrix and their order.

EuclidISTF_Observables_CosmoModel_SpecificationsCase.txt

For example, a Fisher matrix whose file name is

EuclidISTF_GCph_WL_XC_w0wa_flat_optimistic.txt

In addition, the repository contains a set of Python scripts that can be used to plot the contour ellipses and compare them with user-provided matrices.

A.2. How to produce a Fisher matrix using IST:F input

The Euclid Fisher matrices contained in the public repository are obtained following the recipes described in Sect. 3 for the cases analysed in Sects. 4 and 5. Users who wish to validate their own code can build their Fisher matrices following these recipes, using the input files for cosmological quantities that are also included within the repository. These input files have been generated using CAMB, version August 2018. The input files containing the cosmological quantities used to create the Fisher matrices can be found online³¹. The link contains two folders: GC and WL which are described below.

GC folder. The GC folder contains the zip file ISTF_GC.zip. Inside you will find:

• the file fiducial_rate_growth.dat, containing three columns with z, growth rate f(z) – as defined in Eq. (25) – and growth function D(z) – as defined in Eq. (27);

• the subfolder Pk, containing the P_δδ(k, z), as defined in Eq. (28);

• the subfolder Pk-nw containing the P_nw(k, z) that appears in Eq. (83).

Each of the two subfolders contains:

• the fiducial folder, including the P_δδ(k, z) for the fiducial cosmology, with one file for each redshift at which this spectrum is computed. Each file has a suffix z_ii, where ii runs from 00 to 03 corresponding to the four redshift bins used for the Euclid GC_s probe z = {1.0,  1.2,  1.4,  1.65},

• the power spectra used to calculate the derivatives, grouped into folders with names par_step_eps_1p0E-2, with step being the name of the parameter, par being either pl indicating a positive step or mn indicating a negative step. These folders contain, in particular, the power spectra for each step in each cosmological parameter. The steps amplitudes are relative; for each parameter the absolute step will be fidpar(1 + eps) with fidpar the fiducial value of the parameter and eps the value of the relative step.

All the power spectrum files have three columns, being:

WL folder. The WL folder contains:

• the file scaledmeanlum-E2Sa.dat, which contains z (first column) and the ratio ⟨L⟩(z)/L_⋆(z) (second column) described in Eq. (111);

• ISTF_WL_Flat.zip, the matter power spectrum P_δδ(k, z) for the flat case;

• ISTF_WL_NonFlat.zip, the matter power spectrum P_δδ(k, z) for the non-flat case;

Each of these two zip files contains:

• the file pkz-Fiducial.txt, with the fiducial power spectrum;

• the spectra relative to the positive (negative) steps in cosmological parameters, e.g. ns_pl_eps_5p0E-2 indicates a step on n_s in the positive direction of amplitude 5.0 × 10⁻² (the steps amplitudes are relative; for each parameter the absolute step will be fidpar(1 + eps) with fidpar the fiducial value of the parameter and eps the value of the relative step).

Notice that with respect to the GC input files, here many more steps for numerical derivatives are computed; this is because most of the IST:F codes for WL use the derivative scheme described in Sect. 4.3. All the P_δδ(k, z) files have four columns:

How to include an external matrix

In order to compare their own Fisher matrix with the IST:F results, users can use the scripts provided in the following. For that to work, users should add their matrix in the folder corresponding to the case of interest, following the same organisation of the EuclidISTF matrices. The filenames of the new matrices, must follow the convention of the repository, i.e. with a prefix stating the name of the code used, and a suffix containing the observables included in the matrix, the cosmological model used and the specification settings. As an example, a Fisher matrix provided by the user for WL in the optimistic case, for a flat w₀, w_a cosmology should have the path

All_results/optimistic/flat/USERFISHER_WL_w0wa_flat_optimistic.txt

The full description of the settings, with all the corresponding specifications used to obtain the matrices, can be found in Sect. 4. Furthermore, each matrix file provided by users must contain the header listing the cosmological parameters and their order (see any of the EuclidISTF matrices for an example).

A.3. How to compare an external Fisher matrix

Within the main repository³², we also provide a comparison package which consists of:

• a common trunk, comparison.py found in the parent folder.

• three applications (appGC.py, appWL.py, appXC.py) found in PlottingScripts/ handling the different comparisons

• a set of plotting and analysis libraries found in PlottingScripts/.

The python script performing the comparison between an external Fisher matrix and the reference IST:F ones can be launched from the parent folder using The -o and -f arguments are optional.

python comparison.py [-o OBSERVABLES ][-f USERFISHER ].

The -o and -f arguments are optional.

USERFISHER are the labels of the user provided Fisher matrices. Notice that the user’s matrices must be placed in the correct folders of the repository and must follow the filename notation described above, with USERFISHER the prefix of the filename, i.e. the code name. The script will produce plots, similar to the comparison plots of Sect. 4, only for those cases for which an external matrix to compare with is found. In case no USERFISHER argument is provided, a test matrix will be used. This test matrix is contained within the repository only for the cross-correlation comparison cases (see Sect. 4.4).

The OBSERVABLES can be one or more of the cases considered in Sect. 4, i.e.

• GC for spectroscopic galaxy clustering.

• WL.

• XC.

The script can perform one, two or all cases at once. The latter is the default selection if no observable is specified.

Results will be placed in results folders in each of the subcases locations, e.g. in the parent folder one will find. These results folders are automatically generated by the comparison script.

All_results/optimistic/flat/results_XC

These results folders are automatically generated by the comparison script.

As an example, let’s assume that a user wants to compare a WL matrix computed in the optimistic case, for a flat w₀, w_a cosmology. The user will add the matrix as

All_results/optimistic/flat/MyName_WL_w0wa_flat_optimistic.txt

Running the script with the command will compare the new matrix with the Euclid IST:F for this specific case and cosmology, which will be found in

python comparison.py -o WL -f MyName

will compare the new matrix with the Euclid IST:F for this specific case and cosmology, which will be found in

All_results/optimistic/flat/results_WL.

Notice that, as said above, the script does not produce comparison results for those cases in which no external matrix is provided. In the example above, where only the WL, optimistic matrix for w₀, w_a is provided, the script will produce the results only for this specific case also running

python comparison.py -f MyName.

A.4. How to obtain results for the IST:F and external matrices

The python script performing the IST:F ellipses plot can be launched from the parent folder using The -o, -f, -n and -c arguments are optional.

python ellipses.py [-o OBSERVABLES ][-f USERFISHER ] [-n USERNAME ] [-c USERCOLOR].

The -o, -f, -n and -c arguments are optional.

In case no arguments are provided, the script will reproduce Fig. 10 from the paper.

With the option -o a list of OBSERVABLES can be specified, separated by spaces. Available options are:

• All produces a plot like Fig. 10 for different observables;

• XC produces a plot focusing on XC correlations, like Fig. 10;

• GC produces a plot containing only GC_s contours;

• WL produces a plot containing only WL contours.

With the option -f, a USERFISHER can be specified, which is the relative path to the Fisher matrix provided by the user. Notice that the user’s matrices must be placed in the correct folders of the repository and must follow the filename notation described above.

If an external user Fisher matrix is provided, the script will overplot it on top of the other ellipses using a colour and a label name provided by the user. Default colour: black, default name: USER.

The user can change the colour and the label name of its Fisher matrix, with the option -c and -n, respectively.

Ellipses plots will be placed in the following folder, e.g.

All_results/optimistic/flat/results_ellipses.

Additionally to the elliptical contours, the routine also produces a file called: w0waCDM-flat-optimistic_bounds.txt containing the relative errors on all the parameters of interest and a file w0waCDM-flat-optimistic_FoMs.txt containing the FoM for w₀, w_a corresponding to each of the Fisher matrices which were plotted.

The .ini files inside the directory PlottingScripts/ store all the parameters of the plot, and provide a way for the user to gain finer control of the plotting routine, if needed.

All Tables

All Figures