© The Authors 2026

1. Introduction

The weak gravitational lensing of galaxies by the large-scale structure in the Universe, known as cosmic shear, has emerged as a key tool for cosmology given its ability to probe the evolution of matter through cosmic time. This property has allowed for weak lensing analyses (see e.g. Kilbinger 2015 for a review) to be used to constrain the nature of dark energy, proposed as the origin of the observed accelerated expansion of the Universe (Riess et al. 1998, Schmidt et al. 1998, Perlmutter et al. 1999).

Surveys such as Euclid (Laureijs et al. 2011, Euclid Collaboration: Mellier et al. 2025), the forthcoming Nancy Grace Roman Space Telescope (Spergel et al. 2015), and the Vera C. Rubin Observatory Legacy Survey of Space and Time (LSST Science Collaboration 2009), will provide samples of over 1 billion galaxies and are poised to offer an improvement (by at least an order of magnitude) in the precision of cosmic shear measurements over current experiments. These data will be used to investigate the extent to which dark energy can be described as a time-evolving fluid contributing to the mass-energy budget of the Universe. In particular, Euclid will target constraints on a w₀w_aCDM cosmology, where the dark energy equation of state parameter, w, is assumed to take the functional form w = w₀ + (1 − a)w_a (Laureijs et al. 2011), where a is the scale factor.

An intrinsic step with respect to the weak lensing method is the analysis of the 3D radial distribution of matter in the Universe, traced by the galaxies observed in the survey. While a cosmic shear study can employ the ‘3D weak lensing’ analysis technique (Heavens 2003, Kitching et al. 2014, Spurio Mancini et al. 2018) where the shear field is expanded radially in spherical Bessel functions, the common method (i.e., tomography) involves subdividing the observed galaxy population into binned samples at different distances. This technique effectively compresses the 3D cosmological information into a number of slices, through which a 2D two-point analysis of the projected galaxy images is used to recover the radial shear field (Hu 1999).

Comparatively, 3D weak lensing measurements (Castro et al. 2005, Kitching et al. 2007) naturally offer more statistical power for cosmological parameter constraints. However, the information gain is dependent on the chosen parameter(s) and it has been demonstrated (e.g. Hu 2002; Jain & Taylor 2003; Huterer 2002; Schrabback et al. 2010; Benjamin et al. 2013; Taylor et al. 2018b) that the tomographic weak lensing approach recovers the majority of the cosmological information required to make measurements of dark energy with the benefit of a significant reduction in computational cost. Indeed, the current generation of surveys such as the Dark Energy Survey (DES; Abbott et al. 2022), the Kilo Degree Survey (KiDS; Hildebrandt et al. 2017), and the Subaru Hyper Suprime-Cam (HSC; Hikage et al. 2019) have consistently utilised weak lensing tomography to obtain constraints on parameters such as the matter density, Ω_m, the amplitude of matter fluctuations on scales of 8 h⁻¹ Mpc, σ₈, and the dark energy equation of state parameter, w. These analyses provide complementary measurements to constraints derived from alternative probes such as Type Ia supernovae (e.g. Brout et al. 2022), the cosmic microwave background (CMB, e.g. Planck Collaboration VI 2020), and baryon acoustic oscillations (BAO, e.g. Alam et al. 2017).

Notably, the DES collaboration (Abbott et al. 2023) have constrained the dark energy (w₀, w_a) parameters using a tomographic analysis of the 3 × 2 pt signal (the joint data vector of cosmic shear, angular galaxy clustering, and their cross-correlation, see Sect. 3.1). When combining the 3 × 2 pt signal with Planck temperature and E-mode polarisation observables (Planck Collaboration VI 2020), BAO results from eBOSS (Ahumada et al. 2020) and the Pantheon SNe Ia sample (Scolnic et al. 2018), they found $(w_0,w_a)=(-0.95\pm 0.08,-0.4_{-0.3}^{+0.4}),$ which is consistent with a standard Lambda cold dark matter (ΛCDM) model at 1σ.

Comparatively, studies with the Pantheon+ analysis of a larger sample of distinct SNe Ia (see Brout et al. 2022 and references therein) have found $(w_0,w_a)=(-0.{841}_{-0.061}^{+0.066},-0.{65}_{-0.32}^{+0.28})$ when combined with Planck and BAO, which is moderately consistent with the cosmological constant at 2σ. In addition, the Dark Energy Spectroscopic Instrument (DESI Collaboration 2016, 2025) has found support for w₀ > −1, w_a < 0 when combining BAO and Planck CMB data, which is inconsistent with the ΛCDM cosmology at ≳2.5σ. One of the targets of Euclid is to provide a ∼2% and ∼10% measurement on w₀ and w_a, respectively (Laureijs et al. 2011), which will help in further investigating the nature of dark energy and the validity of the standard ΛCDM model.

The structure of this paper is as follows. In Sect. 2, we discuss previous attempts to investigate tomographic binning strategies for weak lensing surveys and highlight the need to determine an appropriate strategy for Euclid. In Sects. 3 and 4, we review the analytical estimators used to measure the 3 × 2 pt signal and present our method to simulate mock Stage IV-like galaxy catalogues. In Sects. 5 and 6, we discuss the inference method that we use to derive constraints on the (w₀, w_a) parameters from the 3 × 2 pt signal, exploring these constraints for different tomographic binning strategies applied to our mock catalogues. In Sect. 7, we explore the effects of catastrophic photometric redshift uncertainties. In Sect. 8, we discuss the overall conclusions of our study to determine an optimum tomographic binning method. Finally, in the appendix, we review the validation testing we conducted for our simulation method and discuss the finer details of the modelling and computational choices that we implemented in this work.

2. Motivations

To construct a tomographic measurement, two analysis choices must be made: the number of bins to split the observed galaxy sample into and how to determine the bin boundaries along the radial (line-of-sight) direction. Each of these choices are governed by both scientific and practical considerations.

With regard to choosing the number of tomographic bins, it is generally expected that a greater number of bins will increase the information captured in a shear analysis, but at a higher computational cost. However, there will be a point at which the information gain starts to saturate, since adding more sample bins increases the Poisson shot noise per bin, which will dominate over the underlying cosmological shear signal of interest (see Sect. 3 for more details). Additionally, the cosmological volume enclosed by a single bin reduces as a larger number of bins is used, resulting in a comparatively lower gain in information, as the lensing signal over a given bin displays less of a change.

Several studies have been conducted to explore this issue and determine the ‘optimum’ tomography for a Stage IV-like weak lensing survey. For example, Taylor et al. (2018a) used a principal component analysis (PCA) to compare the suitability of a range of tomographic binning strategies with respect to the 3D cosmic shear method. For a Euclid-like weak lensing signal, they found that approximately ten tomographic bins with equal populations of galaxies can sufficiently capture the majority (∼97%) of the cosmological information of the 3D approach, but the tomographic method starts to fail at higher redshifts due to information loss and computational cost. Comparatively, a tomographic set-up with bins equally spaced in redshift does not capture as much information as the equipopulated case for a small number of bins; however, it can extract ∼99% of the information with respect to the 3D method for a large number of bins ( ∼ 50). Regardless of binning choice, the computational advantage of tomography is highlighted with respect to a 3D analysis.

This result is consistent with Sipp et al. (2021), who explored different tomographic strategies to assess the statistical properties of a predicted Euclid-like cosmic shear signal. For a range of science metrics, including minimising the uncertainties on individual cosmological parameter constraints and the optimisation of the Euclid dark energy figure of merit (Euclid Collaboration: Blanchard et al. 2020), they found that the optimum tomographic configuration is very close to the equipopulated scenario, with a small number of bins (four to five) sufficient for the cosmological analysis.

In contrast, Kitching et al. (2019) used a self organising map technique (SOM; Kohonen 1990) to define bin boundaries in photometric colour space, while optimising the signal-to-noise ratio (S/N) on the dark energy equation of state parameter, w₀, measured from the cosmic shear power spectrum. For a small number of bins (< 5), the best-performing binning methodology identified by the SOM was found to be a close approximation of the case where the bins are defined to have equal width in redshift.

This issue is further complicated if a ‘3 × 2 pt’ analysis of the data is targeted (see Sect. 3.1), where the shear signal is combined with the tomographic galaxy number density (angular clustering) and galaxy-galaxy lensing (cross-correlation of weak lensing and angular clustering) signals into a joint data vector. The full 3 × 2 pt data vector will offer more constraining power than a shear-only analysis (see e.g. Tutusaus et al. 2020). However, the choice of an appropriate tomographic binning configuration may also depend on which cosmological parameter is being targeted and which part of the data vector is most sensitive to the corresponding science case.

In particular, by considering the photometric redshift estimation techniques for Euclid, Euclid Collaboration: Pocino et al. (2021) found that for a galaxy clustering and galaxy-galaxy lensing analysis, tomographic bins equally spaced in redshift yield a higher dark energy figure of merit (FOM) than equipopulated bins. Furthermore, a relatively large number of bins (13) maximises the amount of cosmological information extracted.

Additionally, Zuntz et al. (2021) have coordinated a tomography challenge in which a range of tomographic binning configurations, guided by machine learning approaches on a test sample, were proposed by the community to target a variety of metrics such as the 3 × 2 pt power spectra S/N and the dark energy FOM. The findings suggest that the optimum tomographic binning set-up used for the analysis is strongly dependent on both the chosen 3 × 2 pt observable and the specific cosmological parameter FOM that is targeted. However, a source galaxy distribution that is split into redshift bins that are equally spaced in fiducial comoving distance was found to be a good general choice.

It is clear that these studies reach a variety of different conclusions and, in addition, they are not all performed for the same target survey set-up. This makes it challenging to use them for choosing an optimum binning scheme for Euclid. Hence, the aim of this work is to present an alternative method to explore a range of 3 × 2 pt tomographic configurations and investigate an optimum tomographic binning approach for the specific survey characteristics expected for the analysis of Euclid’s first data release (DR1).

We construct a cosmological simulation pipeline that generates full Euclid-like cosmic shear catalogues realised on the sky from a given fiducial cosmology, with realistic observational uncertainties such as errors from photometric redshift (‘photo-z’) estimation and shape noise injected into the catalogue observables. From our simulated catalogues we measure the 3 × 2 pt signal and compare constraints on the dark energy equation of state for each tomographic set-up chosen, which will help to further identify the optimum configuration for the cosmological analysis of the 3 × 2 pt observables in the Euclid survey. Our simulation pipeline includes the following stages:

1. The sampling of an observed galaxy population from a predicted Euclid-like galaxy redshift distribution and spatial number density;

2. The calculation of a theoretical 3 × 2 pt signal over a given redshift range and fiducial cosmology;

3. Simulations of the correlated cosmic shear and galaxy density fields on the sky as a set of 2D maps that approximates the 3D evolution of the 3 × 2 pt signal over time;

4. The sampling of galaxy angular positions and weak lensing observables from the density and shear fields, respectively;

5. The injection of observational uncertainties from photo-z estimation and galaxy shape noise into the individual galaxies to generate a mock catalogue mimicking a realistic survey.

We provide a full description of this simulation method and each stage in Sect. 4. We note that the previously discussed tomography studies typically focus on using numerical methods to explore a theoretical shear signal. In comparison, we highlight that our method offers the advantage of imitating a realistic survey observation and measuring the corresponding 3 × 2 pt signal from an underlying cosmology, which provides a more direct and realistic framework with which to investigate the cosmological impact of different tomographic strategies.

3. Background and theory

3.1. Cosmic shear and the 3 × 2 pt signal

Gravitational lensing occurs when the foreground potential of a massive object distorts the light observed from a background source. In the regime of weak lensing (see e.g. Bartelmann & Schneider 2001 for a detailed review), where the distortion effects are small (i.e. at the per cent level), this phenomenon leads to subtle changes to a background galaxy’s size and shape, which can be described using the convergence, κ (spin-0 quantity), and shear, γ (spin-2 quantity), respectively.

For a population of multiple background galaxies, weak lensing leads to correlations in the distortions of individual galaxy shapes. Using tomography, the galaxy sample is split into a set of redshift bins and the cosmic shear signal is then measured using the angular power spectrum under the Limber approximation (Limber 1953, LoVerde & Afshordi 2008) as C_ij(ℓ), where i, j denote the pairs of redshift bins in the sample. Explicitly, for a flat Universe, the tomographic shear power spectrum, $C_{\mathit{ij}}^{\gamma \gamma }(\ell )$, is given by

$$\begin{array}{c}\hfill C_{\mathit{ij}}^{\gamma \gamma }(\ell )={\displaystyle {\int }_0^{{\chi }_{\mathrm{H}}}}\mathrm{d}\chi \frac{q_i^{\kappa }(\chi )q_j^{\kappa }(\chi )}{{\chi }^2}{P}_{\delta }(\frac{\ell +1/2}{\chi },\chi ),\end{array}$$(1)

where we take the matter density parameter at the present day, Ω_m = 0.3, and the Hubble constant, H₀ = 70 km s⁻¹ Mpc⁻¹, throughout this work. Here, the comoving radial coordinate, χ, is integrated up to χ_H, which is the comoving distance to the cosmic horizon. Additionally, a is the scale factor, P_δ is the matter power spectrum, and $q_i^{\kappa }(\chi )$ is the lensing efficiency,

$$\begin{array}{c}\hfill q_i^{\kappa }(\chi )=\frac{3}{2}{\mathrm{\Omega }}_{\mathrm{m}}{\left(\frac{H_0}{c}\right)}^2\frac{\chi }{a(\chi )}{\displaystyle {\int }_{\chi }^{{\chi }_{\mathrm{H}}}}\mathrm{d}\chi \prime \frac{n_i(z)\mathrm{d}z/\mathrm{d}\chi \prime }{{\overline{n}}_i}\frac{\chi \prime -\chi }{\chi \prime },\end{array}$$(2)

which describes the redshift distribution of the background galaxies, n_i(z), weighted by the angular number density,

$$\begin{array}{c}\hfill {\overline{n}}_i={\displaystyle {\int }_0^{\infty }}\mathrm{d}z{n}_i(z),\end{array}$$(3)

where the redshift z = z(χ), and i denotes a given tomographic bin. A full 3 × 2 pt analysis also uses the galaxy angular positions on the sky, δ_g(θ), and the corresponding galaxy clustering radial kernel, q_i^δ_g,

$$\begin{array}{c}\hfill q_i^{{\delta }_g}(k,\chi )=b_i(k,z)n_i(z)\frac{\mathrm{d}z}{\mathrm{d}\chi },\end{array}$$(4)

where b_i(k, z) is the galaxy bias in tomographic bin i, which describes the spatial relation between the underlying matter field and the observed distribution of galaxies as a function of redshift and scale, k, where k = (ℓ+1/2)/χ.

In addition to weak lensing, the 3 × 2 pt signal consists of the galaxy clustering angular power spectrum,

$$\begin{array}{c}\hfill C_{\mathit{ij}}^{{\delta }_g{\delta }_g}(\ell )={\displaystyle {\int }_0^{{\chi }_{\mathrm{H}}}}\mathrm{d}\chi \frac{q_i^{{\delta }_g}(\frac{\ell +1/2}{\chi },\chi )q_j^{{\delta }_g}(\frac{\ell +1/2}{\chi },\chi )}{{\chi }^2}{P}_{\delta }(\frac{\ell +1/2}{\chi },\chi ),\end{array}$$(5)

and the cross correlation of the clustering and shear kernels, the galaxy-galaxy lensing power spectrum, C_ij^δ_gγ(ℓ),

