© The Authors 2024

1. Introduction

The large-scale structure in the Universe has been a key contributor to testing the theories of gravity, the nature of dark energy, and dark matter. It gives us vital information about the growth of structure in the Universe and the history of cosmic expansion. Observations from galaxy survey experiments such as the Sloan Digital Sky Survey (SDSS; Gunn et al. 2006; Strauss et al. 2002), the Wide-field Infrared Survey Explorer (WISE; Schlafly et al. 2019; Wright et al. 2010), the Kilo-Degree Survey (KiDS; Heymans et al. 2021; de Jong et al. 2015), Hyper Suprime-Cam (HSC; Hikage et al. 2019), the Two Micron All Sky Survey (2MASS; Bilicki et al. 2014), and the Dark Energy Survey (DES; Abbott et al. 2018) have been the torchbearers in unveiling the shortcomings of the standard model of cosmology, the ΛCDM model. The upcoming galaxy surveys, including the Vera C. Rubin Observatory Legacy Survey of Space and Time (LSST; Ivezić et al. 2019; LSST Science Collaboration 2009), Euclid (Laureijs et al. 2011), the Nancy Grace Roman Space Telescope (Spergel et al. 2013), the Dark Energy Spectroscopic Instrument (DESI; Dey et al. 2019), and the Spectro-Photometer for the History of the Universe, Epoch of Reionization, and Ices Explorer (SPHEREx; Doré et al. 2014) will play a crucial role in providing a deeper understanding of how our Universe works.

The cosmic microwave background (CMB) gives a view of the early Universe and carries information about the energy components of the Universe. The gravitational lensing of the CMB carries information about the growth of structure at redshift 1 − 3. The CMB lensing signal has been valuable in studying the properties of dark energy and is a complementary probe to galaxy clustering and galaxy weak lensing. Cross correlations between the lensing map of the CMB and tracers of large-scale structure can be used to extract information about the evolution of the large-scale gravitational potential. The importance of cross correlations between CMB lensing convergence and galaxy positions in testing the validity of the standard cosmological model or the ΛCDM model was firmly established with many studies in the past few years (Saraf et al. 2022; Miyatake et al. 2022; Robertson et al. 2021; Krolewski et al. 2021; Darwish et al. 2021 Abbott et al. 2019; Bianchini & Reichardt 2018; Singh et al. 2017; Bianchini et al. 2015, 2016).

Cross-correlation tomography performed with galaxy samples divided into narrow redshift bins also allows us to map the evolution of the cosmological parameters with redshift. Several tomographic cross-correlation studies (e.g. Wang et al. 2023; Yu et al. 2022; White et al. 2022; Pandey et al. 2022; Chang et al. 2023; Sun et al. 2022; Krolewski et al. 2021; Hang et al. 2021; Marques & Bernui 2020; Peacock & Bilicki 2018; Giannantonio et al. 2016) have identified differences in the value of cosmological parameters, such as the σ₈, Ω_m, or the combined S₈ parameter (defined as $S_8\equiv {\sigma }_8\sqrt{{\mathrm{\Omega }}_m/0.3}$), within the ΛCDM model. These low-redshift cross-correlation probes consistently measure a lower value for these parameters than the high-redshift CMB-only measurements from the Planck satellite (Planck Collaboration VI 2020), resulting in the so-called S₈ and σ₈ − Ω_m tensions. Other studies such as Bianchini & Reichardt (2018), Amon et al. (2018), Blake et al. (2016), Giannantonio et al. (2016), and Pullen et al. (2016) reported consistent deviations in the D_g and E_g statistics when they tested the ΛCDM model with different galaxy surveys. The galaxy survey experiments in the future will play a pivotal role in increasing the significance of either tension or agreement in cosmological parameters and in quantifying possible deviations from the standard ΛCDM model.

In this paper, we present a tomographic cross-correlation study through 300 Monte Carlo (MC) simulations of the Vera C. Rubin Observatory Legacy Survey of Space and Time (LSST) galaxy survey and Planck CMB lensing convergence. We include root mean square scatter in redshifts $\frac{\sigma (z)}{1+z}=0.02\text{and}0.05$, but exclude catastrophic redshift errors and photometric calibration errors from our simulations. We superficially create the best possible observations in this way. We show that even for these idealistic observations, the photometric redshift uncertainties in the redshift distributions produce biased estimates of the cross correlation when the redshift bin mismatch of objects due to photometric redshift errors is not properly taken into account. Although a few attempts have been made to account for this effect (Balaguera-Antolínez et al. 2018; Hang et al. 2021), we propose a new fast and efficient method of estimating the tomographic cross correlation in an unbiased way using the scattering matrix formalism introduced by Zhang et al. (2010).

The paper is organised as follows: Sect. 2 presents the simulation setup and the theoretical background, and Sect. 3 describes the method for computing the power spectra, reconstructing true redshift distributions, and estimating the parameters. In Sect. 4, we present the validation results of our simulations and estimation of galaxy auto-power spectra and cross-power spectra between galaxy over-density and CMB lensing convergence. We discuss the correction to the power spectra through the scattering matrix in Sect. 5 and quantify the impact on the linear galaxy bias and on the amplitude of the cross-correlation parameters in Sect. 6. We study the apparent tension on the σ₈ parameter due to the redshift bin mismatch in Sect. 7. Finally, we summarise our results in Sect. 8.

2. Simulations and theory

We used the publicly available code FLASK (Xavier et al. 2016) to generate 300 tomographic MC realisations of correlated log-normal galaxy over-density and CMB lensing convergence fields. The galaxy density followed the LSST photometric redshift distribution profile (Ivezić et al. 2019; LSST Science Collaboration 2009) with a mean redshift at 0.9 and a mean surface number density of 40 arcmin⁻². The simulated CMB lensing convergence field was consistent with Planck observations (Planck Collaboration VIII 2020). The galaxy density field was induced with Poisson noise, and for the CMB convergence, we used the noise power spectrum provided in the Planck 2018 data package¹. Due to computational limitations, the sky area covered in our simulations was 2000 deg². However, the results obtained in this study should remain valid for the planned area of the LSST survey if the errors are appropriately scaled by the fraction of sky coverage. The mock galaxy samples were divided into nine disjoint tomographic bins with redshift intervals of (0.0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.4, 1.8, 2.2, 3.0], as marked by the dashed vertical lines in Fig. 1.

Fig. 1. Effect of photometric redshift scatter on redshift distribution. True LSST redshift distribution (solid blue curve) divided into nine photometric redshift bins (shown by the solid green curves) for Gaussian redshift errors with (a) σ0 = 0.02, (b) σ0 = 0.05, and (c) for modified Lorentzian redshift errors with γ0 = 0.02. The red line denotes the CMB lensing kernel, and the dashed vertical orange lines mark the true redshift bins. The full LSST redshift distribution and CMB lensing kernel are normalised to unit maximum.

The fiducial angular power spectra for each redshift bin we used to generate correlated maps were computed under the Limber approximation (Limber 1953),

$$\begin{array}{c}\hfill C_{\ell }^{\mathit{xy}}={\displaystyle \underset{0}{\overset{\infty }{\int }}}\frac{\mathrm{d}z}{c}\frac{H(z)}{{\chi }^2(z)}{W}^x(z)W^y(z)P(k=\frac{\ell +1/2}{\chi (z)},z),\end{array}$$(1)

where {x, y}∈{κ, g}, κ≡ convergence and g≡ galaxy over-density, c is the speed of light, and P(k, z) is the matter power spectrum generated using the public software CAMB² (Lewis et al. 2000) using the HALOFIT prescription. The kernels W^x(z) connects the observables to the underlying total matter distribution. Assuming a flat Universe (Planck Collaboration VI 2020), the lensing kernel W^κ(z) and the galaxy kernel W^g(z) are expressed as

$$\begin{array}{cc}& W^{\kappa }(z)=\frac{3{\mathrm{\Omega }}_m}{2c}\frac{H_0^2}{H(z)}(1+z)\chi (z)\frac{{\chi }_{\ast }-\chi (z)}{{\chi }_{\ast }}\hfill \end{array}$$(2)

$$\begin{array}{cc}& W^g(z)=b(z)\frac{\mathrm{d}N}{\mathrm{d}z}.\hfill \end{array}$$(3)

Here, H(z) is the Hubble parameter at redshift z, and χ(z) and χ_* are the comoving distances to redshift z and the surface of the last scattering. Ω_m and H₀ are the present-day values of the matter density parameter and Hubble constant, respectively. $\frac{\mathrm{d}N}{\mathrm{d}z}$ stands for the normalised redshift distribution of galaxies, and b(z) is the linear galaxy bias that relates the galaxy over-density to the total underlying matter density (Fry & Gaztanaga 1993). We assumed that the magnification bias does not contribute to the galaxy kernel and leave the study of its effect for future work. The expressions for kernels W^κ(z) and W^g(z) given by Eqs. (2) and (3) are valid in a flat Universe model (i.e. for zero curvature).

