Issue 
A&A
Volume 642, October 2020



Article Number  A191  
Number of page(s)  66  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/202038071  
Published online  20 October 2020 
Euclid preparation
VII. Forecast validation for Euclid cosmological probes
^{1}
Institut de Recherche en Astrophysique et Planétologie (IRAP), Université de Toulouse, CNRS, UPS, CNES, 14 Av. Edouard Belin, 31400 Toulouse, France
^{2}
INFNSezione di Torino, Via P. Giuria 1, 10125 Torino, Italy
^{3}
Dipartimento di Fisica, Universitá degli Studi di Torino, Via P. Giuria 1, 10125 Torino, Italy
^{4}
Dipartimento di Fisica, Universitá degli Studi di Milano “Aldo Pontremoli”, Via Celoria 16, 20133 Milano, Italy
^{5}
INFNSezione di Milano, Via Celoria 16, 20133 Milano, Italy
^{6}
INAFIASF Milano, Via Alfonso Corti 12, 20133 Milano, Italy
^{7}
INAFOsservatorio Astronomico di Roma, Via Frascati 33, 00078 Monte Porzio Catone, Italy
^{8}
AIM, CEA, CNRS, Université ParisSaclay, Université Paris Diderot, Sorbonne Paris Cité, 91191 GifsurYvette, France
^{9}
CEICO, Institute of Physics of the Czech Academy of Sciences, Na Slovance 2, Praha 8, Czech Republic
^{10}
CEA Saclay, DFR/IRFU, Service d’Astrophysique, Bat. 709, 91191 GifsurYvette, France
^{11}
Institut d’Astrophysique de Paris, 98bis Boulevard Arago, 75014 Paris, France
^{12}
Mullard Space Science Laboratory, University College London, Holmbury St Mary, Dorking, Surrey RH5 6NT, UK
^{13}
Université de Genève, Département de Physique Théorique and Centre for Astroparticle Physics, 24 quai ErnestAnsermet, 1211 Genève 4, Switzerland
^{14}
University of California, Berkeley, CA 94720, USA
^{15}
Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth PO1 3FX, UK
^{16}
Institute Lorentz, Leiden University, PO Box 9506, Leiden 2300 RA, The Netherlands
^{17}
School of Physics and Astronomy, Queen Mary University of London, Mile End Road, London E1 4NS, UK
^{18}
Université St Joseph, UR EGFEM, Faculty of Sciences, Beirut, Lebanon
^{19}
Max Planck Institute for Extraterrestrial Physics, Giessenbachstr. 1, 85748 Garching, Germany
^{20}
Departamento de Física, FCFM, Universidad de Chile, Blanco Encalada 2008, Santiago, Chile
^{21}
Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans, s/n, 08193 Barcelona, Spain
^{22}
Institut d’Estudis Espacials de Catalunya (IEEC), 08034 Barcelona, Spain
^{23}
ArgelanderInstitut für Astronomie, Universität Bonn, Auf dem Hügel 71, 53121 Bonn, Germany
^{24}
INAFOsservatorio Astronomico di Brera, Via Brera 28, 20122 Milano, Italy
^{25}
INAFOsservatorio Astronomico di Brera, Via E. Bianchi 46, 23807 Merate, Italy
^{26}
Departamento de Astrofísica, Universidad de La Laguna, 38206 La Laguna, Tenerife, Spain
^{27}
Instituto de Astrofísica de Canarias, Calle Vía Làctea s/n, 38204, San Cristóbal de la Laguna, Tenerife, Spain
^{28}
INAFOsservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via Piero Gobetti, 93/3, 40129 Bologna, Italy
^{29}
Dipartimento di Fisica e Astronomia, Universitá di Bologna, Via Gobetti 93/2, 40129 Bologna, Italy
^{30}
INFNSezione di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy
^{31}
UniversitätsSternwarte München, Fakultät für Physik, LudwigMaximiliansUniversität München, Scheinerstrasse 1, 81679 München, Germany
^{32}
INAFOsservatorio Astronomico di Trieste, Via G. B. Tiepolo 11, 34131 Trieste, Italy
^{33}
INAFOsservatorio Astrofisico di Torino, Via Osservatorio 20, 10025 Pino Torinese (TO), Italy
^{34}
APC, AstroParticule et Cosmologie, Université Paris Diderot, CNRS/IN2P3, CEA/lrfu, Observatoire de Paris, Sorbonne Paris Cité, 10 Rue Alice Domon et Léonie Duquet, 75205 Paris Cedex 13, France
^{35}
Department of Astronomy, University of Geneva, ch. d’Écogia 16, 1290 Versoix, Switzerland
^{36}
INFNSezione di Roma Tre, Via della Vasca Navale 84, 00146 Roma, Italy
^{37}
Department of Mathematics and Physics, Roma Tre University, Via della Vasca Navale 84, 00146 Rome, Italy
^{38}
INAFOsservatorio Astronomico di Capodimonte, Via Moiariello 16, 80131 Napoli, Italy
^{39}
Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Rua das Estrelas, 4150762 Porto, Portugal
^{40}
INFNBologna, Via Irnerio 46, 40126 Bologna, Italy
^{41}
Dipartimento di Fisica e Scienze della Terra, Universitá degli Studi di Ferrara, Via Giuseppe Saragat 1, 44122 Ferrara, Italy
^{42}
INAF, Istituto di Radioastronomia, Via Piero Gobetti 101, 40129 Bologna, Italy
^{43}
Université Côte d’Azur, Observatoire de la Côte d’Azur, CNRS, Laboratoire Lagrange, Nice, France
^{44}
Institut de Física d’Altes Energies IFAE, 08193 Bellaterra, Barcelona, Spain
^{45}
Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências, Universidade de Lisboa, Tapada da Ajuda 1349018, Lisboa, Portugal
^{46}
Department of Physics “E. Pancini”, University Federico II, Via Cinthia 6, 80126 Napoli, Italy
^{47}
INFN section of Naples, Via Cinthia 6, 80126 Napoli, Italy
^{48}
INAFOsservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50125 Firenze, Italy
^{49}
Centre National d’Etudes Spatiales, Toulouse, France
^{50}
Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK
^{51}
University of Nottingham, University Park, Nottingham NG7 2RD, UK
^{52}
ESAC/ESA, Camino Bajo del Castillo, s/n., Urb. Villafranca del Castillo, 28692 Villanueva de la Cañada, Madrid, Spain
^{53}
Université de Lyon, 69622 Lyon, France
^{54}
Université de Lyon 1, Villeurbanne, France
^{55}
CNRS/IN2P3, Institut de Physique Nucléaire de Lyon, Villeurbanne, France
^{56}
Departamento de Física, Faculdade de Ciências, Universidade de Lisboa, Edifício C8, Campo Grande, 1749016 Lisboa, Portugal
^{57}
Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências, Universidade de Lisboa, Campo Grande, 1749016 Lisboa, Portugal
^{58}
AixMarseille Univ, CNRS, CNES, LAM, Marseille, France
^{59}
Department of Physics, Oxford University, Keble Road, Oxford OX1 3RH, UK
^{60}
INFNPadova, Via Marzolo 8, 35131 Padova, Italy
^{61}
Department of Physics & Astronomy, University of Sussex, Brighton BN1 9QH, UK
^{62}
Centre for Extragalactic Astronomy, Department of Physics, Durham University, South Road, Durham DH1 3LE, UK
^{63}
AixMarseille Univ, CNRS/IN2P3, CPPM, Marseille, France
^{64}
Leiden Observatory, Leiden University, Niels Bohrweg 2, 2333 CA Leiden, The Netherlands
^{65}
von Hoerner & Sulger GmbH, SchloßPlatz 8, 68723 Schwetzingen, Germany
^{66}
MaxPlanckInstitut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany
^{67}
Department of Physics and Helsinki Institute of Physics, Gustaf Hällströmin katu 2, 00014 University of Helsinki, Finland
^{68}
Institut de Physique Nucléaire de Lyon, 4, Rue Enrico Fermi, 69622 Villeurbanne Cedex, France
^{69}
Institute of Theoretical Astrophysics, University of Oslo, PO Box 1029 , Blindern 0315, Oslo, Norway
^{70}
Istituto Nazionale di Astrofisica (INAF) – Osservatorio di Astrofisica e Scienza dello Spazio (OAS), Via Gobetti 93/3, 40127 Bologna, Italy
^{71}
Institute for Computational Cosmology, Department of Physics, Durham University, South Road, Durham DH1 3LE, UK
^{72}
INAFOsservatorio Astronomico di Padova, Via dell’Osservatorio 5, 35122 Padova, Italy
^{73}
University of Paris Denis Diderot, University of Paris Sorbonne Cité (PSC), 75205 Paris Cedex 13, France
^{74}
Sorbonne Université, Observatoire de Paris, Université PSL, CNRS, LERMA, 75014 Paris, France
^{75}
Observatoire de Sauverny, Ecole Polytechnique Fédérale de Lausanne, 1290 Versoix, Switzerland
^{76}
INFN, Sezione di Lecce, Via per Arnesano, 193, 73100 Lecce, Italy
^{77}
Department of Mathematics and Physics E. De Giorgi, University of Salento, Via per Arnesano, CPI93, 73100 Lecce, Italy
^{78}
Perimeter Institute for Theoretical Physics, Waterloo, Ontario N2L 2Y5, Canada
^{79}
Department of Physics and Astronomy, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada
^{80}
Centre for Astrophysics, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada
^{81}
Space Science Data Center, Italian Space Agency, Via del Politecnico snc, 00133 Roma, Italy
^{82}
European Space Agency/ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands
^{83}
Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA 91109, USA
^{84}
Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Avenida Complutense 40, 28040 Madrid, Spain
^{85}
Universidad Politécnica de Cartagena, Departamento de Electrónica y Tecnología de Computadoras, 30202 Cartagena, Spain
^{86}
Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700 AV Groningen, The Netherlands
^{87}
IFPU, Institute for Fundamental Physics of the Universe, Via Beirut 2, 34151 Trieste, Italy
^{88}
SISSA, International School for Advanced Studies, Via Bonomea 265, 34136 Trieste TS, Italy
^{89}
INFN, Sezione di Trieste, Via Valerio 2, 34127 Trieste TS, Italy
^{90}
Infrared Processing and Analysis Center, California Institute of Technology, Pasadena, CA 91125, USA
^{91}
CURL & IRMP, Louvain University, 2 Chemin du Cyclotron, 1348 LouvainlaNeuve, Belgium
^{92}
NaXys, University of Namur, Rempart de la Vierge 8, 5000 Namur, Belgium
Received:
2
April
2020
Accepted:
15
July
2020
Aims. The Euclid space telescope will measure the shapes and redshifts of galaxies to reconstruct the expansion history of the Universe and the growth of cosmic structures. The estimation of the expected performance of the experiment, in terms of predicted constraints on cosmological parameters, has so far relied on various individual methodologies and numerical implementations, which were developed for different observational probes and for the combination thereof. In this paper we present validated forecasts, which combine both theoretical and observational ingredients for different cosmological probes. This work is presented to provide the community with reliable numerical codes and methods for Euclid cosmological forecasts.
Methods. We describe in detail the methods adopted for Fisher matrix forecasts, which were applied to galaxy clustering, weak lensing, and the combination thereof. We estimated the required accuracy for Euclid forecasts and outline a methodology for their development. We then compare and improve different numerical implementations, reaching uncertainties on the errors of cosmological parameters that are less than the required precision in all cases. Furthermore, we provide details on the validated implementations, some of which are made publicly available, in different programming languages, together with a reference trainingset of input and output matrices for a set of specific models. These can be used by the reader to validate their own implementations if required.
Results. We present new cosmological forecasts for Euclid. We find that results depend on the specific cosmological model and remaining freedom in each setting, for example flat or nonflat spatial cosmologies, or different cuts at nonlinear scales. The numerical implementations are now reliable for these settings. We present the results for an optimistic and a pessimistic choice for these types of settings. We demonstrate that the impact of crosscorrelations is particularly relevant for models beyond a cosmological constant and may allow us to increase the dark energy figure of merit by at least a factor of three.
Key words: cosmology: observations / cosmological parameters / cosmology: theory
© Euclid Collaboration 2020
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
Euclid^{1} will explore the expansion history of the Universe and the evolution of largescale cosmic structures by measuring shapes and redshifts of galaxies, covering 15 000 deg^{2} 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 threemirror anistigmat and it has two instruments that observe in optical and nearinfrared wavelengths. In the optical, highresolution 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 nearinfrared, 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 nearinfrared instruments share a common fieldofview of 0.53 deg^{2}.
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 welldefined 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, uptodate 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 interscience taskforce for forecasting (IST:F) and it presents the results of the code comparison for forecasts based on Fishermatrix 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 crosscorrelation (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 reassessment 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 nonparametric model was chosen in redshift and in scale, and prior information on the model’s nuisance parameters included the expected performance of spectroscopic galaxygalaxy 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 smallscale (high kmode) matter power spectrum on WL twopoint 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 nonlinear 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 largescale 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^{2} (SKA) is expected to deliver a wealth of radio wavelength data. Its Phase 1 Mid array (SKA1MID) is going to be commissioned on similar timescales to Euclid and is offering a unique opportunity for multiwavelength synergies. As discussed in Kitching et al. (2015), a multitude of crosscorrelation statistics will be provided by comparing Euclid and SKA clustering and weak lensing data, with the additional advantage of the crosscorrelations 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). Crosscorrelation forecasts for an SKA1MID 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 crosscorrelations of shear maps between Euclid and Phase 2 of SKAMID assuming 15 000 deg^{2} 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^{3} (DESI), the Large Synoptic Survey Telescope^{4} (LSST), and the Wide Field Infrared Survey Telescope^{5} (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 largescale 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 eASTROGAM 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 intersurvey 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 Fishermatrix 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 codecomparison 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^{6} 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
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_{0}/a − 1, where a_{0} corresponds to the presentday value of the scale factor, normalised as a_{0} = 1; we convert the t variable into a or z. The Hubble expansion rate, H(z), can be expressed as
with H_{0} ≡ H(z = 0) being the Hubble parameter today, which is commonly written as
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
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^{2}/[a^{2}(z)H^{2}(z)]. With these definitions, Eq. (1) takes the form
where the sum is over all N species considered in the model. We use presentday 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^{2}, which can be redshiftdependent. For the case in which the equation of state parameter is constant in time the energy conservation equation implies
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 massiveneutrinos can be described by a varying equation of state parameter w_{ν}, which matches w_{γ} at early times and w_{c} when they become nonrelativistic. 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 massiveneutrinos as part of the total matter contribution, with Ω_{m, 0} = Ω_{c, 0} + Ω_{b, 0} + Ω_{ν, 0}, since they are nonrelativistic 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 timeindependent 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 redshiftdependent equation of state parameter w_{DE}(z). A commonly used and welltested parameterisation of the time dependence is
where w_{0} 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
Using Eqs. (7) and (9) in Eq. (1), the function E(z) defined in Eq. (3) becomes
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_{0} = −1 and w_{a} = 0, in which case the function E(z) takes the form
2.2. Distance measurements
The comoving distance to an object at redshift z can be computed as
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
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
for a spatially flat universe (K = 0), this becomes
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 nonrelativistic matter on subhorizon scales, can be described by ideal fluid equations (Peebles 1980). Density fluctuations for a given component i are characterised by the density contrast
which quantifies the deviations of the density field ρ_{i}(x, z) around the mean spatial density over space, where x is a threedimensional comoving coordinate at a redshift z. To describe these fluctuations statistically, it is convenient to work in Fourier space by decomposing δ into plane waves,
The power spectrum, P_{i}(k, z), for the generic component i is defined implicitly as
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
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_{0} 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
for a generic component i. In the early Universe during radiation domination, curvature perturbations with comoving scales smaller than the horizon are suppressed, whereas superhorizon 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 matterradiation equality is imprinted on the shape of the transfer function, and hence on the matter power spectrum.
The growth of density fluctuations obeys a secondorder 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
In the ΛCDM scenario with no massiveneutrinos, this equation can be written in terms of the redshift as
where we neglect dark energy perturbations, the prime refers to the derivative with respect to z, and Ω_{m}(z) is given by
The solutions δ_{m}(k, z) of Eq. (22), at late times, are scaleindependent, which motivates the introduction of the growth factor D(z) through
where z_{i} is an arbitrary reference redshift in the matterdominated era. A useful quantity is the growth rate parameter, defined as
Using these definitions^{7}, the growth rate satisfies a firstorder differential equation
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
In ΛCDM the latetime matter growth is scaleindependent, so that the transfer function 𝒯_{m}(k, z) can be split into a scaledependent part T_{m}(k) (normalised so that T_{m} → 1 for k → 0) and the scaleindependent growth factor D(z) introduced above. A convenient way to express the power spectrum defined in Eq. (20) for matter is
with the normalisation constant
where W_{TH}(x) = 3(sin x − x cos x)/x^{3} is the Fourier transform of the tophat filter, and R_{8} = 8 h^{−1} Mpc. The pivot scale k_{0} and the amplitude of the scalar mode A_{s} are absorbed into the normalisation, which is designed to give a desired value of σ_{8}, the root mean square (rms) of presentday linearly evolved density fluctuations in spheres of 8 h^{−1} Mpc, which is given by
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
Since we consider models beyond ΛCDM, we need to specify the behaviour of the dark energy perturbations. We use a minimallycoupled 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 multifield scenarios, Peirone et al. 2017), and to achieve this we use the Parameterised PostFriedmann (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
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})^{8} 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
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 nonvanishing for relativistic species. The second free function Σ_{MG} describes the deflection of light and is defined through
An alternative option is the pair (μ_{MG}, η), with η defined as
which reduces to^{9}
for negligible shear, as happens at the (low) redshifts of interest for Euclid. Equivalently, combining Eqs. (33) and (35), one has
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:
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 nonlinear 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 semianalytical 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
with
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 massiveneutrinos 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.
Fig. 1.
Function μ_{ν}(k). It represents the scale dependent correction to f(z) in Eqs. (39) and (40), evaluated here at ∑_{i}m_{ν, i} = 0.06 and 0.15 eV. 
Given the current upper limits on the total neutrino mass (PalanqueDelabrouille et al. 2020), the growth suppression is of the order of 0.6% at maximum, and mainly affects scales k > 0.1 h Mpc^{−1}. For simplicity, we ignore this subpercent effect in our forecasts, assuming f to be scale independent for values ∑_{i}m_{ν, i} = 0.06, 0.15 eV tested, and computed in the massless limit. For ∑_{i}m_{ν, i} = 0.06 eV solar oscillation experiments constrain neutrinos to be in a normal hierarchy, and for ∑_{i}m_{ν, 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 redshiftdependent one, 𝒯(k; z), while keeping the scaleindependent linear growth factor D(z) as given in the absence of massiveneutrino freestreaming (Takada et al. 2006; Takada 2006; Eisenstein & Hu 1997; see Eq. 25 of Eisenstein & Hu 1997)^{10}. 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 massiveneutrino 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 (JelicCizmek et al. 2020). Furthermore, magnification can also be important for the crosscorrelations between weak lensing and galaxy clustering of photometricallydetected 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 bestfit 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_{0} = −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;

σ_{8}, 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 nonflat 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_{0} 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
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
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.
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 parameterindependent. In both cases, one of the two terms in Eq. (43) is nonvanishing. 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^{11},
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
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 offdiagonal elements in the parameter covariance matrix,
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. Reinverting 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 loglikelihood function is then locally close to the function of a hyperdimensional ellipsoid, defined by the Fisher matrix of Eq. (42), and the projection on any twodimensional subspace is simply an ellipse. The parametric form of the confidence level ellipses is defined by the semiminor and semimajor 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 ycomponent to the xcomponent of the larger eigenvector (as given by Eq. (40.72) of the Particle Data Group’s Review of Particle Properties section on statistics^{12}). They are:
where A is a constant factor defined as A^{2} = 2.3, 6.18, 11.8 for twoparameter 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
The above formula was utilised in the original dark energy FoM definition from Albrecht et al. (2006), which considers the w_{0}, w_{a} parameterization defined in Sect. 2.1:
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 ddimensional vector p of random variables as
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 crosscorrelations 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:
which is 1 if parameters are fully uncorrelated. Offdiagonal nonzero 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
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 crosscovariance 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)
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:
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.
Parameter values of our fiducial cosmological model, corresponding to those of Planck Collaboration XIV (2016), both in the baseline ΛCDM case and in the extensions considered.
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 simulationderived 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.
Specifications for the spectroscopic galaxy redshift survey.
Expected number density of observed Hα emitters for the Euclid spectroscopic survey.
Specifications for the Euclid photometric weak lensing survey.
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
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^{−4}.
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,
and
where for a given α, β pair where α ≤ 2 and similarly for β^{13}. Since
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
where the correlation coefficient ; we also note that from Eq. (46). From this, one has
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
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^{−4}, which is a reasonable aim for a good numerical implementation.
We note that the FoM accuracy is simply given by
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
For many parameter pairs, for example (w_{0}, 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
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^{−4}. In other words, one certainly does not need to push near machine precision. We note here that for twosided or threepoint finite differences for evaluation of derivatives with a parameter stepsize ϵ, the errors go as ϵ^{2}, while fivepoint differences go as ϵ^{4}.
3.2. Recipe for spectroscopic galaxy clustering
This section describes the forecasting procedure recommended for our galaxy clustering observable – the full, anisotropic, and redshiftdependent 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:
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 latetime nonlinearities and redshiftspace distortions (RSDs), which damp the baryon acoustic oscillation (BAO) signal in the galaxy power spectrum, and cause the socalled fingersofGod 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 redshiftdependent set of parameters is also specific to galaxy clustering. Anisotropies in the power spectrum with respect to the lineofsight 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σ_{8}(z) ≡ f(z)σ_{8}(z) defined in Sect. 2.3; note that here we generalise Eq. (30) to be redshiftdependent by changing the redshift argument to the power spectrum in that integral. The additional parameters bσ_{8}(z) ≡ b(z)σ_{8}(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σ_{8}(z) obtained at each redshift provide us with darkenergymodelindependent 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σ_{8}(z)], P_{s}(z), σ_{p}(z_{i}), and σ_{v}(z_{i}), we transform the constraints on the remaining set of redshiftdependent 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}, σ_{8}, ∑m_{ν}} with or without the extensions {Ω_{DE, 0}, w_{0}, 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 wavevector k and the lineofsight direction , and z is the redshift of the density field^{15}. 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 nonlinear, regime of perturbation theory.
There are five main effects that we need to model:

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

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 lineofsight.

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,

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

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 nonlinearities 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:
where P_{δδ}(k; z) (see Eq. (28)) is the linear matter power spectrum, which, when multiplied by the square of the redshiftdependent effective bias of the galaxy sample, b^{2}(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 noncosmological contribution due to the lineofsight component of the peculiar velocity
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 redshiftspace galaxy power spectra (here labelled as “zs” for redshiftspace), is given by Kaiser (1987)
where σ_{8}(z) is defined in Eq. (30), but generalised to be redshiftdependent. Written in this way, Eq. (72) illustrates that the quantities controlling the anisotropies in the two dimensional power spectrum are the combinations b(z)σ_{8}(z) and f(z)σ_{8}(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
where the factor F_{z}(k, μ; z) is given by
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),
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^{−1}(z) (Wang et al. 2013).
As described before, P_{s}(z) is a scaleindependent 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 shotnoise term as a free parameter (see e.g. GilMarí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 wavevector k in the directions parallel and perpendicular to the lineofsight, k_{∥} and k_{⊥}, as
where the coefficients q_{⊥} and q_{∥} are given by the ratio of the angulardiameter distance and the Hubble parameter in the true and reference cosmologies as
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
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)
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
where all k = k(k_{ref}, μ_{ref}) and μ = μ(μ_{ref}).
3.2.2. Nonlinear 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 nonlinear scales into the analysis. However, this requires modelling nonlinear 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 nonlinear 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 nonlinearities, 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
However, we keep them separate to allow for some theoretical uncertainty in our practical model of nonlinearities.
The first effect of nonlinearity is that the description of RSD requires an additional factor that we model as a Lorentzian, which accounts for the fingerofGod (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 socalled “dewiggled” power spectrum, P_{dw}(k, μ; z), (see Wang et al. 2013). Therefore we model the true nonlinear underlying redshiftspace power spectrum as (cf. Eq. (72) above)
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
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 “nowiggle” 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
is the nonlinear damping factor of the BAO signal derived by Eisenstein et al. (2007). The parameter σ_{v}(z) controls the strength of the nonlinear damping of the BAO signal in all directions in threedimensions.
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
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 nonlinear 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 nonlinear 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 nonlinear 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 nonlinear model for the observed power spectrum of Hα emitters. We can now write down the full observed, nonlinear anisotropic galaxy power spectrum as
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 redshiftdependent 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 redshiftdependent 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_{0} and w_{a}. On the other hand, it is the redshiftdependent parameter fσ_{8}(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 σ_{8}(z) = σ_{8}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^{16}. 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,
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
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
Here, P_{obs}(k, μ; z) is the full observed twodimensional 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:
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:
We note here that for the redshiftdependent 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
since
The sum in Eq. (92) then reduces to a single term with i = j:
This includes both the offdiagonal terms containing the covariance between different redshift bins of the same redshiftdependent parameter, and the covariance between different zbins of different parameters, leaving only the nonzero covariances between redshiftdependent parameters within the same redshift bin. In other words,
if α(z_{i}) and β(z_{j}) are redshiftdependent parameters, and in the following we abbreviate these subscripts to ij.
Explicitly, the modelindependent 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 zindependent parameters falls to 4). Schematically, let us assume that the Fisher matrix for one redshift bin is
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
where we must sum over the contributions to the shape parameter elements from each redshift, but treat the redshift bins as creating separate redshiftdependent 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:
(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σ_{8}(z_{i}), bσ_{8}(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 nonlinear effects. The full list of parameters, {θ_{i}}, of the initial Fisher matrix then includes redshiftindependent parameters, {θ_{zi}}, and redshiftdependent ones {θ_{zd}}: {θ_{i}}={θ_{zd}}+{θ_{zi}}. These parameters are listed in Eq. (69) and recalled below for convenience:

5 × N_{z} redshiftdependent parameters, {θ_{zd}}={ln D_{A}(z_{i}), ln H(z_{i}), ln[fσ_{8}(z_{i})], ln[bσ_{8}(z_{i})], P_{s}(z_{i})}

4 cosmological redshiftindependent parameters and 2 nonlinear 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σ_{8}(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 GConly 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_{0}, w_{a}, σ_{8}, γ}.
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
and so
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 DEindependent 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})σ_{8}(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_{0} and w_{a}. We therefore first marginalise out the observational nuisance parameters, bσ_{8}(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} redshiftdependent 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 largescale 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 tracefree 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 largescale structure has a zero mean because we assume isotropy and homogeneity of the Universe, but its twopoint correlation function and its power spectrum (as well as higherorder correlations) contain cosmological information that probes both the background evolution of the Universe and the growth of cosmic structure; these twopoint 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 twopoint 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 sphericalharmonic measurement of the twopoint 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 threedimensional shear field. This redshift binning is known as “tomography” and is required in order to achieve highprecision 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_{0}, w_{a}, σ_{8}, γ} 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 precomputation, 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 largescale 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 smallscale 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^{−1}) 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 broadband 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 nonlinear 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 largescale structure along the lineofsight. For an individual galaxy we express this, to linear order, as
where γ is the cosmological shear and ϵ^{I} is the intrinsic (unlensed) ellipticity. This ellipticity is a spin2 quantity with zero mean over a large survey area and a nonzero twopoint 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 twopoint 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
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
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 nonflat 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 nonLimber approximated case to hyperspherical 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 twopoint 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 twopoint correlation function of Eq. (102), which results in four terms
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 intrinsicintrinsic (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 linearalignment model:
The weight function for IA can be conveniently written as
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
where the function ℱ_{IA}(z) reads
We refer to this as the extended nonlinear alignment (eNLA) model. For ℱ_{IA}(z) = 1, the model reduces to the socalled nonlinear alignment (NLA) model widely used in the literature (Bridle & King 2007). In Eq. (111) ⟨L⟩(z) and L_{⋆}(z) are the redshiftdependent 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
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^{17} (Laureijs et al. 2011; Smail et al. 1994):
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 equipopulated bins, which are
with and . We note that, since we use 10 equipopulated 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
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 selfcalibrated. These are nevertheless not entirely flexible and can be inferred independently, with some errors, for instance with clusteringbased 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 selfcalibrated 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 nonzero for autocorrelation (intrabin) power spectra, but is zero for crosscorrelation (interbin) power spectra, because ellipticities of galaxies at different redshifts should be uncorrelated. This term can be written as
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 flatsky and Limber approximations, we can write the total observed tomographic shear angular power spectrum as
which, in an expanded form reads
where the first three terms give the contributions from the cosmological (theoretical) shear, cosmological shearIA, 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
having defined the following kernel functions
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 nonzero only over the redshift range where n_{i}(z) and n_{j}(z) are both nonvanishing.
3.3.2. Nonlinear scales for cosmic shear
The standard cosmic shear power spectrum is sensitive to kmodes in the range 0 < k ≤ 7 h Mpc^{−1} (Taylor et al. 2018a), which therefore necessitates a modelling of highk behaviour of the matter power spectrum. Whilst it is possible to exclude some highk 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 kmode, a method known as “kcut 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 nonlinear power spectrum can decrease the FoM by up to 40%, depending on the choices of model and prior. We therefore need to model the highk 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 nonlinear, highk, 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 nonlinear 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, nonlinear 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 socalled 1 and 2halo terms. The 1halo 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 2halo term describes correlations between distinct dark matter haloes, dominates on scales larger than the ones where the 1halo 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 2halo terms are computed. Additionally, both methods depend on parameters that are fitted to reproduce the matter power spectrum measured from stateoftheart (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 CDMdominated universe) to models with constant equation of state dark energy. The 1 and 2halo 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 highresolution Nbody simulations. The refined HF method provides an accurate prediction of the nonlinear matter power spectrum over a wide range of wavenumbers k ≤ 30 h Mpc^{−1} at redshifts 0 ≤ z ≤ 10: the precision is 5% for k ≤ 1 h Mpc^{−1} at 0 ≤ z ≤ 10 and 10% for 1 h Mpc^{−1} ≤ k ≤ 10 h Mpc^{−1} at 0 ≤ z ≤ 3. Since massive neutrinos are included in the analysis, we consider the standard HF method with the extension for massiveneutrinos 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 nonstandard scenarios including massiveneutrinos and some modified gravity theories. It takes into account the halo mass function and the halo density profile for the 1halo 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 Nbody simulations. The resulting nonlinear matter power spectrum is accurate to ∼5% up to k ≤ 10h Mpc^{−1} 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 nonlinear recipes available at the time of this work, which also include massiveneutrino effects.

Bird. Bird et al. (2012) modified the nonlinear HF version from Smith et al. (2003) to account for massive neutrino effects in the nonlinear evolution of the total matter power spectrum. This is a fitting to Nbody 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 massiveneutrino effects at the nonlinear level. Also in this case, Nbody simulations, which include a massiveneutrino 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 massiveneutrinos is not adopted in this work.

TakaBird. As reported in the documentation of the CAMB webpage^{18}, on March 2014 modified massiveneutrino parameters were implemented in the nonlinear 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 massiveneutrinos are included. The fitting parameters that account for nonlinear corrections in the presence of massiveneutrinos 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_{0}, 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 TakahashiBird, TakaBird, power spectrum is used to correct for nonlinearities 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
In the presence of a nonstandard 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
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
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
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 nonGaussian contributions (Takada & Hu 2013; Cooray & Hu 2001; Hamilton et al. 2006; Hu & Kravtsov 2003; Kayo et al. 2013). The nonGaussian terms involve the convergence trispectrum, and there are large theoretical uncertainties on how to model this quantity, both in the quasilinear and in the nonlinear 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 signaltonoise 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 nonGaussian 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 nonGaussian modelling, we decide not to implement full nonGaussian 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 signaltonoise is not directly linked to constraints on parameters beyond a global amplitude; parameter sensitivities are generally scaledependent (Copeland et al. 2018) and therefore respond differently to smallℓ cuts and the correlation introduced by nonGaussian covariances. However, the amplitude of the impact of nonGaussian 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 righthand 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
where C^{ϵϵ}(ℓ) is defined in Eq. (118). No crosscorrelation 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
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:
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_{0}, w_{a}, h}, whilst in the second term the derivatives of the power spectrum may be further decomposed into
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_{0}, 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 crosscorrelation 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 crosscorrelation (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 crosscorrelation of photometric galaxy clustering and cosmic shear. Both latter probes are quantified as 2D projected observables, which allows us to compute the crossobservables and crosscovariance matrix in a simple, unified scheme. The crosscorrelation 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 crosscorrelation 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
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
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
where we have defined the overall weight function , which includes intrinsic alignment contributions
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
and represent, respectively, the cosmic shear power spectrum, the sheargalaxy power spectrum and the galaxygalaxy power spectrum (using photometric galaxy clustering).
3.4.2. Nonlinear scales for probe combination
For the modelling of nonlinear scales in the cosmic shear, photometric galaxy clustering, and crosscorrelations 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 crosscorrelations, 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_{0} 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 crosscorrelations 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 crosscorrelation, 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
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 fourthorder covariance as
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
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 signaltonoise 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 nonGaussian 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 crosscorrelation, is given in the case of the fourth order covariance by
and in the case of the secondorder covariance as
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 autocorrelation, crosscorrelation between galaxy clustering and cosmic shear galaxy clustering autocorrelation. The indices in the sum ij and mn run over all unique pairs of tomographic bins (i ≤ j, m ≤ n) for the cosmic shear autocorrelation and the galaxy clustering autocorrelation, while all pairs of tomographic bins are considered to take into account all the crosscorrelations 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
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 crosscorrelation 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 crosscorrelations. 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 nonlinear 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 nonlinear corrections and how we project our Fisher matrices to derive constraints on commonly used dark energy and gravity parameters, such as w_{0}, 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 crosscorrelation 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.
Summary of the codes involved in the comparison.
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 higherorder 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 nonlinear 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 σ_{8} 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^{19} 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 nonlinear 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 postprocessing. 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 crosschecks 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 nonlinear 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 YahiaCherif 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 nearestneighbour, the down nearestneighbour and the rightdown nearestneighbour); 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 crosscorrelations. 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^{−1} (computed by default via onestep row reduction) by comparing it to the Moore–Penrose pseudoinverse and to the inverse obtained via eigendecomposition (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, 0}h^{2}, Ω_{c, 0}h^{2}, h, n_{s}, log10^{10}A_{s}, (w_{0}, 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^{20} (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 σ_{8}, 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 nonlinear 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 (3point and 5point stencils) with a userdefined 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 crosscorrelation 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 semianalytical 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^{(GCph, WL)}: whether the code has been validated to forecast crosscorrelations between galaxy clustering and weak lensing, both from the photometric survey;

curvature: whether the code has been validated to forecast nonflat 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 codecomparison 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 nonzero 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 twodimensional galaxy power spectrum, P_{obs}(k_{ref}, μ_{ref}; z), in Eq. (87). This is achieved using the full nonlinear 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σ_{8}(z) and bσ_{8}(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 σ_{8}(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 σ_{8}(z).

The pure matter power spectrum normalised by , . Although these redshiftindependent 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 “nowiggles” 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} = 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 redshiftindependent shape parameters and five redshiftdependent 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 nonlinear model, we add two parameters to the set and hence two rows and columns to the Fisher matrix. We then project these initial “modelindependent” matrices into the final, modelspecific parameter spaces following Sect. 3.2.5.
A note about derivatives. Derivatives of log_{10}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 threepoint rules:
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 nonlinear parameters, we also use 3point derivatives to calculate the following:
We use different approaches for the derivatives with respect to redshiftdependent 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σ_{8}(z) and ln bσ_{8}(z) we use the 3point rule:
For ln H and ln D_{A} we choose 5point 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^{21}. The derivatives of the galaxy power spectrum with respect to the Hubble parameter and the angular diameter distance are:
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
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 debugging. We now describe each of these settings in detail. A summary of all settings is presented in Table 7.
Summary of the cases considered in our code comparison.
Linear setting. In this case we only consider the observed power spectrum for scales k < k_{max} = 0.25 h Mpc^{−1} for all redshift bins and use the lineartheory prediction for P_{obs}(k, μ). This is achieved by fixing the nonlinear parameters σ_{p} = σ_{v} = 0 in the nonlinear 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 nonflat 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 redshiftdependent 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 nonlinear model for the observed power spectrum of Eq. (87) and consider scales k < k_{max} = 0.25 h Mpc^{−1}. 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 nonlinearities^{22}. 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 nonlinear 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
Here, we have assumed that σ_{p} can be written as σ_{p}(z_{mean})D(z)/D(z_{mean}). Furthermore σ_{8}(z) can also be written in terms of a σ_{8}(z_{mean}), as it was the case for σ_{p}(z) and σ_{v}(z), and Eq. (152) then reads:
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
whereas the damping term remains as
We calculate the derivatives with respect to the two new nonlinear parameters using the same threepoint 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^{−1}. 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^{−1}; however, we keep σ_{p} = σ_{v} fixed to their reference nonzero values, computed using Eq. (81). This means that nonlinear 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 nonflat 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 :
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
Now the Jacobian is given by a (4 + 3 × N_{z})×8 matrix. The function σ_{8}(z) can be expressed as
and consequently fσ_{8}(z) is given by
where f(z) is the numerical solution of the firstorder 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, γ:
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.
PRnflat 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σ_{8}(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σ_{8}(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 nonflat w_{0} + w_{a} + γ cosmology), which is equal for all the four cases listed above^{23}. The Jacobian is explictly
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 redshiftdependent 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 σ_{8}. 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, topleft 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σ_{8}(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%.
Fig. 2.
Comparison of the errors for the intermediate case with k_{max} = 0.30 h Mpc^{−1}. The points indicate deviation of each code from the median of the codes, the light grey area indicates the maximum deviation from the median of the marginalised errors, while the dark grey area indicates the same but for the unmarginalised errors. Topleft panel: four shape parameters. The rest of the panels show redshiftdependent parameters. Specifically, the topright panel is for the angular diameter distance logarithm, ln D_{A}, the bottomleft for the Hubble function, ln H, and the bottomright for the growth of structure, ln fσ_{8}. As explained in the text, the values of these parameters in each redshift bin are treated as separate, independent parameters. However, to indicate that the quantities would be in reality expected to evolve in redshift, we connect the dots with dotted lines in the marginalised case. In the unmarginalised case, we don’t plot lines, since the shaded area connects the points already. The nonlinear parameters σ_{p} and σ_{v}, as well as ln bσ_{8} and P_{s}, have been marginalised over. 
Finally in Fig. 3 we also show the relative difference of the full marginalised projected parameters among six GC codes, for the nonflat ΛCDM + w_{0} + w_{a} + γ case, using the nonlinear recipe up to k_{max} = 0.30 h Mpc^{−1}, corresponding to the optimistic setting. The agreement between the outcomes are at most of the order of a few percent.
Fig. 3.
Comparison for the projected parameters. Comparison of the errors for the fullymarginalised final parameters among six GC codes, for the nonflat ΛCDM + w_{0} + w_{a} + γ case, using the nonlinear recipe up to k_{max} = 0.30 h Mpc^{−1}, corresponding to the optimistic case described in the text. 
4.3. Weaklensing code comparison
We now describe the code comparison and validation of the cosmicshear weaklensing Fishermatrix 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 weaklensing 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 equispaced 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
with λ_{k} = λ_{min} + (k − 1)Δλ, and Δλ = (λ_{max} − λ_{min})/N_{ℓ}. The width of the bin is
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 “semianalytical” 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 semianalytical recipe, namely derivatives of the kernel and of the matter and IA power spectra. We address them separately in the following discussion.
Semianalytical 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 nonflat 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^{24}
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)^{2} 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
With regard to the comoving distance, we give below the logarithmic derivatives of interest. For the general nonflat case, these are
where . Again, when the flatness constraint is enforced, the derivative with respect to Ω_{m, 0} includes an extra term and becomes
while all the others remain the same, provided Ω_{DE, 0} = 1 − Ω_{m, 0}.
Semianalytical 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:
where we recall that the index a ∈ {γγ, Iγ, II}. Let us consider first the lensing kernel (a = γγ). We then obtain
with
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 shearIA kernel (i.e. a = Iγ), it is convenient to first split it as
with
The nonnull derivatives can then be written as
Finally, the derivatives of the IA–IA kernel (a = II) are
where the first row holds for derivatives with respect to parameters Ω_{m, 0}, Ω_{DE, 0}, w_{0}, and w_{a}, and ∂_{μ} is a shorthand notation for ∂/∂p_{μ}.
Semianalytical 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_{0}, 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 kz 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
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_{0}, w_{a}, h, n_{s}, σ_{8}}:
We note that the second term in Eqs. (193) and (194) is only present for Ω_{m, 0}, Ω_{DE, 0}, w_{0}, w_{a}, and γ, since the growth factor D(z) depends on these parameters only.
The derivatives with respect to the IA parameters read
for the δI term, and
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 nonlinear 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 nonflat models, as well as including modelling of IA systematics (eNLA model described in Sect. 3.3.1), and implemention of nonstandard parameters such as w, w_{a}, γ discussed in Eq. (32).
Summary of a selection of the cases investigated for the weaklensing code comparison, ordered according to the number of parameters involved.
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 semianalytical 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 semianalytical 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 cosmicshear 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.
Fig. 4.
Comparison of the errors on the marginalised (light grey) and unmarginalised (dark grey) cosmological parameters. Top panels: we compare seven WL codes for the flat w_{0}, w_{a} cosmology, using the “TakaBird” recipe for the nonlinear power spectrum and a maximum multipole of ℓ_{max} = 1500 (pessimistic cut, top left panel) or ℓ_{max} = 5000 (optimistic cut, top right panel). This is case 3 in Table 8. Bottom panels: we compare three WL codes (the ones available for a nonflat geometry that include γ), for a γ, w_{0}, w_{a} cosmology, using the “TakaBird” nonlinear recipe for the power spectrum and a maximum multipole of ℓ_{max} = 1500 (pessimistic cut, bottom left panel) and a maximum multipole of ℓ_{max} = 5000 (optimistic cut, bottom right panel). This is case 6 in Table 8. 
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 crosscorrelation 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 crosscorrelation 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 ϑ = {Ω_{b}h^{2}, Ω_{c}h^{2}, h, n_{s}, A_{s}, w_{0}, w_{a}} as the primary set, and then projects on the chosen set of derived parameters. For CosmicFish, CosmoSIS, and TotallySAF we use the threepoint 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 semianalytical 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 nonlinear power spectrum the “TakaBird” modified HaloFit (HF) prescription of Takahashi et al. (2012), with additional corrections for massiveneutrinos 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 massiveneutrino and two massless neutrinos) and we include intrinsic alignment in the weaklensing 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 crosscorrelation of the two probes; we further consider optimistic and pessimistic cuts when adding spectroscopic galaxy clustering as follows.
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 nonlinear 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 nonlinear 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%.
Fig. 5.
Percentage difference in the errors on the final parameters among CosmicFish, CosmoSIS, STAFF, and TotallySAF for the GC_{ph}+WL+XC combination (top panels) and the GC_{s}+GC_{ph}+WL+XC combination (bottom panels). We assume here a flat (w_{0}, w_{a}) cosmology. This comparison is done both in the pessimistic and optimistic configurations (left and right panels, respectively); we note, however, that for the top panels, since no spectroscopic galaxy clustering is included, there is no need to perform a cut in z (as the pessimistic setting would require). The light grey band correspond to the marginal errors on the parameters, while the dark grey band indicates the errors when we do not marginalise over the other parameters. 
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_{10}k,log_{10}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.
we note that k depends on H(z) via Eq. (78). This is a semianalytical 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 semianalytical 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^{−1}. Moreover, we found that the semianalytical derivative is not stable when varying the kbinning, 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 nonlinear corrections to P_{obs}(k, μ; z) in Eq. (87).
For a generic function f(x)^{25}, the threepoint stencil derivative of f with respect to x at x_{0} is given by
with ε ≪ 1 an increment, usually an arbitrary choice. The threepoint stencil approximates the derivative up to order ε^{2} errors. To be more precise, another choice is represented by the fivepoint stencil
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 illconditioned, 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 singleprobe results can be stepdependent, but they become stepindependent 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 redshiftdependent 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.
kbinning 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_{10}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^{26}.
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_{0}w_{a}flat 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.
Fig. 6.
Left panel: FoM as a function of the steps for the shape parameters, right panel: FoM as a function of the steps for the ln D_{A}(z) parameter. Both results are for spectroscopic GC_{s} only. 
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 lookup table. We found that the lensing kernel of Eq. (120), , is a wellbehaved 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 nonzero 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_{10}k, z, log_{10}P_{A}(k, z)} over the range 0.0005 ≤ k/Mpc^{−1} ≤ 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^{−1} and P_{A}(k, z) in Mpc^{3} in order to avoid unnecessary complications when taking derivatives with respect to the dimensionless Hubble constant h.
Numerical versus semianalytic 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
According to this relation, one can either: (i) first compute the integral and then the derivatives, for instance through finitedifference 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 “semianalytic” 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 semianalytic 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 cosmicshear 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^{−1}. 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, 0}r(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 xline in parameter space in, say, 15 x_{i} points around the fiducial value x_{0} (including x_{0}). 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 xrange 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 semianalytically. FisherMathica uses this method for the matter power spectrum derivative, while the derivatives are computed with standard three and fivepoint stencil methods. On the other hand, CosmicFish, CosmoSIS, and TotallySAF rely on threepoint 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 fivepoint 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_{0}, w_{a}, σ_{8}} 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 offdiagonal elements of the Fisher matrix – turn out to be quite sensitive to the choice of ε and to whether a three or fivepoint 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)^{27}. 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 illconditionned^{28}. In the case of weak lensing and galaxy clustering the condition number is of the order of 2 × 10^{6}. In such a case, in principle, a precision better than 0.5 × 10^{−7} 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 Gaussiandistributed random percentage of all the elements of the WL Fisher matrix (using the TakaBird nonlinear prescription) according to the approach proposed by YahiaCherif 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^{−5}. The 65% confidence level is reached at 10^{−4}. 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^{−4}.
The above results have implications on the accuracy required for the derivatives entering the Fisher matrix. Standard methods of derivation are the Npoint 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 stepsize below 10^{−5} using CAMB and 10^{−4} 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 nonlinear 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^{−1} 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}, σ_{8}}. 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 σ_{8} 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 Fishermatrix codes for the cosmological parameters of interest (baseline and extensions), whose fiducial parameters have been illustrated in Table 1:
We do this for the following scenarios:

ΛCDM within a spatially flat (and nonflat) geometry;

w_{0}, w_{a} within a spatially flat (and nonflat) geometry;

w_{0}, w_{a}, γ within a spatially flat (and nonflat) geometry.
In the ΛCDM case Ω_{DE, 0} ≡ Ω_{Λ}, the equationofstate parameter and its derivative are fixed to w_{0} = −1 and w_{a} = 0, and the modified gravity parameter γ is fixed to γ = 6/11^{29}. 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 crosscorrelation terms only between the photometric galaxy clustering and the weaklensing 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 crosscorrelation of the two probes; and (b) a stronger (weaker) cut when adding GC_{s}. For the pessimistic case, the WL nonlinear 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 crosscorrelation analysis is limited to the lowest of these multipoles, with . For GC_{s} we reach a value of k_{max} = 0.25 h Mpc^{−1}. The optimistic case extends the regime of GC_{s} to k_{max} = 0.30 h Mpc^{−1}, 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 nonGaussian 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 nonlinear regime up to ℓ_{max} = 5000. We can summarise this as follows :
For GC_{s}, we have further tested two more settings, as they can provide a valuable reference for comparison:
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 nonlinear 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 (nonflat) cosmology. Relative marginalised errors are also detailed in Table 9.
Fig. 7.
Marginalised 1σ errors on cosmological parameters (from the Fisher matrix for an indicative code), relative to their corresponding fiducial value for a flat (upper panels) and nonflat (lower panels) spatial geometry, in ΛCDM. Left (right) panels: pessimistic (optimistic) settings, as described in the text. Each histogram refers to different observational probes. For both a spatially flat and a spatially nonflat cosmology, we show results for GC_{s}, WL, GC_{s}+WL, GC_{s}+WL+GC_{ph}, and GC_{s}+WL+GC_{ph}+XC. 
Marginalised 1σ errors on cosmological parameters (from the Fisher matrix for an indicative code), relative to their corresponding fiducial value for a flat spatial geometry, in ΛCDM.
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 nonflat 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 crosscorrelation 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 nonflat Λ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_{0}, w_{a}
We now move to the case in which the equation of state of a darkenergy 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_{0}, 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 (nonflat) 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, crosscorrelations have a significant impact on the estimation of the errors, and help to further tighten uncertainties, in particular on w_{0}.
Fig. 8.
Marginalised 1σ errors on cosmological parameters, relative to their corresponding fiducial values for a flat (upper panels) and nonflat (lower panels) spatial geometry, in the (w_{0}, w_{a}) cosmology. Left (right) panels: pessimistic (optimistic) settings, as described in the text. The histogram refers to different observational probes. For both a spatially flat and a spatially nonflat cosmology, we show results for GC_{s}, WL, GC_{s}+WL, GC_{s}+WL+GC_{ph}, and GC_{s}+WL+GC_{ph}+XC. For w_{a} we show the absolute error, since a relative error is not possible for a fiducial value of 0. 
For this set of cosmologies, we can also estimate the FoM for w_{0}, 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_{0}, 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.
Figure of merit (FoM) values for the w_{0} and w_{a} parameters in flat and nonflat cosmologies for GC_{s}, GC_{ph}, WL, and different combinations of these observables, including crosscorrelations between photometric galaxy clustering and weak lensing.
The combination of all three probes (without crosscorrelations) reaches a FoM of 122 (376) in the pessimistic (optimistic) settings, for a spatially flat cosmology. Including crosscorrelations 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 crosscorrelations). This demonstrates the importance, for a w_{0}, w_{a} scenario, of including crosscorrelations 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 crosscorrelation. 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 nondiagonal terms are reduced, i.e. correlations among the corresponding parameters are removed; this in turn allows for better constraints on those parameters.
Fig. 9.
Correlation matrix P defined in Eq. (50), obtained from the covariance matrix in a (w_{0}, w_{a}) flat cosmology, for an optimistic setting. Left panel: photometric survey combination, without crosscorrelation. Right panel: crosscorrelation term. The inclusion of crosscorrelation between GC_{ph} and WL significantly reduces the correlation among the darkenergy parameters w_{0}, w_{a} and the standard parameters n_{s}, h, and Ω_{b, 0}. The FoC reduces from 1907 in the case with no crosscorrelations to 824 in the GC_{ph}+WL+ case. This, in turn, leads to reduced uncertainties on all 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 nonGaussian 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 nonGaussian 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 nonflat 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_{0}, 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_{0}, 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.
Fig. 10.
Fishermatrixmarginalised contours for the Euclid space mission for a w_{0}, w_{a} flat cosmology, for GC_{s} (purple), WL (blue), their combination (orange), and with the addition of GC_{ph} and its crosscorrelation with WL (yellow). All combinations correspond to an optimistic setting, as defined in Eq. (206). 
Fig. 11.
Same as in Fig. 10, focusing on the w_{0}, w_{a} plane (left panel) and the Ω_{m, 0}, σ_{8} plane (right panel). 
Furthermore, in Fig. 12 we show the photometric survey alone (blue) including GC_{ph} and WL, and the impact of their crosscorrelations (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.
Fig. 12.
Fishermatrixmarginalised contours for the Euclid space mission for a w_{0}, w_{a} cosmology, for a flat cosmology and with an optimistic setting, of the photometric survey GC_{ph} + WL with (dark red) and without (blue) their crosscorrelation term XC. The combination in yellow is the same as in Fig. 10 and includes all probes. 
A comparison of yellow and red contours shows that the inclusion of crosscorrelations 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_{0}, w_{a}, Ω_{m, 0}, and σ_{8}, 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_{0}, 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_{0}, w_{a}, γ scenario the relative marginal errors for GC_{s}, WL, and their combination, are shown in Fig. 13 for flat and nonflat 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.
Fig. 13.
Marginalised 1σ errors on cosmological parameters, relative to their corresponding fiducial values for a flat (upper panels) and nonflat (lower panels) spatial geometry, in (w_{0}, w_{a}, γ) cosmology. Left (right) panels: pessimistic (optimistic) settings, as described in the text. The histogram refers to different observational probes: GC_{s} (purple), WL (blue), GC_{s}+WL (orange); none of the available codes allow us to include GC_{ph} or crosscorrelation, while also allowing γ to vary. For w_{a} we show the absolute error, since a relative error is not possible for a fiducial value of 0. 
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 nonlinear settings, since this is essential to grasp any information from this probe.
Marginalised 1σ errors on cosmological parameters, relative to their corresponding fiducial value for a flat spatial geometry, in a w_{0}, w_{a}, γ cosmology.
Marginalised 1σ errors on cosmological parameters, relative to their corresponding fiducial value for a flat spatial geometry, in a w_{0}, w_{a}, γ cosmology.
Finally, the combination of all probes, including crosscorrelations, 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 nonflat cases. For the other set of parameters, i.e. w_{0}w_{a} and w_{0}w_{a}γ, for both flat and nonflat cases, the errors decrease by about 25−30%, with the exception of γ for which we gain 15% information only.
Marginalised 1σ errors on cosmological parameters, relative to their corresponding fiducial value for a flat spatial geometry, in a w_{0}, w_{a} cosmology.
Finally, comparing the intermediate with the optimistic setting, i.e. keeping k_{max} = 0.3 h Mpc^{−1} but assuming perfect knowledge on the nonlinear parameters σ_{p} and σ_{v}, we gain information over all the cosmological parameters. For example, considering the projection case w_{0}, w_{a} flat, the errors on the dark energy parameters decrease by about 100% and 150%, for w_{0} and w_{a}, respectively. If curvature is also considered in the analysis, the errors on w_{0} and w_{a} decrease by 28% only.
6. Conclusions
In this paper we have presented the validation of Fishermatrix forecasts for the Euclid satellite, as a result of activity carried out by the interscience 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 crosscorrelation 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 preexisting 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^{−4}, 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 codecomparison 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 nonlinear contribution, and how to take different cuts depending on the emerging presence of nonGaussian 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 nonlinear 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_{0}, w_{a}) scenario in which the equation of state of dark energy is changing with time; the case (w_{0}, 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 nonflat geometry. For all models we present a series of histogram bar plots (Fig. 7 for ΛCDM and Fig. 8 for w_{0}, w_{a}; Fig. 13 for w_{0}, w_{a}, γ) that allow marginalised 1σ errors on cosmological parameters to be visualised for different probe combinations.
In the (w_{0}, 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 nonflat cosmology is allowed.
We have demonstrated that the impact of crosscorrelations in redshift among GC_{ph} and WL is particularly relevant in models beyond ΛCDM: XC improves the FoM for a (w_{0}, 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 simulationderived 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.
We note that throughout, as is common in the community, we use interchangeable arguments of functions where the arguments are also functionally related. For example D(z) or D(a) where a ≡ 1/(1 + z). In all cases it should be clear from context of the equation why the function is expressed in the manner that it is.
Technically the covariance is the inverse of the Hessian matrix of the likelihood, and the Fisher matrix is the expectation of the Hessian; furthermore the inverse of the Fisher matrix does not give the expectation of the covariance, but yields an upper limit to the errors. Therefore we should use another symbol other than C. However for the sake of clarity of expression we use the notation as stated here.
In more than two dimensions, the numerator in Eq. (60) is replaced by the cofactor of element F_{αβ}.
In practice, observers measure the power spectrum in redshift bins, so z here becomes the effective redshift of the bin, which for simplicity we just approximate here to be the bin’s central redshift (see Linder et al. 2014, for a discussion of the effective redshift).
Although other more realistic choices are possible, we do not explore them here. We are, however, confident that the choice of n(z) does not impact the results of the comparison significantly (Kitching et al. 2019), that is to say, if two codes agree for a given n(z), they still agree if a different n(z) is adopted.
Here we use k < k_{max} = 0.25 h Mpc^{−1} as “pessimistic”, but one may argue that it is not pessimistic or even optimistic, depending on ones confidence in modelling of the nonlinear power spectrum. For example the 1halo term is a 30% correction to linear theory at this scale. If the model beyond linear scales is insufficient then this is an issue of potential bias in the results; however, this should still enable us to capture the expected uncertainties in fitted parameters correctly.
Acknowledgments
Authors’ contributions: All authors have significantly contributed to this publication and/or made multiyear essential contributions to the Euclid project which made this work possible. The paperspecific 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 YahiaCherif, Elisabetta Majerotto, Victoria Yankelevich. Weak Lensing Task: Vincenzo F. Cardone, Santiago Casas, Stefano Camera, Matteo Martinelli, Isaac Tutusaus, Thomas Kitching, Martin Kilbinger. Crosscorrelations 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.euclidec.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 Multiwavelength 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 YahiaCherif acknowledge support from french space agency CNES. Sebastien Clesse aknowledges support from the Belgian Fund for Research F.R.SFNRS. Martin Kunz and Fabien Lacasa acknowledge financial support from the Swiss National Science Foundation. Eric Linder acknowledges support by NASA ROSES grant 12EUCLID120004. 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 DITP 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 ESP201789838C31R, 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 BonnCologne 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
 Abazajian, K. N., Adshead, P., Ahmed, Z., et al. 2016, ArXiv eprints [arXiv:1610.02743] [Google Scholar]
 Abbott, T. M. C., Abdalla, F. B., Alarcon, A., et al. 2018a, Phys. Rev. D, 98, 043526 [Google Scholar]
 Abbott, T. M. C., Abdalla, F. B., Annis, J., et al. 2018b, MNRAS, 480, 3879 [NASA ADS] [CrossRef] [Google Scholar]
 Albrecht, A., Bernstein, G., Cahn, R., et al. 2006, ArXiv eprints [arXiv:astroph/0609591] [Google Scholar]
 Alsing, J., Kirk, D., Heavens, A., & Jaffe, A. H. 2015, MNRAS, 452, 1202 [CrossRef] [Google Scholar]
 Amendola, L., Kunz, M., & Sapone, D. 2008, JCAP, 0804, 013 [Google Scholar]
 Amendola, L., Pettorino, V., Quercellini, C., & Vollmer, A. 2012, Phys. Rev. D, 85, 103008 [NASA ADS] [CrossRef] [Google Scholar]
 Amendola, L., Appleby, S., Avgoustidis, A., et al. 2018, Living Rev. Relativ., 21, 2 [Google Scholar]
 Bacon, D. J., et al. 2020, PASA, 37, e007 [CrossRef] [Google Scholar]
 Baldauf, T., Mirbabayi, M., Simonović, M., & Zaldarriaga, M. 2015, Phys. Rev. D, 92, 043514 [NASA ADS] [CrossRef] [Google Scholar]
 Ballinger, W. E., Peacock, J. A., & Heavens, A. F. 1996, MNRAS, 282, 877 [NASA ADS] [CrossRef] [Google Scholar]
 Baxter, E. J., Omori, Y., Chang, C., et al. 2019, Phys. Rev. D, 99, 023508 [NASA ADS] [CrossRef] [Google Scholar]
 Bernardeau, F., van de Rijt, N., & Vernizzi, F. 2012, Phys. Rev. D, 85, 063509 [CrossRef] [Google Scholar]
 Beutler, F., Seo, H.J., Saito, S., et al. 2017, MNRAS, 466, 2242 [NASA ADS] [CrossRef] [Google Scholar]
 Bianchi, D., Guzzo, L., Branchini, E., et al. 2012, MNRAS, 427, 2420 [NASA ADS] [CrossRef] [Google Scholar]
 Bird, S., Viel, M., & Haehnelt, M. G. 2012, MNRAS, 420, 2551 [NASA ADS] [CrossRef] [Google Scholar]
 Blas, D., Lesgourgues, J., & Tram, T. 2011, JCAP, 7, 034 [NASA ADS] [CrossRef] [Google Scholar]
 Bonaldi, A., Harrison, I., Camera, S., & Brown, M. L. 2016, MNRAS, 463, 3686 [CrossRef] [Google Scholar]
 Bridle, S., & King, L. 2007, New J. Phys., 9, 444 [NASA ADS] [CrossRef] [Google Scholar]
 Bueno Belloso, A., GarciaBellido, J., & Sapone, D. 2011, JCAP, 1110, 010 [CrossRef] [Google Scholar]
 Bunn, E. F. 1995, Ph.D. Thesis, University of California, Berkeley, USA [Google Scholar]
 Camera, S., & Nishizawa, A. 2013, Phys. Rev. Lett., 110, 151103 [CrossRef] [Google Scholar]
 Camera, S., Santos, M. G., Bacon, D. J., et al. 2012, MNRAS, 427, 2079 [CrossRef] [Google Scholar]
 Camera, S., Fornasa, M., Fornengo, N., & Regis, M. 2013, ApJ, 771, L5 [CrossRef] [Google Scholar]
 Camera, S., Fornasa, M., Fornengo, N., & Regis, M. 2015a, JCAP, 1506, 029 [CrossRef] [Google Scholar]
 Camera, S., Maartens, R., & Santos, M. G. 2015b, MNRAS, 451, L80 [CrossRef] [Google Scholar]
 Camera, S., Santos, M. G., & Maartens, R. 2015c, MNRAS, 448, 1035 [CrossRef] [Google Scholar]
 Camera, S., Harrison, I., Bonaldi, A., & Brown, M. L. 2017, MNRAS, 464, 4747 [CrossRef] [Google Scholar]
 Camera, S., Fonseca, J., Maartens, R., & Santos, M. G. 2018, MNRAS, 481, 1251 [Google Scholar]
 Carbone, C., Mangilli, A., & Verde, L. 2011a, JCAP, 9, 028 [CrossRef] [Google Scholar]
 Carbone, C., Verde, L., Wang, Y., & Cimatti, A. 2011b, JCAP, 3, 030 [CrossRef] [Google Scholar]
 Carbone, C., Fedeli, C., Moscardini, L., & Cimatti, A. 2012, JCAP, 3, 023 [NASA ADS] [CrossRef] [Google Scholar]
 Cardona, W., Durrer, R., Kunz, M., & Montanari, F. 2016, Phys. Rev. D, 94, 043007 [NASA ADS] [CrossRef] [Google Scholar]
 Carron, J. 2013, A&A, 551, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Casas, S., Amendola, L., Baldi, M., Pettorino, V., & Vollmer, A. 2016, JCAP, 1, 045 [CrossRef] [Google Scholar]
 Casas, S., Kunz, M., Martinelli, M., & Pettorino, V. 2017, Phys. Dark Univ., 18, 73 [CrossRef] [Google Scholar]
 Casas, S., Pauly, M., & Rubio, J. 2018, Phys. Rev. D, 97, 043520 [CrossRef] [Google Scholar]
 Castro, P. G., Heavens, A. F., & Kitching, T. D. 2005, Phys. Rev. D, 72, 023516 [NASA ADS] [CrossRef] [Google Scholar]
 Chisari, N., Codis, S., Laigle, C., et al. 2015, MNRAS, 454, 2736 [NASA ADS] [CrossRef] [Google Scholar]
 Chisari, N. E., Alonso, D., Krause, E., et al. 2019, ApJS, 242, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Coe, D. 2009, ArXiv eprints [arXiv:0906.4123] [Google Scholar]
 Cooray, A., & Hu, W. 2001, ApJ, 554, 56 [NASA ADS] [CrossRef] [Google Scholar]
 Cooray, A., & Sheth, R. 2002, Phys. Rep., 372, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Copeland, D., Taylor, A., & Hall, A. 2018, MNRAS, 480, 2247 [CrossRef] [Google Scholar]
 De Angelis, A., Tatischeff, V., Grenier, I. A., et al. 2018, J. High Energy Astrophys., 19, 1 [NASA ADS] [CrossRef] [Google Scholar]
 de Putter, R., Doré, O., & Takada, M. 2013, ArXiv eprints [arXiv:1308.6070] [Google Scholar]
 Desjacques, V., Jeong, D., & Schmidt, F. 2018, Phys. Rep., 733, 1 [NASA ADS] [CrossRef] [Google Scholar]
 di Porto, C., Amendola, L., & Branchini, E. 2012, MNRAS, 419, 985 [NASA ADS] [CrossRef] [Google Scholar]
 Duncan, C. A. J., Heymans, C., Heavens, A. F., & Joachimi, B. 2016, MNRAS, 457, 764 [CrossRef] [Google Scholar]
 Eisenstein, D. J., & Hu, W. 1997, ApJ, 511, 5 [Google Scholar]
 Eisenstein, D. J., & Hu, W. 1998, ApJ, 496, 605 [NASA ADS] [CrossRef] [Google Scholar]
 Eisenstein, D. J., Seo, H.J., & White, M. 2007, ApJ, 664, 660 [NASA ADS] [CrossRef] [Google Scholar]
 Fang, W., Hu, W., & Lewis, A. 2008, Phys. Rev. D, 78, 087303 [NASA ADS] [CrossRef] [Google Scholar]
 Fonseca, J., Camera, S., Santos, M., & Maartens, R. 2015, ApJ, 812, L22 [CrossRef] [Google Scholar]
 Gatti, M., Vielzeuf, P., Davis, C., et al. 2018, MNRAS, 477, 1664 [CrossRef] [Google Scholar]
 Ghosh, B., Durrer, R., & Sellentin, E. 2018, JCAP, 2018, 008 [CrossRef] [Google Scholar]
 Giannantonio, T., Porciani, C., Carron, J., Amara, A., & Pillepich, A. 2012, MNRAS, 422, 2854 [Google Scholar]
 GilMarín, H., Percival, W. J., Brownstein, J. R., et al. 2016, MNRAS, 460, 4188 [Google Scholar]
 Grieb, J. N., Sánchez, A. G., SalazarAlbornoz, S., et al. 2017, MNRAS, 467, 2085 [NASA ADS] [Google Scholar]
 Hamann, J., Hannestad, S., & Wong, Y. Y. Y. 2012, JCAP, 11, 052 [Google Scholar]
 Hamilton, A. J. S. 1998, in The Evolving Universe, ed. D. Hamilton, Astrophys. Space Sci. Lib., 231, 185 [Google Scholar]
 Hamilton, A. J. S., Rimes, C. D., & Scoccimarro, R. 2006, MNRAS, 371, 1188 [Google Scholar]
 Harrison, I., Camera, S., Zuntz, J., & Brown, M. L. 2016, MNRAS, 463, 3674 [Google Scholar]
 Haynsworth, E. V. 1968, BMN 20 [Google Scholar]
 Heitmann, K., Lawrence, E., Kwan, J., Habib, S., & Higdon, D. 2014, ApJ, 780, 111 [Google Scholar]
 Hilbert, S., Xu, D., Schneider, P., et al. 2017, MNRAS, 468, 790 [Google Scholar]
 Ho, S., Cuesta, A., Seo, H.J., et al. 2012, ApJ, 761, 14 [Google Scholar]
 Hojjati, A., Pogosian, L., & Zhao, G.B. 2011, JCAP, 1108, 005 [Google Scholar]
 Hollenstein, L., Sapone, D., Crittenden, R., & Schaefer, B. M. 2009, JCAP, 0904, 012 [Google Scholar]
 Howlett, C., Lewis, A., Hall, A., & Challinor, A. 2012, JCAP, 1204, 027 [Google Scholar]
 Hu, W. 1999, ApJ, 522, L21 [Google Scholar]
 Hu, W., & Kravtsov, A. V. 2003, ApJ, 584, 702 [Google Scholar]
 Hu, W., & Sawicki, I. 2007, Phys. Rev. D, 76, 104043 [Google Scholar]
 Hu, B., Raveri, M., Frusciante, N., & Silvestri, A. 2014, Phys. Rev. D, 89, 103530 [Google Scholar]
 Jain, B., Spergel, D., Bean, R., et al. 2015, ArXiv eprints [arXiv:1501.07897] [Google Scholar]
 JelicCizmek, G., Lepori, F., Bonvin, C., & Durrer, R. 2020, ArXiv eprints [arXiv:2004.12981] [Google Scholar]
 Jimenez, R., Kitching, T., PeñaGaray, C., & Verde, L. 2010, JCAP, 2010, 035 [Google Scholar]
 Joachimi, B., Cacciato, M., Kitching, T. D., et al. 2015, SSR, 193, 1 [Google Scholar]
 Kaiser, N. 1987, MNRAS, 227, 1 [Google Scholar]
 Kaiser, N. 1992, ApJ, 388, 272 [Google Scholar]
 Kayo, I., Takada, M., & Jain, B. 2013, MNRAS, 429, 344 [Google Scholar]
 Kiakotou, A., Elgarøy, Ø., & Lahav, O. 2008, Phys. Rev. D, 77, 063005 [Google Scholar]
 Kiessling, A., Cacciato, M., Joachimi, B., et al. 2015, SSR, 193, 67 [Google Scholar]
 Kilbinger, M. 2015, Rep. Prog. Phys., 78, 086901 [Google Scholar]
 Kilbinger, M., Heymans, C., Asgari, M., et al. 2017, MNRAS, 472, 2126 [Google Scholar]
 Kirk, D., Brown, M. L., Hoekstra, H., et al. 2015, SSR, 193, 139 [Google Scholar]
 Kitching, T. D., & Amara, A. 2009, MNRAS, 398, 2134 [Google Scholar]
 Kitching, T. D., Amara, A., Abdalla, F. B., Joachimi, B., & Refregier, A. 2009, MNRAS, 399, 2107 [Google Scholar]
 Kitching, T., Bacon, D., Brown, M., et al. 2015, PoS, AASKA14, 146 [Google Scholar]
 Kitching, T. D., Alsing, J., Heavens, A. F., et al. 2017, MNRAS, 469, 2737 [Google Scholar]
 Kitching, T. D., Taylor, P. L., Capak, P., Masters, D., & Hoekstra, H. 2019, Phys. Rev. D, 99, 063536 [Google Scholar]
 Lacasa, F. 2018, A&A, 615, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lacasa, F., & Rosenfeld, R. 2016, JCAP, 8, 005 [Google Scholar]
 Lahav, O., Lilje, P. B., Primack, J. R., & Rees, M. J. 1991, MNRAS, 251, 128 [Google Scholar]
 Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, ArXiv eprints [arXiv:1110.3193] [Google Scholar]
 Lemos, P., Challinor, A., & Efstathiou, G. 2017, JCAP, 5, 014 [Google Scholar]
 Lesgourgues, J. 2011a, ArXiv eprints [arXiv:1104.2932] [Google Scholar]
 Lesgourgues, J. 2011b, ArXiv eprints [arXiv:1104.2934] [Google Scholar]
 Lesgourgues, J., & Tram, T. 2011, JCAP, 9, 032 [Google Scholar]
 Lewis, A., Challinor, A., & Lasenby, A. 2000, ApJ, 538, 473 [Google Scholar]
 Linder, E. V. 2005, Phys. Rev. D, 72, 043529 [Google Scholar]
 Linder, E. V., & Jenkins, A. 2003, MNRAS, 346, 573 [Google Scholar]
 Linder, E. V., Oh, M., Okumura, T., Sabiu, C. G., & Song, Y.S. 2014, Phys. Rev. D, 89, 063525 [Google Scholar]
 Majerotto, E., Guzzo, L., Samushia, L., et al. 2012, MNRAS, 424, 1392 [Google Scholar]
 Majerotto, E., Sapone, D., & Schäfer, B. M. 2016, MNRAS, 456, 109 [Google Scholar]
 Markovič, K., Percival, W. J., Scodeggio, M., et al. 2017, MNRAS, 467, 3677 [Google Scholar]
 Massey, R., Rhodes, J., Refregier, A., et al. 2004, AJ, 127, 3089 [Google Scholar]
 Mead, A. J., Heymans, C., Lombriser, L., et al. 2016, MNRAS, 459, 1468 [Google Scholar]
 Montanari, F., & Durrer, R. 2015, JCAP, 2015, 070 [Google Scholar]
 Mueller, E.M., Percival, W., Linder, E., et al. 2018, MNRAS, 475, 2122 [Google Scholar]
 Olivari, L. C., Dickinson, C., Battye, R. A., et al. 2018, MNRAS, 473, 4242 [Google Scholar]
 PalanqueDelabrouille, N., Yèche, C., Schöneberg, N., et al. 2020, JCAP, 04, 038 [Google Scholar]
 Park, Y., et al. 2016, Phys. Rev. D, 94, 063533 [Google Scholar]
 Peebles, P. J. E. 1980, The LargeScale Structure of the Universe (Princeton: Princeton University Press) [Google Scholar]
 Peirone, S., Martinelli, M., Raveri, M., & Silvestri, A. 2017, Phys. Rev. D, 96, 063524 [Google Scholar]
 Planck Collaboration XIV. 2016, A&A, 594, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VI. 2020, A&A, 641, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pogosian, L., & Silvestri, A. 2016, Phys. Rev. D, 94, 104014 [Google Scholar]
 Pourtsidou, A., Bacon, D., & Crittenden, R. 2017, MNRAS, 470, 4251 [Google Scholar]
 Pozzetti, L., Hirata, C. M., Geach, J. E., et al. 2016, A&A, 590, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2007, Numerical Recipes: The Art of Scientific Computing, 3rd edn. (New York, NY, USA: Cambridge University Press) [Google Scholar]
 Racca, G. D., Laureijs, R., Stagnaro, L., et al. 2016, in Proc. SPIE, SPIE Conf. Ser., 9904, 99040O [Google Scholar]
 Rassat, A., Amara, A., Amendola, L., et al. 2008, ArXiv eprints [arXiv:0810.0003] [Google Scholar]
 Ratra, B., & Peebles, P. J. E. 1988, Phys. Rev. D, 37, 3406 [Google Scholar]
 Raveri, M., Hu, B., Frusciante, N., & Silvestri, A. 2014, Phys. Rev. D, 90, 043513 [Google Scholar]
 Raveri, M., Martinelli, M., Zhao, G., & Wang, Y. 2016a, ArXiv eprints [arXiv:1606.06273] [Google Scholar]
 Raveri, M., Martinelli, M., Zhao, G., & Wang, Y. 2016b, ArXiv eprints [arXiv:1607.01005] [Google Scholar]
 Rhodes, J., Nichol, R. C., Aubourg, É., et al. 2017, ApJS, 233, 21 [Google Scholar]
 Rimes, C. D., & Hamilton, A. J. S. 2005, MNRAS, 360, L82 [NASA ADS] [Google Scholar]
 Sapone, D., & Majerotto, E. 2012, Phys. Rev. D, 85, 123529 [Google Scholar]
 Sapone, D., Kunz, M., & Kunz, M. 2009, Phys. Rev. D, 80, 083519 [Google Scholar]
 Sapone, D., Kunz, M., & Amendola, L. 2010, Phys. Rev. D, 82, 103535 [Google Scholar]
 Sapone, D., Majerotto, E., Kunz, M., & Garilli, B. 2013, Phys. Rev. D, 88, 043503 [Google Scholar]
 Sapone, D., Majerotto, E., & Nesseris, S. 2014, Phys. Rev. D, 90, 023012 [Google Scholar]
 Sato, M., Hamana, T., Takahashi, R., et al. 2009, ApJ, 701, 945 [Google Scholar]
 Scoccimarro, R., & Frieman, J. 1996, ApJS, 105, 37 [Google Scholar]
 Semboloni, E., Hoekstra, H., Schaye, J., van Daalen, M. P., & McCarthy, I. G. 2011, MNRAS, 417, 2020 [Google Scholar]
 Seo, H.J., & Eisenstein, D. J. 2003, ApJ, 598, 720 [Google Scholar]
 Seo, H.J., & Eisenstein, D. J. 2007, ApJ, 665, 14 [Google Scholar]
 Smail, I., Ellis, R. S., & Fitchett, M. J. 1994, MNRAS, 270, 245 [Google Scholar]
 Smith, R. E., Peacock, J. A., Jenkins, A., et al. 2003, MNRAS, 341, 1311 [Google Scholar]
 Spurio Mancini, A., Reischke, R., Pettorino, V., Schäfer, B. M., & Zumalacárregui, M. 2018, MNRAS, 480, 3725 [Google Scholar]
 Takada, M. 2006, Phys. Rev. D, 74, 043505 [Google Scholar]
 Takada, M., & Hu, W. 2013, Phys. Rev. D, 87, 123504 [NASA ADS] [CrossRef] [Google Scholar]
 Takada, M., & Jain, B. 2009, MNRAS, 395, 2065 [NASA ADS] [CrossRef] [Google Scholar]
 Takada, M., Komatsu, E., & Futamase, T. 2006, Phys. Rev. D, 73, 083520 [CrossRef] [Google Scholar]
 Takahashi, R., Sato, M., Nishimichi, T., Taruya, A., & Oguri, M. 2012, ApJ, 761, 152 [NASA ADS] [CrossRef] [Google Scholar]
 Tanidis, K., Camera, S., & Parkinson, D. 2020, MNRAS, 491, 4869 [CrossRef] [Google Scholar]
 Taylor, P. L., Kitching, T. D., & McEwen, J. D. 2018a, Phys. Rev. D, 98, 043532 [NASA ADS] [CrossRef] [Google Scholar]
 Taylor, P. L., Kitching, T. D., McEwen, J. D., & Tram, T. 2018b, Phys. Rev. D, 98, 023522 [NASA ADS] [CrossRef] [Google Scholar]
 Taylor, P. L., Bernardeau, F., & Kitching, T. D. 2018c, Phys. Rev. D, 98, 083514 [NASA ADS] [CrossRef] [Google Scholar]
 Tegmark, M. 1997, Phys. Rev. Lett., 79, 3806 [NASA ADS] [CrossRef] [Google Scholar]
 Tegmark, M., Taylor, A., & Heavens, A. 1997, ApJ, 480, 22 [NASA ADS] [CrossRef] [Google Scholar]
 Tenneti, A., Mandelbaum, R., & Di Matteo, T. 2016, MNRAS, 462, 2668 [NASA ADS] [CrossRef] [Google Scholar]
 Thiele, L., Duncan, C. A. J., & Alonso, D. 2020, MNRAS, 491, 1746 [NASA ADS] [CrossRef] [Google Scholar]
 Tutusaus, I., Lamine, B., Blanchard, A., et al. 2016, Phys. Rev. D, 94, 123515 [NASA ADS] [CrossRef] [Google Scholar]
 Tutusaus, I., Martinelli, M., Cardone, V. F., et al. 2020, A&A, in press, https://doi.org/10.1051/00046361/202038313 [Google Scholar]
 Vogeley, M. S., & Szalay, A. S. 1996, ApJ, 465, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, Y. 2008, Phys. Rev. D, 77, 123525 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, Y. 2010, Mod. Phys. Lett. A, 25, 3093 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, Y. 2012, MNRAS, 423, 3631 [CrossRef] [Google Scholar]
 Wang, Y., Percival, W., Cimatti, A., et al. 2010, MNRAS, 409, 737 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, Y., Chuang, C.H., & Hirata, C. M. 2013, MNRAS, 430, 2446 [CrossRef] [Google Scholar]
 Weltman, A., Bull, P., Camera, S., et al. 2020, PASA, 37, e002 [CrossRef] [Google Scholar]
 Wetterich, C. 1988, Nucl. Phys. B, 302, 668 [NASA ADS] [CrossRef] [Google Scholar]
 YahiaCherif, S., Blanchard, A., Camera, S., et al. 2020, A&A, in press, https://doi.org/10.1051/00046361/201937312 [Google Scholar]
 Yankelevich, V., & Porciani, C. 2019, MNRAS, 483, 2078 [CrossRef] [Google Scholar]
 Zhang, F. 2005, in The Schur Complement and its Applications (New York: Springer), Numerical Methods and Algorithms, 4 [CrossRef] [Google Scholar]
 Zhao, G.B., Pogosian, L., Silvestri, A., & Zylberberg, J. 2009, Phys. Rev. D, 79, 083513 [NASA ADS] [CrossRef] [Google Scholar]
 Zuntz, J., Paterno, M., Jennings, E., et al. 2015, Astron. Comput., 12, 45 [NASA ADS] [CrossRef] [Google Scholar]
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^{30}. 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 nonflat 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_{0}, 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 userprovided 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^{31}. 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 Pknw 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_1p0E2, 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 scaledmeanlumE2Sa.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 nonflat case;
Each of these two zip files contains:

the file pkzFiducial.txt, with the fiducial power spectrum;

the spectra relative to the positive (negative) steps in cosmological parameters, e.g. ns_pl_eps_5p0E2 indicates a step on n_{s} in the positive direction of amplitude 5.0 × 10^{−2} (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_{0}, 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^{32}, 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 crosscorrelation 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_{0}, 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_{0}, 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: w0waCDMflatoptimistic_bounds.txt containing the relative errors on all the parameters of interest and a file w0waCDMflatoptimistic_FoMs.txt containing the FoM for w_{0}, 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
Parameter values of our fiducial cosmological model, corresponding to those of Planck Collaboration XIV (2016), both in the baseline ΛCDM case and in the extensions considered.
Expected number density of observed Hα emitters for the Euclid spectroscopic survey.
Summary of a selection of the cases investigated for the weaklensing code comparison, ordered according to the number of parameters involved.
Marginalised 1σ errors on cosmological parameters (from the Fisher matrix for an indicative code), relative to their corresponding fiducial value for a flat spatial geometry, in ΛCDM.
Figure of merit (FoM) values for the w_{0} and w_{a} parameters in flat and nonflat cosmologies for GC_{s}, GC_{ph}, WL, and different combinations of these observables, including crosscorrelations between photometric galaxy clustering and weak lensing.
Marginalised 1σ errors on cosmological parameters, relative to their corresponding fiducial value for a flat spatial geometry, in a w_{0}, w_{a}, γ cosmology.
Marginalised 1σ errors on cosmological parameters, relative to their corresponding fiducial value for a flat spatial geometry, in a w_{0}, w_{a}, γ cosmology.
Marginalised 1σ errors on cosmological parameters, relative to their corresponding fiducial value for a flat spatial geometry, in a w_{0}, w_{a} cosmology.
All Figures
Fig. 1.
Function μ_{ν}(k). It represents the scale dependent correction to f(z) in Eqs. (39) and (40), evaluated here at ∑_{i}m_{ν, i} = 0.06 and 0.15 eV. 

In the text 
Fig. 2.
Comparison of the errors for the intermediate case with k_{max} = 0.30 h Mpc^{−1}. The points indicate deviation of each code from the median of the codes, the light grey area indicates the maximum deviation from the median of the marginalised errors, while the dark grey area indicates the same but for the unmarginalised errors. Topleft panel: four shape parameters. The rest of the panels show redshiftdependent parameters. Specifically, the topright panel is for the angular diameter distance logarithm, ln D_{A}, the bottomleft for the Hubble function, ln H, and the bottomright for the growth of structure, ln fσ_{8}. As explained in the text, the values of these parameters in each redshift bin are treated as separate, independent parameters. However, to indicate that the quantities would be in reality expected to evolve in redshift, we connect the dots with dotted lines in the marginalised case. In the unmarginalised case, we don’t plot lines, since the shaded area connects the points already. The nonlinear parameters σ_{p} and σ_{v}, as well as ln bσ_{8} and P_{s}, have been marginalised over. 

In the text 
Fig. 3.
Comparison for the projected parameters. Comparison of the errors for the fullymarginalised final parameters among six GC codes, for the nonflat ΛCDM + w_{0} + w_{a} + γ case, using the nonlinear recipe up to k_{max} = 0.30 h Mpc^{−1}, corresponding to the optimistic case described in the text. 

In the text 
Fig. 4.
Comparison of the errors on the marginalised (light grey) and unmarginalised (dark grey) cosmological parameters. Top panels: we compare seven WL codes for the flat w_{0}, w_{a} cosmology, using the “TakaBird” recipe for the nonlinear power spectrum and a maximum multipole of ℓ_{max} = 1500 (pessimistic cut, top left panel) or ℓ_{max} = 5000 (optimistic cut, top right panel). This is case 3 in Table 8. Bottom panels: we compare three WL codes (the ones available for a nonflat geometry that include γ), for a γ, w_{0}, w_{a} cosmology, using the “TakaBird” nonlinear recipe for the power spectrum and a maximum multipole of ℓ_{max} = 1500 (pessimistic cut, bottom left panel) and a maximum multipole of ℓ_{max} = 5000 (optimistic cut, bottom right panel). This is case 6 in Table 8. 

In the text 
Fig. 5.
Percentage difference in the errors on the final parameters among CosmicFish, CosmoSIS, STAFF, and TotallySAF for the GC_{ph}+WL+XC combination (top panels) and the GC_{s}+GC_{ph}+WL+XC combination (bottom panels). We assume here a flat (w_{0}, w_{a}) cosmology. This comparison is done both in the pessimistic and optimistic configurations (left and right panels, respectively); we note, however, that for the top panels, since no spectroscopic galaxy clustering is included, there is no need to perform a cut in z (as the pessimistic setting would require). The light grey band correspond to the marginal errors on the parameters, while the dark grey band indicates the errors when we do not marginalise over the other parameters. 

In the text 
Fig. 6.
Left panel: FoM as a function of the steps for the shape parameters, right panel: FoM as a function of the steps for the ln D_{A}(z) parameter. Both results are for spectroscopic GC_{s} only. 

In the text 
Fig. 7.
Marginalised 1σ errors on cosmological parameters (from the Fisher matrix for an indicative code), relative to their corresponding fiducial value for a flat (upper panels) and nonflat (lower panels) spatial geometry, in ΛCDM. Left (right) panels: pessimistic (optimistic) settings, as described in the text. Each histogram refers to different observational probes. For both a spatially flat and a spatially nonflat cosmology, we show results for GC_{s}, WL, GC_{s}+WL, GC_{s}+WL+GC_{ph}, and GC_{s}+WL+GC_{ph}+XC. 

In the text 
Fig. 8.
Marginalised 1σ errors on cosmological parameters, relative to their corresponding fiducial values for a flat (upper panels) and nonflat (lower panels) spatial geometry, in the (w_{0}, w_{a}) cosmology. Left (right) panels: pessimistic (optimistic) settings, as described in the text. The histogram refers to different observational probes. For both a spatially flat and a spatially nonflat cosmology, we show results for GC_{s}, WL, GC_{s}+WL, GC_{s}+WL+GC_{ph}, and GC_{s}+WL+GC_{ph}+XC. For w_{a} we show the absolute error, since a relative error is not possible for a fiducial value of 0. 

In the text 
Fig. 9.
Correlation matrix P defined in Eq. (50), obtained from the covariance matrix in a (w_{0}, w_{a}) flat cosmology, for an optimistic setting. Left panel: photometric survey combination, without crosscorrelation. Right panel: crosscorrelation term. The inclusion of crosscorrelation between GC_{ph} and WL significantly reduces the correlation among the darkenergy parameters w_{0}, w_{a} and the standard parameters n_{s}, h, and Ω_{b, 0}. The FoC reduces from 1907 in the case with no crosscorrelations to 824 in the GC_{ph}+WL+ case. This, in turn, leads to reduced uncertainties on all parameters. 

In the text 
Fig. 10.
Fishermatrixmarginalised contours for the Euclid space mission for a w_{0}, w_{a} flat cosmology, for GC_{s} (purple), WL (blue), their combination (orange), and with the addition of GC_{ph} and its crosscorrelation with WL (yellow). All combinations correspond to an optimistic setting, as defined in Eq. (206). 

In the text 
Fig. 11.
Same as in Fig. 10, focusing on the w_{0}, w_{a} plane (left panel) and the Ω_{m, 0}, σ_{8} plane (right panel). 

In the text 
Fig. 12.
Fishermatrixmarginalised contours for the Euclid space mission for a w_{0}, w_{a} cosmology, for a flat cosmology and with an optimistic setting, of the photometric survey GC_{ph} + WL with (dark red) and without (blue) their crosscorrelation term XC. The combination in yellow is the same as in Fig. 10 and includes all probes. 

In the text 
Fig. 13.
Marginalised 1σ errors on cosmological parameters, relative to their corresponding fiducial values for a flat (upper panels) and nonflat (lower panels) spatial geometry, in (w_{0}, w_{a}, γ) cosmology. Left (right) panels: pessimistic (optimistic) settings, as described in the text. The histogram refers to different observational probes: GC_{s} (purple), WL (blue), GC_{s}+WL (orange); none of the available codes allow us to include GC_{ph} or crosscorrelation, while also allowing γ to vary. For w_{a} we show the absolute error, since a relative error is not possible for a fiducial value of 0. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.