$$\begin{array}{c}\hfill C_{\mathit{ij}}^{{\delta }_g\gamma }(\ell )={\displaystyle {\int }_0^{{\chi }_{\mathrm{H}}}}\mathrm{d}\chi \frac{q_i^{{\delta }_g}(\frac{\ell +1/2}{\chi },\chi )q_j^{\kappa }(\chi )}{{\chi }^2}{P}_{\delta }(\frac{\ell +1/2}{\chi },\chi ).\end{array}$$(6)

We note that these expressions for the 3 × 2 pt signal do not include the effects of intrinsic alignments, which we did not consider in this work.

3.2. Estimators for the 3 × 2 pt signal

3.2.1. Cosmic shear estimators

The shear, γ(θ), which describes how weak lensing ‘stretches’ and ‘rotates’ an image, is a spin-2 field. It can be expressed as a complex quantity,

$$\begin{array}{c}\hfill \gamma (\mathit{\theta })={\gamma }_1(\mathit{\theta })+\mathrm{i}{\gamma }_2(\mathit{\theta }).\end{array}$$(7)

On the full sky, the shear signal can be expressed in spherical harmonic space via the _±2Y_ℓ, m(θ) spin-2 harmonics defined on the sphere,

$$\begin{array}{c}\hfill {\gamma }_1(\mathit{\theta })\pm \mathrm{i}{\gamma }_2(\mathit{\theta })={\displaystyle \sum _{\ell ,m}}{}_{\pm 2}{{\stackrel{\sim }{\gamma }}_{\ell ,m}}_{\pm 2}{Y}_{\ell ,m}(\mathit{\theta }),\end{array}$$(8)

where the ${\stackrel{\sim }{\gamma }}_{\ell ,m}$ coefficients can be expressed using the gradient E- and curl B-modes analogous to the Q and U Stokes parameters for an electromagnetic field (Crittenden et al. 2002, Kamionkowski et al. 1998, Stebbins 1996, Kaiser 1992),

$$\begin{array}{c}\hfill {}_{\pm 2}{\stackrel{\sim }{\gamma }}_{\ell ,m}=E_{\ell ,m}\pm \mathrm{i}{B}_{\ell ,m}={\displaystyle \int {({\gamma }_1\pm \mathrm{i}{\gamma }_2)}_{\pm 2}}{Y}_{\ell ,m}^{\ast }(\mathit{\theta })\mathrm{d}\mathit{\theta },\end{array}$$(9)

and the coefficients for each mode are given by

$$\begin{array}{cc}\hfill E_{\ell ,m}& =\frac{1}{2}{\displaystyle \int [\gamma {(\mathit{\theta })}_{+2}{Y}_{\ell ,m}^{\ast }(\mathit{\theta })+{\gamma }^{\ast }{(\mathit{\theta })}_{-2}{Y}_{\ell ,m}^{\ast }(\mathit{\theta })]}\mathrm{d}\mathit{\theta },\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill B_{\ell ,m}& =-\frac{\mathrm{i}}{2}{\displaystyle \int [\gamma {(\mathit{\theta })}_{+2}{Y}_{\ell ,m}^{\ast }(\mathit{\theta })-{\gamma }^{\ast }{(\mathit{\theta })}_{-2}{Y}_{\ell ,m}^{\ast }(\mathit{\theta })]}\mathrm{d}\mathit{\theta },\hfill \end{array}$$(11)

where each of these integrals are done over the full sphere. From the harmonic coefficients, we define the E-mode, B-mode, and EB angular power spectra as

$$\begin{array}{cc}\hfill ⟨E_{\ell ,m}{E}_{{\ell }^{\prime },m^{\prime }}^{\ast }⟩& ={\delta }_{\ell {\ell }^{\prime }}{\delta }_{m{m}^{\prime }}{C}^{\mathit{EE}}(\ell ),\hfill \end{array}$$(12)

$$\begin{array}{cc}\hfill ⟨B_{\ell ,m}{B}_{\ell \prime ,m\prime }^{\ast }⟩& ={\delta }_{\ell \ell \prime }{\delta }_{mm\prime }{C}^{\mathit{BB}}(\ell ),\hfill \end{array}$$(13)

$$\begin{array}{cc}\hfill ⟨E_{\ell ,m}{B}_{\ell \prime ,m\prime }^{\ast }⟩& ={\delta }_{\ell \ell \prime }{\delta }_{mm\prime }{C}^{\mathit{EB}}(\ell ),\hfill \end{array}$$(14)

where δ_XY is the Kronecker delta and the triangular brackets denote an ensemble average over many realisations. Since we are only observing one Universe and are limited by the number of m modes available, each C(ℓ) value is distributed among 2ℓ+1 degrees of freedom. Consequently, we can only measure an estimate of the true power spectrum, $\widehat{C}(\ell )$, which, for the given tomographic bins i, j, we define as

$$\begin{array}{c}\hfill {\widehat{C}}_{\mathit{ij}}^{\alpha \beta }(\ell )=\frac{1}{2\ell +1}{\displaystyle \sum _{m=-\ell }^{+\ell }}{\alpha }_{\ell ,m}^{(i)}{\beta }_{\ell ,m}^{(j)\ast },\end{array}$$(15)

where the indices α, β denote the modes for the shear power spectra, $[{\widehat{C}}_{\mathit{ij}}^{\mathit{EE}}(\ell ),{\widehat{C}}_{\mathit{ij}}^{\mathit{BB}}(\ell ),{\widehat{C}}_{\mathit{ij}}^{\mathit{EB}}(\ell )]$. In the weak lensing regime, assuming the Born approximation, it has been shown that to leading order, the shear field contains negligible power arising from the B-mode curl component (Hilbert et al. 2009, Krause & Hirata 2010), and the coefficient B_ℓ, m is taken as zero (e.g. Giocoli et al. 2016).

Therefore, the shear signal is typically fully characterised by the E-mode, which we use to estimate the theoretical shear power spectrum in this work (i.e. $C_{\mathit{ij}}^{\gamma \gamma }(\ell )=⟨{\widehat{C}}_{\mathit{ij}}^{\mathit{EE}}(\ell )⟩$). The B-mode power spectrum can be used as an important null test for systematics. Finally, observed shear power spectra will typically contain a noise component. For this work, we consider the shape noise, associated with the intrinsic elliptical shape of galaxies (see also Sect. 4.6.2). This effect is random and uncorrelated between tomographic bins, so the final measured shear signal, ${\widehat{C}}_{\mathit{ij}}^{EE,\mathrm{obs}}(\ell )$, contains a contribution from the shape noise in the autocorrelation power spectra,

$$\begin{array}{c}\hfill {\widehat{C}}_{\mathit{ij}}^{EE,\mathrm{obs}}(\ell )={\widehat{C}}_{\mathit{ij}}^{\mathit{EE}}(\ell )+N_i^{ϵ}(\ell ){\delta }_{\mathit{ij}},\end{array}$$(16)

where $N_i^{ϵ}(\ell )$ is the power spectrum of the shape noise in the autocorrelation tomographic power spectra.

3.2.2. Galaxy clustering estimators

The angular galaxy clustering signal is estimated using the normalised galaxy number count on a pixelised sky, which traces the underlying matter distribution. For the tomographic bin, i, we formulate the real-space overdensity estimator,

$$\begin{array}{c}\hfill {\widehat{\delta }}_g^{(i)}(\mathit{\theta })=\frac{n_i^{\mathrm{p}}(\mathit{\theta })-{\overline{n}}_i^{\mathrm{p}}}{{\overline{n}}_i^{\mathrm{p}}},\end{array}$$(17)

where $n_i^{\mathrm{p}}(\mathit{\theta })$ is the per-pixel integer galaxy count in the redshift bin and ${\overline{n}}_i^{\mathrm{p}}$ is the average number of galaxies observed per pixel in the same bin.

Analogously to the shear field, the measured galaxy overdensity field can be represented in spherical harmonic space via the overdensity coefficients, d_ℓ, m,

$$\begin{array}{c}\hfill {\widehat{\delta }}_g^{(i)}(\mathit{\theta })={\displaystyle \sum _{\ell ,m}}{d}_{\ell ,m}^{(i)}{Y}_{\ell ,m}(\mathit{\theta }).\end{array}$$(18)

From the d_ℓ, m coefficients, we can measure the angular power spectrum of the estimated density field,

$$\begin{array}{c}\hfill D_{\mathit{ij}}(\ell )=\frac{1}{2\ell +1}{\displaystyle \sum _{m=-\ell }^{+\ell }}{d}_{\ell ,m}^{(i)}{d}_{\ell ,m}^{(j)\ast }.\end{array}$$(19)

The measured power spectrum D_ij(ℓ) contains contributions from both the underlying galaxy clustering signal and the Poisson noise associated with the discretised number counts in Eq. (17). Explicitly, the Poisson noise contributes to the auto-power spectra in a tomographic analysis, while the power spectrum of the Poisson noise in a tomographic bin i is given by (e.g. Loureiro et al. 2019, Nicola et al. 2020)

$$\begin{array}{c}\hfill N_i^{\mathrm{P}}(\ell )={\mathrm{\Omega }}_{\mathrm{p}}\frac{\overline{w}}{{\overline{n}}_i^{\mathrm{p}}},\end{array}$$(20)

where Ω_p is the angular size in steradians of the pixels that form the pixelised overdensity map, and $\overline{w}$ is the average value of the mask over all pixels on the sky. For this work, we used a binary mask where the pixels were set to 1 inside the observed footprint and 0 outside the footprint. The Poisson noise and underlying galaxy clustering signal are uncorrelated, so the angular power spectrum of the estimated density field is simply the sum of the two components:

$$\begin{array}{c}\hfill D_{\mathit{ij}}(\ell )=C_{\mathit{ij}}^{{\delta }_g{\delta }_g}(\ell )+N_i^{\mathrm{P}}(\ell ){\delta }_{\mathit{ij}}.\end{array}$$(21)

3.2.3. Galaxy-galaxy lensing estimators

To measure the galaxy-galaxy lensing signal, we can estimate the angular power spectrum from the harmonic coefficients of the density and shear fields. By taking B_ℓ, m = 0, the galaxy-galaxy lensing power spectrum is characterised by the cross-correlation of the density field and the shear E-mode,

$$\begin{array}{c}\hfill {\widehat{C}}_{\mathit{ij}}^{{\delta }_g\gamma }(\ell )={\widehat{C}}_{\mathit{ij}}^{{\delta }_{g}E}(\ell )=\frac{1}{2\ell +1}{\displaystyle \sum _{m=-\ell }^{+\ell }}{d}_{\ell ,m}^{(i)}{E}_{\ell ,m}^{(j)\ast },\end{array}$$(22)

for the given tomographic bins i, j.

3.3. Estimators on the cut sky

The above framework is defined for the case of a full sky being observed. In reality there will be partial sky coverage due to unobserved regions, or bright stars, meaning that a survey mask, W(θ), is applied to the shear and overdensity fields, which leads to mode mixing in the derived power spectra. These ‘pseudo’ power spectra measured on the cut sky, ${\stackrel{\sim }{C}}_{\ell }$, can be related to the true underlying full-sky elements, C_ℓ, by a mode-mixing matrix, M_ℓℓ′ (see e.g. Peebles 1973, Brown et al. 2005, and Hikage et al. 2011), which contains the coupled harmonic space information about the survey mask,

$$\begin{array}{c}\hfill ⟨{\stackrel{\sim }{C}}_{\ell }^{(i,j)}⟩={\displaystyle \sum _{\ell \prime }}{M}_{\ell \ell \prime }^{(i,j)}{C}_{\ell \prime }^{(i,j)}.\end{array}$$(23)

Here, the angular brackets denote an expectation value, and the tomographic dependence of the power spectra are introduced using the superscripts (i, j). In addition, C_ℓ denotes the elements of a column vector representing C(ℓ) evaluated at integer ℓ, which we will adopt as a shorthand description of the power spectrum for the remainder of this work.

For our analysis, we worked with the Pseudo-C_ℓ power spectrum of the signal measured directly from the catalogue data, ${\stackrel{\sim }{S}}_{\ell }$. The expected value for this measured signal is calculated as the sum of the cosmological signal on the masked sky and the noise,

$$\begin{array}{cc}\hfill ⟨{\stackrel{\sim }{S}}_{\ell }^{(i,j)}⟩& =⟨{\stackrel{\sim }{C}}_{\ell }^{(i,j)}⟩+⟨N_{\ell }^{(i,j)}⟩\hfill \\ \hfill & =\left({\displaystyle \sum _{\ell \prime }}{M}_{\ell \ell \prime }^{(i,j)}{C}_{\ell \prime }^{(i,j)}\right)+⟨N_{\ell }^{(i,j)}⟩.\hfill \end{array}$$(24)

Here, the noise term $N_{\ell }^{(i,j)}$ arises due to shape noise for the cosmic shear field and Poisson noise for the galaxy clustering field, and is only non-zero for the auto-power spectra. A full description of the matrix M_ℓℓ′ for each 3 × 2 pt component is given in Appendix A. Throughout our work we measure the 3 × 2 pt Pseudo-C_ℓ power spectra using NaMaster (Alonso et al. 2019).

3.4. Bandpowers

Following the technique of Hivon et al. (2002), it is often convenient to bin the Pseudo-C_ℓ measured on the cut sky in ℓ space. This gives a reduction in the errors on the power spectrum estimator and allows the mixing matrix (Eq. (23)) to remain invertible in the case of large and complex sky cuts. For a set of bins, denoted by the index b, with bin boundaries $\{{\ell }_{\mathrm{low}}^b<{\ell }_{\mathrm{high}}^b<{\ell }_{\mathrm{low}}^{b+1}<{\ell }_{\mathrm{high}}^{b+1}<\dots \}$ we take a binning operator, P_bℓ, defined as