We used a redshift-dependent model of the galaxy bias (Solarz et al. 2015; Moscardini et al. 1998; Fry 1996),

$$\begin{array}{c}\hfill b(z)=1+\frac{b_0-1}{D(z)},\end{array}$$(4)

where we took b₀ ≡ b(z = 0) = 1.3, and D(z) is the linear growth function normalised to unity at z = 0

$$\begin{array}{c}\hfill D(z)=exp\{-{\displaystyle \underset{0}{\overset{z}{\int }}}\frac{{[{\mathrm{\Omega }}_m(z^{\prime })]}^{\gamma }}{1+z^{\prime }}\mathrm{d}{z}^{\prime }\},\end{array}$$(5)

where γ = 0.55 is the growth index for General Relativity (Linder 2005).

For every simulation, the code FLASK produces CMB convergence map and galaxy number count maps for each tomographic bin along with the catalogue of galaxy redshifts. We call these maps true datasets. We generated photometric redshifts, z_p, for galaxies in every simulation by drawing positive random values from a Gaussian distribution 𝒩(z_t, σ₀(1 + z_t)) (hereafter, Case-I) and from modified Lorentzian function L(z_p; z_t, γ₀(1 + z_t),a) (hereafter, Case-II) given by

$$\begin{array}{c}\hfill L(z_p;z_t,{\gamma }_0(1+z_t),a)\propto [1+\frac{1}{2a}(\frac{z_p-z_t}{{\gamma }_0(1+z_t)}{)}^2{]}^{-a}.\end{array}$$(6)

We adopted two different values of σ₀, 0.02 and 0.05, for the Gaussian distribution, and we took γ₀ = 0.02,  a = 2 for the modified Lorentzian distribution to study the dependence of our results on the strength of the redshift scatter. The galaxies were again divided into nine tomographic bins based on their photometric redshifts. In Fig. 1, we show the true LSST redshift distribution (solid blue curve) divided into nine photometric redshift bins (solid green lines) and the CMB lensing kernel (solid red curve) in the redshift range 0 ≤ z ≤ 3. The green lines represent how the disjoint true redshift bins transform after the photometric redshift errors are introduced. The dashed vertical orange lines mark the boundaries of true redshift bins.

We built galaxy over-density maps from photometric number count maps (hereafter, photometric datasets) with a HEALPix³ (Górski et al. 2005) resolution parameter N_side = 1024 using

$$\begin{array}{c}\hfill g(\widehat{n})=\frac{n(\widehat{n})-\overline{n}}{\overline{n}},\end{array}$$(7)

where $n(\widehat{n})$ is the number of galaxies at an angular position $\widehat{n}$, and $\overline{n}$ is the mean number of galaxies per pixel. In Table 1, we present based on one realisation the comparison of the mean number of objects per pixel (for N_side = 1024) and median redshift between the true and photometric datasets.

For the simulations and analyses presented in this paper, we adopted the flat ΛCDM cosmology with the best-fit Planck + WP + highL + lensing cosmological parameters as described in Planck Collaboration VI (2020). Here, WP refers to WMAP polarisation data at low multipoles, highL are the high-resolution CMB data from the Atacama Cosmology Telescope (ACT), and South Pole Telescope (SPT) and lensing refer to the inclusion of Planck CMB lensing data in the parameter likelihood.

3. Method

In this section, we outline the method for computing power spectra from maps and for estimating the true redshift distribution from the photometric redshift distribution.

3.1. Estimating the power spectra

We extracted the full sky power spectra from partial sky power spectra for every tomographic bin using the pseudo-C_ℓ method based on the MASTER algorithm (Hivon et al. 2002), taking into account mode coupling induced by incomplete sky coverage and pixelisation effects. We estimated the full sky power spectra in linearly spaced multipole bins with Δℓ = 30 between 50 ≤ ℓ ≤ 1500. The noise-subtracted mean full sky power spectrum over N_sim realisations was computed as (Saraf et al. 2022)

$$\begin{array}{c}\hfill {\overline{C}}_L^{\mathit{xy}}\equiv ⟨{\widehat{C}}_L^{\mathit{xy}}⟩=\frac{1}{N_{\text{sim}}}{\displaystyle \sum _{i=1}^{N_{\text{sim}}}}{\widehat{C}}_L^{xy,i}-{⟨N_L^{\mathit{xy}}⟩}_{\text{MC}},\end{array}$$(8)

where ${\widehat{C}}_L^{x\mathit{y},i}$ represents the full sky power spectrum estimate for the ith simulation, and ${⟨N_L^{x\mathit{y}}⟩}_{\text{MC}}$ is the average noise power spectrum from the Monte Carlo simulations. The errors associated with the mean power spectrum were computed from the diagonal of the power spectrum covariance matrix ${\text{Cov}}_{L{L}^{\prime }}^{xy}$ as

$$\begin{array}{c}\hfill \mathrm{\Delta }{\overline{C}}_L^{\mathit{xy}}=(\frac{{\text{Cov}}_{\mathit{LL}}^{\mathit{xy}}}{N_{\text{sim}}}{)}^{1/2},\end{array}$$(9)

where

$$\begin{array}{c}\hfill {\text{Cov}}_{L{L}^{\prime }}^{\mathit{xy}}=\frac{1}{N_{\text{sim}}-1}{\displaystyle \sum _{i=1}^M}({\widehat{C}}_L^{xy,i}-{\overline{C}}_L^{\mathit{xy}})({\widehat{C}}_{L^{\prime }}^{xy,i}-{\overline{C}}_{L^{\prime }}^{\mathit{xy}}).\end{array}$$(10)

3.2. Estimating the true redshift distribution

The true redshift distribution $\frac{\mathrm{d}N(z_t)}{\mathrm{d}{z}_t}$ was estimated from the observed photometric redshift distribution $\frac{\mathrm{d}N(z_p)}{\mathrm{d}{z}_p}$ and some quantification of errors on the photometric redshifts (from a cross validation with some spectroscopic survey or posteriors from machine-learning methods). These errors are often expressed by conditional probabilities p(z_p − z_t|z_t) and p(z_t − z_p|z_p), which we call photometric redshift error distributions. The true and photometric redshift distributions are then related through (Sheth & Rossi 2010)

$$\begin{array}{c}\hfill \frac{\mathrm{d}N(z_t,z_p)}{\mathrm{d}{z}_t\mathrm{d}{z}_p}=\frac{\mathrm{d}N(z_t)}{\mathrm{d}{z}_t}p(z_p-z_t|z_t)=\frac{\mathrm{d}N(z_p)}{\mathrm{d}{z}_p}p(z_t-z_p|z_p).\end{array}$$(11)

Thus, depending on whether we have estimates of p(z_p − z_t|z_t) or p(z_t − z_p|z_p), the method for estimating the true redshift distribution is called deconvolution or convolution, respectively.

3.2.1. Convolution method

When p(z_t − z_p|z_p) is known, we can estimate $\frac{\mathrm{d}N(z_t)}{\mathrm{d}{z}_t}$ for each tomographic bin i using the convolution

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}^i(z_t)}{\mathrm{d}{z}_t}={\displaystyle \int \mathrm{d}}{z}_p\frac{\mathrm{d}N(z_p)}{\mathrm{d}{z}_p}{W}^i(z_p)p^i(z_t-z_p|z_p),\end{array}$$(12)

where $\frac{\mathrm{d}N(z_p)}{\mathrm{d}{z}_p}$ is the observed photometric redshift distribution of the objects, and Wⁱ(z_p) is the window function defining the ith redshift bin given by a step function,

$$\begin{array}{c}\hfill W^i(z)=\{\begin{array}{cc}1,\hfill & \text{if}{z}_{\text{min}}^i\le z<z_{\text{min}}^{i+1}\hfill \\ 0,\hfill & \text{Otherwise}.\hfill \end{array}\end{array}$$(13)

Generally, pⁱ(z_t − z_p|z_p) is fitted with parametric functions such as a Gaussian with an assumed zero mean (Sun et al. 2022; Marques & Bernui 2020) or a modified Lorentzian (Hang et al. 2021; Peacock & Bilicki 2018). Sheth & Rossi (2010) showed that the quantity p(z_t − z_p|z_p) is biased and is not centred on zero. In our study, we fit the error distributions p(z_t − z_p|z_p) with a modified Lorentzian function (for Case-II) and a sum of three Gaussians (for Case-I),

$$\begin{array}{c}\hfill \mathcal{N}(x)={\displaystyle \sum _{i=1}^3}{A}_{i}exp[\frac{{(x-{\mu }_i)}^2}{2{\sigma }_i^2}],\end{array}$$(14)

where A, μ, σ control the amplitude, mean, and width of the individual Gaussians. The sum of the Gaussians can account for the bias in p(z_t − z_p|z_p) and for other characteristic features of the error distributions, such as non-Gaussian wings and a higher peak in the centre. We also checked the fit with a higher number of Gaussians, but did not find any improvement in the quality of fit beyond three Gaussians. For each of the tomographic bins, we fit for A_i, μ_i, σ_i, and estimated the true redshift distribution using Eq. (12).

3.2.2. Deconvolution method

The true redshift distribution $\frac{\mathrm{d}N(z_t)}{\mathrm{d}{z}_t}$ can be estimated by a deconvolution method when p(z_p − z_t|z_t) is known,

$$\begin{array}{c}\hfill \frac{\mathrm{d}N(z_p)}{\mathrm{d}{z}_p}={\displaystyle \int \mathrm{d}}{z}_t\frac{\mathrm{d}N(z_t)}{\mathrm{d}{z}_t}p(z_p-z_t|z_t).\end{array}$$(15)

We fit p(z_p − z_t|z_t) with a single Gaussian for Case-I and with a modified Lorentzian function for Case-II. We find the mean to be consistent with zero, in agreement with the unbiased nature of p(z_p − z_t|z_t) (Sheth & Rossi 2010). For Case-I, we also fit the error distribution with a sum of Gaussians to find no significant improvement in the fit quality. Padmanabhan et al. (2005) proposed a deconvolution method based on Tikhonov regularisation, which lacks a general method for quantifying the impact of the penalty function on the reconstruction of $\frac{\mathrm{d}N(z_t)}{\mathrm{d}{z}_t}$. We used a different approach for the deconvolution. It is based on the convolution theorem and kernel-based regularisation (Meister 2009). The true redshift distribution in our approach was estimated as

$$\begin{array}{c}\hfill \frac{\mathrm{d}N(z_t)}{\mathrm{d}{z}_t}={\mathcal{F}}^{-1}\left[\frac{\mathcal{F}\left[\frac{\mathrm{d}N(z_p)}{\mathrm{d}{z}_p}\right]}{\mathcal{F}[p(z_p-z_t|z_t)]}\right],\end{array}$$(16)

where ℱ and ℱ⁻¹ represent the Fourier and inverse Fourier transforms, respectively. We show the performance of our deconvolution method through a toy example in Appendix A.

We estimated the true redshift distribution for the entire redshift range, 0 < z ≤ 3. Then, the true redshift distribution for each tomographic bin, i, can be expressed as

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}^i}{\mathrm{d}{z}_t}=\frac{\mathrm{d}N}{\mathrm{d}{z}_t}{W}^i(z_t),\end{array}$$(17)

where Wⁱ(z_t) is the window function given by Eq. (13).

The corresponding photometric redshift distribution for every tomographic bin i is given by

$$\begin{array}{c}\hfill \frac{\mathrm{d}{N}^i(z_p)}{\mathrm{d}{z}_p}={\displaystyle \int \mathrm{d}}{z}_t\frac{\mathrm{d}{N}^i(z_t)}{\mathrm{d}{z}_t}{p}^i(z_p-z_t|z_t),\end{array}$$(18)

where pⁱ(z_p − z_t|z_t) is the error distribution for bin i. Equation (15) does not follow convolution strictly near z = 0 because negative redshifts are unphysical. The reconstructed true redshift distribution is therefore inaccurate close to redshift z = 0, and we expect these inaccuracies to affect the first two tomographic bins to some extent.

3.3. Galaxy bias and cross-correlation amplitude

For every tomographic bin, we estimated two parameters: the linear galaxy bias b, and the amplitude of the cross-power spectrum A. We assumed the linear galaxy bias to be constant in each tomographic bin. This is a fairly good assumption for redshift bins, which are narrow relative to the redshift dependence of the bias. The amplitude of the cross-power spectrum acts as a rescaling of the observed cross-power spectrum to the fiducial theoretical power spectrum. From an unbiased estimation of the cross-power spectrum, we expect A = 1 in the ΛCDM cosmology.

The galaxy auto-power spectrum scales as b², while the cross-power spectrum depends on the product b × A. To break this degeneracy, we performed a maximum likelihood estimation on the joint data vector ${\widehat{C}}_L=({\widehat{C}}_L^{\kappa g},{\widehat{C}}_L^{\mathit{gg}})$, using the likelihood function