$$\begin{array}{c}\hfill P_{b\ell }=\{\begin{array}{cc}\frac{1}{2\pi }\frac{\ell (\ell +1)}{{\ell }_{\mathrm{low}}^{b+1}-{\ell }_{\mathrm{low}}^b},\hfill & \mathrm{if}2\le {\ell }_{\mathrm{low}}^b\le \ell <{\ell }_{\mathrm{low}}^{b+1}\hfill \\ 0,\hfill & \mathrm{otherwise}.\hfill \end{array}\end{array}$$(25)

The elements of the binned Pseudo-C_ℓ power spectrum are termed ‘bandpowers’ and calculated as

$$\begin{array}{c}\hfill {\stackrel{\sim }{\mathcal{C}}}_b={\displaystyle \sum _{\ell }}{P}_{b\ell }{\stackrel{\sim }{C}}_{\ell }.\end{array}$$(26)

We note that our approach to estimate angular power spectra is similar to the analysis that is expected to be applied to the Euclid DR1 sample (Euclid Collaboration: Tessore et al. 2025).

4. Methods: Simulating mock catalogues

In this section, we describe the modelling process and present the simulation pipeline with which we construct Stage IV-like mock cosmic shear catalogues¹.

For a given fiducial cosmology and background galaxy redshift distribution, we first generated a theoretical prediction for the full-sky 3 × 2 pt power spectra, $[C_{\mathit{ij}}^{\gamma \gamma }(\ell ),C_{\mathit{ij}}^{{\delta }_g{\delta }_g}(\ell ),C_{\mathit{ij}}^{{\delta }_g\gamma }(\ell )]$, at multiple, finely spaced redshift intervals. From the 3 × 2 pt data vectors, we generated a realisation on the sky of the galaxy clustering, δ_g(θ), and shear, [γ₁(θ),γ₂(θ)], fields at each point in redshift space covered by our chosen galaxy sample. In this way, we approximated the full 3D cosmological information with a set of 2D sky maps that represents the behaviour of the 3 × 2 pt observable fields as the background Universe evolves over time and redshift.

Finally, we performed a random sampling routine and Poisson sampled the galaxy clustering fields at each redshift to populate the sky with galaxies that trace the underlying matter distribution. We constructed a mock catalogue by then assigning the weak lensing observables [κ(θ),γ₁(θ),γ₂(θ)] from the field values of the cosmic shear maps at the positions of each pixel hosting a galaxy.

This simulation pipeline consisted of a number of distinct stages of calculation that were executed in series. We describe each step in detail below. In Fig. 1 we present a flowchart diagram showing the workflow of our simulation to produce the final mock catalogues.

Fig. 1. Flowchart describing the structure of our simulation to generate Stage IV-like mock galaxy catalogues from 2D realisations on the sky of the 3 × 2 pt signal that traces an underlying redshift distribution and fiducial cosmology. The red dashed and green dotted sections represent the simulation stages that are executed using CosmoSIS and FLASK, respectively.

4.1. Generation of the n(z) redshift sample

The first step is the construction of a galaxy population that follows a given n(z) redshift distribution. Following Sipp et al. (2021), we take a redshift probability distribution, p(z), that has the functional form typical of a Stage IV survey,

$$\begin{array}{c}\hfill p(z)\propto {\left(\frac{z}{z_0}\right)}^{2}exp[-{\left(\frac{z}{z_0}\right)}^{1.5}],\end{array}$$(27)

where z₀ is the characteristic redshift constant chosen such that the median redshift of the sample satisfies $z_{\mathrm{median}}=\sqrt{2}{z}_0$. We expect that this chosen form for p(z) will be a good approximation for the Euclid DR1 galaxy sample. We note that the measured cosmological constraints would change if the real galaxy distribution is observed to be significantly different.

Following the definition of the redshift probability distribution, we then set the simulated redshift range, [z_min, z_max], and a precision, dz, to which galaxy redshifts are generated. A raw galaxy sample, n(z), is then created by performing a sampling routine to draw a chosen total number of galaxies, N_gal, that traces the p(z) redshift probability distribution.

We highlight that the functional form for p(z), the total observed galaxy density, N_gal, and the redshift range and precision, [z_min, z_max, dz], are left as free choices in our simulation that can arbitrarily be changed to match the specifications of a chosen survey.

4.2. Calculation of theoretical 3 × 2 pt data vector for each galaxy redshift

Following the simulation of a galaxy population, n(z), we aim to generate 2D realisations on the sky of the galaxy clustering and weak lensing fields for each redshift that we sample galaxies at. To do this, we first need a prediction for the 3 × 2 pt power spectra at each redshift. This set of power spectra is constructed by converting the raw n(z) sample into a ‘binned’ population and creating a table where every column contains the galaxy number count at each specific, discrete redshift that has been sampled in the p(z).

This process allows us to retain the correlations in the 3 × 2 pt fields between the different, finely spaced redshift slices. In Fig. 2 we show the Euclid DR1-like redshift distribution we consider for our work and demonstrate the effects of binning and normalising the population. This ultimately allows us to generate a set of 2D slices of the 3 × 2 pt signal that approximates the full 3D cosmological information.

Fig. 2. Illustration of the method used to simulate a mock weak lensing survey. 2D maps of the correlated 3 × 2 pt fields are generated at finely sampled points in redshift. These 2D maps are used to approximate the full 3D cosmological information of the 3 × 2 pt signal. From the overdensity fields, we Poisson sample a galaxy population that traces the underlying n(z) distribution. We then assign the correlated weak lensing observables to each galaxy from the shear field values at the galaxy’s angular position on the sky at a given redshift. We show an early approximation of the Euclid DR1 footprint at the bottom of the figure which we have used to create our Euclid-like simulations. The ‘observed’ region is shown in yellow. Note: the actual Euclid DR1 footprint will be significantly different to that shown here.

In addition to the ‘binned’ n(z), we require a matter power spectrum, P_δ(k), which consists of both a linear and non-linear component, and a galaxy bias model, b(k, z), which represents the scale- and redshift-dependent relation between the distribution of galaxies and the underlying matter field that they trace. For the purposes of this work we assume a linear model and take the galaxy bias as a global constant, b. While it is beyond the scope of this study, we note that our framework has the flexibility to incorporate both scale- and redshift-dependence in the bias model, and the effects of these choices could be investigated in a future study.

To construct the final theoretical 3 × 2 pt power spectra at each redshift, we use CosmoSIS (Zuntz et al. 2015), a parameter estimation code which specialises in the joint modelling of cosmological power spectra and exploration of parameter constraints for cosmic shear and galaxy clustering studies. Explicitly, we run a CosmoSIS pipeline, which executes the following modules:

• Code for Anisotropies in the Microwave Background (CAMB, Lewis et al. 2000), used to model the linear matter power spectrum, which describes the evolution of density perturbations in the early Universe.

• halofit_takahashi (Smith et al. 2003, Takahashi et al. 2012), used to model the matter power spectrum in the non-linear regime by scaling the linear component using fitting functions derived from simulations.

• constant_bias, used to define the functional form of the galaxy bias for modelling the galaxy clustering and galaxy-galaxy lensing signals. We made the assumption that on the large scales we focus on in this work, the galaxy bias can be sufficiently described using a constant.

• load_nz, used to process and normalise the ‘binned’ n(z) galaxy population based on the user-defined redshift boundaries of each column. The resulting n(z) table products can then be used to generate the 3 × 2 pt signal.

• project_2d, used for calculation of the theoretical 3 × 2 pt data vector by evaluating the Limber approximation to project the 3D line-of-sight information into 2D angular power spectra.

For the background cosmology in the CosmoSIS pipeline and our overall simulation, we worked with a flat w₀w_a Universe with Ω_m = 0.3 and H₀ = 70 km s⁻¹ Mpc⁻¹. For the matter power spectrum generated using the CAMB module in CosmoSIS, we take the latest best fit Planck values (see Planck Collaboration VI 2020 and references therein):

• optical depth to reionisation, τ = 0.05;

• scalar spectral index, n_s = 0.96;

• primordial amplitude of power, A_s = 2.1 × 10⁹.

For the constant galaxy bias model, we chose a value b = 0.4, which is chosen to avoid numerical issues encountered when generating the shear and clustering fields on the sky (see Sect. 4.4 for a further discussion of this point). Lastly, we chose to work with the fiducial values (w₀, w_a) = (−1, 0). We emphasise that the values chosen for the cosmological parameters listed in this section are free variables that can arbitrarily be changed in our simulation pipeline.

4.3. Generation of weak lensing and clustering maps

To simulate the cosmological 3 × 2 pt signal on the sky, we generated a set of full-sky maps of the observables contained in the 3 × 2 pt power spectra calculated at each redshift sampled from the n(z) distribution.

For the galaxy clustering component, we require the observed galaxy overdensity field, δ_g(θ), a scalar quantity, while for weak lensing, there are three observable quantities on the sky: the scalar convergence field, κ(θ), describing the lensing-induced magnification of a source; and the components of the spin-2 shear field, γ₁(θ) and γ₂(θ).

These fields can be generated on the sky using HEALPix (Górski et al. 2005), a pixelisation technique which offers the ability to represent data on the sphere. To generate the pixelised maps from the 3 × 2 pt data vector, we used the Full-sky Lognormal Astro-fields Simulation Kit (FLASK; Xavier et al. 2016), which creates realisations of an arbitrary set of correlated fields on the sky from their given tomographic power spectra via a Cholesky decomposition.

While FLASK has the capacity to generate both lognormal and Gaussian-distributed realisations of the chosen fields, for the purposes of this work, we chose to consider only Gaussian fields, on the premise that complete information of a Gaussian field is directly captured in the power spectrum.

4.4. Construction of mock catalogue

The FLASK routine yields a 2D HEALPix map of the galaxy clustering, convergence, and cosmic shear fields at each redshift, z_s, which is sampled in the p(z) distribution. To construct a mock catalogue, we first place galaxies at angular positions on the sky by Poisson-sampling the galaxy density field, δ_g^z_s(θ), at each redshift, z_s. We convert this into an observed galaxy population by Poisson-sampling the cell value to generate a pixelised integer number count of galaxies on the sky at the redshift slice, z_s,

$$\begin{array}{c}\hfill n_{\mathrm{gal}}^{z_s}(\mathit{\theta })=\mathrm{Poisson}\left\{{\overline{n}}_{\mathrm{gal}}^{z_s}[1+{\delta }_g(\mathit{\theta })]\right\},\end{array}$$(28)

where ${\overline{n}}_{\mathrm{gal}}^{z_s}$ is the average number of galaxies per pixel on the sky at redshift z_s, which is defined from the sampled galaxy distribution, n(z), in Eq. (27). In principle, this may vary on the sky due to effects such as the variable depth of a survey, but we did not consider such effects in this work.

To execute a Poisson sampling of the density field, the minimum value of the field must satisfy δ_g(θ)≥ − 1. Due to the nature of Gaussian realisations on the sky, there can be random pixels generated that have the unphysical value δ_g < −1. To reduce the impact of this effect, we took a redshift-independent galaxy bias value of b = 0.4, which reduces the frequency of these random occurrences. We note that such values of b < 1, while less realistic, are a choice taken by other studies in the literature (e.g. Tessore et al. 2023) that have investigated Gaussian realisations of the 3 × 2 pt fields on the sky. Quantitatively, we would expect that (to first order) this choice will not affect the cosmic shear measurements, since the bias term does not feature in the weak lensing kernel (Eqs. (1), (2)). For the angular clustering and galaxy-galaxy lensing signals, the bias term acts to scale the amplitude of the power spectrum. Since we used a constant bias over all angular scales and redshift, we expect that the choice of b would change the absolute values of constraints on (w₀, w_a) measured from our simulations, but would not alter the relative measurements that compare the constraints between different tomographic binning schemes.

In addition, we rejected the redshift range z < 0.3, in which the evolution of galaxy clustering is more likely to produce such unphysical pixels. However, we note that at this redshift range there is also the least amount of information contained in the cosmic shear signal, since there is comparatively little foreground structure to induce the weak lensing effect.

We assigned an angular (RA, Dec) position to each galaxy by placing them at the centre of the pixel they are Poisson-sampled in. The final galaxy ‘observation’ to create a mock shear catalogue is then a straightforward assignment process. Each galaxy is described with a data vector of its (RA, Dec, z) values. We then append to this observed data the value of the cosmic shear [κ(θ),γ₁(θ),γ₂(θ)] fields at the specific sky position and redshift that the galaxy is found at. By repeating this procedure for each galaxy in our sampled population, we can produce a final cosmic shear catalogue, which:

• consists of a galaxy population that traces a survey-specified number density and predicted redshift probability distribution, p(z);

• statistically represents on the sky the information contained in the theoretical 3 × 2 pt power spectra calculated from a chosen underlying cosmology.

4.5. Pipeline overview

The overall workflow of our pipeline is presented in Fig. 1. Our simulation-based approach offers an approximate, but relatively rapid technique for exploring the impact of different analysis choices on the precision of cosmological constraints that can be achieved with Euclid. While a full N-body simulation to create a catalogue would provide the most complete astrophysical information, it is extremely computationally costly to generate and would typically yield only a few realisations. Correspondingly, introducing or removing any additional astrophysical or observational features in the simulation to test their impact on observables and cosmological measurements would be complex and inefficient.

As discussed in Sect. 2, previous studies on the choice of 3 × 2 pt tomography have typically focussed only on exploring a pre-defined signal, without generating a mock observation on the sky. Our method offers an alternative to such approaches. Since we are replicating a real analysis by measuring the 3 × 2 pt signal from the simulated data, we have a unique capability to introduce complex observational and astrophysical effects at the catalogue level, and rapidly explore their effects on measurements of cosmological parameters. We note that by using the map-based simulation approach, we can practically generate and work with many realisations of the Universe to further improve the statistical power of the analysis. To demonstrate our method, we focussed on two sources of uncertainty: photo-z errors and shape noise. We discuss them in the following sections and emphasise that a comprehensive suite of systematics can be explored in a follow-up investigation(s).

4.6. Simulated noise and errors

The basic execution of our pipeline contains intrinsic Poisson noise in the galaxy clustering signal, which has the noise power spectrum given in Eqs. (20), (21). In addition to the Poisson noise, we can also increase the realism of our simulation by introducing redshift uncertainties and shape noise in the galaxies, which are included as an optional feature in our simulation.

4.6.1. Redshift uncertainties

For a realistic survey, we require either spectroscopic or photometric analysis of individual galaxies to measure their redshift. While the spectroscopic measurement offers a very precise estimate of the galaxy’s true redshift, it is unfeasible to yield a spectroscopic redshift for every galaxy in a survey. Comparatively, there are much larger uncertainties when deriving redshifts via the photometric method, which relies on fitting the galaxy’s observed spectral energy distribution (see e.g. Euclid Collaboration: Ilbert et al. 2021 for an overview of redshift measurement techniques for Euclid).

When constructing the final survey catalogues in our pipeline (Sect. 4.4) we can optionally introduce uncertainty in the galaxy redshifts associated with such measurement techniques. Explicitly, we expect that, to first order, the photometric uncertainty on a measured redshift will be Gaussian-distributed with standard deviation, σ_z, dependent on the redshift,

$$\begin{array}{c}\hfill {\sigma }_z={\sigma }_z^{\mathrm{phot}}(1+z),\end{array}$$(29)

where ${\sigma }_z^{\mathrm{phot}}$ is a constant. To introduce this effect into our pipeline, for each galaxy at redshift z, we convert the ‘true’ redshift of the galaxy in the catalogue to a ‘measured’ redshift by drawing an estimated value from a Gaussian distribution with standard deviation σ_z and mean z.

In addition to the Gaussian-distributed uncertainty in the photo-zs, there will be some catastrophic outliers due to limitations in the photometric fitting procedure. Explicitly, the measured SED of the galaxy can be fit to a specific template (see e.g. Euclid Collaboration: Desprez et al. 2020 for an overview). However, confusion between given pairs of spectral lines can lead to a measured galaxy redshift that is catastrophically inaccurate. In the technique of Jouvel et al. (2011), we can model a ‘catastrophic’ redshift estimate using the expression,

$$\begin{array}{c}\hfill z_{\mathrm{cata}}=(1+z)\frac{{\lambda }_{\mathrm{break}-\mathrm{rf}}}{{\lambda }_{\mathrm{break}-\mathrm{cata}}}-1,\end{array}$$(30)

where λ_{break − rf} is the true rest-frame wavelength of a break feature in an SED model, and λ_{break − cata} is the wavelength of the break feature that is predicted in the fitting procedure. As a first-order approximation, in the redshift range z ≤ 2.5, Jouvel et al. (2011) suggest that the Lyman-α line, Lyman-break, Balmer-break and 4000-Å break (D4000), and permutations between pairs of these four features, will contribute to the catastrophic redshift estimation. Within our simulation, we include the option to model this effect by converting the true redshifts of a chosen percentage of the galaxy population into a catastrophic measurement based on a given set of spectral line pairs.

4.6.2. Shape noise

A further significant source of noise in a weak lensing analysis is the shape noise. Galaxies have an intrinsic ellipticity, ϵ^int, which acts as an irreducible source of confusion when estimating the weak lensing-induced shape distortion, γ. Under the assumption that galaxies have no preferred ellipticity on average, we can model the shape noise to be Gaussian-distributed with standard deviation σ(ϵ^int), known as the total ellipticity dispersion.

We introduce this as an optional effect in our simulation by drawing shape noise error quantities from a Gaussian distribution with standard deviation σ(ϵ^int). This error is then attributed to each galaxy by adding the component $\sigma ({ϵ}^{\mathrm{int}})/\sqrt{2}$ to each of the γ₁(θ) and γ₂(θ) galaxy observable values in the catalogue.

Finally, when measuring tomographic shear power spectra from the catalogues, we model the uncertainty due to this shape noise using Eq. (16) for the cut sky signal, with the following expression for the shape noise contribution to the E-mode power spectrum,

$$\begin{array}{c}\hfill {\stackrel{\sim }{N}}_i^{ϵ}(\ell )=\frac{{[\sigma ({ϵ}^{\mathrm{int}})/\sqrt{2}]}^2}{1/N_i^{\mathrm{P}}},\end{array}$$(31)

where $N_i^{\mathrm{P}}$ is the Poisson noise power spectrum associated with the tomographic bin i (Eq. (20)). In this way, we ensure that the shape noise only contributes to the auto-correlation power spectra.

5. Inference routine

5.1. Grid-based Gaussian likelihood

The central aim of this work is to explore the 3 × 2 pt signal for different tomographic configurations, and identify the optimum tomographic binning strategy that provides the best constraints for the time evolving dark energy equation of state. To derive parameter constraints, we perform a grid-based inference routine using a Gaussian likelihood in which we vary the (w₀, w_a) parameters in a w₀w_aCDM cosmology (see Sect. 1) while all other cosmological parameters are held at the fixed fiducial values presented in Sect. 4.2. An in-depth discussion on the validity of the Gaussian likelihood and the description of the grid-based framework adopted is given in Upham et al. (2021). Explicitly, the multivariate Gaussian likelihood is given by

$$\begin{array}{c}\hfill \mathcal{L}=\frac{1}{{(2\pi )}^{k/2}{|\mathsf{C}|}^{1/2}}exp\left[{(\mathit{D}-\mathit{M})}^{\mathrm{T}}{\mathsf{C}}^{-1}(\mathit{D}-\mathit{M})\right],\end{array}$$(32)

where D is the (mock) observed data vector of length k, M is a given theory data vector, and C is the covariance matrix for the data vector D. The 2D posterior distribution, p(Θ | D), of the model parameters, Θ, is then given by Bayes’ theorem,

$$\begin{array}{c}\hfill p(\mathbf{\Theta }|\mathit{D})\propto \pi (\mathbf{\Theta })\mathcal{L}(\mathit{D}|\mathbf{\Theta }),\end{array}$$(33)

where ℒ(D | Θ) is the likelihood and π(Θ) is the prior knowledge. For this work, we assumed a flat prior in (w₀, w_a).

5.2. Numerical covariance matrix

A key ingredient of the Bayesian method is the covariance matrix, C. For our work, we use a numerical covariance matrix in which we estimate the errors using a large number of realisations. After measuring the bandpowers of the 3 × 2 pt Pseudo-C_ℓ data vector, ${\stackrel{\sim }{C}}_b$ (Eq. (26)) for each realisation, the elements of the covariance matrix are then constructed from the scatter of the realisations with respect to the mean,

$$\begin{array}{c}\hfill {\mathrm{C}}_{bb\prime }=⟨({\stackrel{\sim }{C}}_b-{\stackrel{\sim }{C}}_{b,av})({\stackrel{\sim }{C}}_{b\prime }-{\stackrel{\sim }{C}}_{b\prime ,av})⟩,\end{array}$$(34)

where ${\stackrel{\sim }{C}}_{b,av}$ is the mean recovered bandpower in a band, b, and the angled brackets denote an average over all realisations. Then, C_bb′ are the elements of the covariance matrix, C.

While the numerical covariance matrix provides an estimate of the true errors intrinsic in our simulation method, it will also contain numerical noise due to the finite number of simulations used. By definition, a greater number of simulations leads to more suppression of this noise, until the limit of infinite realisations where the noise will average to zero and leave only the unbiased, ‘true’ signal in the covariance matrix. This unbiased estimate of the covariance will contain information on both the cosmological 3 × 2 pt signal, and the effects of Poisson noise and shape noise in the measurement.

To derive unbiased constraints using the numerical covariance matrix, we apply the Hartlap correction (Hartlap et al. 2007) to the inverse of the covariance matrix, C⁻¹,

$$\begin{array}{c}\hfill {\mathsf{C}}^{-1}\to \frac{N_S-N_D-2}{N_S-1}{\mathsf{C}}^{-1},\end{array}$$(35)

where N_S is the number of simulations used, and N_D is the size of the data vector. For the Hartlap correction to be applicable, we require the condition N_S > N_D + 2. The largest data vector considered in this work will be a full 3 × 2 pt analysis with ten tomographic bins and ten bandpowers, which gives N_D = 2100. Hence, we generated 3000 realisations of our simulated catalogues to define the numerical covariance matrix.

6. Impact of tomography on dark energy constraints

We now present constraints on the (w₀, w_a) parameters measured from the 3000 realisations of our simulated 3 × 2 pt catalogue data using the inference method described in Sect. 5, for a range of different tomographic binning configurations. For a Euclid-like survey, we take a sample of 3 × 10⁸ galaxies drawn in the range 0.3 ≤ z < 2.7 from the p(z) defined in Eq. (27) using the constant z₀ = 0.636 to yield a sample with a median redshift z_med = 0.9, matching that predicted for Euclid (Laureijs et al. 2011). We place these galaxies within an approximation of the Euclid DR1 footprint (∼2600 deg², presented in Fig. 2) to achieve the target galaxy number density, 30 gal/arcmin², on the sky. We use these survey characteristics for all results presented in this work.

For generating the 3 × 2 pt fields on the sky, we use a HEALPix grid with resolution N_side = 1024. This allows us to generate catalogues from power spectra that are evaluated up to angular scales of ℓ = 2000. Figure 2 shows the n(z) model that we use alongside examples of the generated 3 × 2 pt fields and the approximate Euclid DR1 mask that we use.

We sample galaxies to a redshift resolution dz = 0.1, which we note is slightly more pessimistic than the target accuracy for the Euclid photo-z estimation for weak lensing science, σ_z/(1 + z)≤0.05. However, due to our simulation methodology, we find that this is the limiting precision that allows for a consistent recovery of the fiducial cosmology from the 2-point statistics measured from the catalogues. In Appendix B we present our validation methodology to investigate the accuracy of our simulation method and demonstrate the ability of our pipeline to self-consistently reproduce the underlying 3 × 2 pt signal predicted by the input cosmology.

Following the procedure detailed in Sect. 3.2, we measured the 3 × 2 pt Pseudo-C_ℓ bandpowers of the data in the catalogues to derive tomographic constraints on dark energy, meaning that the theoretical model for the full-sky power spectra needs to be convolved with the mixing matrices associated with the DR1 mask. In addition, a model for the noise components (see Eq. (24)) is required for a direct comparison. We measure cosmic shear power spectra in the multipole range 100 ≤ ℓ ≤ 1500, and power spectra for the angular clustering and galaxy-galaxy lensing components in the range 100 ≤ ℓ ≤ 600, to isolate the scales where the assumption of the constant galaxy bias remains well-motivated. We evaluate all power spectra using ten bandpowers binned logarithmically in angular scale.

By fixing all other parameters in the likelihood analysis, we find that the inference process is extremely finely tuned, and even a ≪1% systematic effect in the 3 × 2 pt bandpowers (see also Appendices B.1 and B.2) leads to a biased recovery of the fiducial parameter values, (w₀, w_a) = (−1, 0). Hence, to examine the effects of tomographic binning on the true cosmology, we choose to work with the fiducial data vector but use a numerical covariance matrix, which we derive from our simulated measurements (see Sect. 5.2). This ensures that we will be probing the fluctuations about the true cosmology, but we do not expect that the size of the errorbars and the (w₀, w_a) contours will change, since we retain the realistic simulated noise that is observed in the 3 × 2 pt data from our catalogues.

To evaluate the performance of a given tomographic binning configuration, we plot the areas enclosed within the (w₀, w_a) posterior contours as a function of the number of tomographic bins used. We consider three binning cases: bins equipopulated with galaxies; bins equally spaced in redshift; and bins equally spaced in fiducial comoving distance. We consider constraints for both the ‘no-noise’ simulation and the ‘realistic’ set-up including Gaussian shape noise and Gaussian photo-z estimation uncertainty. For the ‘noisy’ simulation, we consider an intrinsic shape noise in the source galaxies that is described by a Gaussian parameterised by a constant, σ(ϵ^int) = 0.3 (e.g. Euclid Collaboration: Paykari et al. 2020), and inject this uncertainty into the shear field values assigned to each galaxy in the catalogue using the method described in Sect. 4.6.2.

To model the photo-z uncertainty, we use Gaussian distributed errors with the redshift-dependent deviation discussed in Sect. 4.6.1, where we take the constant ${\sigma }_z^{\mathrm{phot}}=0.05$ to match the target uncertainty for the Euclid analysis as defined in Laureijs et al. (2011). While we model such uncertainties in the redshift and galaxy shapes, in this section we consider only unbiased redshift and shear estimates, which we note would need to be controlled in a precision 3 × 2 pt analysis. While a comprehensive study of these issues is beyond the scope of this work, in Sect. 7 we consider the effects of contamination by catastrophic redshift errors.

For both the no-noise and noisy cases, we investigate: a 1 × 2 pt cosmic shear only; a 1 × 2 pt angular clustering only; and lastly a full 3 × 2 pt analysis. For the noisy simulation measurements, we additionally plot the 1 and 2σ posterior contours in the (w₀, w_a) plane for a selection of equipopulated binning constraints, which we find are similar in shape to the equivalent contours in the no-noise measurements.

6.1. Signal-only simulations

6.1.1. 1 × 2 pt cosmic shear only analysis, no-noise simulation

For a catalogue simulation in the absence of shape noise and photo-z uncertainty, we measure the tomographic Pseudo-C_ℓ power spectra of the cosmic shear signal only, and plot in Fig. 3 the area enclosed in the (w₀, w_a) plane by the 1 and 2σ contours as a function of the number of redshift bins used for the analysis. We present these results for the three binning choices: equipopulated (solid line); equal redshift width (dotted line); and bins equally spaced in fiducial comoving distance (dashed line).

Fig. 3. Areas enclosed by the 1σ (Blue) and 2σ contours in the (w0, wa) plane, for different numbers of redshift bins used in a tomographic analysis of the cosmic shear signal measured from 3000 realisations of our simulation. We show in circular markers joined with solid lines the areas measured for the equipopulated binning choice; in square markers joined with dashed lines the equal comoving distance bins; and in triangular markers joined with dotted lines the equal redshift width binning choice. Since the cosmic shear component alone is relatively weakly constraining, the 1 bin measurement does not yield a closed contour in (w0, wa) within the ranges of the parameter grid. Hence, the data point for this case represents a lower bound of the true value, which we represent by using a vertical arrow.

We find that, comparatively, the equal redshift width bins cover the smallest area in (w₀, w_a) for a given number of redshift bins, indicating that this binning choice would yield the optimum constraints on the dark energy equation of state. We note that for the particular case of a single bin the shear signal is only weakly constraining and does not enclose a 2σ area within the prior boundaries of the parameter space. However, increasing the analysis to include just two bins provides sufficient constraining power to yield closed contours and place constraints on the dark energy parameters.

For a two-bin analysis, the difference in the contour area between each choice is at a considerable level of ∼15–25%, but once three or four bins are considered for the analysis, the relative gain or loss between the binning choices reduces to the level of < 10%, on average. Throughout this work, we will define the percentage change in the 2σ contour area between any pair(s) of binning set-ups as the relative ‘information gain or loss’ associated with the binning choices. Following the definition of the dark energy figure of merit (FOM; Albrecht et al. 2006), whereby the FOM is taken as the inverse of the area of the 2σ (w₀, w_a) contour, our measured comparison between the areas for different binning choices similarly represents the level of improvement or degradation in the dark energy FOM value.

We note that for the equal redshift width and the equal comoving distance binning choices, the maximum number of bins that is used for the analysis is lower (four and five, respectively) than the maximum number we use for the equipopulated case. This is due to the fact that as a greater number of bins is used, there are fewer galaxies that trace the sky in a given bin, which leads to an undersampling of the observed field where pixels in the HEALPix map are left unfilled. For the case of a single realisation, this feature could be interpreted as a part of the mask itself and convolved with the mixing matrix formalism (see Sect. 3.3). However, for the purposes of this work in which we consider multiple thousands of realisations, we highlight a crucial complication in that the numbers and locations of the unfilled pixels change per realisation, and then per tomographic bin, and per binning choice within each realisation.

Hence, a new coupling matrix would need to be generated for each permutation and combination of these variations since the cosmological signal is different in each case. The simulation would be both computationally impractical, and would require a fundamentally different scientific analysis that would deserve a separate, independent investigation beyond this work. Regardless, a key point we emphasise is that as the number of redshift bins used for the tomographic analysis increases, the more similar the different binning choices become, until the limit of an infinite number of bins at which point the binning choices are identical.

From Fig. 3 we find that the area of the 1σ (w₀, w_a) contours for the equipopulated binning choice starts to converge to a consistent value by six to seven bins. We find that only for no-noise measurements, there is additional uncertainty in the contour areas for relatively large numbers of tomographic bins, due to the redshift sampling precision, but we do not expect that this impacts any of the overall conclusions presented here and in the remainder of this work.

In addition, we note that while the (w₀, w_a) constraints derived from the numerical covariance matrix are unbiased, there will be uncertainty in measurements of the contour areas due to excess numerical noise. Indeed, for a real measurement, any statistical gain in the errors on (w₀, w_a) that can be extracted by continuously introducing more bins will require increasing accuracy in estimating the covariance matrix.

6.1.2. 1 × 2 pt angular clustering only analysis, no-noise simulation

We present in this section the equivalent analysis of Sect. 6.1.1 for the clustering component of the 3 × 2 pt signal, showing in Fig. 4 the areas of (w₀, w_a) contours derived for the three different binning choices as a function of the number of bins used for a tomographic Pseudo-C_ℓ power spectrum analysis.

Fig. 4. Contour areas of the 1 and 2σ constraints on (w0, wa), measured from the tomographic angular clustering component of our no-noise simulation. We plot the contour areas measured as a function of the number of redshift bins used in a tomographic analysis, considering the equipopulated, equal comoving distance, equally spaced in redshift binning strategies. The 1 bin measurement does not sufficiently constrain either parameter within the prior volume, hence we include a vertical arrow to denote that this data point is a lower bound.

Consistent with the cosmic shear measurements in Sect. 6.1.1, we find that for a clustering-only analysis in the absence of noise, the use of two bins is required to sufficiently place firm constraints on the (w₀, w_a) parameters. Comparatively, our measurements demonstrate clearly that the equipopulated tomographic bins provide the optimum constraints on (w₀, w_a), which is in direct contrast with the results in Fig. 3 that suggest that the equipopulated binning is the worst choice for the shear-only measurement.

Moreover, in the range of three to five tomographic bins, we find that there is a degradation in the clustering-only constraints on (w₀, w_a) of ∼20% for the equal comoving distance bins, and ∼40% for the equal redshift width bins, with respect to the equipopulated choice. This difference between the binning choices is a factor of ∼2 greater than the difference measured for the no-noise shear-only analysis.

Regarding the number of tomographic bins considered, it is clear that the rate of decrease of the area enclosed by the (w₀, w_a) contours converges to a constant value by six to seven bins, indicating that the information gain starts to saturate beyond this point. For the equipopulated binning choice we calculate that by using six bins the 2σ (w₀, w_a) contour decreases by ∼80% compared to a two-bin analysis, and for every extra tomographic bin introduced beyond seven bins, there will only be a ∼1–2% further improvement in the parameter constraints with respect to the two-bin case.

6.1.3. Full 3 × 2 pt analysis, no-noise simulation

For the no-noise simulations, we lastly considered a full 3 × 2 pt analysis. We plot the results of the measurements on (w₀, w_a) in Fig. 5, demonstrating a comparison in the posterior constraints between the different binning choices. In contrast with measurements of the cosmic shear or angular clustering components alone, we find that the full 3 × 2 pt data vector has sufficient constraining power to yield closed contours in (w₀, w_a) for a single redshift bin. This effectively demonstrates the joint power that the weak lensing and galaxy clustering signals have in investigating the effect of dark energy on the growth of structure in the Universe, since they probe different redshift-dependent information of the matter field along a line of sight (see also Tutusaus et al. 2020).

Fig. 5. Areas enclosed by the 1 and 2σ contours of the posterior constraints on (w0, wa), derived from measurements of the full tomographic 3 × 2 pt signal from our ‘no-noise’ mock catalogues. We plot the contour areas as a function of the number of redshift bins used in the tomographic analysis for the three different binning choices considered in this work.

By increasing the number of tomographic bins, we find that the equipopulated binning choice provides the smallest errors, followed by the equal comoving distance bins and the equal redshift width bins – the same order of preference as for the clustering only analysis, which suggests that the tomographic behaviour of the clustering signal is the dominant component of the 3 × 2 pt data vector. The degradation in the (w₀, w_a) constraints with respect to the equipopulated choice for a given number of bins is ≲10% for the equal comoving distance and ∼10–15% for the equal redshift width choices, which is smaller in magnitude than for the clustering analysis. However, these magnitudes are consistent with the relative performances of the different binning choices in the shear-only analysis at small numbers of bins, which highlights the conclusion that the ‘optimum’ binning choice is likely to be dependent on which observable is being targeted in the 3 × 2 pt measurement.

We find that the contour areas, for the 3 × 2 pt case, converge by around six to eight bins. This is in good agreement with measurements of the shear- and clustering-only constraints (Figs. 3, 4). Beyond this number of bins, we find that the information gain in constraints on dark energy, for any and all fields of the 3 × 2 pt signal, is relatively small in general. In this regime, we expect that for a real analysis the marginal statistical gains on the dark energy parameter constraints will likely be outweighed by the increased demand for accuracy of the systematics modelling, the likelihood method, and the estimation of the covariance matrix.

The results for the noise-free ‘limiting-case’ simulations are summarised in Table 1, alongside the equivalent conclusions for a realistic noisy set-up, which we discuss in detail next.

6.2. Realistic set-up including photo-z and shape noise

We go on to present constraints on the (w₀, w_a) parameters derived from tomographic measurements of 3000 realisations of our catalogue simulation pipeline, in the presence of survey noise arising from Gaussian uncertainty in the photo-z redshift estimation of galaxies and Gaussian shape noise. The characteristics of the noise are chosen to represent the level of uncertainty that is expected to be achieved in the Euclid DR1 survey (see further details in Sect. 4.6). As presented in the no-noise analysis in Sect. 6.1 we examine each of the cosmic shear, angular clustering and full 3 × 2 pt signals, and measure the 1 and 2σ constraints across different binning choices as a function of the number of tomographic bins used for the analysis.

6.2.1. 1 × 2 pt cosmic shear only analysis, noisy simulation

In Fig. 6 we plot the 2D posterior constraints on the time evolving dark energy parameters (w₀, w_a) from noisy tomographic cosmic shear measurements using one, two, three, five, and ten equipopulated bins. We find that in the presence of realistic noise, the cosmic shear is relatively weakly constraining and the constraints on both parameters do not change significantly when increasing the number of tomographic bins up to ten. Additionally, we note that within the prior boundaries of the grid-based likelihood approach, only the 1σ contour yields an enclosed area, while the 2σ contour coverage in the parameter space is extended such that it is computationally intractable to constrain. However, since we are sampling and fully characterising the parameter volume for the 1σ contour for all binning numbers considered, we do not believe that the true behaviour of the 2σ contours would yield different results or conclusions.

Fig. 6. 1 and 2σ constraints on the (w0, wa) parameters derived from measurements of the tomographic cosmic shear signal in our noisy simulation that includes contributions from Gaussian-distributed photo-z errors and Gaussian-distributed shape noise. We present the constraints for different numbers of equipopulated bins but find that the signal has relatively weak constraining power that does not significantly improve as the number of bins increases. Hence, the contours corresponding to each different number of bins lie on top of each other.

In Fig. 7, we measured the area of the 1σ contours in the (w₀, w_a) plane for each binning choice. Here, we find that tomographic bins equally spaced in fiducial comoving distance produce the best constraints, followed by the equal redshift width bins and finally the equipopulated bins. This represents a different order of preference compared to the no-noise set of cosmic shear measurements (Fig. 3), for which equal redshift bins give the optimum constraints. However, the relative difference between the binning choices is at a minimal level of < 5% and indicates that the gain that can be leveraged by choosing a specific binning choice for the cosmic shear is relatively small, regardless of the presence of realistic noise.

Fig. 7. The 1σ contour areas of constraints in (w0, wa) from measurements of the tomographic cosmic shear signal in our noisy simulation. We consider three different binning strategies and find that the 1 bin measurement does not have sufficient constraining power to generate closed contours within the prior ranges of the parameter grid. Hence, we plot the unbounded area measured within the prior ranges, and indicate that this is a lower bound of the true value by using a vertical arrow.

In line with the 2D contour plot, we also find only a small reduction in the (w₀, w_a) area going from two to ten equipopulated bins of ∼10%, which is a considerably lower gain than in the no-noise simulation, which demonstrates a ≳99% gain between the same number of equipopulated bins. The improvement in the (w₀, w_a) constraints shows a step-like behaviour between six to ten bins, which we attribute to the finite pixel size of the parameter grid. As the difference in the enclosed (w₀, w_a) contour areas is minimal in this range, we conclude that the convergence in the dark energy constraints for the noisy shear measurements is achieved at roughly seven to eight tomographic bins.

The simulated uncertainties from the shape noise in an individual galaxy’s shear estimate, and the photo-z uncertainties, which result in overlap in the redshift distributions of different tomographic bins, will both lead to a decrease in the constraining power. It is not necessarily straightforward to determine which one contributes more to the overall degradation here. One could explore simulations in which there is shape noise included but no photo-z uncertainty, and vice versa, and conduct a similar analysis to investigate how the (w₀, w_a) constraints behave in each scenario. In particular, the former case would be more analogous to a spectroscopic weak lensing survey, in which the redshifts are precisely known and tomographic bins do not overlap. However, the fundamental survey characteristics such as the redshift range, survey footprint, and number density on the sky would change, and fundamentally require different tuning in our simulation parameters. It would then be difficult to make a direct comparison with the results presented here.

Additionally, the reverse scenario, in which photo-z uncertainties are included but shape noise is not present in the simulation, is unachievable in a real observation since observed galaxies have an intrinsic shape. As such, we consider that both of these intermediate simulation set-ups are beyond the scope of this work, but would nonetheless be interesting to explore in a future study.

6.2.2. 1 × 2 pt angular clustering only analysis, noisy simulation

In this section, we explore the tomographic constraints made on (w₀, w_a) by the clustering component of the 3 × 2 pt signal in the noisy simulation. In Fig. 8 we show the 2D contours measured from varying the number of redshift bins used in the equipopulated binning choice, and in Fig. 9 we plot the areas enclosed by the 1 and 2σ contours in the (w₀, w_a) plane as a function of the number of tomographic bins, for all three binning choices considered.

Fig. 8. 1 and 2σ posterior constraints on (w0, wa) measured from the tomographic angular clustering signal in our mock catalogues which include simulated noise from Gaussian photo-z uncertainties and shape noise. We show the contours for different numbers of equipopulated bins used for the tomographic analysis.

Fig. 9. Contour areas of the 1 and 2σ constraints on (w0, wa) for measurements of the tomographic angular clustering signal in our noisy simulation. We show the contour areas measured as a function of the number of redshift bins used, for the three different tomographic binning choices considered in this work.

By inspecting Fig. 8, we find in particular that the shapes of the contours in the (w₀, w_a) plane for the clustering analysis are drastically different than the equivalent contour shapes derived from the cosmic shear analysis in Fig. 6. This property clearly highlights the power of the 3 × 2 pt combined probe to break degeneracies in the dark energy parameters and provide highly precise measurements of the equation of state.

Collectively, the behaviour of (w₀, w_a) posterior areas measured as a function of the number of tomographic redshift bins is remarkably consistent with the no-noise simulation. In particular, the relative order of preference of the binning choices is identical across the noisy and no-noise analyses, whereby the equipopulated bins provide the optimum constraints, followed by the equal comoving distance bins and finally the equal redshift bins.

Moreover, we find that simply using two bins is sufficient, for all binning choices in the noisy simulations, to make measurements of (w₀, w_a) which indicates that the presence of the photo-z uncertainty and shape noise does not significantly degrade the constraining power of the cosmological signal. In an absolute sense, we find that compared to the no-noise measurements, the contour areas increase by ∼5–10% on average for the equipopulated case, and ∼5% for both the equal comoving distance and equal redshift binning choices. We note that for the equipopulated choice, the increase is only due to the degradation of the cosmological signal by the photo-z uncertainties, as the shape noise does not affect the angular galaxy positions and the Poisson noise is the same in both the no-noise and noisy simulations.

Most importantly, we highlight that, with and without the realistic noise included in the simulation, the improvement in the (w₀, w_a) constraints starts to saturate at six or seven tomographic bins, which is demonstrated by the convergence and flattening of the distributions corresponding to the 1 and 2σ measurements in Fig. 9.

Consequently, we conclude that for an analysis of the angular clustering only, the optimum tomographic binning recommended by these results is the equipopulated choice. Here, the choice of about six or seven bins is likely to be sufficient to achieve the optimum performance in the noisy case.

6.2.3. Full 3 × 2 pt analysis, noisy simulation

Finally, we consider the effects of the tomographic binning choices on the (w₀, w_a) parameters for a full 3 × 2 pt analysis in the presence of noise, the target experimental and analytical set-up with which to compare to the aims of the Euclid DR1 survey. In Fig. 10, we first plot the 1 and 2σ (w₀, w_a) contours for tomographic measurements of the 3 × 2 pt signal using one, two, three, five, and ten equipopulated bins, and in Fig. 11 we compare the areas enclosed by such contours across different binning strategies.

Fig. 10. 1 and 2σ constraints on (w0, wa) from tomographic Pseudo-Cℓ measurements of the full 3 × 2 pt signal in our noisy simulation. We show constraints measured for different numbers of equipopulated bins used for the tomographic analysis. The cross marks the fiducial cosmology (w0, wa) = (−1, 0) used for the simulation.

Fig. 11. Contour areas enclosed by the 1 and 2σ constraints on the dark energy (w0, wa) parameters, measured for the full 3 × 2 pt signal in the presence of Gaussian shape noise and photo-z uncertainty. We vary the number of redshift bins used for each of the equipopulated, equally spaced in fiducial comoving distance, and equal redshift width binning strategies.

While a single tomographic bin used for analysis is only weakly constraining for measurements of the noisy cosmic shear and clustering signals only (Figs. 6, 8), we find that by considering the joint information contained in the 3 × 2 pt signal, the single bin can fully constrain the (w₀, w_a) parameters and identify a time evolving nature of the dark energy equation of state. Consequently, this result suggests that the single-bin analysis could be used as an important first-order measurement to investigate whether the photo-z uncertainties and shape noise are correctly accounted for in the modelling.

With respect to the different binning choices in Fig. 11, the noisy 3 × 2 pt analysis is similarly consistent with the no-noise results and demonstrates that the equipopulated bins yield measurements on (w₀, w_a) that have the least error, followed by the equal comoving distance bins and then the equal redshift width bins. Since the same order of preference is found in the clustering-only investigation in Sect. 6.2.2, and we find that the different binning choices only make per cent level changes to the (w₀, w_a) contour areas for the shear only analysis in Sect. 6.2.1, we conclude that the performance of the different binning choices in the full 3 × 2 pt signal is likely to be dominated by the response of the angular clustering component to the tomographic binning.

For the noisy measurements, the increase in the (w₀, w_a) areas between each binning choice in successive order of preference is more consistent with ∼10–15% at low numbers of bins (two-four), compared to ∼5–10% for the no-noise analysis, which demonstrates that in this range the specific binning choice becomes more important in deriving optimum constraints on dark energy with increasing realism in the simulation.

However, we find that beyond six or seven tomographic redshift bins used for the analysis, the (w₀, w_a) areas corresponding to the equipopulated constraints start to plateau, demonstrating the same saturation in the information gain present in the noisy angular clustering measurements in Fig. 11. In Table 1, we summarise our findings on the behaviour of the (w₀, w_a) contours for the different tomographic measurements of the noisy cosmic shear-only, angular clustering-only, and full 3 × 2 pt signals. Also reported in Table 1 are the limiting case no-noise results, for all 3 × 2 pt probes.

Empirically, we find that the areas of the 1 and 2σ contours in the (w₀, w_a) plane, measured as a function of the number of equipopulated tomographic redshift bins used, are well described by a 1/x function of the form y = a/x + b. We use the SciPy (Virtanen et al. 2020) curve_fit routine to derive (a, b) = (0.00401, 0.00094),(0.0105, 0.00263) for each of the 1 and 2σ distributions, respectively.

We present these functions fit to the areas of the equipopulated constraints in Fig. 12, which demonstrates empirically how the rate of improvement in the (w₀, w_a) areas starts to decrease as the number of bins used in the tomographic analysis increases. Using the fitting functions as a baseline, we find that compared to a one-bin analysis, the areas of the (w₀, w_a) contours decrease by 69% when using seven bins, and 72% when using ten bins, which indicates only a further 3% gain in information between the use of seven and ten bins.

Fig. 12. Areas of the 1 and 2σ constraints on (w0, wa), measured from our noisy catalogue simulations, plotted in circular markers as a function of the number of equipopulated bins used for the tomographic analysis. Black solid lines show the plot of the best fit functions of the form y = a/x + b that we fit to these data sets, evaluated using the SciPy curve_fit routine.

By extending the fitting functions to explore the theoretical behaviour of tomographic analyses that use more than ten bins (which we find we cannot practically measure with our method due to the presence and randomness of untraced pixels on the sky), we can estimate the area of the (w₀, w_a) contours for the particular case of 13 equipopulated bins. We note this use of 13 bins was found in Euclid Collaboration: Pocino et al. (2021) to improve the dark energy FOM by 15% for a 2 × 2 pt analysis of galaxy clustering and galaxy-galaxy lensing, based on photo-z studies using the Euclid Collaboration Flagship Simulation (Euclid Collaboration: Castander et al. 2025).

Using our extended fitting functions in Fig. 12, we find that the (w₀, w_a) contour areas for 13 bins decrease by 17% and 7% compared to analyses using seven and ten tomographic bins, respectively. Additionally, by measuring the total improvement on the (w₀, w_a) constraints with respect to the single-bin analysis, we find that this corresponds to only a further 5% gain in the precision on the parameters when increasing from seven to 13 bins in the measurement, and a 2% level gain when going from ten to 13 bins.

Overall, these results indicate that increasing the number of tomographic bins beyond ten leads to improvements in the (w₀, w_a) constraints that are less substantial than those found by Euclid Collaboration: Pocino et al. (2021). However, we note that each of these results has been determined using different methods and with different assumed survey characteristics.

In particular, the latter study considers a 2 × 2 pt analysis and focuses on the configuration and cosmological information that is expected to be obtained from the complete Euclid survey, rather than just the first data release (DR1) that our analysis has focussed on. Consequently, the underlying n(z), spatial density of galaxies on the sky, and scale cuts taken for the 3 × 2 pt signal considered in Euclid Collaboration: Pocino et al. (2021) will be different than those considered in this study. Additionally, while we have considered an analytical method to simulate galaxy redshifts and their uncertainties, Euclid Collaboration: Pocino et al. (2021) adopted machine learning methods to estimate the photometric redshifts.

We note for our analysis that there are still errors associated with the data points due to the fact we can only work with a finite number of realisations, which in an ideal case would be included in the derivation of the fitting function and its parameters. As such, the relative statistical gains and comparisons we present here should be considered only as a guide and not absolute results. In principle, a quantification of the errors could be achieved for example by repeating the entire analysis to measure the set of (w₀, w_a) contour areas over multiple batches of 3000 realisations, and then examining the deviation across all batches. However, this is unfortunately computationally expensive and impractical to achieve in this study.

Nonetheless, we assert that the errors will be equivalent in magnitude for the different binning choices considered, and hence conclude that for a full 3 × 2 pt tomographic analysis of our simulated Euclid DR1-like catalogues, in the presence of realistic noise, the optimum binning choice is the equipopulated case.

While there perhaps could still be merit in using ten tomographic bins to extract a few more per cent precision in (w₀, w_a), we note that the use of more bins would increase the computational expense. Moreover, it is likely that any potential benefits are only realisable, within the framework of a theoretical study, such as the one presented here, where one has full control and knowledge of the uncertainties that are introduced.

Comparatively, for a real Stage IV analysis, it is likely that this extra information probed by increasing the number of bins beyond about seven will be difficult to leverage, especially in light of the additional model and inference complexity that we do not consider here. Such effects include the shear bias in the shape measurements of galaxies (e.g. Hirata & Seljak 2003, Jansen et al. 2024, Euclid Collaboration: Congedo et al. 2024), the effects of magnification in the galaxy sample and further high order terms (e.g. Duncan et al. 2022, Deshpande et al. 2020, Euclid Collaboration: Deshpande et al. 2024), or the effects of marginalisation over nuisance parameters, such as the galaxy bias, or other cosmological parameters within a given prior range.

We note additionally the impact of intrinsic alignments on 3 × 2 pt cosmology (e.g. Troxel & Ishak 2015), which we do not model in this study. It has been proposed that the systematic effects of intrinsic alignments can be quantified by the use of tomography (King & Schneider 2002), and we highlight that such an investigation, combined with the effects of further systematics not considered in this study, can be examined in later work.

Indeed, Euclid Collaboration: Pocino et al. (2021) consider marginalising over both intrinsic alignments and a constant galaxy bias in each tomographic bin for a 2 × 2 pt forecast of photometric galaxy clustering and galaxy-galaxy lensing with Euclid. By targeting the dark energy figure of merit, they report that the marginalisation procedure motivates the use of fewer tomographic bins compared to conditional constraints, since less cosmological information can be extracted.

In addition, they found that for the marginalised 2 × 2 pt constraints, there is negligible difference between the equal redshift width and equipopulated binning choices. Comparatively, there is a ∼30% difference between these binning choices for conditional constraints for a clustering-only analysis using 13 tomographic bins. This result indicates that the choice of tomographic binning strategy for each 3 × 2 pt probe may have a complex dependence on different systematics that may be incorporated in the analysis. A detailed study covering different models for a range of systematics, in which their effect on tomography is considered both individually and collectively, is a study that is beyond the scope of this paper, but will be an important investigation for future work.

Lastly, we note that the true errors in the photo-z estimates, the shear estimates and/or the measured power spectra may not be purely Gaussian (e.g. Takada & Jain 2009). We explore an example of this scenario in the following section (Sect. 7) where we consider the effects of catastrophic photo-z uncertainties in the estimation of the galaxy redshifts and the subsequent impact on the inference on the dark energy parameters.

7. Catastrophic photo-z uncertainties

In this section, we demonstrate the flexibility of our simulation method to incorporate further astrophysical or survey-like effects by injecting catastrophic photo-z uncertainties into our catalogue simulation and then propagating the resulting biased measurements into cosmological constraints on (w₀, w_a). Photo-z estimation codes can misidentify features in the spectral energy distribution (SED) of a galaxy, leading to a measured redshift that significantly fails to recover the true redshift of the galaxy. For the purposes of the 3 × 2 pt signal, the misestimated galaxy redshift will lead to an incorrect interpretation of the structure along the line of sight that has induced the shear in the galaxy, which will limit our ability to make accurate measurements of dark energy and its possible time-evolving behaviour.

To achieve its science goals, Euclid will target a 5% level of contamination by catastrophic photo-z errors (Laureijs et al. 2011) in the real analysis. Hence, we introduced the effects of this systematic into our simulation by randomly selecting 5% of our galaxy sample in the noisy set-up and distorting their estimated redshift measurement in the catalogue to a catastrophic estimation using Eq. (30). The relation between the true redshift and the catastrophic redshift, z_cata, modelled here is governed by the ratio between the rest-frame wavelength of a break feature in the SED and the wavelength that is predicted by the photo-z estimation technique. As a first-order investigation, we considered catastrophic errors that result from the rest-frame Lyman-α line being confused with the Balmer and D4000 breaks, i.e. taking the pairs [λ_{break − rf}, λ_{break − cata}]=[Ly-α, Balmer-break],  [Ly-α, D4000-break] and splitting the 5% contamination equally between the two pairs of misidentified features. We finally model the resulting catastrophic photo-z estimate as a Gaussian-distributed value with mean z_cata and standard deviation σ_z, cata = 0.1.

Following this procedure, we measured the tomographic Pseudo-C_ℓ power spectra from these biased mock catalogues using the standard technique adopted in this work, as described in Sect. 3.2. The presence of the catastrophic photo-z errors will induce a change, ΔC_ℓ, in the measured power spectra compared to the standard noisy case. To determine the impact of this bias on the dark energy (w₀, w_a) constraints, we: 1) added the measured ΔC_ℓ onto the true, theoretical 3 × 2 pt data vector of the fiducial cosmology; and 2) performed the grid-based Gaussian likelihood analysis on this distorted 3 × 2 pt signal using an underlying tomographic redshift distribution n(z) and a numerical covariance matrix that have both been derived from the simulations that include Gaussian noise. As such, the final constraints on (w₀, w_a) that we measure will be directly comparable to the Gaussian noise-like results of Sect. 6.2.3. Any difference will be due to the 5% contamination by catastrophic photo-zs, where we do not retain knowledge of which individual galaxies are affected, and the use of a covariance matrix that does not include sufficient information on the catastrophic photo-zs. As such, we emphasise that we do not attempt to mitigate against the presence of the catastrophic outliers to derive our dark energy constraints, which we expect to be pessimistic. Comparatively, a real analysis could attempt to correct for these uncertainties in the assumed n(z) distribution, or through nuisance parameters to marginalise over in the inference.

In Fig. 13, we show the 1, 2, and 3σ constraints on (w₀, w_a) for three different tomographic analyses of the full 3 × 2 pt signal: a single redshift bin; five equipopulated bins, and ten equipopulated bins, comparing in blue the results from the Gaussian noise-like simulations in Sect. 6.2.3 and in red the results from the catastrophic photo-z analysis averaged over 3000 realisations.

Fig. 13. Constraints on the dark energy (w0, wa) parameters measured from: simulations that include only Gaussian shape noise and Gaussian photo-z uncertainties (blue) and simulations that include Gaussian shape noise, Gaussian photo-z uncertainties, and a 5% level of catastrophic photo-z errors (outliers) that result from the Lyman-α line in the galaxy SED being misidentified as the Balmer and D4000 break features (red). From left to right, we plot the constraints derived from tomographic analyses of the full 3 × 2 pt signal using a single redshift bin, five equipopulated bins, and ten equipopulated bins.

We find that the measurements across the three different tomographic binning choices all consistently yield the constraints (w₀ > −1,  w_a < 0). These parameters are strongly biased from the underlying ΛCDM cosmology and demonstrate that the presence of a 5% level of contamination from catastrophic photo-zs, if entirely unaccounted for in the modelling, would lead us to infer significantly incorrect conclusions on a time-evolving behaviour of the dark energy equation of state. Indeed, such effects of catastrophic photo-z uncertainties for weak lensing surveys have been extensively studied in the literature, e.g. Hearin et al. (2010), Ma et al. (2006), and Sun et al. (2009), for which our results act to support and consolidate.