$$\begin{array}{cc}& \mathcal{L}({\widehat{C}}_L|b,A)=\frac{1}{\sqrt{{(2\pi )}^{N_L}det({\text{Cov}}_{L{L}^{\prime }})}}\times \hfill \\ \hfill & \times \text{exp}\}-\frac{1}{2}[{\widehat{C}}_L-C_L(b,A)]{({\text{Cov}}_{L{L}^{\prime }})}^{-1}[{\widehat{C}}_{L^{\prime }}-C_{L^{\prime }}(b,A)]\{,\hfill \end{array}$$(19)

where C_L(b, A) is the joint theoretical power spectrum template, which is defined as $C_L(b,A)=(A{C}_L^{\kappa g}(b),C_L^{\mathit{gg}}(b))$, and the covariance matrix is given as

$$\begin{array}{c}\hfill {\text{Cov}}_{L{L}^{\prime }}=\left[\begin{array}{cc}{\text{Cov}}_{L{L}^{\prime }}^{\kappa g,\kappa g}& {\text{Cov}}_{L{L}^{\prime }}^{\kappa g,gg}\\ \\ {\text{Cov}}_{L{L}^{\prime }}^{\kappa g,gg}& {\text{Cov}}_{L{L}^{\prime }}^{gg,gg}.\end{array}\right]\end{array}$$(20)

The individual elements of Cov_LL′ are approximated by the expression (Saraf et al. 2022)

$$\begin{array}{cc}& {\text{Cov}}_{L{L}^{\prime }}^{AB,CD}=\frac{1}{(2{\ell }_{L^{\prime }}+1)\mathrm{\Delta }\ell f_{\text{sky}}}[\sqrt{C_L^{\mathit{AC}}{C}_{L^{\prime }}^{\mathit{AC}}{C}_L^{\mathit{BD}}{C}_{L^{\prime }}^{\mathit{BD}}}\hfill \\ \hfill & +\sqrt{C_L^{\mathit{AD}}{C}_{L^{\prime }}^{\mathit{AD}}{C}_L^{\mathit{BC}}{C}_{L^{\prime }}^{\mathit{BC}}}]{\delta }_{L{L}^{\prime }},\hfill \end{array}$$(21)

where {A, B, C, D}∈{κ, g}, Δℓ is the multipole binwidth, and we assumed a common sky coverage fraction between the CMB convergence and the galaxy density fields. Equation (20) was used to estimate the parameters from a single realisation. In Sect. 6, we use the average power spectra from 300 simulations to estimate the parameters. We then divide the covariance matrix Cov_LL′ by the total number of simulations N_sim to obtain the covariance matrix for the average.

We used flat priors b ∈ [0, 10] and A ∈ [ − 5, 5] to estimate the parameters, and the remaining cosmological parameters were kept constant with values from our fiducial background cosmology described in Sect. 2. To effectively sample the parameter space, we used the publicly available software package EMCEE (Foreman-Mackey et al. 2013). The best-fit value of the parameters are medians of their posterior distributions, with ±1σ errors being the 16th and 84th percentile, respectively.

4. Results

In this section, we present the results of estimating the true redshift distribution and power spectra for datasets simulated by the code FLASK. We searched for any systematics in our estimates from the datasets with photometric redshift errors. The results of the tests for datasets without errors are presented in Appendix B.

4.1. Estimating the true redshift distribution

The estimation of the true redshift distribution requires fitting the error function (i.e. either p(z_t − z_p|z_p) or p(z_p − z_t|z_t)) with parametric functions. We used a single Gaussian for p(z_p − z_t|z_t) and the sum of three Gaussian to fit p(z_t − z_p|z_p) for Case-I. For Case-II, we used a single modified Lorentzian function to fit both p(z_t − z_p|z_p) and p(z_p − z_t|z_t). Figures 2a and 2b show the quality of fit to the error functions with σ₀ = 0.02 and compares the single Gaussian with the three Gaussian fit of p(z_t − z_p|z_p) for Case-I. The sum of three Gaussians provides a significantly better fit by accurately capturing the non-Gaussian wings of the error function. Figure 2c shows the fit made to the error distribution p(z_p − z_t|z_t) with a single modified Lorentzian function (Case-II). The uncertainties in these fits remain within 0.25%, and we do not expect these sub-percent uncertainties to bias the power spectra or the parameter estimates.

Fig. 2. Modelling the redshift error distribution with parametric functions. (a) Fit made to the error distribution p(zp − zt|zt) with σ0 = 0.02 (solid red line) using a single Gaussian (blue dashed line). (b) Fit made to the error distribution p(zt − zp|zp) with σ0 = 0.02 (solid red line) using a single Gaussian (dotted orange line), and sum of three Gaussians (dashed blue line). (c) Fit made to the error distribution p(zp − zt|zt) with γ0 = 0.02 (solid red line) using a single modified Lorentzian function (dashed blue line).

In Fig. 3, we show the relative error, averaged over 300 simulations, of the redshift distribution reconstructed using the deconvolution method for the Gaussian error distribution with σ₀ = 0.02, 0.05 and for the modified Lorentzian error distribution with γ₀ = 0.02. The reconstructed redshift distribution is within ∼2% for the entire redshift range. The strongest deviations occur near boundaries and are possibly due to sharp cuts in the redshift distribution at z = 0 and z = 3 (Sect. 3.2.2).

Fig. 3. Average relative difference between the redshift distributions reconstructed using the deconvolution method and the true distribution from 300 simulations for (a) Gaussian error distribution with σ0 = 0.02, (b) a Gaussian error distribution with σ0 = 0.05, and (c) a modified Lorentzian error distribution with γ0 = 0.02. The shaded region represents the 1σ deviations in reconstruction between 300 simulations.

4.2. Power spectra from photometric datasets

We present the power spectra for simulations with photometric redshift errors estimated using the method described in Sect. 3.1. The noise-subtracted average galaxy auto-power spectra $({\overline{C}}_{\ell }^{\mathit{gg}})$ and cross-power spectra between CMB lensing and galaxy over-density $({\overline{C}}_{\ell }^{\kappa g})$ extracted from photometric datasets for three tomographic bins are shown in Fig. 4. The blue circles and green squares show the reconstructed power spectra for Case-I with σ₀ = 0.02 and σ₀ = 0.05, respectively, and the orange stars show the average power spectra for Case-II with γ₀ = 0.02. The red lines (solid, dashed, and dotted) are the theoretical power spectra computed using the redshift distributions estimated by the more common convolution method (Sect. 3.2.1) for σ₀ = 0.02, σ₀ = 0.05 and γ₀ = 0.02, respectively. We also present the relative difference between the extracted average power spectra and their theoretical expectations. The relative difference between estimated and theoretical power spectra for the other tomographic bins are shown in Appendix C.

Fig. 4. Comparison of estimated power spectra for different photometric redshift errors. Top: average galaxy auto-power spectrum. Bottom: cross-power spectrum reconstructed from 300 realisations of photometric datasets for three tomographic bins. The blue, green, and orange symbols show the reconstructed power spectra for Case-I (σ0 = 0.02,  σ0 = 0.05) and Case-II (γ0 = 0.02), respectively. The solid, dashed, and dotted red lines are the corresponding theoretical power spectra computed with redshift distributions estimated by the convolution method using Eq. (12). The error bars are computed from the covariance matrix of simulations using Eq. (9).

The estimated galaxy auto-power spectra are smaller than expectations in every tomographic bin, with offsets varying between 2 − 15% for σ₀ = 0.02, between 15 − 40% for σ₀ = 0.05, and between 20 − 40% for γ₀ = 0.02. We note that even with a smaller redshift scatter of γ₀ = 0.02, we obtain large offsets in the galaxy auto-power spectra, which is due to the broad tails of the modified Lorentzian error distribution. The cross-power spectra show comparatively smaller biases of < 5% in every tomographic bin for both Case-I and Case-II. We find similar offsets when the true redshift distributions (and hence, the theoretical power spectra) were computed using the deconvolution method (following Eq. (18)). This shows that the offsets in the power spectra are not related to the method of choice. Furthermore, since the level of offsets in the power spectra depends on the strength of photometric redshift scatter (Case-I) and type of redshift error distribution (Case-II), this confirms that the origin of these offsets is rooted in the leakage of objects from one redshift bin to the next due to photometric redshift errors. These deviations will also affect the parameter estimation from the power spectra.

5. Correction for leakage through scattering matrix

As a result of the errors in photometric redshifts, a fraction of galaxies observed in a given photometric redshift bin come from other redshift bins. The leakage of objects across redshift bins changes the strength of the correlation in a tomographic analysis and the results in the non-zero correlation between different redshift bins. In this section, we attempt to counter the redshift bin mismatch of objects through a scattering matrix. When we divide galaxies into n tomographic bins, then Zhang et al. (2010) showed that noise-subtracted galaxy auto-power spectra between the ith and jth photometric bins, $C_{\mathit{ij}}^{gg,\text{ph}}$, are related to the noise-subtracted galaxy auto-power spectra from the kth true redshift bin $C_{\mathit{kk}}^{gg,\text{tr}}$ by

$$\begin{array}{c}\hfill C_{\mathit{ij}}^{gg,\text{ph}}(\ell )={\displaystyle \sum _k}{P}_{\mathit{ki}}{P}_{\mathit{kj}}{C}_{\mathit{kk}}^{gg,\text{tr}}(\ell ),\end{array}$$(22)

when there are no cross correlations between true redshift bins. Equation (22) has a generalised form for the case when true redshift bins have non-zero correlations, but using disjoint true redshift bins significantly reduces the complexity. The elements of the scattering matrix P_ij are defined as the ratio $N_{i\to j}/N_j^{\text{ph}}$, where N_i → j is the number of galaxies coming from the ith true redshift bin to the jth photometric bin, and $N_j^{\text{ph}}$ is the total number of galaxies in the jth photometric bin. This definition also produces a natural normalisation ${\displaystyle \sum _i}{P}_{\mathit{ij}}=1$. A similar relation for the cross-power spectra between the galaxy over-density in redshift bin i and the CMB lensing convergence can be obtained as

$$\begin{array}{c}\hfill C_i^{\kappa g,\text{ph}}(\ell )={\displaystyle \sum _k}{P}_{\mathit{ki}}{C}_{\mathit{kk}}^{\kappa g,\text{tr}}(\ell ).\end{array}$$(23)

When we collect P_ij as elements of the scattering matrix P, then we can compactly write

$$\begin{array}{cc}& C^{gg,\text{ph}}=P^{\text{T}}{C}^{gg,\text{tr}}P\hfill \end{array}$$(24)

$$\begin{array}{cc}& C^{\kappa g,\text{ph}}=P^{\text{T}}{C}^{\kappa g,\text{tr}},\hfill \end{array}$$(25)

where P^T denotes the transpose of matrix P. Equations (24) and (25) show that the redistribution of galaxies across redshift bins due to photometric redshift errors results in a non-trivial relation between the photometric and true power spectra, weighted by the elements of the scattering matrix. Thus, to properly mitigate the effects of leakage on power spectra, a precise estimation of the scattering matrix is necessary.

Zhang et al. (2017) proposed an algorithm for solving problems similar to Eq. (24) based on the non-negative matrix factorization (NMF) method. This algorithm simultaneously approximates the matrices P and C^gg, tr. However, the NMF method proves to be computationally challenging for cases with many tomographic and multipole bins. Here, we propose an alternative method for the fast and efficient computation of the scattering matrix based on the true and photometric redshift distributions. We first estimate the true redshift distribution, $\frac{\mathrm{d}N(z_t)}{\mathrm{d}{z}_t}$, for the entire redshift range following Eq. (16), and then, we use Eq. (18) to compute the redshift distribution $\frac{\mathrm{d}{N}^i(z_p)}{\mathrm{d}{z}_p}$ for every tomographic bin i. The elements of the scattering matrix P_ij can then be computed directly with the relation

$$\begin{array}{c}\hfill P_{\mathit{ij}}=\frac{\underset{z_{min}^j}{\overset{z_{min}^{j+1}}{\int }}\mathrm{d}{z}_p\frac{\mathrm{d}{N}^i}{\mathrm{d}{z}_p}}{\underset{z_{min}^j}{\overset{z_{min}^{j+1}}{\int }}\mathrm{d}{z}_p\frac{\mathrm{d}N}{\mathrm{d}{z}_p}},\end{array}$$(26)

where $\frac{\mathrm{d}N}{\mathrm{d}{z}_p}$ is the observed photometric redshift distribution of the galaxies, and $z_{min}^j$ is the lower limit of the jth redshift bin. Our method of computing the elements of the scattering matrix is significantly faster than the NMF method and is only subject to an accurate estimation of the error distribution and the true redshift distribution.

Figure 5 shows our estimation of the scattering matrix with our proposed method using redshift distributions for the three cases of σ₀ = 0.02, σ₀ = 0.05, and γ₀ = 0.02. The average value of the scattering matrix ⟨P⟩ and its standard deviation σ(P) averaged over 300 simulations is shown in the top and middle panels of Fig. 5, respectively. We note that the scattering matrix elements corresponding to the first true redshift bin have the strongest standard deviation. This behaviour is expected as the objects near z = 0 do not strictly follow convolution, as discussed in Sect. 3.2.2. The accuracy of the estimated scattering matrix can be verified using the true scattering matrix P^True, which is computed based on exactly counting the number of objects that move from one redshift bin to the next in the catalogue generated by the code FLASK. In the bottom panel of Fig. 5, we show the difference between the scattering matrix computed with our method and P^True, averaged over 300 realisations. We find that |P − P^True|< 0.006 for all elements of the scattering matrix. The strongest differences occur in the first and last tomographic bins, that is, near the boundaries of the redshift range we simulated in our analysis. Hence, the overall precision and accuracy in the estimation of the scattering matrix is found to be good and can be used to correct the redshift bin mismatch for the power spectra.

Fig. 5. Performance of the scattering matrix computed with the method described in Sect. 5. Top: average scattering matrix computed from 300 realisations. Middle: its standard deviation. Bottom: difference between the scattering matrix computed from our method and the true scattering matrix computed from simulated catalogues. Left column: scattering matrix computed for σ0 = 0.02. Middle column: scattering matrix computed for σ0 = 0.05. Right column: scattering matrix computed for γ0 = 0.02.

Given an estimate of the true redshift distribution deconvoluted from the observed photometric redshift distribution and scattering matrix, the effect of the redshift bin mismatch on the objects can be corrected for in two ways: either by transforming the true theoretical power spectra C^th,tr into C^th,ph using Eqs. (24) and (25) and comparing them to the estimated photometric power spectra ${\widehat{C}}^{\text{ph}}$, or by inverting Eqs. (24) and (25) to transform the estimated photometric power spectra ${\widehat{C}}^{\text{ph}}$ into true power spectra ${\widehat{C}}^{\text{tr}}$ and comparing them to the theoretical true power spectra C^th,tr. We used the former approach for the figures that compare the power spectra and the latter approach to estimate the parameters.

In Figs. 6–8, we compare for the three tomographic bins shown in Fig. 4 the noise-subtracted average estimated photometric power spectra with the theoretical power spectra corrected for the bin mismatch leakage using Eqs. (24) and (25). The power spectra for other tomographic bins are presented in Appendix D. The theoretical power spectra after correction for leakage agree completely with the estimated power spectra in all bins, except for the first and last tomographic bin. The disparity in the first and last bins results directly from the inaccuracy of the convolution method near the lower and upper bounds of the redshift distribution considered in the analysis. Nevertheless, we note that even for these tomographic bins, the agreement with corresponding theoretical power spectra improves.

Fig. 6. Comparison of estimated angular power spectra with theoretical power spectra after correction for redshift bin mismatch. Top: average galaxy auto-power spectrum. Bottom: cross-power spectrum reconstructed from 300 realisations of photometric datasets with σ0 = 0.02, shown for three tomographic bins. The red line represents the theoretical power spectrum corrected for the redshift bin mismatch using Eqs. (24) and (25). The error bars are computed from the covariance matrix of simulations using Eq. (9).

Fig. 7. Comparison of estimated angular power spectra with theoretical power spectra after correction for redshift bin mismatch. Top: average galaxy auto-power spectrum. Bottom: cross-power spectrum reconstructed from 300 realisations of photometric datasets with σ0 = 0.05, shown for three tomographic bins. The red line represents the theoretical power spectrum corrected for the redshift bin mismatch using Eqs. (24) and (25). The error bars are computed from the covariance matrix of simulations using Eq. (9).

Fig. 8. Comparison of estimated angular power spectra with theoretical power spectra after correction for redshift bin mismatch. Top: average galaxy auto-power spectrum. Bottom: cross-power spectrum reconstructed from 300 realisations of photometric datasets with a modified Lorentzian error distribution, shown for three tomographic bins. The red line represents the theoretical power spectrum corrected for the redshift bin mismatch using Eqs. (24) and (25). The error bars are computed from the covariance matrix of simulations using Eq. (9).

6. Parameter estimation

In previous sections, we observed that the power spectra in every tomographic bin become biased by the leakage of objects across redshift bins, which can be corrected for by an accurate estimation of the scattering matrix. In this section, we study the impact of the leakage on the estimation of the redshift-dependent linear galaxy bias b and on the amplitude of the cross correlation A from tomographic bins, estimated using the maximum likelihood estimation method discussed in Sect. 3.3.

Before accounting for the leakage, we estimated the linear galaxy bias and cross-correlation amplitude for every tomographic bin using the average galaxy power spectra and the average cross-power spectra estimated from the photometric datasets. The theoretical power spectrum templates for a tomographic bin i were computed using Eq. (1) with the redshift distributions given by Eq. (12). To estimate the parameters b and A after correcting for the redshift bin mismatch, we transformed the extracted photometric power spectra (${\overline{C}}^{gg,\text{ph}}$ and ${\overline{C}}^{\kappa g,\text{ph}}$) into true power spectra by inverting Eqs. (24) and (25). The theoretical power spectrum templates for the likelihood estimation after correction for leakage were computed by substituting Eq. (17) for the redshift distribution in Eq. (1). It is important to note that the photometric power spectra in the tomographic analysis are a combination of true power spectra as represented in Eqs. (24) and (25). Thus, the linear galaxy bias in a photometric redshift bin is also a combination of the linear galaxy bias from the true redshift bins. The parameters can also be directly estimated over the estimated photometric power spectra by properly defining the covariance matrix in the likelihood function. However, transforming the estimated photometric power spectra into true power spectra to estimate the parameter avoids complexities based on defining the covariance matrix and also reduces the computation time.

In Figs. 9–11, we compare the posterior distribution of the parameters estimated from the average power spectra with σ₀ = 0.02, σ₀ = 0.05, and γ₀ = 0.02 for three tomographic bins before and after leakage correction. The lighter and darker shaded contour represents the 68% and 95% confidence intervals. The red lines mark the true values of the galaxy bias and cross-correlation amplitude. The best-fit values of the linear galaxy bias and cross-correlation amplitude from all tomographic bins with ±1σ errors, estimated before and after correction for leakage, are quoted in Tables 2–4 for σ₀ = 0.02, σ₀ = 0.05, and γ₀ = 0.02. Column b_true contains the true values of bias for every tomographic bin. The linear galaxy bias for every tomographic bin estimated from photometric datasets is smaller than their expected value for both Case-I and Case-II, whereas the amplitude of cross correlation is consistently higher than the expected value of unity. However, both parameters become consistent with their expected values after we corrected for the redshift bin mismatch of the objects. In Fig. 12 we show the variations in the redshift evolution of the linear galaxy bias parameter due to the effects from the redshift bin mismatch of the objects for Case-I and Case-II. The dashed black line marks the fiducial evolution of the galaxy bias used in our simulations. The linear galaxy bias shows marginal deviations from its true values after correction with the scattering matrix, and the amplitude of the cross correlation is perfectly consistent with its expected value of unity within 1σ.

Fig. 9. Parameter posteriors obtained from the maximum likelihood estimation performed over the average power spectra from 300 simulations with σ0 = 0.02. The top and bottom panels show the posteriors before and after correction for leakage, respectively. The three vertical dashed lines are the median value of the posterior and ±1σ errors. The 68% and 95% confidence contours are shown in darker and lighter shades, respectively. The red lines are the true values of parameters used for simulations.

Fig. 10. Parameter posteriors obtained from the maximum likelihood estimation performed over the average power spectra from 300 simulations with σ0 = 0.05. The top and bottom panels show the posteriors before and after correction for leakage, respectively. The three vertical dashed lines are the median value of the posterior and ±1σ errors. The 68% and 95% confidence contours are shown in darker and lighter shades, respectively. The red lines are the true values of the parameters used for simulations.

Fig. 11. Parameter posteriors obtained from the maximum likelihood estimation performed over the average power spectra from 300 simulations with a modified Lorentzian error distribution. The top and bottom panels show the posteriors before and after the correction for leakage, respectively. The three vertical dashed lines are the median value of the posterior and ±1σ errors. The 68% and 95% confidence contours are shown in darker and lighter shades, respectively. The red lines are the true values of the parameters used for simulations.

Fig. 12. Comparison of the linear galaxy bias parameter estimated from the average of 300 simulations before (blue circles) and after (red squares) the correction for the redshift bin mismatch of objects, with (a) the Gaussian error distribution with σ0 = 0.02, (b) the Gaussian error distribution with σ0 = 0.05, and (c) the modified Lorentzian error distribution with γ0 = 0.02. The error bars on the parameters correspond to a single realisation.

In Fig. 13, we show the relative difference between the mean and the fiducial value (in terms of the standard deviation for a single realisation) of the linear galaxy bias and the amplitude of the cross correlation for σ₀ = 0.02 (left column), σ₀ = 0.05 (middle column), and γ₀ = 0.02 (right column). The blue circles and red squares represent the parameter estimates before and after the correction for leakage, respectively. The error bars on the data points correspond to the average estimated power spectra. The top panel of Fig. 13 shows the relative difference for the linear galaxy bias parameter. When the scatter of objects across redshift bins is not properly accounted for, the galaxy bias can deviate between 5 − 30σ when σ₀ = 0.02, and by 25 − 120σ when σ₀ = 0.05 and γ₀ = 0.02. These large deviations on the linear galaxy bias are visible because the errors from the likelihood estimation (quoted in Tables 2–4) are between 0.15 − 0.5% for a single realisation. This shows that the estimates for the galaxy bias are very tightly constrained. The bottom panel of Fig. 13 shows the difference values for the amplitude of the cross correlation, which can deviate by up to ∼1.2σ with σ₀ = 0.02, by up to ∼3.5σ when σ₀ = 0.05, and by up to ∼3.0σ with γ₀ = 0.02. As clearly conveyed by Fig. 13, the parameters of the linear galaxy bias and of the cross-correlation amplitude become consistent with their expected values after we corrected for the effect of the redshift bin mismatch of the objects through our scattering matrix formalism.

Fig. 13. Comparison of the relative difference for linear galaxy bias (top) and for the cross-correlation amplitude estimated from the average power spectra (bottom) before (blue circles) and after (red squares) the correction for the redshift bin mismatch with (a) a Gaussian error distribution with σ0 = 0.02, (b) a Gaussian error distribution with σ0 = 0.05, and (c) a modified Lorentzian error distribution with γ0 = 0.02. The error bars on the data points correspond to the average spectra, and the standard deviations σ(b) and σ(A) correspond to a single realisation.

7. Note on estimating σ₈

We showed in Sect. 6 that the scatter of objects across redshift bins can lead to a biased parameter estimation and can thus alter our inferences about the cosmological model. In this section, we estimate the effects of leakage on the σ₈ parameter. Peacock & Bilicki (2018) proposed a method for computing the σ₈ parameter from the cross-correlation amplitude using the relation

$$\begin{array}{c}\hfill {\sigma }_8(z)=A(z){\sigma }_{8,0}D(z),\end{array}$$(27)

where σ_8, 0 is the value of the σ₈ parameter at redshift z = 0, and D(z) is the linear growth function given by Eq. (5). We computed the value of σ_8, 0 for our assumed background cosmology using the software CAMB. In Fig. 14, we show the impact of the scattering of objects on the σ₈ parameter computed using Eq. (27) for Case-I (in the left and middle panels) and Case-II (right panel). The dashed black lines are the fiducial evolution of the σ₈ parameter with redshift. We measure a higher than expected value of σ₈ up to ∼1σ for σ₀ = 0.02 before correcting for the redshift bin mismatch of the objects. With σ₀ = 0.05 and γ₀ = 0.02, the σ₈ parameter is biased between 1 − 3σ. As expected, the biases on the amplitude of the cross correlation are directly reflected in the σ₈ parameter. Because of the broad wings of the modified Lorentzian error distribution, the deviation in the σ₈ parameter even for a smaller redshift scatter γ₀ = 0.02 is as strong as for a Gaussian error distribution with a larger redshift scatter, that is, σ₀ = 0.05. Hence, a precise modelling of the error distribution is necessary to quantify the effects of the redshift scatter. We obtain an unbiased estimate of the σ₈ parameter after leakage correction, and thus, it becomes crucial to correct for the mismatch between different redshift bins in a tomographic analysis to obtain unbiased parameter estimates.

Fig. 14. Comparison of the σ8 parameter estimated from the average power spectra of 300 simulations before and after the correction for leakage for (a) a Gaussian error distribution with σ0 = 0.02, (b) a Gaussian error distribution with σ0 = 0.05, and (c) a modified Lorentzian error distribution with γ0 = 0.02. The dashed black line is the fiducial evolution of the σ8 parameter used during simulations. The error bars on the parameters correspond to a single realisation.

8. Summary and conclusions

We presented the tomographic study of cross correlations by performing simulations of the Planck CMB lensing convergence and galaxy density field mimicking properties of the LSST photometric survey. We used the code FLASK to simulate log-normal fields and divided the galaxies into nine redshift bins. We considered Gaussian photometric redshift errors with a standard deviation of σ(z) = 0.02(1 + z) and σ(z) = 0.05(1 + z) and a modified Lorentzian photometric redshift error with a deviation of 0.02(1 + z), but we did not include catastrophic redshift errors or photometric calibration errors, which we keep for future studies. In this way, we generated an ideal observational scenario in our simulations that was free from other systematics, which is crucial for demonstrating the importance of the redshift bin mismatch of the objects. We computed the galaxy auto-power spectrum and the cross-power spectrum between the galaxy over-density and the CMB convergence fields, and we used these power spectra to estimate two parameters, the redshift-dependent linear galaxy bias b, and amplitude of the cross correlation A, employing the maximum likelihood method.

We estimated the true redshift distribution from a simulated photometric redshift distribution by the convolution method (Sect. 3.2.1). In addition, we also estimated the true redshift distribution with the deconvolution method (Sect. 3.2.2). The most important quantity for an accurate recovery of the true redshift distribution is the precise estimation of the error distributions p(z_t − z_p|z_p) for the convolution and p(z_p − z_t|z_t) for the deconvolution methods. We estimated the error functions with sub-percent accuracy (< 0.25%) by fitting a single Gaussian to p(z_p − z_t|z_t) and summing of three Gaussians to p(z_t − z_p|z_p) (shown in Fig. 2) for the Gaussian error distributions. We find that the sum of the Gaussians accurately captures the peculiarities of the error function p(z_t − z_p|z_p) such as non-Gaussian tails and higher peaks in the centre. For the modified Lorentzian error distribution, we fit both p(z_p − z_t|z_t) and p(z_t − z_p|z_p) with a single modified Lorentzian function.

The galaxy auto-power spectra measured from photometric datasets are found to be consistently smaller in every bin than their fiducial predictions. The offsets vary between 2 − 15% for simulations with σ(z) = 0.02(1 + z) and between 15 − 40% for σ(z) = 0.05(1 + z) as well as for the modified Lorentzian error distribution. The measured cross-power spectra are also biased with smaller deviations (< 5%) for both cases. The measured power spectra are inconsistent with their expectations because the objects scatter from one redshift bin to the next due to the photometric redshift errors. This conclusion is consistent with the fact that the deviations are larger in the case of photometric redshift scatter σ(z) = 0.05(1 + z). The offsets in the power spectra also vary when we changed from a Gaussian error distribution to a modified Lorentzian error distribution due to the broad wings of the modified Lorentzian function.

To alleviate the differences in the power spectra, we implemented the scattering matrix approach introduced by Zhang et al. (2010) to counter the effect of the redshift bin mismatch of the objects. The scattering matrix describes the fraction of objects in a photometric redshift bin that comes from different true redshift bins. The power spectra in photometric redshift bins then transform as a linear combination of the power spectra from different true redshift bins (Eqs. (24) and (25)). Zhang et al. (2017) proposed an algorithm based on the non-negative matrix factorisation method to solve similar numerical problems, but this method is computationally challenging for the many data points in the power spectra and the number of tomographic bins. To circumvent these challenges, we proposed an alternative method for a fast and accurate computation of the scattering matrix based on the reconstruction of the true redshift distribution by the deconvolution method (see Sect. 5). We showed in Fig. 5 that our new method for computing the scattering matrix is robust and only proves inefficient in the first tomographic bin because of a cut in the redshift distribution at boundary z = 0. With a precise estimation of the scattering matrix, we corrected the theoretical power spectra for the tomographic bins to compare with the estimated galaxy power spectra from simulated photometric datasets. Figures 6 and 7 showed that the scattering matrix method makes the estimated power spectra consistent with the leakage-corrected theoretical power spectra.

We quantified the impact of the redshift bin mismatch of the objects on the estimation of the linear galaxy bias and on the amplitude of the cross correlation. To estimate the parameters after the correction for leakage, we transformed the estimated photometric power spectra into estimated true power spectra by inverting Eqs. (24) and (25) (as described in Sect. 6). The best-fit values of these parameters before and after the correction for leakage are quoted in Tables 2 and 3 for Case-I (σ(z) = 0.02(1 + z) and σ(z) = 0.05(1 + z)) and in Table 4 for Case-II. Without accounting for the leakage, we estimate lower values for the linear galaxy bias by 5 − 30σ when σ₀ = 0.02, and by 25 − 110σ when σ₀ = 0.05 and γ₀ = 0.02. The amplitude of the cross correlation, on the other hand, is estimated to be higher than its fiducial value of unity up to ∼1.2σ with σ₀ = 0.02, ∼3.5σ when σ₀ = 0.05, and up to ∼3.0σ when γ₀ = 0.02. It is important to note here that the estimates of a lower galaxy bias and higher amplitude are not to be generalised for every tomographic analysis. The offsets of the power spectra and parameters estimated from photometric datasets in a tomographic study strongly depend on the photometric redshift error distributions and on the redshift distribution of the objects. After correcting for leakage by using the scattering matrix, both parameters are very well constrained with their expected values.

The amplitude of the cross correlation is an indicator of the validity of the background cosmological model. Thus, when the bias resulting from photometric redshift errors is not corrected for, incorrect inferences are inevitable when testing cosmological models with tomographic analyses. Other estimators frequently used to test the cosmological models, such as the D_G (Giannantonio et al. 2016) or E_G (Pullen et al. 2016; Zhang et al. 2007) statistics, also employ the ratio of the cross-power spectra to the galaxy auto-power spectra, and they are therefore prone to similar systematics. We studied the relation between the amplitude of the cross correlation and the more familiar σ₈ parameter in Sect. 7. The σ₈ parameter deviates by up to ∼1σ when σ(z) = 0.02(1 + z), and up to ∼3σ with σ(z) = 0.05(1 + z) and with a modified Lorentzian error distribution. We showed that the offsets in the amplitude that arise due to the scatter of the objects are synonymous with the deviations in the σ₈ parameter. With next-generation galaxy surveys such as the Vera C. Rubin Observatory Legacy Survey of Space and Time (LSST; Ivezić et al. 2019; LSST Science Collaboration 2009), Euclid (Laureijs et al. 2011), and the Dark Energy Spectroscopic Instrument (DESI; Dey et al. 2019), the tomographic approach will emerge as a powerful tool for placing stringent constraints on the validity of cosmological models. We therefore propose that the scattering matrix approach developed and presented in this paper be strictly used for future tomographic studies.

Acknowledgments

The authors would like to thank Agnieszka Pollo and Maciej Bilicki for their valuable comments and discussion. CSS thanks Deepika Bollimpalli and Swayamtrupta Panda for discussions on the deconvolution method. The work has been supported by the Polish Ministry of Science and Higher Education grant DIR/WK/2018/12. The authors acknowledge the use of CAMB, HEALPix, EMCEE and FLASK software packages.

References

Appendix A: Validation of the deconvolution method

In this section, we present the performance of the deconvolution method as described in section 3.2.2 using a toy example. We generated a fiducial redshift distribution by random sampling from the LSST photometric redshift distribution profile (Ivezić et al. 2019; LSST Science Collaboration 2009) with a mean redshift 0.9, which we call the true distribution. We convolved the true distribution with a Gaussian distribution with μ = 0.1, σ = 0.02 and call the result the observed distribution. In the left column of Fig. A.1, we show the true and observed distributions by solid blue and red lines, respectively. We show the robustness of our deconvolution method by applying it on unsmoothed distributions. In the right column of Fig. A.1, we compare the true distribution (solid blue line) with the distribution recovered using our deconvolution method (dashed red line). The recovered distribution agrees well with the true redshift distribution, which validates our deconvolution method for reconstructing the true redshift distribution.

Fig. A.1. Performance of the deconvolution method for estimating the true redshift distribution. Left: True (solid blue line) and observed (solid red line) distributions for our toy example. Right: Comparison of the true redshift distribution (solid blue line) with the distribution recovered from the deconvolution method (dashed red line). The recovered distribution agrees excellently with the true redshift distribution.

Appendix B: Power spectra from true datasets

It is important to check for any systematics that may arise when preparing true datasets with the code FLASK. These systematics, if significant, also affect any inferences made from photometric datasets. In Fig. B.1 we show for every tomographic bin the relative difference between noise-subtracted average galaxy auto-power spectra estimated before adding photometric redshift errors and their theoretical expectation with error bars computed from Eq. (9). Fig. B.2 and B.3 present the relative difference for the cross-power spectrum and CMB convergence auto-power spectrum, respectively. The power spectra are consistent with their theoretical expectations in all tomographic bins, and thus, the FLASK simulations are free from any internal systematics.

Fig. B.1. Relative errors on the average galaxy auto-power spectrum reconstructed from 300 simulations generated by the code FLASK without adding photometric redshift errors. The error bars are computed from the covariance matrix of simulations using Eq. (9).

Fig. B.2. Relative errors on the average cross-power spectrum between galaxy over-density and CMB lensing reconstructed from 300 simulations without adding photometric redshift errors. The error bars are computed from the covariance matrix of simulations using Eq. (9).

Fig. B.3. Relative errors on the average CMB convergence auto-power spectrum reconstructed from 300 simulations. The error bars are computed from the covariance matrix of simulations using Eq. (9).

Appendix C: Power spectra from photometric datasets

Fig. C.1. Average galaxy auto-power spectrum reconstructed from 300 simulations after adding photometric redshift errors (complementary to Fig. 4). The error bars are computed from the covariance matrix of simulations using Eq. (9).

Fig. C.2. Average cross-power spectrum between galaxy over-density and CMB convergence reconstructed from 300 simulations after adding photometric redshift errors (complementary to Fig. 4). The error bars are computed from the covariance matrix of simulations using Eq. (9).

Appendix D: Power spectra after correction for leakage

Fig. D.1. Average galaxy auto-power spectrum reconstructed from 300 simulations of photometric datasets with σ0 = 0.02, computed through the scattering matrix (complementary to Fig. 6). The error bars are computed from the covariance matrix of simulations using Eq. (9).

Fig. D.2. Average cross-power spectrum between galaxy over-density and CMB convergence reconstructed from 300 simulations of photometric datasets with σ0 = 0.02, shown for three tomographic bins computed through the scattering matrix (complementary to Fig. 6). The error bars are computed from the covariance matrix of simulations using Eq. (9).

Fig. D.3. Average galaxy auto-power spectrum reconstructed from 300 simulations of photometric datasets with σ0 = 0.05, computed through the scattering matrix (complementary to Fig. 7). The error bars are computed from the covariance matrix of simulations using Eq. (9).

Fig. D.4. Average cross-power spectrum between galaxy over-density and CMB convergence reconstructed from 300 simulations of photometric datasets with σ0 = 0.05, shown for three tomographic bins computed through the scattering matrix (complementary to Fig. 7). The error bars are computed from the covariance matrix of simulations using Eq. (9).

Fig. D.5. Average galaxy auto-power spectrum reconstructed from 300 simulations of photometric datasets with a modified Lorentzian error distribution, computed through the scattering matrix (complementary to Fig. 8). The error bars are computed from the covariance matrix of simulations using Eq. (9).

Fig. D.6. Average cross-power spectrum between galaxy over-density and CMB convergence reconstructed from 300 simulations of photometric datasets with a modified Lorentzian error distribution, shown for three tomographic bins computed through the scattering matrix (complementary to Fig. 8). The error bars are computed from the covariance matrix of simulations using Eq. (9).

All Tables

All Figures

Fig. 1. Effect of photometric redshift scatter on redshift distribution. True LSST redshift distribution (solid blue curve) divided into nine photometric redshift bins (shown by the solid green curves) for Gaussian redshift errors with (a) σ0 = 0.02, (b) σ0 = 0.05, and (c) for modified Lorentzian redshift errors with γ0 = 0.02. The red line denotes the CMB lensing kernel, and the dashed vertical orange lines mark the true redshift bins. The full LSST redshift distribution and CMB lensing kernel are normalised to unit maximum.

In the text

Fig. 2. Modelling the redshift error distribution with parametric functions. (a) Fit made to the error distribution p(zp − zt|zt) with σ0 = 0.02 (solid red line) using a single Gaussian (blue dashed line). (b) Fit made to the error distribution p(zt − zp|zp) with σ0 = 0.02 (solid red line) using a single Gaussian (dotted orange line), and sum of three Gaussians (dashed blue line). (c) Fit made to the error distribution p(zp − zt|zt) with γ0 = 0.02 (solid red line) using a single modified Lorentzian function (dashed blue line).

In the text

Fig. 3. Average relative difference between the redshift distributions reconstructed using the deconvolution method and the true distribution from 300 simulations for (a) Gaussian error distribution with σ0 = 0.02, (b) a Gaussian error distribution with σ0 = 0.05, and (c) a modified Lorentzian error distribution with γ0 = 0.02. The shaded region represents the 1σ deviations in reconstruction between 300 simulations.

In the text

Fig. 4. Comparison of estimated power spectra for different photometric redshift errors. Top: average galaxy auto-power spectrum. Bottom: cross-power spectrum reconstructed from 300 realisations of photometric datasets for three tomographic bins. The blue, green, and orange symbols show the reconstructed power spectra for Case-I (σ0 = 0.02,  σ0 = 0.05) and Case-II (γ0 = 0.02), respectively. The solid, dashed, and dotted red lines are the corresponding theoretical power spectra computed with redshift distributions estimated by the convolution method using Eq. (12). The error bars are computed from the covariance matrix of simulations using Eq. (9).

In the text

Fig. 5. Performance of the scattering matrix computed with the method described in Sect. 5. Top: average scattering matrix computed from 300 realisations. Middle: its standard deviation. Bottom: difference between the scattering matrix computed from our method and the true scattering matrix computed from simulated catalogues. Left column: scattering matrix computed for σ0 = 0.02. Middle column: scattering matrix computed for σ0 = 0.05. Right column: scattering matrix computed for γ0 = 0.02.

In the text

Fig. 6. Comparison of estimated angular power spectra with theoretical power spectra after correction for redshift bin mismatch. Top: average galaxy auto-power spectrum. Bottom: cross-power spectrum reconstructed from 300 realisations of photometric datasets with σ0 = 0.02, shown for three tomographic bins. The red line represents the theoretical power spectrum corrected for the redshift bin mismatch using Eqs. (24) and (25). The error bars are computed from the covariance matrix of simulations using Eq. (9).

In the text

Fig. 7. Comparison of estimated angular power spectra with theoretical power spectra after correction for redshift bin mismatch. Top: average galaxy auto-power spectrum. Bottom: cross-power spectrum reconstructed from 300 realisations of photometric datasets with σ0 = 0.05, shown for three tomographic bins. The red line represents the theoretical power spectrum corrected for the redshift bin mismatch using Eqs. (24) and (25). The error bars are computed from the covariance matrix of simulations using Eq. (9).

In the text

Fig. 8. Comparison of estimated angular power spectra with theoretical power spectra after correction for redshift bin mismatch. Top: average galaxy auto-power spectrum. Bottom: cross-power spectrum reconstructed from 300 realisations of photometric datasets with a modified Lorentzian error distribution, shown for three tomographic bins. The red line represents the theoretical power spectrum corrected for the redshift bin mismatch using Eqs. (24) and (25). The error bars are computed from the covariance matrix of simulations using Eq. (9).

In the text

Fig. 9. Parameter posteriors obtained from the maximum likelihood estimation performed over the average power spectra from 300 simulations with σ0 = 0.02. The top and bottom panels show the posteriors before and after correction for leakage, respectively. The three vertical dashed lines are the median value of the posterior and ±1σ errors. The 68% and 95% confidence contours are shown in darker and lighter shades, respectively. The red lines are the true values of parameters used for simulations.

In the text

Fig. 10. Parameter posteriors obtained from the maximum likelihood estimation performed over the average power spectra from 300 simulations with σ0 = 0.05. The top and bottom panels show the posteriors before and after correction for leakage, respectively. The three vertical dashed lines are the median value of the posterior and ±1σ errors. The 68% and 95% confidence contours are shown in darker and lighter shades, respectively. The red lines are the true values of the parameters used for simulations.

In the text

Fig. 11. Parameter posteriors obtained from the maximum likelihood estimation performed over the average power spectra from 300 simulations with a modified Lorentzian error distribution. The top and bottom panels show the posteriors before and after the correction for leakage, respectively. The three vertical dashed lines are the median value of the posterior and ±1σ errors. The 68% and 95% confidence contours are shown in darker and lighter shades, respectively. The red lines are the true values of the parameters used for simulations.

In the text

Fig. 12. Comparison of the linear galaxy bias parameter estimated from the average of 300 simulations before (blue circles) and after (red squares) the correction for the redshift bin mismatch of objects, with (a) the Gaussian error distribution with σ0 = 0.02, (b) the Gaussian error distribution with σ0 = 0.05, and (c) the modified Lorentzian error distribution with γ0 = 0.02. The error bars on the parameters correspond to a single realisation.

In the text

Fig. 13. Comparison of the relative difference for linear galaxy bias (top) and for the cross-correlation amplitude estimated from the average power spectra (bottom) before (blue circles) and after (red squares) the correction for the redshift bin mismatch with (a) a Gaussian error distribution with σ0 = 0.02, (b) a Gaussian error distribution with σ0 = 0.05, and (c) a modified Lorentzian error distribution with γ0 = 0.02. The error bars on the data points correspond to the average spectra, and the standard deviations σ(b) and σ(A) correspond to a single realisation.

In the text

Fig. 14. Comparison of the σ8 parameter estimated from the average power spectra of 300 simulations before and after the correction for leakage for (a) a Gaussian error distribution with σ0 = 0.02, (b) a Gaussian error distribution with σ0 = 0.05, and (c) a modified Lorentzian error distribution with γ0 = 0.02. The dashed black line is the fiducial evolution of the σ8 parameter used during simulations. The error bars on the parameters correspond to a single realisation.

In the text

Fig. A.1. Performance of the deconvolution method for estimating the true redshift distribution. Left: True (solid blue line) and observed (solid red line) distributions for our toy example. Right: Comparison of the true redshift distribution (solid blue line) with the distribution recovered from the deconvolution method (dashed red line). The recovered distribution agrees excellently with the true redshift distribution.

In the text

	Fig. B.1. Relative errors on the average galaxy auto-power spectrum reconstructed from 300 simulations generated by the code FLASK without adding photometric redshift errors. The error bars are computed from the covariance matrix of simulations using Eq. (9).
In the text

	Fig. B.2. Relative errors on the average cross-power spectrum between galaxy over-density and CMB lensing reconstructed from 300 simulations without adding photometric redshift errors. The error bars are computed from the covariance matrix of simulations using Eq. (9).
In the text

	Fig. B.3. Relative errors on the average CMB convergence auto-power spectrum reconstructed from 300 simulations. The error bars are computed from the covariance matrix of simulations using Eq. (9).
In the text

	Fig. C.1. Average galaxy auto-power spectrum reconstructed from 300 simulations after adding photometric redshift errors (complementary to Fig. 4). The error bars are computed from the covariance matrix of simulations using Eq. (9).
In the text

	Fig. C.2. Average cross-power spectrum between galaxy over-density and CMB convergence reconstructed from 300 simulations after adding photometric redshift errors (complementary to Fig. 4). The error bars are computed from the covariance matrix of simulations using Eq. (9).
In the text

	Fig. D.1. Average galaxy auto-power spectrum reconstructed from 300 simulations of photometric datasets with σ0 = 0.02, computed through the scattering matrix (complementary to Fig. 6). The error bars are computed from the covariance matrix of simulations using Eq. (9).
In the text

	Fig. D.2. Average cross-power spectrum between galaxy over-density and CMB convergence reconstructed from 300 simulations of photometric datasets with σ0 = 0.02, shown for three tomographic bins computed through the scattering matrix (complementary to Fig. 6). The error bars are computed from the covariance matrix of simulations using Eq. (9).
In the text

	Fig. D.3. Average galaxy auto-power spectrum reconstructed from 300 simulations of photometric datasets with σ0 = 0.05, computed through the scattering matrix (complementary to Fig. 7). The error bars are computed from the covariance matrix of simulations using Eq. (9).
In the text

	Fig. D.4. Average cross-power spectrum between galaxy over-density and CMB convergence reconstructed from 300 simulations of photometric datasets with σ0 = 0.05, shown for three tomographic bins computed through the scattering matrix (complementary to Fig. 7). The error bars are computed from the covariance matrix of simulations using Eq. (9).
In the text

	Fig. D.5. Average galaxy auto-power spectrum reconstructed from 300 simulations of photometric datasets with a modified Lorentzian error distribution, computed through the scattering matrix (complementary to Fig. 8). The error bars are computed from the covariance matrix of simulations using Eq. (9).
In the text

	Fig. D.6. Average cross-power spectrum between galaxy over-density and CMB convergence reconstructed from 300 simulations of photometric datasets with a modified Lorentzian error distribution, shown for three tomographic bins computed through the scattering matrix (complementary to Fig. 8). The error bars are computed from the covariance matrix of simulations using Eq. (9).
In the text