Since the constraining power of the 3 × 2 pt signal increases with the number of redshift bins used for the tomographic analysis, the statistical disagreement with the true underlying cosmology also increases with the number of tomographic bins. Explicitly, for a single-bin analysis we find that if the catastrophic outliers are not modelled in the 3 × 2 pt data vector, the power spectrum measurements are in tension with ΛCDM at ∼2σ, while for five and ten bins, the tension is at a ≳3σ level. It is clear that the dark energy (w₀, w_a) parameters are highly sensitive to the catastrophic photo-z uncertainties. Since we do not attempt to correct for these errors to derive the constraints presented here, these results aim to quantify the extent to which these uncertainties would need to be minimised or mitigated in the real analysis to determine the true nature of dark energy.

In addition to the biased recovery of the true cosmology, we find that the presence of the catastrophic photo-z errors leads to unusual structures in the 2D posterior. For our standard no-noise and noisy simulations (for which the photo-z uncertainties are purely Gaussian), we expect that the 3 × 2 pt measurements will generate smooth, elliptical contours, as demonstrated in Fig. 10. In contrast, for the five- and ten-bin analyses we see that the contours are less elliptical in shape and have sharper features in the catastrophic photo-z analysis. We attribute this behaviour to the ΔC_ℓ change in the 3 × 2 pt power spectra that the catastrophic redshift errors induce. The C_ℓ bias is non-trivial and random across all redshifts, and hence the space of cosmological models in (w₀, w_a) that the resulting power spectra are compatible with is complex and possibly discontinuous.

We conclude that the introduction of the 5% contamination by catastrophic photo-z uncertainties can severely compromise the attempt to measure accurate cosmological constraints from the 3 × 2 pt signal if not adequately mitigated. In particular, the constraints yield a best fit model that is in tension with ΛCDM at ≳3σ for analyses using 5 or more tomographic redshift bins, and the posterior distributions in (w₀, w_a) demonstrate sharp features.

Additionally, we note that catastrophic photo-z uncertainties may be caused by further effects in addition to the confusion between the Lyman-α line, and the Balmer and D4000 breaks that we have considered throughout this work. For example, ‘uniform’ catastrophic outliers may be caused when the light from a galaxy is contaminated by light from a nearby source on the sky which is at a different redshift (see e.g. Hearin et al. 2010). A full investigation into the nature of the bias induced by further causes of catastrophic photo-z errors is beyond the scope of this work. However, we emphasise that our simulation framework has the flexibility to incorporate a range of different distributions for the catastrophic photo-z uncertainty that could be examined in a future study.

These results emphasise that an accurate modelling of the catastrophic photo-z errors would need to be accounted for in the theoretical model of the cosmological signal, the covariance matrix, and the likelihood distribution. Further methods to mitigate against catastrophic uncertainties include marginalisation over a nuisance parameter (Bernstein & Huterer 2010) capturing catastrophic errors; or selectively removing galaxies from the redshift distribution (Hearin et al. 2010), for example by monitoring the posterior likelihood distribution of galaxy redshift estimates (Nishizawa et al. 2010).

8. Discussion and conclusions

Investigations attempting to determine an ‘optimum’ binning configuration for a tomographic 3 × 2 pt analysis, mostly using machine learning or numerical methods (Taylor et al. 2018a, Kitching et al. 2019, Sipp et al. 2021, Euclid Collaboration: Pocino et al. 2021, Zuntz et al. 2021, see Sect. 2), have returned inconsistent conclusions on both the binning choice and number of bins that ought to be targeted. In this work, we present an alternative method to explore the tomographic binning with the aim to mimic a real survey and subsequent cosmological analysis in a robust and comprehensive manner. We have developed a simulation pipeline that generates multiple realisations of mock galaxy catalogues for a Stage IV-like survey, which are constructed by generating a set of correlated 3 × 2 pt 2D maps on the sky, sampling a target redshift range, to model the full 3D cosmological information in the Universe.

We simulated 3000 realisations of our mock catalogues and measured the 3 × 2 pt Pseudo-C_ℓ power spectra, evaluated over an approximation of the Euclid DR1 footprint, for a range of tomographic configurations. We propagated these measurements into a grid-based Gaussian likelihood analysis to place constraints on the dark energy equation of state parameters, (w₀, w_a), and derive the optimum binning choice by comparing the size of the contours in (w₀, w_a) across all tomographic binning configurations considered.

For a limiting case ‘no-noise’ simulation, where we did not consider any effects of redshift uncertainties due to photo-z estimation techniques or shape noise associated with the observed shear of a galaxy, we found that the optimum binning choice for a full 3 × 2 pt analysis is given by redshift bins that are equally populated with galaxies. For the angular galaxy clustering signal, the optimum binning choice is also the equipopulated case, which is likely to be the dominant component of the 3 × 2 pt data vector. However, for an analysis of the cosmic shear signal only, bins that are equally spaced in redshift yield the best constraints on (w₀, w_a), which indicates that the optimum binning choice is dependent on which observable is targeted.

For a more realistic simulation in which we include Gaussian shape noise and Gaussian photo-z uncertainties in the mock catalogues, we find that the equipopulated bins remain the best choice for measuring (w₀, w_a) for the full 3 × 2 pt and angular clustering only analyses. However, the cosmic shear signal itself gives the best constraints when bins that are equally spaced in fiducial comoving distance are used, where the fiducial cosmology we work with is the standard ΛCDM.

We note that for a different underlying cosmology, the redshift bin boundaries for the equal comoving distance binning choice would change. However, we emphasise that as greater numbers of tomographic bins are used for the analysis, the derived constraints on (w₀, w_a) for each binning choice would converge together regardless of the fiducial cosmology. Hence, we do not expect that our results would change significantly in this range.

Collectively, we find that the relative gain or loss between the different binning choices is least significant for the cosmic shear component, in particular, at the ∼5% level when noise is included in the simulation. This result suggests that the cosmic shear signal is less sensitive to the tomographic binning choice and there may not be significant benefit in targeting a given binning choice in the real analysis.

In comparison, we find that the difference in size of the (w₀, w_a) contour areas, between the different binning choices in successive order of preference, is ∼10% on average for the full 3 × 2 pt and clustering-only measurements in the no-noise set-up. For the realistic noise simulation, this difference is slightly more amplified at ∼10–15%. We conclude for these signals that the gain is considerable and it provides clear motivation to choose the equipopulated binning for a real Stage IV analysis.

With respect to the number of redshift bins chosen for tomography, we expect that the rate of information gain will decrease as the number of bins increases, until a point at which the area of the (w₀, w_a) contours will converge to a consistent value that does not continue to improve when arbitrarily adding more bins for the tomographic analysis. For each of the cosmic shear and clustering components individually, we find that this convergence is reached at around six to eight bins for the 1 and 2σ contours across both the no-noise and noisy simulations. For each of these signals, the minimum number of bins required to fully constrain the dark energy parameters is two, and using this measurement as a baseline, we expect that for both components in the noisy set-up, there is only a further 𝒪(1%) level improvement in the (w₀, w_a) constraints when increasing the number of tomographic bins successively from six to ten.

In terms of the full 3 × 2 pt analysis, we find a convergence in the (w₀, w_a) areas at the use of seven or eight tomographic redshift bins for the no-noise measurements. For the noisy simulation, the behaviour is similar and empirically. We find that the rate at which the (w₀, w_a) contour areas decrease (as a function of the number of equipopulated tomographic redshift bins) is well described by a 1/x function of the form y = a/x + b, where we derived (a, b) = (0.00401, 0.00094),(0.0105, 0.00263) for the 1 and 2σ trends, respectively. Using this function to theoretically predict the contour size for the 13 bin analysis that has been proposed for Euclid (e.g. Euclid Collaboration: Pocino et al. 2021), we measured the total gain in information with respect to the single redshift bin constraints, finding that there is only a 5% further improvement in measurements of (w₀, w_a) when going from seven to 13 bins.

We note that for higher numbers of redshift bins, the 3 × 2 pt data vector becomes larger and the information gain on (w₀, w_a) becomes increasingly limited by the demand for accuracy in the covariance matrix. This is manifested as an increase in uncertainty due to excess random noise in the numerical covariance matrix used in this work; however, the trade-off would also be important to consider in a real analysis, particularly where an analytical covariance may be used.

Furthermore, there will be additional uncertainties present in a real survey, such as from the shear calibration, redshift-dependent astrophysical systematics (e.g. intrinsic alignments), and the nature of non-Gaussianities in the likelihood or the errors. The accuracy to which these are accounted for will then limit whether the information gain can actually be realised and, hence, the extent to which arbitrarily increasing the number of redshift bins would be preferred. In principle, the flexibility of our simulation approach offers the capability to include the effects of these additional uncertainties in our mock catalogues and then propagate the resulting power spectra into cosmological constraints on (w₀, w_a). We note that such additional effects may lead to a convergence in the (w₀, w_a) areas at an even lower number of redshift bins or yield a different order of preference for the binning choice; however such a study is beyond the scope of this paper.

Nevertheless, for the 3 × 2 pt and angular clustering analyses such additional effects would have to be considerable since the relative difference between the different binning choices is reasonably significant at ∼10–15%. Hence, we reassert that the best performing binning strategy, namely, the equipopulated case, is likely to be a good choice for these observables in the DR1 analysis for Euclid. Furthermore, in the limit that the observation is dominated by shape noise and photo-z uncertainty, we believe that any given binning choice is unlikely to result in significant degradation in the constraints achievable using the cosmic shear signal on (w₀, w_a).

Additionally, we note that throughout this work we have considered a fixed scale range of 100 ≤ ℓ ≤ 600 for angular clustering and galaxy-galaxy lensing and 100 ≤ ℓ ≤ 1500 for cosmic shear when deriving dark energy constraints from tomographic Pseudo-C_ℓ power spectra. For the real analysis of the Euclid sample, we expect that the upper limit is likely to be a conservative estimate for both probes. The effect of changing these angular scale cuts on conclusions of the optimum tomographic binning strategy could represent an informative extension of this study. However, such an investigation would require a careful consideration of a range of further modelling and analytical choices to adopt for the 3 × 2 pt measurements. In particular, as we probe increasingly smaller scales by including larger ℓ modes, the assumptions of the constant galaxy bias or even the Gaussianity of the underlying 3 × 2 pt fields are likely to break down.

Furthermore, there is a considerable range of freedom in how a scale cut could be applied to each component of the 3 × 2 pt signal and it is not trivial to predict the degeneracy in the behaviour of the (w₀, w_a) constraints when different combinations of the probes have different cuts applied. We further note that the practical and computational limitations of the simulation-based method will be exacerbated as smaller scales are introduced. Notably, the smaller scales will increase the presence of the unphysical pixels realised on the sky (those with δ_g < −1 which mask out observable regions of the mock survey, see also Sect. 4.3); furthermore, it may increase the shot noise in the numerical covariance matrix if the 3 × 2 pt power spectra need to be resolved to a greater number of angular bandpowers to capture the small scale information. Hence, we have considered a thorough examination into the impact of the chosen angular scales on tomography as an important study, but this is beyond the scope of this work. However, as a first-order investigation, we have explored extending the assumption of the linear galaxy bias for the clustering and galaxy-galaxy lensing signals over the range 100 ≤ ℓ ≤ 1500 to match the cosmic shear measurements. Ultimately, we found that this does not change any of our overall conclusions.

Lastly, we have considered the effects of catastrophic photo-z uncertainties for our dark energy constraints, where the errors associated with the redshift estimation of galaxies is no longer Gaussian. We injected a 5% level of contamination by catastrophic photo-z errors, randomly distributed across all redshifts, caused by the Lyman-α line being confused with the Balmer and D4000 breaks.

If such a contamination is unaccounted for in the modelling, we find that a full 3 × 2 pt analysis using a single tomographic bin would return measurements of (w₀, w_a) that are in tension with the fiducial ΛCDM cosmology at ∼2σ. For a five- or ten-bin analysis with equipopulated bins, we find that the (w₀, w_a) constraints are in considerable tension with ΛCDM at ≳3σ and there are sharp features in the posterior distribution. This result reaffirms the necessity to minimise the presence of these catastrophic redshift uncertainties, and to develop techniques that model their effects in the cosmological signal, the covariance, and even the nature of the likelihood distribution.

Acknowledgments

We thank colleagues in the Euclid Consortium for useful comments and feedback. We thank the anonymous referee for their comments. We also thank Duncan Austin, Tom Harvey and Erik Rosenberg (JBCA) for their helpful discussions around modelling approaches for photo-z outliers and the numerical covariance matrix. JHWW acknowledges support in the form of a PhD studentship from the UK Science and Technology Facilities Council (STFC) and the University of Manchester Centre for Doctoral Training. MLB and CAJD acknowledge funding from the STFC (grant number ST/X001229/1). The Euclid Consortium acknowledges the European Space Agency and a number of agencies and institutes that have supported the development of Euclid, in particular the Agenzia Spaziale Italiana, the Austrian Forschungsförderungsgesellschaft funded through BMIMI, the Belgian Science Policy, the Canadian Euclid Consortium, the Deutsches Zentrum für Luft- und Raumfahrt, the DTU Space and the Niels Bohr Institute in Denmark, the French Centre National d’Etudes Spatiales, the Fundação para a Ciência e a Tecnologia, the Hungarian Academy of Sciences, the Ministerio de Ciencia, Innovación y Universidades, the National Aeronautics and Space Administration, the National Astronomical Observatory of Japan, the Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the Research Council of Finland, the Romanian Space Agency, the Swiss Space Office (SSO) at the State Secretariat for Education, Research, and Innovation (SERI), and the United Kingdom Space Agency. A complete and detailed list is available on the Euclid web site (www.euclid-ec.org/consortium/community/).

References

Appendix A: 3×2pt pseudo-C_ℓ estimators

In this section, we summarise the expressions for the tomographic Pseudo-C_ℓ power spectra on the cut sky for a 3×2pt analysis (as presented in Sect. 3.3). For a spin-0 clustering field, the partial sky overdensity, ${\stackrel{\sim }{\delta }}_g(\mathit{\theta })$, is related to the full-sky field by the W(θ) window function,

$$\begin{array}{c}\hfill {\stackrel{\sim }{\delta }}_g(\mathit{\theta })=W(\mathit{\theta }){\delta }_g(\mathit{\theta }),\end{array}$$(A.1)

and the spherical harmonic coefficients of the cut-sky field, ${\stackrel{\sim }{d}}_{\ell ,m}$, are related to the full-sky coefficients according to

$$\begin{array}{c}\hfill {\stackrel{\sim }{d}}_{\ell ,m}={\displaystyle \sum _{{\ell }^{\prime },m^{\prime }}}{}_0{W}_{\ell {\ell }^{\prime }}^{m{m}^{\prime }}{d}_{{\ell }^{\prime },m^{\prime }},\end{array}$$(A.2)

where ₀W_ℓℓ′^mm′ is the harmonic space window function characterised by the mask. For a general spin-s field, the spin-weighted harmonic window function is given in terms of the spin-weighted _sY_ℓ, m(θ) polynomials and integrating over the full sphere:

$$\begin{array}{c}\hfill {}_s{W}_{\ell {\ell }^{\mathit{prime}}}^{m{m}^{\mathit{prime}}}=\int \mathrm{d}\mathit{\theta }{}_s{Y}_{{\ell }^{\mathit{prime}},m^{\mathit{prime}}}(\mathit{\theta })W(\mathit{\theta }){}_s{Y}_{\ell ,m}^{\ast }(\mathit{\theta }).\end{array}$$(A.3)

This framework can be extended to relate the observed cut-sky shear field components, ${\stackrel{\sim }{\gamma }}_{1,2}(\mathit{\theta })$, to their full-sky counterparts, γ_1, 2(θ), via a spin-2 harmonic decomposition,

$$\begin{array}{cc}\hfill {\stackrel{\sim }{\gamma }}_1(\mathit{\theta })\pm \mathrm{i}{\stackrel{\sim }{\gamma }}_2(\mathit{\theta })& =W(\mathit{\theta })[{\gamma }_1(\mathit{\theta })\pm \mathrm{i}{\gamma }_2(\mathit{\theta })]\hfill \\ \hfill & ={\displaystyle \sum _{\ell ,m}}{({\stackrel{\sim }{E}}_{\ell ,m}\pm \mathrm{i}{\stackrel{\sim }{B}}_{\ell ,m})}_{\pm 2}{Y}_{\ell ,m}(\mathit{\theta }).\hfill \end{array}$$(A.4)

Following Lewis et al. (2001), the cut sky harmonic coefficients $({\stackrel{\sim }{E}}_{\ell ,m},{\stackrel{\sim }{B}}_{\ell ,m})$ are related to the spin-weighted harmonic window functions according to:

$$\begin{array}{cc}\hfill {\stackrel{\sim }{E}}_{\ell ,m}& ={\displaystyle \sum _{{\ell }^{\prime },m^{\prime }}}(E_{{\ell }^{\prime },m^{\prime }}{W}_{\ell {\ell }^{\prime }m{m}^{\prime }}^{+}+B_{{\ell }^{\prime },m^{\prime }}{W}_{\ell {\ell }^{\prime }m{m}^{\prime }}^{-}),\hfill \end{array}$$(A.5)

$$\begin{array}{cc}\hfill {\stackrel{\sim }{B}}_{\ell ,m}& ={\displaystyle \sum _{{\ell }^{\prime },m^{\prime }}}(B_{{\ell }^{\prime },m^{\prime }}{W}_{\ell {\ell }^{\prime }m{m}^{\prime }}^{+}-E_{{\ell }^{\prime },m^{\prime }}{W}_{\ell {\ell }^{\prime }m{m}^{\prime }}^{-}),\hfill \end{array}$$(A.6)

where we have defined the compound matrices,

$$\begin{array}{cc}\hfill W_{\ell {\ell }^{\prime }m{m}^{\prime }}^{+}& =\frac{1}{2}\left({}_2{W}_{\ell {\ell }^{\prime }}^{m{m}^{\prime }}{+}_{-2}{W}_{\ell {\ell }^{\prime }}^{m{m}^{\prime }}\right),\hfill \end{array}$$(A.7)

$$\begin{array}{cc}\hfill W_{\ell {\ell }^{\prime }m{m}^{\prime }}^{-}& =\frac{\mathrm{i}}{2}\left({}_2{W}_{\ell {\ell }^{\prime }}^{m{m}^{\prime }}{-}_{-2}{W}_{\ell {\ell }^{\prime }}^{m{m}^{\prime }}\right),\hfill \end{array}$$(A.8)

and where Eq. (A.3) is evaluated for s = ±2.

The cut-sky Pseudo-C_ℓ power spectra are related to the full-sky signal by a mixing matrix, M_ℓℓ′ (see e.g. the appendix in Brown et al. 2005 for a derivation and discussion). Explicitly, for each 3×2pt component, the Pseudo-C_ℓ power spectra are given by the following matrix expressions.

• Angular clustering,

  $$\begin{array}{c}\hfill ⟨{\stackrel{\sim }{C}}_{\ell }^{{\delta }_g(i){\delta }_g(j)}⟩={\displaystyle \sum _{{\ell }^{\prime }}}{W}_{\ell {\ell }^{\prime }}^{00}(i,j)C_{{\ell }^{\prime }}^{{\delta }_g(i){\delta }_g(j)}.\end{array}$$(A.9)

• Cosmic shear,

  $$\begin{array}{cc}\hfill & \left(\begin{array}{c}⟨{\stackrel{\sim }{C}}_{\ell }^{E(i)E(j)}⟩\\ ⟨{\stackrel{\sim }{C}}_{\ell }^{E(i)B(j)}⟩\\ ⟨{\stackrel{\sim }{C}}_{\ell }^{B(i)B(j)}⟩\end{array}\right)={\displaystyle \sum _{{\ell }^{\prime }}}[0cm0.9cm\hfill \\ \hfill & \left(\begin{array}{ccc}{W}_{\ell {\ell }^{\prime }}^{++}(i,j)& W_{\ell {\ell }^{\prime }}^{-+}(i,j)+W_{\ell {\ell }^{\prime }}^{+-}(i,j)& W_{\ell {\ell }^{\prime }}^{--}(i,j)\\ -W_{\ell {\ell }^{\prime }}^{+-}(i,j)& W_{\ell {\ell }^{\prime }}^{++}(i,j)-W_{\ell {\ell }^{\prime }}^{--}(i,j)& W_{\ell {\ell }^{\prime }}^{-+}(i,j)\\ W_{\ell {\ell }^{\prime }}^{--}(i,j)& -W_{\ell {\ell }^{\prime }}^{-+}(i,j)-W_{\ell {\ell }^{\prime }}^{+-}(i,j)& W_{\ell {\ell }^{\prime }}^{++}(i,j)\\ \end{array}\right),\hfill \\ \hfill & \left(\begin{array}{c}{C}_{{\ell }^{\prime }}^{E(i)E(j)}\\ C_{{\ell }^{\prime }}^{E(i)B(j)}\\ C_{{\ell }^{\prime }}^{B(i)B(j)}\end{array}\right)0cm0.9cm].\hfill \end{array}$$(A.10)

• Galaxy-galaxy lensing

  $$\begin{array}{c}\hfill \left(\begin{array}{c}⟨{\stackrel{\sim }{C}}_{\ell }^{{\delta }_g(i)E(j)}⟩\\ ⟨{\stackrel{\sim }{C}}_{\ell }^{{\delta }_g(i)B(j)}⟩\end{array}\right)={\displaystyle \sum _{{\ell }^{\prime }}}\left(\begin{array}{cc}{W}_{\ell {\ell }^{\prime }}^{0+}(i,j)& W_{\ell {\ell }^{\prime }}^{0-}(i,j)\\ -W_{\ell {\ell }^{\prime }}^{0-}(i,j)& W_{\ell {\ell }^{\prime }}^{0-}(i,j)\end{array}\right)\left(\begin{array}{c}{C}_{{\ell }^{\prime }}^{{\delta }_g(i)E(j)}\\ C_{{\ell }^{\prime }}^{{\delta }_g(i)B(j)}\end{array}\right).\end{array}$$(A.11)

Here, we have defined the tomographic mixing functions,

$$\begin{array}{c}\hfill W_{\ell {\ell }^{\prime }}^{\mathit{MN}}(i,j)=\frac{1}{2\ell +1}{\displaystyle \sum _{m{m}^{\prime }}}{W}_{\ell {\ell }^{\prime }m{m}^{\prime }}^{M,i}{\left(W_{\ell {\ell }^{\prime }m{m}^{\prime }}^{N,j}\right)}^{\ast },\end{array}$$(A.12)

for M, N = (0, +, −). Here, the W_{ℓℓ′mm′}^+/− are as defined in Eq. (A.8), and $W_{\ell \ell \prime mm\prime }^0{=}_0{W}_{\ell \ell \prime }^{mm\prime }$ from Eq. (A.2).

Appendix B: Results and validation

In this section, we present the validation framework we used to investigate the accuracy of our simulation method to reproduce a target tomographic 3×2pt signal for both a no-noise and noisy set-up. We measured the 3×2pt Pseudo-C_ℓ power spectra over 3000 realisations of our mock catalogues, and demonstrate that the sample generated by our pipeline can self-consistently reproduce the underlying power spectra predicted by the input cosmology to an encouraging accuracy. For our validation framework, we adopted the simulation set-up parameters that are defined at the start of Sect. 6, and measured the Pseudo-C_ℓ power spectra using ten log-spaced bandpowers, taking the range 100 ≤ ℓ ≤ 600 for angular galaxy clustering and galaxy-galaxy lensing, and 100 ≤ ℓ ≤ 1500 for cosmic shear.

B.1. No-noise simulation

For our first-order validation demonstration, we present results for a ‘no-noise’ simulation which includes no systematic contribution from redshift uncertainty or shape noise. In Fig. B.1 we present the binned n(z) for this validation tomography and the corresponding measurements of the Pseudo-C_ℓ bandpowers of the data for each of the 3×2pt components. Alongside, we plot the theoretical prediction, constructed using Eqs. (24–A.12), and the normalised (fractional) residuals for each spectrum, which we denote as Δ_f. We calculated the errorbars using the standard deviation of the data points measured from the 3000 realisations.

Fig. B.1. Top-left: Binned n(z) for a tomographic 3×2pt analysis of our mock catalogues using three equipopulated bins, with no contribution from shape noise or photo-z uncertainty in the simulation. Top-right, bottom-left, bottom-right: Measured 3×2pt Pseudo-Cℓ power spectra (cosmic shear, angular clustering, galaxy-galaxy lensing, respectively) measured over 3000 realisations for this chosen tomography, shown in crosses. We bin the power spectra into ten log-spaced bandpowers and measure the signal in the range 100 ≤ ℓ ≤ 1500 for cosmic shear, and 100 ≤ ℓ ≤ 600 for galaxy clustering and galaxy-galaxy lensing. The solid black lines show the input model power spectra. Beneath each individual power spectrum we plot the normalised residuals of the data, defined as the fractional quantity ($C_{\ell }^{\mathrm{Measured}}/C_{\ell }^{\mathrm{Theory}})-1$, which we denote as Δf on the y-axis label. The vertical double-headed arrows in the off-diagonal subplots indicate that the fractional residuals, and their errors, for these spectra are extremely large due to the division by an extremely small signal.

In general, we find that for the ‘no-noise’ set-up, the 3×2pt measurements from our catalogues show very good agreement with the fiducial model. Collectively, the residuals of the 3×2pt data vector demonstrate that our method can accurately reconstruct the underlying 3×2pt signal to a sub-per cent level of accuracy, which we highlight as a very encouraging performance of our simulation approach.

Comparatively, the galaxy clustering power spectra show that the measured signal in the auto-correlation bins shows a minimal ( ∼ 0.5%) under-recovery with respect to the theoretical prediction. This behaviour is not present in the shear or galaxy-galaxy lensing power spectrum, which both show a more random scatter about the model. Additionally, we note that for a no-noise simulation, in which source galaxy redshifts are estimated with complete precision, there is no overlap between the redshift distributions of different tomographic bins. Hence, the cross correlation power spectra of the galaxy clustering signal, and the galaxy-galaxy lensing cross spectra, in which a foreground (low redshift) shear field is correlated with a background clustering field, each have extremely low signal. The corresponding normalised residuals are therefore extremely large due to the small number statistics associated with the very small quantities.

B.2. Realistic, ‘noisy’ realisations with shape noise and photo-z uncertainty

In this section, we present results of the 3×2pt Pseudo-C_ℓ power spectra measured from our simulated mock catalogues that now include contributions from shape noise and photometric redshift estimation uncertainty. We modelled the uncertainties arising from each of these sources using Gaussian-distributed errors, taking σ(ϵ^int) = 0.3 and ${\sigma }_z^{\mathrm{phot}}=0.05$ as defined in Sect. 6.

In Fig. B.2, we present the bandpowers of the 3×2pt Pseudo-C_ℓ power spectra measured over 3000 simulations using the same validation framework of three bins equipopulated with galaxies as considered in the no-noise case in Appendix B.1. We show the normalised residuals under each power spectrum component in Fig. B.2, and present the binned n(z) in the upper left corner.

Fig. B.2. Top-left: Binned n(z) for a three equipopulated bin tomographic set-up, measured from our mock catalogues that include realistic Gaussian shape noise and photo-z uncertainties in the simulation. Top-right, bottom-left, bottom-right: Cosmic shear, angular clustering, and galaxy-galaxy lensing components (respectively) of the 3×2pt data vector. We show in crosses the Pseudo-Cℓ power spectra measured over 3000 realisations of our simulation for this chosen tomography. We bin the multipoles into ten log-spaced bandpowers and measure the signal in the range 100 ≤ ℓ ≤ 1500 for cosmic shear, and 100 ≤ ℓ ≤ 600 for galaxy clustering and galaxy-galaxy lensing. We plot the theoretical prediction in black alongside the measured data, and underneath each individual power spectrum we plot the normalised residuals, which we define as the fractional quantity ($C_{\ell }^{\mathrm{Measured}}/C_{\ell }^{\mathrm{Theory}})-1$, denoted as Δf on the y-axis label.

We find that the 3×2pt power spectra measured from the simulated noisy catalogues are in good general agreement with the theoretical prediction. In the presence of Gaussian shape noise and photo-z errors, the shear and galaxy-galaxy lensing signals both show a similar level of recovery of the fiducial model, at the ∼1% level. This is consistent with that seen in the noise-free 3×2pt simulations. However, for the galaxy clustering, the measured auto-correlation power spectra are consistently ∼2% lower than the prediction from the theoretical model. This is significantly larger than the ∼0.5% deficit found in the corresponding no-noise simulations. Furthermore, the cross-correlation power spectra between bins are in more considerable disagreement, with the recovered power spectra being 5–6% in excess of the theoretical prediction. While the absolute values of these cross spectra are over a magnitude smaller than the auto-correlation signals, this 5–6% level discrepancy in the cross-correlation clustering power spectra represents the limiting precision of our simulation approach.

We attribute the above discrepancies to a small (but non-negligible) fundamental inconsistency between our simulation process and our model predictions. The former is an attempt to simulate the full 3D distribution and evolution of matter in the Universe, approximating it as a stack of 2D fields on the sky that have been evaluated using a tomographic projection over finely binned, discretised points in redshift. The resulting catalogues are then binned into much broader redshift bins for the cosmological power spectrum analysis. In contrast, the theoretical model directly predicts the signal in these much broader redshift bins, again using a tomographic projection but now performing the Limber integration over a much broader redshift kernel (see Eqs. 1–6). Consequently, our simulation and modelling approaches contain different levels of approximation. This is the fundamental reason for the small discrepancies described above. Furthermore, one expects that the disagreement would be exacerbated in the presence of photometric redshift errors, which is what we have found.

All Tables

All Figures

	Fig. 1. Flowchart describing the structure of our simulation to generate Stage IV-like mock galaxy catalogues from 2D realisations on the sky of the 3 × 2 pt signal that traces an underlying redshift distribution and fiducial cosmology. The red dashed and green dotted sections represent the simulation stages that are executed using CosmoSIS and FLASK, respectively.
In the text

Fig. 2. Illustration of the method used to simulate a mock weak lensing survey. 2D maps of the correlated 3 × 2 pt fields are generated at finely sampled points in redshift. These 2D maps are used to approximate the full 3D cosmological information of the 3 × 2 pt signal. From the overdensity fields, we Poisson sample a galaxy population that traces the underlying n(z) distribution. We then assign the correlated weak lensing observables to each galaxy from the shear field values at the galaxy’s angular position on the sky at a given redshift. We show an early approximation of the Euclid DR1 footprint at the bottom of the figure which we have used to create our Euclid-like simulations. The ‘observed’ region is shown in yellow. Note: the actual Euclid DR1 footprint will be significantly different to that shown here.

In the text

Fig. 3. Areas enclosed by the 1σ (Blue) and 2σ contours in the (w0, wa) plane, for different numbers of redshift bins used in a tomographic analysis of the cosmic shear signal measured from 3000 realisations of our simulation. We show in circular markers joined with solid lines the areas measured for the equipopulated binning choice; in square markers joined with dashed lines the equal comoving distance bins; and in triangular markers joined with dotted lines the equal redshift width binning choice. Since the cosmic shear component alone is relatively weakly constraining, the 1 bin measurement does not yield a closed contour in (w0, wa) within the ranges of the parameter grid. Hence, the data point for this case represents a lower bound of the true value, which we represent by using a vertical arrow.

In the text

Fig. 4. Contour areas of the 1 and 2σ constraints on (w0, wa), measured from the tomographic angular clustering component of our no-noise simulation. We plot the contour areas measured as a function of the number of redshift bins used in a tomographic analysis, considering the equipopulated, equal comoving distance, equally spaced in redshift binning strategies. The 1 bin measurement does not sufficiently constrain either parameter within the prior volume, hence we include a vertical arrow to denote that this data point is a lower bound.

In the text

	Fig. 5. Areas enclosed by the 1 and 2σ contours of the posterior constraints on (w0, wa), derived from measurements of the full tomographic 3 × 2 pt signal from our ‘no-noise’ mock catalogues. We plot the contour areas as a function of the number of redshift bins used in the tomographic analysis for the three different binning choices considered in this work.
In the text

Fig. 6. 1 and 2σ constraints on the (w0, wa) parameters derived from measurements of the tomographic cosmic shear signal in our noisy simulation that includes contributions from Gaussian-distributed photo-z errors and Gaussian-distributed shape noise. We present the constraints for different numbers of equipopulated bins but find that the signal has relatively weak constraining power that does not significantly improve as the number of bins increases. Hence, the contours corresponding to each different number of bins lie on top of each other.

In the text

Fig. 7. The 1σ contour areas of constraints in (w0, wa) from measurements of the tomographic cosmic shear signal in our noisy simulation. We consider three different binning strategies and find that the 1 bin measurement does not have sufficient constraining power to generate closed contours within the prior ranges of the parameter grid. Hence, we plot the unbounded area measured within the prior ranges, and indicate that this is a lower bound of the true value by using a vertical arrow.

In the text

	Fig. 8. 1 and 2σ posterior constraints on (w0, wa) measured from the tomographic angular clustering signal in our mock catalogues which include simulated noise from Gaussian photo-z uncertainties and shape noise. We show the contours for different numbers of equipopulated bins used for the tomographic analysis.
In the text

	Fig. 9. Contour areas of the 1 and 2σ constraints on (w0, wa) for measurements of the tomographic angular clustering signal in our noisy simulation. We show the contour areas measured as a function of the number of redshift bins used, for the three different tomographic binning choices considered in this work.
In the text

	Fig. 10. 1 and 2σ constraints on (w0, wa) from tomographic Pseudo-Cℓ measurements of the full 3 × 2 pt signal in our noisy simulation. We show constraints measured for different numbers of equipopulated bins used for the tomographic analysis. The cross marks the fiducial cosmology (w0, wa) = (−1, 0) used for the simulation.
In the text

	Fig. 11. Contour areas enclosed by the 1 and 2σ constraints on the dark energy (w0, wa) parameters, measured for the full 3 × 2 pt signal in the presence of Gaussian shape noise and photo-z uncertainty. We vary the number of redshift bins used for each of the equipopulated, equally spaced in fiducial comoving distance, and equal redshift width binning strategies.
In the text

	Fig. 12. Areas of the 1 and 2σ constraints on (w0, wa), measured from our noisy catalogue simulations, plotted in circular markers as a function of the number of equipopulated bins used for the tomographic analysis. Black solid lines show the plot of the best fit functions of the form y = a/x + b that we fit to these data sets, evaluated using the SciPy curve_fit routine.
In the text

Fig. 13. Constraints on the dark energy (w0, wa) parameters measured from: simulations that include only Gaussian shape noise and Gaussian photo-z uncertainties (blue) and simulations that include Gaussian shape noise, Gaussian photo-z uncertainties, and a 5% level of catastrophic photo-z errors (outliers) that result from the Lyman-α line in the galaxy SED being misidentified as the Balmer and D4000 break features (red). From left to right, we plot the constraints derived from tomographic analyses of the full 3 × 2 pt signal using a single redshift bin, five equipopulated bins, and ten equipopulated bins.

In the text

Fig. B.1. Top-left: Binned n(z) for a tomographic 3×2pt analysis of our mock catalogues using three equipopulated bins, with no contribution from shape noise or photo-z uncertainty in the simulation. Top-right, bottom-left, bottom-right: Measured 3×2pt Pseudo-Cℓ power spectra (cosmic shear, angular clustering, galaxy-galaxy lensing, respectively) measured over 3000 realisations for this chosen tomography, shown in crosses. We bin the power spectra into ten log-spaced bandpowers and measure the signal in the range 100 ≤ ℓ ≤ 1500 for cosmic shear, and 100 ≤ ℓ ≤ 600 for galaxy clustering and galaxy-galaxy lensing. The solid black lines show the input model power spectra. Beneath each individual power spectrum we plot the normalised residuals of the data, defined as the fractional quantity ($C_{\ell }^{\mathrm{Measured}}/C_{\ell }^{\mathrm{Theory}})-1$, which we denote as Δf on the y-axis label. The vertical double-headed arrows in the off-diagonal subplots indicate that the fractional residuals, and their errors, for these spectra are extremely large due to the division by an extremely small signal.

In the text

Fig. B.2. Top-left: Binned n(z) for a three equipopulated bin tomographic set-up, measured from our mock catalogues that include realistic Gaussian shape noise and photo-z uncertainties in the simulation. Top-right, bottom-left, bottom-right: Cosmic shear, angular clustering, and galaxy-galaxy lensing components (respectively) of the 3×2pt data vector. We show in crosses the Pseudo-Cℓ power spectra measured over 3000 realisations of our simulation for this chosen tomography. We bin the multipoles into ten log-spaced bandpowers and measure the signal in the range 100 ≤ ℓ ≤ 1500 for cosmic shear, and 100 ≤ ℓ ≤ 600 for galaxy clustering and galaxy-galaxy lensing. We plot the theoretical prediction in black alongside the measured data, and underneath each individual power spectrum we plot the normalised residuals, which we define as the fractional quantity ($C_{\ell }^{\mathrm{Measured}}/C_{\ell }^{\mathrm{Theory}})-1$, denoted as Δf on the y-axis label.

In the text