© ESO 2020

1. Introduction

The discovery of cosmic shear has provided us with a new and powerful cosmological tool to empirically test the standard model of cosmology and to determine its parameters. Contrary to the analysis of the cosmic microwave background (CMB, e.g., by Planck Collaboration VI 2018), cosmic shear is more sensitive to the properties of the low-redshift large-scale structure and, thus, provides an excellent consistency check for the standard model. Current cosmic shear surveys are particularly sensitive to the parameter $S_8={\sigma }_8\sqrt{{\mathrm{\Omega }}_{\mathrm{m}}/0.3}$, where σ₈ characterizes the normalization of the matter power spectrum and Ω_m is the matter density parameter. Constraints on S₈ from the three current major cosmic shear results are all consistent with the CMB analysis by Planck Collaboration VI (2018). It is interesting to note, however, that they all favor values that are slightly lower than the Planck constraints of S₈ = 0.830 ± 0.013. Hikage et al. (2019) report $S_8=0.{800}_{-0.028}^{+0.029}$ from an analysis of the Subaru Hyper Suprime-Cam survey, Hildebrandt et al. (2020, hereafter H20) obtained $0.{737}_{-0.036}^{+0.040}$ from KiDS+VIKING data, and Troxel et al. (2018) constrain S₈ = 0.782 ± 0.027 using the Dark Energy Survey (DES). Combined analyses of DES and KiDS data (Joudaki et al. 2019; Asgari et al. 2020) result in a ∼3σ tension with the CMB value for S₈. If this tension is not the manifestation of an unaccounted systematic effect, in either the cosmic shear surveys (Mandelbaum 2018) or the Planck mission (Addison et al. 2016), it certainly merits attention. It could be interpreted as a sign of new physics exemplified by massive neutrinos (Battye & Moss 2014), time-varying dark energy, or modified gravity (Planck Collaboration XIV 2016) and coupling within the dark sector (Kumar et al. 2019). It could also, however, prove to be a simple statistical coincidence.

For current cosmic shear surveys, the estimated systematic error is becoming comparable in magnitude to the statistical error, implying that for next-generation surveys, a significant reduction of systematic errors is necessary. With surveys like the Large Synoptic Survey Telescope (LSST, Ivezic et al. 2008) and Euclid (Laureijs et al. 2011) soon to begin, systematic effects in gravitational lensing have received a large amount of attention (see Mandelbaum 2018, and references therein).

In this paper, we focus on systematic effects induced by variation in survey depth that is so far unaccounted for in cosmic shear analyses (Vale et al. 2004). For a survey with a fixed exposure time, varying atmospheric conditions, dithering strategies and galactic extinction all contribute to an inhomogeneous limiting magnitude as a function of sky position. In order to assess the impact of variable depth for current and future surveys, we build an analytical model for the effect based on the survey specifications of the Kilo-Degree Survey (KiDS, Kuijken et al. 2015). To the first order, the depth variation in KiDS can be modeled by a piece-wise constant depth function, which varies between each 1 deg² square pointing. KiDS object detection is defined in the r-band as these images were chosen to be significantly deeper in comparison to the other optical and near-infra red filters. We, therefore, quantify survey depth with the limiting r-band magnitude, as defined in de Jong et al. (2017). We defer the study of multi-band variable depth and its impact on photometric redshift accuracy for a future work.

This work is complementary to the analysis of Guzik & Bernstein (2005), who investigate the effect of a general, position-dependent multiplicative shear bias on the shear power spectrum. In principle, the varying depth of the source galaxy sample that we study here could be recast as a varying effective shear bias. The inclusion of an inhomogeneous distribution of source galaxies has also been explored using mock catalogs of the Subaru Hyper Suprime-Cam Survey (Shirasaki et al. 2019), with a focus on resulting estimation of the cosmic shear covariance matrix.

In Sect. 2, we will introduce two simple toy models to understand this effect and analyze the impact on the cosmic shear power spectrum. In Sect. 3, we will estimate the effect on the shear correlation functions ξ_± using a semi-analytic model. We will present our results in Sect. 4. In Sect. 5, we will discuss our results and comment on the impact of our used simplifications. In the appendices, we present the full derivation of our model for finite field surveys. We assume the standard weak gravitational lensing formalism, a summary of which can be found in Bartelmann & Schneider (2001).

2. Simple, analytic toy models

For our first analysis, we assume that all the matter between the sources and observer is concentrated in a single lens plane of distance D_d from the observer. If we then distribute sources at varying distances D_s, two effects become apparent: firstly, the lensing efficiency D_ds/D_s varies, where D_ds is the distance between the lens plane and the respective source. Secondly, and more importantly, for a more distant source, more matter is concentrated between the source and the observer, leading to a stronger shear signal.

Assuming that the depth and, thus, the source redshift population, only varies between pointings of the camera, an observer will measure a shear signal that is modified by a step-like depth-function, γ^obs(θ) = W(θ)γ(θ), where W is proportional to the mean of the lensing efficiency D_ds/D_s of one pointing and γ denotes the shear that this pointing would experience if it were of the average depth. We can parametrize W as W(θ) = 1 + w(θ). This implies that ⟨w(θ)⟩ = 0 holds, where ⟨ ⋅ ⟩ denotes the average over all pointings.

2.1. Modeling the power spectrum

In our first model, we describe the impact of varying depth on the power spectrum, following the simplifications described above. In accordance with the definition of the shear power spectrum

$$\begin{array}{c}\hfill \langle \widehat{\gamma }(\mathit{\ell }){\widehat{\gamma }}^{\ast }({\mathit{\ell }}^{\mathbf{\prime }})\rangle ={(2\pi )}^2\delta (\mathit{\ell }-{\mathit{\ell }}^{\mathbf{\prime }})P(|\mathit{\ell }|),\end{array}$$(1)

where $\widehat{\gamma }$ denotes the Fourier transform of γ, we define the observed power spectrum via

$$\begin{array}{c}\hfill P^{\mathrm{obs}}(\mathit{\ell })\equiv \frac{1}{{(2\pi )}^2}{\displaystyle \int {\mathrm{d}}^2{\ell }^{\prime }\langle {\widehat{\gamma }}^{\mathrm{obs}}(\mathit{\ell }){\widehat{\gamma }}^{\mathrm{obs}\ast }({\mathit{\ell }}^{\mathbf{\prime }})\rangle .}\end{array}$$(2)

We note that due to the depth-function, both the assumptions of homogeneity and isotropy break down, which means that we can neither assume isotropy in the power spectrum, nor can we assume that $⟨{\widehat{\gamma }}^{\mathrm{obs}}(\mathit{\ell }){\widehat{\gamma }}^{\mathrm{obs}\ast }(\mathit{\ell }\mathbf{\prime })⟩$ vanishes for ℓ ≠ ℓ′. This estimator provides a natural extension to the definition of the regular power spectrum and, in the case of a homogeneous depth distribution, is reduced back to the original estimator. To model a constant depth on each individual pointing, α, we can choose random variables, w_α, that only need to satisfy ⟨w_α⟩ = 0. As we assume an infinite number of pointings, α can assume any two-dimensional integer value ℤ² and we can parametrize w(θ) as

$$\begin{array}{c}\hfill w(\mathit{\theta })={\displaystyle \sum _{\mathit{\alpha }\in {\mathbb{Z}}^2}{w}_{\mathit{\alpha }}\mathrm{\Xi }(\mathit{\theta }-L\mathit{\alpha }),}\end{array}$$(3)

with the box-function

$$\begin{array}{c}\hfill \mathrm{\Xi }(\mathit{\theta })=\{\begin{array}{c}1\mathit{\theta }\in {[-\frac{L}{2},\frac{L}{2}]}^2\hfill \\ 0\mathrm{else}\hfill \end{array},\end{array}$$(4)

where L is the sidelength of one pointing. Following the calculations in Appendix A.1, we derive

$$\begin{array}{c}\hfill P^{\mathrm{obs}}(\mathit{\ell })=P(\mathit{\ell })+⟨w^2⟩{\displaystyle \int \frac{{\mathrm{d}}^2{\mathit{\ell }}^{\mathbf{\prime }}}{{(2\pi )}^2}\widehat{\mathrm{\Xi }}(\mathit{\ell }-{\mathit{\ell }}^{\mathbf{\prime }})P({\mathit{\ell }}^{\mathbf{\prime }}).}\end{array}$$(5)

Here we have denoted $\langle w^2\rangle \equiv \langle w_{\mathit{\alpha }}^2\rangle$ as the dispersion of the depth-function, since the statistical properties of this function do not depend on the pointing α. The Fourier transform of the box function, $\widehat{\mathrm{\Xi }}$, is a 2-dimensional sinc-function (see Appendix A.1). The observed power spectrum, P^obs, is thus composed of the original power spectrum P(ℓ) from Eq. (1), plus a convolution of the power spectrum with a sinc-function, scaling with the variance of the function w(θ).

2.2. Modeling the shear correlation functions

Measures that are more convenient for the inference of cosmological information from observational data are the shear correlation functions ξ_±, which are defined as

$$\begin{array}{c}\hfill {\xi }_{\pm }(\theta )=⟨{\gamma }_{\mathrm{t}}{\gamma }_{\mathrm{t}}⟩(\theta )\pm ⟨{\gamma }_{\times }{\gamma }_{\times }⟩(\theta ).\end{array}$$(6)

Here, γ_t and γ_×, denote the tangential and cross-component of the shear for a galaxy pair with respect to their relative orientation (see Schneider et al. 2002a). The shear correlation functions are the prime estimators to quantify a cosmic-shear signal, since it is simple to include a weighting of the shear measurements into the correlation functions and, contrary to the power spectrum, one does not have to worry about the shape of the survey footprint or masked regions, or model the noise contribution. For this analysis, we follow the assumption that a deeper pointing shows a stronger shear signal γ^obs(θ) = W(θ)γ(θ) as described above. This assumption implies that a higher redshift just increases the amplitude of the shear signal, but as can clearly be seen by inspecting shear correlation functions of different redshift distributions, the change of the signal is extremely scale-dependent and not just a multiplication with a constant factor. In other words, not just the average shear changes as a function of redshift, but also its entire two-point statistics. However, this should serve as a reasonable first approximation for small variations in mean source redshift. Additionally, we assume that a greater depth does not only lead to a stronger average shear, but also to a higher galaxy number density, implying a correlation between those two quantities.

We denote by Nⁱ(θ) the average weighted number of galaxies¹ per pointing in redshift bin i and by Wⁱ(θ) the weighting of average shear. The observed correlation functions ${\xi }_{\pm }^{ij,\text{obs}}(\theta )$ now change from uniform depth, ${\xi }_{\pm }^{ij,\text{uni}}(\theta )$, via

$$\begin{array}{c}\hfill {\xi }_{\pm }^{ij,\text{obs}}(\theta )=\frac{\langle N^i({\mathit{\theta }}^{\mathbf{\prime }})N^j({\mathit{\theta }}^{\mathbf{\prime }}+\mathit{\theta }){\gamma }_{\text{t}}^{i,\text{obs}}({\mathit{\theta }}^{\mathbf{\prime }}){\gamma }_{\text{t}}^{j,\text{obs}}({\mathit{\theta }}^{\mathbf{\prime }}+\mathit{\theta })\rangle }{\langle N^i({\mathit{\theta }}^{\mathbf{\prime }})N^j({\mathit{\theta }}^{\mathbf{\prime }}+\mathit{\theta })\rangle }\\ \pm \frac{\langle N^i({\mathit{\theta }}^{\mathbf{\prime }})N^j({\mathit{\theta }}^{\mathbf{\prime }}+\mathit{\theta }){\gamma }_{\times }^{i,\text{obs}}({\mathit{\theta }}^{\mathbf{\prime }}){\gamma }_{\times }^{j,\text{obs}}({\mathit{\theta }}^{\mathbf{\prime }}+\mathit{\theta })\rangle }{\langle N^i({\mathit{\theta }}^{\mathbf{\prime }})N^j({\mathit{\theta }}^{\mathbf{\prime }}+\mathit{\theta })\rangle }\hfill \end{array}$$(7)

$$=\frac{\langle N^i({\mathit{\theta }}^{\mathbf{\prime }})N^j({\mathit{\theta }}^{\mathbf{\prime }}+\mathit{\theta })W^i({\mathit{\theta }}^{\mathbf{\prime }})W^j({\mathit{\theta }}^{\mathbf{\prime }}+\mathit{\theta })\rangle }{\langle N^i({\mathit{\theta }}^{\mathbf{\prime }})N^j({\mathit{\theta }}^{\mathbf{\prime }}+\mathit{\theta })\rangle }{\xi }_{\pm }^{ij,\text{uni}}(\theta ),$$(8)

where the average, ⟨ ⋅ ⟩, represents both an ensemble average as well as an average over the position θ′. Assuming that the depth of different pointings is uncorrelated, the only important property of a galaxy pair is whether or not they lie in the same pointing. We denote the probability that a random galaxy pair of separation θ lies in the same pointing by E(θ). This function is depicted in Fig. 1, and an analytic expression is derived in Appendix A.2.

Fig. 1. Probability E(θ) that a random pair of galaxies of distance θ lie in the same 1 deg2 pointing.

To compute the modified shear correlation functions, we parametrize the number densities Nⁱ(θ) = ⟨Nⁱ⟩[1 + nⁱ(θ)] and the weight Wⁱ(θ) = 1 + wⁱ(θ) and, as in Eq. (4), interpret nⁱ(θ) as a function with average ⟨nⁱ⟩ = 0 that is constant on each pointing. We can see that ⟨nⁱ(θ′)n^j(θ′ + θ)⟩ = E(θ)⟨nⁱ(θ′)n^j(θ′)⟩ = E(θ)⟨nⁱn^j⟩ holds and compute:

$$\begin{array}{cc}& \frac{\langle N^i({\mathit{\theta }}^{\mathbf{\prime }})N^j({\mathit{\theta }}^{\mathbf{\prime }}+\mathit{\theta })W^i({\mathit{\theta }}^{\mathbf{\prime }})W^j({\mathit{\theta }}^{\mathbf{\prime }}+\mathit{\theta })\rangle }{⟨N^i⟩⟨N^j⟩}\hfill \\ \hfill & =1+\langle n^i{w}^i\rangle +\langle n^j{w}^j\rangle +E(\theta )[\langle n^i{n}^j\rangle +\langle n^i{w}^j\rangle +\langle n^j{w}^i\rangle \hfill \\ \hfill & +\langle w^i{w}^j\rangle +\langle n^i{n}^j{w}^i\rangle +\langle n^i{n}^j{w}^j\rangle +\langle n^i{w}^i{w}^j\rangle +\langle n^j{w}^i{w}^j\rangle \hfill \\ \hfill & +\langle n^i{n}^j{w}^i{w}^j\rangle ].\hfill \end{array}$$(9)

Ignoring correlations higher than second order in nⁱ and wⁱ,² and performing the same calculation for the denominator of Eq. (8) we find

$$\begin{array}{cc}\hfill {\xi }_{\pm }^{ij,\mathrm{obs}}(\theta )& =[1+\langle n^i{w}^i\rangle +\langle n^j{w}^j\rangle +E(\theta )(\langle n^i{n}^j\rangle +\langle n^i{w}^j\rangle \hfill \\ \hfill & +\langle n^j{w}^i\rangle +\langle w^i{w}^j\rangle )]{[1+E(\theta )\langle n^i{n}^j\rangle ]}^{-1}{\xi }_{\pm }^{ij,\mathrm{uni}}(\theta ).\hfill \end{array}$$(10)

A model correlation function for a cosmic shear survey is usually calculated by taking the average redshift distribution of a redshift bin, weighted by the number density. Ignoring that the depth is correlated on scales of one pointing (here at $\theta \le \sqrt{2}°$) is equivalent to setting E(θ) ≡ 0. We note that there is still a correlation between N and W for the same galaxy. Performing the same calculations as above, this yields a relation between the correlation function of uniform depth, ${\xi }_{\pm }^{ij,\text{uni}}$, and the one that is usually modeled, ${\xi }_{\pm }^{ij}$:

$$\begin{array}{c}\hfill {\xi }_{\pm }^{\mathit{ij}}(\theta )=(1+\langle n^i{w}^i\rangle +\langle n^j{w}^j\rangle ){\xi }_{\pm }^{ij,\mathrm{uni}}(\theta ).\end{array}$$(11)

When an observer now calculates the model correlation functions ${\xi }_{\pm }^{ij}$ without accounting for varying depth between pointings, the ratio between modeled and observed correlation functions becomes:

$$\begin{array}{cc}\hfill \frac{{\xi }_{\pm }^{\mathit{ij}}(\theta )}{{\xi }_{\pm }^{ij,\mathrm{obs}}(\theta )}\approx & [1+\langle n^i{w}^i\rangle +\langle n^j{w}^j\rangle +E(\theta )\langle n^i{n}^j\rangle ]\hfill \\ \hfill & \times [1+\langle n^i{w}^i\rangle +\langle n^j{w}^j\rangle +E(\theta )(\langle n^i{n}^j\rangle +\langle n^i{w}^j\rangle \hfill \\ \hfill & {+\langle n^j{w}^i\rangle +\langle w^i{w}^j\rangle )]}^{-1}.\hfill \end{array}$$(12)

It is interesting to note that ${\xi }_{\pm }^{ij}={\xi }_{\pm }^{ij,\text{obs}}$ holds wherever E(θ) = 0, so we expect the observed and the modeled correlation functions to be equivalent on scales where the depth is uncorrelated. One thing left to determine is how to define the weight-function W(θ). For this, we will refer the reader to the beginning of Sect. 4.

3. A semi-analytic model

The previously derived analytic model describes, how varying depth between pointings modifies the correlation function due to the correlation between number density and the average redshift of source galaxies. While this model serves as an intuitive first approximation, it completely ignores any effects from the large scale structure (LSS) between the closest and the most distant galaxy. Therefore, we do not expect this model to yield accurate, quantitative results for cosmic shear surveys.

Below we derive a more sophisticated model that includes the effects of the LSS. While it is computationally more expensive, it improves the accuracy of the model for cosmic shear surveys, which are sensitive to the exact redshift distributions of sources as well as the underlying cosmology.

An inspection of KiDS-data showed that the redshift distribution of sources is highly correlated with the limiting magnitude in the r-band. We thus chose to separate the survey into ten quantiles, sorted by r-band depth, that is, if a pointing had a shallower depth than 90% of the other pointings, it would belong to the first quantile, and so on. For each quantile m and each tomographic redshift bin i we can extract a weighted number of galaxies $N_m^i$ and a source redshift distribution $p_m^i(z)$ following the direct spectroscopic calibration method of H20. In Fig. 2, the average redshift and weighted number of galaxies are plotted for each quantile of each redshift bin, whereas a selection of source redshift distributions is depicted in Fig. 3. A table of the limiting magnitudes for each quantile can be found in Table C.1.

Fig. 2. Weighted number of galaxies N and average redshift ⟨z⟩ in the KiDS+VIKING-450 survey (KV450, Wright et al. 2019) in pointings of different depth for each of the five tomographic bins used in H20. Each color corresponds to one redshift bin of H20. A single point represents one quantile of the respective redshift bin, where the fainter points denote pointings of shallower depth.

Fig. 3. Source redshift distributions $p_m^i(z)$ for a selection of very shallow pointings (blue), average pointings (yellow) and very deep pointings (green). The percentage points in the legend denote to which quantile a pointing belongs, when all are ordered by their depth.

Given the two comoving distance probability distributions of sources, ${\mathcal{L}}_m^i(\chi )$ and ${\mathcal{L}}_n^j(\chi )$, we can compute the shear correlation functions from the underlying matter power spectrum, P_δ(k, χ), via (Kaiser 1992)

$$\begin{array}{cc}& {\xi }_{\pm ,mn}^{\mathit{ij}}(\theta )={\displaystyle {\int }_0^{\infty }\frac{\mathrm{d}\ell \ell }{2\pi }{J}_{0,4}(\ell \theta )P_{\kappa ,mn}^{\mathit{ij}}(\ell ),}\hfill \end{array}$$(13)

$$\begin{array}{cc}& P_{\kappa ,mn}^{\mathit{ij}}(\ell )=\frac{9H_0^4{\mathrm{\Omega }}_{\mathrm{m}}^2}{4c^4}{\displaystyle {\int }_0^{{\chi }_{\mathrm{H}}}\mathrm{d}\chi \frac{g_m^i(\chi )g_n^j(\chi )}{a^2(\chi )}{P}_{\delta }(\frac{\ell }{f_K(\chi )},\chi ),}\hfill \end{array}$$(14)

$$\begin{array}{cc}& g_m^i(\chi )={\displaystyle {\int }_{\chi }^{{\chi }_{\mathrm{H}}}\mathrm{d}{\chi }^{\prime }{\mathcal{L}}_m^i({\chi }^{\prime })\frac{f_K({\chi }^{\prime }-\chi )}{f_K({\chi }^{\prime })}·}\hfill \end{array}$$(15)

Here, J_n denotes the nth order Bessel Functions, f_K(χ) is the comoving angular diameter distance and χ_H is the comoving distance to the horizon. The parameters H₀ and c denote the Hubble constant and the speed of light.

Using Eq. (13), we can compute the model correlation functions, ${\xi }_{\pm ,mn}^{\mathit{ij}}(\theta )$, for each pair of quantiles m, n and redshift bins i, j.³ When measuring the shear correlation functions of a survey, we take the weighted average of tangential and cross shears of all pairs of galaxies (see Hildebrandt et al. 2017). If, for a single pair of galaxies, one galaxy lies in the mth quantile of redshift bin i and the second one lies in the nth quantile of redshift bin j, then their contribution to the observed correlation functions is, on average, ${\xi }_{\pm ,mn}^{\mathit{ij}}(\theta )$. This means that, if we know each of those single correlation functions, we can reconstruct the total correlation functions via a weighted average of the single functions. Formally, we define

$$\begin{array}{c}\hfill {\xi }_{\pm }^{ij,\mathrm{obs}}(\theta )=\frac{{\sum }_{m,n}{\mathcal{P}}_{\mathit{mn}}^{\mathit{ij}}(\theta ){\xi }_{\pm ,mn}^{\mathit{ij}}(\theta )}{{\sum }_{m,n}{\mathcal{P}}_{\mathit{mn}}^{\mathit{ij}}(\theta )},\end{array}$$(16)

where ${\mathcal{P}}_{\mathit{mn}}^{\mathit{ij}}$ is a weighting of the correlation functions, which has to be proportional to the probability that a galaxy pair of separation θ comes from quantiles m and n. In this analysis, we will assume an uncorrelated distribution of depth and neglect boundary effects as well as the sample variance of the depth-distribution between pointings. We will later discuss the validity of these assumptions as well as possible mitigation strategies.

To calculate ${\mathcal{P}}_{\mathit{mn}}^{\mathit{ij}}(\theta )$, we imagine two arbitrary (infinitesimally small) surface elements d²θ₁ and d²θ₂ of separation θ on the sky. For the case m ≠ n, we know that the two galaxies contributing to ${\mathcal{P}}_{\mathit{mn}}^{\mathit{ij}}(\theta )$ have to lie in different pointings, else they would automatically be in the same quantile. The probability that the surface elements are within different pointings is [1 − E(θ)]. Furthermore, the first element d²θ₁ has to lie in quantile m, the probability of which is 1/10. The pointing of the second element d²θ₂ has to be of quantile n; the probability of that is also equal to 1/10. The probability that a galaxy pair populates those surface elements is proportional to the weighted number of galaxies $N_m^i$ and $N_n^j$. We get for n ≠ m:

$$\begin{array}{c}\hfill {\mathcal{P}}_{\mathit{mn}}^{\mathit{ij}}(\theta )=[1-E(\theta )]\frac{1}{100}{N}_m^i{N}_n^j.\end{array}$$(17)

For the calculation of ${\mathcal{P}}_{\mathit{mm}}^{\mathit{ij}}(\theta )$, we have to account for a different possibility: In case that the galaxies lie in the same pointing, they automatically are in the same quantile. We therefore obtain

$$\begin{array}{c}\hfill {\mathcal{P}}_{\mathit{mn}}^{\mathit{ij}}(\theta )=E(\theta )\frac{1}{10}{N}_m^i{N}_m^j{\delta }_{\mathit{mn}}+[1-E(\theta )]\frac{1}{100}{N}_m^i{N}_n^j,\end{array}$$(18)

where δ_mn denotes the Kronecker delta. Inserting this into Eq. (16), we compute

$$\begin{array}{cc}\hfill {\xi }_{\pm ,mn}^{ij,\mathrm{obs}}(\theta )=& \frac{1}{C}{\displaystyle \sum _{m=1}^{10}{N}_m^i}\{E(\theta )N_m^j{\xi }_{\pm ,mm}^{\mathit{ij}}(\theta )\hfill \\ \hfill & +\frac{[1-E(\theta )]}{10}{\displaystyle \sum _{n=1}^{10}{N}_n^j{\xi }_{\pm ,mn}^{\mathit{ij}}(\theta )}\},\hfill \end{array}$$(19)

with the normalization

$$\begin{array}{c}\hfill C={\displaystyle \sum _{m=1}^{10}{N}_m^i[E(\theta )N_m^j+\frac{[1-E(\theta )]}{10}\sum _{n=1}^{10}{N}_n^j].}\end{array}$$(20)

A mathematically more rigorous derivation of this function can be found in Appendix A.3.

Computing this for all five redshift bins of the KV450-survey, forces us to calculate and co-add 1275 correlation functions⁴. Since the variation in depth is a relatively small effect, even tiny numerical errors can add up, skewing the calculations. Additionally, calculating 10³ correlation functions is computationally expensive. However, if we examine Eq. (15), we see that the comoving distance distribution of sources enters linearly. This, in turn, implies that in Eqs. (14) and (13), both source distance distributions enter linearly, meaning that, instead of adding correlation functions, we can add their respective redshift distributions and compute the correlation functions of that. In particular, we can define the combined number of galaxies Nⁱ and average comoving distance probability distribution, ℒⁱ(χ), of tomographic bin i as

$$\begin{array}{c}\hfill N^i\equiv {\displaystyle \sum _m{N}_m^i,}{\mathcal{L}}^i(\chi )=\frac{{\sum }_m{N}_m^i{\mathcal{L}}_m^i(\chi )}{{\sum }_m{N}_m^i}·\end{array}$$(21)

Defining ${\xi }_{\pm }^{ij}$ as the correlation functions between the average comoving distance distributions ℒⁱ(χ) and ℒ^j(χ), we find:

$$\begin{array}{cc}\hfill {\displaystyle \sum _{m,n}{N}_m^i{N}_n^j{\xi }_{\pm ,mn}^{\mathit{ij}}(\theta )=}& \frac{9H_0^4{\mathrm{\Omega }}_{\mathrm{m}}^2}{4c^4}\sum _{m,n}{N}_m^i{N}_n^j{\int }_0^{\infty }\frac{\mathrm{d}\ell \ell }{2\pi }{J}_{0,4}(\ell \theta )\hfill \\ \hfill & \times {\displaystyle {\int }_0^{{\chi }_{\mathrm{H}}}\frac{\mathrm{d}\chi }{a^2(\chi )}{P}_{\delta }(\frac{\ell }{f_K(\chi )},\chi ){\int }_{\chi }^{{\chi }_{\mathrm{H}}}\mathrm{d}{\chi }^{\prime }{\mathcal{L}}_m^i({\chi }^{\prime })\frac{f_K({\chi }^{\prime }-\chi )}{f_K({\chi }^{\prime })}}\hfill \\ \hfill & \times {\displaystyle {\int }_{\chi }^{{\chi }_{\mathrm{H}}}\mathrm{d}{\chi }^{″}{\mathcal{L}}_n^j({\chi }^{″})\frac{f_K({\chi }^{″}-\chi )}{f_K({\chi }^{″})}}\hfill \\ \hfill =& N^i{N}^j\frac{9H_0^4{\mathrm{\Omega }}_{\mathrm{m}}^2}{4c^4}{\displaystyle {\int }_0^{\infty }\frac{\mathrm{d}\ell \ell }{2\pi }{J}_{0,4}(\ell \theta ){\int }_0^{{\chi }_{\mathrm{H}}}\frac{\mathrm{d}\chi }{a^2(\chi )}{P}_{\delta }(\frac{\ell }{f_K(\chi )},\chi )}\hfill \\ \hfill & \times {\displaystyle {\int }_{\chi }^{{\chi }_{\mathrm{H}}}\mathrm{d}{\chi }^{\prime }\frac{{\sum }_m{N}_m^i{\mathcal{L}}_m^i({\chi }^{\prime })}{N^i}\frac{f_K({\chi }^{\prime }-\chi )}{f_K({\chi }^{\prime })}}\hfill \\ \hfill & \times {\displaystyle {\int }_{\chi }^{{\chi }_{\mathrm{H}}}\mathrm{d}{\chi }^{″}\frac{{\sum }_n{N}_n^j{\mathcal{L}}_n^j({\chi }^{″})}{N^j}\frac{f_K({\chi }^{″}-\chi )}{f_K({\chi }^{″})}}\hfill \\ \hfill =& N^i{N}^j{\xi }_{\pm }^{\mathit{ij}}(\theta ).\hfill \end{array}$$(22)

Consequently, we can apply this to Eq. (19), yielding

$$\begin{array}{cc}\hfill {\xi }_{\pm }^{ij,\mathrm{obs}}(\theta )=& \frac{1}{C}\{E(\theta ){\displaystyle \left[\sum _{m=1}^{10}{N}_m^i{N}_m^j{\xi }_{\pm ,mm}^{\mathit{ij}}(\theta )\right]}\hfill \\ \hfill & +\frac{[1-E(\theta )]}{10}{\xi }_{\pm }^{\mathit{ij}}(\theta )N^i{N}^j\}.\hfill \end{array}$$(23)

For each pair of redshift bins we, thus, only have to compute eleven correlation functions, which reduces the number of functions to compute from 1275 to 165.

4. Results

We compare the analytic and semi-analytic models for a variable-depth cosmic shear measurement in a KiDS-like survey. While the application of the semi-analytic method is straightforward, for the analytic method we need to decide how to estimate the weight function W from the given redshift data. Following the separation of a survey into quantiles as in Sect. 3, we define W(θ) ≡ W_n whenever θ is in a pointing of quantile n. For the determination of W_n we test two approaches: As a first method, following Van Waerbeke et al. (2006), Bernardeau et al. (1997), we estimate

$$\begin{array}{c}\hfill W_n\propto {⟨z⟩}_n^{0.85},\end{array}$$(24)

where ⟨z⟩_n is the average redshift of quantile n. As a second method, we define

$$\begin{array}{c}\hfill W_n\propto \sqrt{⟨{\gamma }_{\mathrm{t}}{\gamma }_{\mathrm{t}}⟩({\theta }_{\mathrm{ref}})+⟨{\gamma }_{\times }{\gamma }_{\times }⟩({\theta }_{\mathrm{ref}})}=\sqrt{{\xi }_{+,nn}^{\mathit{ij}}({\theta }_{\mathrm{ref}})},\end{array}$$(25)

where the ${\xi }_{+,nn}^{\mathit{ij}}({\theta }_{\mathrm{ref}})$ denotes the model correlation function defined in Sect. 3, evaluated at a characteristic scale θ_ref, that needs to be chosen.

While the first method suffers from the fact that the power-law index only holds for sources of redshifts 1 ≲ z ≲ 2, the second method is sensitive to the angular range θ_ref, at which the shear correlation functions are evaluated, which is fairly arbitrary. For θ_ref ≈ 11′, which is roughly in the logarithmic middle between the range of the correlation functions, $[0\stackrel{\prime }{.}5,300\prime ]$, the two calibration methods agree. The choice of other values for θ_ref leads to a different amplitude of the change ${\xi }_{\pm }/{\xi }_{\pm }^{\text{obs}}$, but does not affect its shape. Generally, a smaller θ_ref leads to a stronger effect, in particular, the highest amplitude of the change is at ${\theta }_{\mathrm{ref}}=0\stackrel{\prime }{.}5$.

4.1. Effect on the shear correlation functions

In this section, we calculate both the analytic (Eq. (12)) and semi-analytic (Eq. (23)) models for the shear correlation function from a KiDS-like variable depth survey. We adopt the tomographic bins defined in H20 and their resulting best-fit cosmological parameters to present. In Fig. 4, the ratio between our models for the observed correlation functions ${\xi }_{\pm }^{\text{obs}}$, and the standard theoretical prediction that assumes uniform depth are shown. We find that the level of variation in the depth of the KiDS survey increases the amplitude of the observed shear correlation function, on sub-pointing scales, by up to 5% relative to the uniform-depth case.

Fig. 4. Ratio of correlation functions measured for a uniform depth survey, ξ±, and a KiDS-like variable depth survey, ${\xi }_{\pm }^{\text{obs}}$, cross-correlating five tomographic bins (as denoted in the upper left corner of each panel). The upper left triangle depicts the ratios of ξ+, whereas the lower right triangle depicts the ratios of ξ−. Results from mock KiDS-like data (red) can be compared to analytic models from Sect. 2 (average redshift weighting, green), and the semi-analytic model from Sect. 3 (blue solid). Mock data is limited to the angular regime which is not significantly impacted by resolution effects. As the mocks only take galaxies with z <  2 into account, they slightly underestimate the effect. Applying the same redshift-cutoff to the semi-analytical model (blue dashed) yields a near-perfect agreement on sub-pointing scales. Therefore, the seemingly better agreement between the mocks and the analytic method is purely coincidental. The models adopt the best-fit cosmology of H20.

We compare our models to mock KiDS-like data created using a modified version of the Full-sky Lognormal Astro-fields Simulation Kit (FLASK, Xavier et al. 2016; Joachimi et al., in prep.; Lin et al., in prep.). Using FLASK lognormal fields, we generate galaxy mocks with coherent clustering and lensing signals. Adopting a linear relation between the limiting r-band magnitude and the effective number density, fit to the KiDS data, we imprint the variable depth of the full KiDS-1000 footprint (Kuijken et al. 2019) on a HEALPIX⁵ grid of Nside = 4096. Our resolution choice represents a compromise between minimising the mock computation time and maximising the accuracy of the recovered shear signal at small angular scales. With Nside = 4096, ξ₊ is accurate to 7% above ∼1 arcmin, and ξ₋ is accurate to 10% above ∼6 arcmin. For each HEALPIX pixel in the KiDS-1000 data, the limiting r-band magnitude defines the effective number density of sources, and the average source redshift distribution for each tomographic redshift bin (see Fig. 3). In the uniform depth case, redshifts and number densities are sampled from the average of these tomographic sets. The ratio of the lensing signals from the two mocks is computed, averaged over 2000 shape-noise-free realizations, and is shown in Fig. 4 (red).

The results can be seen in Fig. 4. We observe that for high-redshift bins, all methods yield consistent results. For low-redshift bins, there are discrepancies between the different models. However, the average-redshift weighted analytic model, as shown here, is only valid for high redshifts, whereas the auto-correlation-ξ₊-weighted model (not shown) is entirely dependent on the choice of θ_ref. For θ_ref = 11′, both analytic models agree very well for all redshift bins. For ${\theta }_{\mathrm{ref}}=0\stackrel{\prime }{.}5$, the auto-correlation weighted method agrees with the semi-analytic one pretty well for ξ₊. As ξ₋ is affected much stronger by this effect⁶, the analytic method is not able to trace this change for any choice of θ_ref. Furthermore we note, that the effect seems to be strongest in the first redshift bin, which is not surprising, as there the average redshift between pointings varies the most (compare Fig. 2).

The simulations and the models seem to be in relatively good agreement, but there are some differences. It is noticeable that in the simulations, the value ${\xi }_{\pm }/{\xi }_{\pm }^{\text{obs}}$ consistently stays below unity at large scales, which can be attributed to the fact that the depth of different pointings is not completely uncorrelated, as was assumed in the models.

An additional difference between the models and simulations is, that in the models we neglect boundary effects and the sample-variance of the depth-distribution between pointings, whereas the simulations were performed with the KiDS-1000 footprint. In Appendix B, we develop a model to extract the correction ${\xi }_{\pm }^{\mathit{ij}}/{\xi }_{\pm }^{ij,\mathrm{obs}}$ for a specific survey footprint. With this model, we can estimate the impact of a correlated distribution of depth, the sample variance of the depth-distribution and boundary effects. We find that for a square footprint of 450 deg² or 1000 deg² with an uncorrelated depth-distribution, finite field effects are negligible.

In general, the semi-analytic model predicts a stronger effect than the mocks. This is due to the fact that the mocks are subject to a redshift-cutoff at z = 2, meaning that they do not take the high-redshift tail of galaxies into account. This is of particular importance in the first redshift bin, where this feature is especially pronounced (compare Fig. 3). When the same cutoff is applied to the models, the agreement on sub-pointing scales is striking. In particular this means that the mocks slightly underestimate the effect of varying depth.

4.2. Impact on cosmological parameter constraints

As the next step, we assess how the observational depth variations propagate to cosmological parameters inferred from ${\xi }_{\pm }^{ij,\text{obs}}$ compared to ${\xi }_{\pm }^{ij}$. For this test, we choose a fiducial cosmology, Φ, and determine the relative change in Ω_m and σ₈ compared to a reference setup with uniform depth. All other cosmological parameters are kept fixed. First, we compute the reference correlation functions, ${\xi }_{\pm }^{ij}(\theta ,\mathbf{\Phi })$, for each pair of redshift bins i, j using NICAEA as described in Sect. 3. Then we derive the observed correlation functions, ${\xi }_{\pm }^{ij,\text{obs}}(\theta ,\mathbf{\Phi })$, from Eq. (23). Using the Markov-chain Monte Carlo sampler EMCEE, we sample correlation functions ${\xi }_{\pm }^{ij}(\theta ,{\mathbf{\Phi }}^{\prime })$ for different cosmologies Φ′ and find the likelihood distribution, given the data vector ${\xi }_{\pm }^{ij,\text{obs}}(\theta ,\mathbf{\Phi })$ and the covariance-matrix computed in H20. This yields an estimate of the shift in Ω_m and σ₈ introduced by varying depth.

As can be seen in Fig. 5, the impact of varying depth is insignificant compared to the uncertainties for a KiDS-450-like survey. To get a rough estimate for the impact on future surveys, we divide our covariance-matrix by a factor of 30 to model a KiDS-15 000-like survey, which approximately accounts for the increased survey area of LSST and Euclid with respect to KiDS-450. Here the impact on Ω_m, σ₈ and S₈ is significant at the level of approximately 1σ. As our modified covariance-matrix does not account for the factor of ∼4 expected increase in galaxy number density for LSST and Euclid, we note that this is likely to be a lower estimate for the significance of the effect. Even though Euclid is a space-based mission and, therefore, will not suffer from variable atmospheric effects, the key photometric redshift measurement uses data from several ground-based surveys, including KiDS. Placing a selection criteria on redshift estimation success, will therefore lead to depth variations in the source galaxy sample. While the data from LSST will be practically free of variations in depth after 10 years of observations, the first few years data will include significant depth variation. The impact may be even stronger than the KiDS-like analysis presented here as the multi-band KiDS depth variation was minimised using seeing-dependent data acquisition. This is in contrast to the seeing-agnostic multi-band cadence of LSST.

Fig. 5. Recovered cosmological parameters for a variable-depth (Observed) and uniform-depth (Reference) KiDS-450-like (left) and KiDS-15 000-like survey (right). The KiDS-450-like figure was computed using the covariance matrix of H20. For the KiDS-15 000-like survey, we divided the covariance matrix of H20 by 30. This approximately accounts for the increased survey area of next-generation experiments, but does not factor in the increased number density and higher redshifts. Hence, this exercise provides a rough indication of the significance of varying depth effects in stage IV surveys. Both figures were computed using a fiducial cosmology of Ωm = 0.25 and σ8 = 0.85.

Calculating the correction ${\xi }_{\pm }^{ij}/{\xi }_{\pm }^{ij,\text{obs}}$ for varying values of Ω_m and σ₈, reveals a nontrivial dependence on the cosmology, which can be seen in Fig. C.1. For various combinations of Ω_m and σ₈ within the 95% confidence limit of KV450, we report a variation in ${\xi }_{\pm }/{\xi }_{\pm }^{\text{obs}}$ of a few percentage points on small scales.

4.3. Variable depth contribution to B-modes

To check for remaining systematics, a weak lensing signal can be divided into two components, the so-called E- and B-modes (Crittenden et al. 2002; Schneider et al. 2002b). To leading order, B-modes cannot be created by astrophysical phenomena and are thus an excellent test for remaining systematics. Direct E- and B-mode decomposition for cosmic shear surveys can be provided by Complete Orthogonal Sets of E- and B-mode Integrals (COSEBIs, Schneider et al. 2010, hereafter S10), as they can easily be applied to real data. We note that the non-existence of B-modes does not necessarily imply that the sample is free of remaining systematics. To estimate the B-modes created by this effect, we extract the COSEBIs from the correlation functions, ${\xi }_{\pm }^{\text{obs}}$, that have been modified under the semi-analytic model, and from a reference set of correlation functions ξ_±. To be most sensitive to the effect of varying depth, we choose logarithmic COSEBIs with an angular range of ${\theta }_{\min }=0\stackrel{\prime }{.}5$ to θ_max = 72′. As the B-modes of the reference correlation functions are zero, they are a good test for numerical errors in the calculations. Motivated by the discussion in Asgari et al. (2017), we calculated the COSEBIs from correlation functions binned in 400 000 linear bins, which yielded neglible numerical errors. We report a consistent B-mode behaviour across all redshift-bins, which can be seen in Fig. 6. However, we compare the B-modes to the ones measured by Asgari et al. (2020) in KV450, which were consistent with zero. Since the B-modes created by varying depth are smaller than these by a factor of 50, we conclude that this effect cannot create measurable B-modes in the KV450 survey. It should be noted that the difference in E-modes is as large as the B-modes, which suggests that any significant change in the cosmological parameters due to varying depth will also yield a significant detection of B-modes. Additionally, the created pattern is very characteristic, which makes it easy to recognize in a B-mode analysis of an actual survey (see Asgari et al. 2019).

Fig. 6. Difference in the E-modes (top left) and B-modes (bottom right) between the reference and the observed correlation functions. For comparison: Scaled total E-modes of the reference correlation function En and scaled B-modes measured in the KV450 survey $B_n^{\text{KV450}}$. All E- and B-modes were calculated using the logarithmic COSEBIs in S10 for an angular range of ${\theta }_{\min }=0\stackrel{\prime }{.}5$, θmax = 72′.

5. Discussion

With our semi-analytic model we describe the impact of varying depth in ground-based cosmic shear surveys. During our analysis we made several simplifications, which we discuss below.

In the most general terms, we analyze the effects of a position-dependent selection function on cosmic shear surveys. In our analysis, this selection function was governed by the KiDS r-band depth of a pointing. This neglects a number of other effects: The depth in different bands and the seeing of a pointing will also modify the number densities and redshift distributions on the scale of a pointing (although those variations are also correlated with r-band depth and thus at least partly accounted for), whereas dithering strategies as well as imperfections in the telescope and CCD cause modifications on sub-pointing scales. However, several tests showed that these effects are subdominant compared to the variations caused by the r-band depth.

We assumed an uncorrelated distribution of the depth-function and neglected boundary effects as well as the sample variance of the depth-distribution between pointings. While the boundary effects arising from a finite survey footprint have a small impact on the shape of the function E(θ)⁷, the governing factor is the sample variance of the depth-distribution. We assumed that the probability for any pointing to be in quantile n is exactly the expectation value, namely 1/10. While this would be true for an infinitely large survey with an uncorrelated distribution of the depth-function, it does not necessarily hold for a real survey. However, our analysis in Appendix B suggests that these effects are not significant for the KV450 survey. In the models, we also assumed an uncorrelated distribution of the depth-function. As can be seen in Fig. 4, this approximation introduces a small error when compared to the simulations.

In our MCMC runs, we did not account for degeneracies with other cosmological parameters or observational effects. In particular, intrinsic alignments and baryon feedback also modify the correlation functions primarily on small scales, so they are likely to be degenerate with the effect of varying depth (Troxel & Ishak 2015). In an MCMC run that includes these nuisance parameters, we suspect that the parameters for intrinsic alignments and baryon feedback change to mitigate this effect, so that the impact on cosmological parameters will be smaller than in our results.

Despite these repercussions, we are confident to say that the effects of varying depth are not significant for the KV450 survey. In particular, this means that a varying depth cannot explain the tension between observations of the low-redshift Universe and results from analysis of the CMB.

For next-generation surveys like Euclid and LSST, we have demonstrated that if variable depth is unaccounted for in the analysis, the resulting bias is likely to be significant. A detailed LSST and Euclid study that uses a realistic variation of depth should therefore be conducted. If these studies reach a similar conclusion, variable-depth bias could be circumvented in likelihood analyses by including a cosmology-dependent correction for this effect using the semi-analytical model presented in this paper (see Fig. C.1).

Although the analytic model (Eq. (12)) does not fully describe the effect, it can be used as a valid approximation to estimate the importance of varying depth for an arbitrary survey: As can be seen in Fig. 2, in the KiDS-survey the characteristic changes both in redshift and number density for the third redshift bin are about 0.1. Setting ⟨nⁱn^j⟩=⟨wⁱw^j⟩=⟨nⁱw^j⟩ = 0.01 yields ${\xi }_{\pm }^{ij}/{\xi }_{\pm }^{ij,\text{obs}}(0)\approx 0.97$, which roughly agrees with the actual results (see Fig. 4). For a survey with only half the variation in depth, one gets ${\xi }_{\pm }^{ij}/{\xi }_{\pm }^{ij,\text{obs}}(0)\approx 0.993$. This method allows us to estimate a threshold for an acceptable variation of depth, given a required precision for the shear correlation functions.

Additionally, it is interesting to note that E(θ) is the azimuthal average of the function E(θ) derived in Appendix A.2, which is not isotropic. Therefore, it would be possible to observe a direction-dependent correlation function ${\xi }_{\pm }^{ij,\text{obs}}(\mathit{\theta })$ in future surveys. An anisotropy in the observed correlation function could be a sign for the influence of varying depth.

The variations in depth will also affect the covariance of a survey, both because they modify the signal and because they introduce an additional term of sample variance in terms of the distribution of depth. This effect will be investigated in a forthcoming publication (Joachimi et al., in prep.; Lin et al., in prep.).

Acknowledgments

We thank the anonymous referee for very constructive and helpful comments. The results in this paper are based on observations made with ESO Telescopes at the La Silla Paranal Observatory under programme IDs 177.A-3016, 177.A-3017, 177.A-3018 and 179.A-2004, and on data products produced by the KiDS consortium. The KiDS production team acknowledges support from: Deutsche Forschungsgemeinschaft, ERC, NOVA and NWO-M grants; Target; the University of Padova, and the University Federico II (Naples). We acknowledge support from the European Research Council under grant numbers 770935 (HH,JLvdB) and 647112 (CH,MA,CL). SH acknowledges support from the German Research Foundation (DFG SCHN 342/13). H. Hildebrandt is supported by a Heisenberg grant of the Deutsche Forschungsgemeinschaft (Hi 1495/5-1). CH acknowledges support from the Max Planck Society and the Alexander von Humboldt Foundation in the framework of the Max Planck-Humboldt Research Award endowed by the Federal Ministry of Education and Research. KK acknowledges support by the Alexander von Humboldt Foundation. TT acknowledges funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 797794. Some of the results in this paper have been derived using the HEALPix (Górski et al. 2005) package. This research has made use of NASA’s Astrophysics Data System and adstex (https://github.com/yymao/adstex). Author contributions: all authors contributed to the development and writing of this paper. The authorship list is given in two groups: the lead authors (SH, PS, HHi) followed by key contributors to the scientific analysis in alphabetical order.

References

Appendix A: Detailed calculations

A.1 Calculation of the power spectrum

Here we perform the calculation for the observed power spectrum P^obs(ℓ). For this, we assume an infinitely large field in order to perform our integration over ℝ². In reality, finite field effects would play a role. We begin with the calculation of the correlation for the Fourier transformed shear:

$$\begin{array}{cc}& \langle {\widehat{\gamma }}^{\mathrm{obs}}(\mathit{\ell }){\widehat{\gamma }}^{\mathrm{obs}\ast }({\mathit{\ell }}^{\mathbf{\prime }})\rangle \hfill \\ \hfill & ={\displaystyle \langle \int {\mathrm{d}}^2\theta \int {\mathrm{d}}^2{\theta }^{\prime }W(\mathit{\theta })W({\mathit{\theta }}^{\mathbf{\prime }})\gamma (\mathit{\theta }){\gamma }^{\ast }({\mathit{\theta }}^{\mathbf{\prime }})exp(i\mathit{\ell }\mathit{\theta }-i{\mathit{\ell }}^{\mathbf{\prime }}{\mathit{\theta }}^{\mathbf{\prime }})\rangle }\hfill \\ \hfill & ={\displaystyle \langle \int {\mathrm{d}}^2\theta \int {\mathrm{d}}^2{\theta }^{\prime }W(\mathit{\theta })W({\mathit{\theta }}^{\mathbf{\prime }})exp(i\mathit{\ell }\mathit{\theta }-i{\mathit{\ell }}^{\mathbf{\prime }}{\mathit{\theta }}^{\mathbf{\prime }})\int \frac{{\mathrm{d}}^{2}k}{{(2\pi )}^2}\int \frac{{\mathrm{d}}^2{k}^{\prime }}{{(2\pi )}^2}\widehat{\gamma }(\mathit{k}){\widehat{\gamma }}^{\ast }({\mathit{k}}^{\mathbf{\prime }})exp(-i\mathit{k}\mathit{\theta }+i{\mathit{k}}^{\mathbf{\prime }}{\mathit{\theta }}^{\mathbf{\prime }})\rangle }\hfill \\ \hfill & ={\displaystyle \langle \int {\mathrm{d}}^2\theta \int {\mathrm{d}}^2{\theta }^{\prime }\int \frac{{\mathrm{d}}^{2}k}{{(2\pi )}^2}\int \frac{{\mathrm{d}}^2{k}^{\prime }}{{(2\pi )}^2}P(\mathit{k}){(2\pi )}^2\delta (\mathit{k}-{\mathit{k}}^{\mathbf{\prime }})exp[i(\mathit{\ell }\mathit{\theta }-{\mathit{\ell }}^{\mathbf{\prime }}{\mathit{\theta }}^{\mathbf{\prime }}-\mathit{k}\mathit{\theta }+{\mathit{k}}^{\mathbf{\prime }}{\mathit{\theta }}^{\mathbf{\prime }})]W(\mathit{\theta })W({\mathit{\theta }}^{\mathbf{\prime }})\rangle }\hfill \\ \hfill & ={\displaystyle \langle \int \frac{{\mathrm{d}}^{2}k}{{(2\pi )}^2}P(\mathit{k})\int {\mathrm{d}}^2\theta W(\mathit{\theta })exp[i\mathit{\theta }(\mathit{\ell }-\mathit{k})]\int {\mathrm{d}}^2{\theta }^{\prime }W({\mathit{\theta }}^{\mathbf{\prime }})exp[-i{\mathit{\theta }}^{\mathbf{\prime }}({\mathit{\ell }}^{\mathbf{\prime }}-\mathit{k})]\rangle }\hfill \\ \hfill & ={\displaystyle \langle \int \frac{{\mathrm{d}}^{2}k}{{(2\pi )}^2}P(\mathit{k})\widehat{W}(\mathit{\ell }-\mathit{k}){\widehat{W}}^{\ast }({\mathit{\ell }}^{\mathbf{\prime }}-\mathit{k})\rangle .}\hfill \end{array}$$(A.1)

It is important to keep in mind, that the ensemble averages of the weight function are independent of the ensemble averages of the shear values, meaning ⟨W(θ)γ(θ)⟩ = ⟨W(θ)⟩⟨γ(θ)⟩. We can define W(θ) = 1 + w(θ) with ⟨w(θ)⟩ = 0, which leads to the expression

$$\begin{array}{cc}& \langle {\widehat{\gamma }}^{\mathrm{obs}}(\mathit{\ell }){\widehat{\gamma }}^{\mathrm{obs}\ast }({\mathit{\ell }}^{\mathbf{\prime }})\rangle \hfill \\ \hfill & ={\displaystyle \langle \int \frac{{\mathrm{d}}^{2}k}{{(2\pi )}^2}P(\mathit{k})\{{(2\pi )}^4\delta (\mathit{\ell }-\mathit{k})\delta ({\mathit{\ell }}^{\mathbf{\prime }}-\mathit{k})+{(2\pi )}^2[\widehat{w}(\mathit{\ell }-\mathit{k})\delta ({\mathit{\ell }}^{\mathbf{\prime }}-\mathit{k})+{\widehat{w}}^{\ast }({\mathit{\ell }}^{\mathbf{\prime }}-\mathit{k})\delta (\mathit{\ell }-\mathit{k})]+\widehat{w}(\mathit{\ell }-\mathit{k})\widehat{w}({\mathit{\ell }}^{\mathbf{\prime }}-\mathit{k})\}\rangle }\hfill \\ \hfill & ={(2\pi )}^2\delta (\mathit{\ell }-{\mathit{\ell }}^{\mathbf{\prime }})P(\mathit{\ell })+[⟨\widehat{w}(\mathit{\ell }-{\mathit{\ell }}^{\mathbf{\prime }})⟩P({\mathit{\ell }}^{\mathbf{\prime }})+⟨{\widehat{w}}^{\ast }({\mathit{\ell }}^{\mathbf{\prime }}-\mathit{\ell })⟩P(\mathit{\ell })]+{\displaystyle \langle \int \frac{{\mathrm{d}}^{2}k}{{(2\pi )}^2}\widehat{w}(\mathit{\ell }-\mathit{k}){\widehat{w}}^{\ast }({\mathit{\ell }}^{\mathbf{\prime }}-\mathit{k})P(\mathit{k})\rangle }\hfill \\ \hfill & ={(2\pi )}^2\delta (\mathit{\ell }-{\mathit{\ell }}^{\mathbf{\prime }})P(\mathit{\ell })+{\displaystyle \langle \int \frac{{\mathrm{d}}^{2}k}{{(2\pi )}^2}\widehat{w}(\mathit{\ell }-\mathit{k}){\widehat{w}}^{\ast }({\mathit{\ell }}^{\mathbf{\prime }}-\mathit{k})P(\mathit{k})\rangle ,}\hfill \end{array}$$(A.2)

where in the final step we have used that the average $⟨\widehat{\mathit{w}}(\mathit{\ell })⟩$ vanishes. Up until now, we have not specified our weight-function w. We parametrize it as

$$\begin{array}{c}\hfill w(\mathit{\theta })={\displaystyle \sum _{\mathit{\alpha }\in {\mathbb{Z}}^2}{w}_{\mathit{\alpha }}\mathrm{\Xi }(\mathit{\theta }-L\mathit{\alpha }),\mathrm{with}\mathrm{the}\mathrm{box}\text{-}\mathrm{function}\mathrm{\Xi }(\mathit{\theta })=\{\begin{array}{c}1\mathit{\theta }\in {[-\frac{L}{2},\frac{L}{2}]}^2\hfill \\ 0\mathrm{else}\hfill \end{array}.}\end{array}$$(A.3)

Here, the w_α are random variables, drawn from the random distribution describing the survey depths. For the Fourier transform we compute:

$$\begin{array}{c}\hfill \widehat{w}(\mathit{\ell })={\displaystyle \sum _{\mathit{\alpha }\in {\mathbb{Z}}^2}{w}_{\mathit{\alpha }}exp(-iL\mathit{\ell }·\mathit{\alpha })\widehat{\mathrm{\Xi }}(\mathit{\ell }),}\end{array}$$(A.4)

where

$$\begin{array}{c}\hfill \widehat{\mathrm{\Xi }}(\mathit{\ell })=\frac{4sin\left(\frac{L{\ell }_1}{2}\right)sin\left(\frac{L{\ell }_2}{2}\right)}{{\ell }_1{\ell }_2},\end{array}$$(A.5)

is a 2-dimensional sinc function. Assuming an uncorrelated weight-distribution (⟨w_αw_β⟩ = 0 for α ≠ β) and setting $\langle w^2\rangle \equiv \langle w_{\mathit{\alpha }}^2\rangle$ for each α, we get

$$\begin{array}{cc}& {\displaystyle \langle \int \frac{{\mathrm{d}}^{2}k}{{(2\pi )}^2}\widehat{w}(\mathit{\ell }-\mathit{k}){\widehat{w}}^{\ast }({\mathit{\ell }}^{\mathbf{\prime }}-\mathit{k})P(\mathit{k})\rangle }\hfill \\ \hfill & ={\displaystyle \langle \int \frac{{\mathrm{d}}^{2}k}{{(2\pi )}^2}\sum _{\mathit{\alpha },\mathit{\beta }}{w}_{\mathit{\alpha }}{w}_{\mathit{\beta }}exp[-iL(\mathit{\ell }-\mathit{k})·\mathit{\alpha }]\widehat{\mathrm{\Xi }}(\mathit{\ell }-\mathit{k})exp[iL({\mathit{\ell }}^{\mathbf{\prime }}-\mathit{k})·\mathit{\beta }]{\widehat{\mathrm{\Xi }}}^{\ast }({\mathit{\ell }}^{\mathbf{\prime }}-\mathit{k})P(\mathit{k})\rangle }\hfill \\ \hfill & ={\displaystyle \int \frac{{\mathrm{d}}^{2}k}{{(2\pi )}^2}\sum _{\mathit{\alpha }}⟨w^2⟩exp[-iL(\mathit{\ell }-\mathit{k})·\mathit{\alpha }+iL({\mathit{\ell }}^{\mathbf{\prime }}-\mathit{k})·\mathit{\alpha }]\widehat{\mathrm{\Xi }}(\mathit{\ell }-\mathit{k}){\widehat{\mathrm{\Xi }}}^{\ast }({\mathit{\ell }}^{\mathbf{\prime }}-\mathit{k})P(\mathit{k}).}\hfill \end{array}$$(A.6)

Using this result, we can obtain the observed power spectrum

$$\begin{array}{c}\hfill P^{\mathrm{obs}}(\mathit{\ell })=\frac{1}{{(2\pi )}^2}{\displaystyle \int {\mathrm{d}}^2{\ell }^{\prime }\langle {\widehat{\gamma }}^{\mathrm{obs}}(\mathit{\ell }){\widehat{\gamma }}^{\mathrm{obs}\ast }({\mathit{\ell }}^{\mathbf{\prime }})\rangle ,}\end{array}$$(A.7)

by performing the ℓ′-integration in (A.2):

$$\begin{array}{cc}\hfill P^{\mathrm{obs}}(\mathit{\ell })=& P(\mathit{\ell })+{\displaystyle \int \frac{{\mathrm{d}}^2{\mathit{\ell }}^{\mathbf{\prime }}}{{(2\pi )}^2}\int \frac{{\mathrm{d}}^2\mathit{k}}{{(2\pi )}^2}\sum _{\mathit{\alpha }}⟨w^2⟩exp[-iL(\mathit{\ell }-\mathit{k})·\mathit{\alpha }+iL({\mathit{\ell }}^{\mathbf{\prime }}-\mathit{k})·\mathit{\alpha }]\widehat{\mathrm{\Xi }}(\mathit{\ell }-\mathit{k})\widehat{\mathrm{\Xi }}({\mathit{\ell }}^{\mathbf{\prime }}-\mathit{k})P(\mathit{k})}\hfill \\ \hfill =& P(\mathit{\ell })+{\displaystyle \int \frac{{\mathrm{d}}^2\mathit{k}}{{(2\pi )}^2}\sum _{\mathit{\alpha }}⟨w^2⟩exp[-iL(\mathit{\ell }-\mathit{k})·\mathit{\alpha }]\widehat{\mathrm{\Xi }}(\mathit{\ell }-\mathit{k})P(k)\int \frac{{\mathrm{d}}^2\mathit{\ell }}{{(2\pi )}^2}{\widehat{\mathrm{\Xi }}}^{\ast }({\mathit{\ell }}^{\mathbf{\prime }}-\mathit{k})exp[iL({\mathit{\ell }}^{\mathbf{\prime }}-\mathit{k})·\mathit{\alpha }]}\hfill \\ \hfill =& P(\mathit{\ell })+⟨w^2⟩{\displaystyle \int \frac{{\mathrm{d}}^2\mathit{k}}{{(2\pi )}^2}\widehat{\mathrm{\Xi }}(\mathit{\ell }-\mathit{k})P(\mathit{k})\sum _{\mathit{\alpha }}exp[-iL(\mathit{\ell }-\mathit{k})·\mathit{\alpha }]\mathrm{\Xi }(L\mathit{\alpha })}\hfill \\ \hfill =& P(\mathit{\ell })+⟨w^2⟩{\displaystyle \int \frac{{\mathrm{d}}^2\mathit{k}}{{(2\pi )}^2}\widehat{\mathrm{\Xi }}(\mathit{\ell }-\mathit{k})P(\mathit{k}),}\hfill \end{array}$$(A.8)

which is a convolution of the power spectrum and the 2-dimensional sinc function (A.5). We note that due to the statistical inhomogeneity of the field, many usually adapted conventions fail. In particular, ⟨γ(θ)γ^*(θ′)⟩ does not only depend on the separation vector θ′ − θ, but also on the position θ. For example, the Fourier transform of the observed power spectrum yields ⟨γ(0)γ^*(θ)⟩, but not the shear correlation function.

A.2 The function E(θ)

When computing the shear correlation between a pair of galaxies, it is of central importance whether or not those two galaxies lie in the same pointing. We want to model the probability that a pair of galaxies with separation θ lies in the same pointing by the function E(θ), which we will derive here: Given one square field of length L (in our case L = 60′) and a separation vector θ = (θ₁, θ₂), without loss of generality we can assume θ₁, θ₂ ≥ 0. The dashed square in Fig. A.1 represents all possible positions that the first galaxy can take, such that the second galaxy is still within the same pointing. The volume of this square equals

$$\begin{array}{c}\hfill V(|\mathit{\theta }|,\varphi )=[L-|\mathit{\theta }|cos(\varphi )][L-|\mathit{\theta }|sin(\varphi )],\end{array}$$(A.9)

Fig. A.1. Graphic representation on how to obtain the function E(θ). For a separation vector θ, the dashed square represents the area of galaxies that have their partner in the same pointing.

where ϕ represents the polar angle of the vector θ. The function E(θ) then simply equals V(|θ|,ϕ)/L². To exclude negative volumes (which could occur when |θ| > 1 holds), we need to add the Heaviside step function ℋ:

$$\begin{array}{cc}\hfill E(\mathit{\theta })=& [1-\frac{|\mathit{\theta }|}{L}cos(\varphi )][1-\frac{|\mathit{\theta }|}{L}sin(\varphi )]\mathcal{H}[1-\frac{|\mathit{\theta }|}{L}cos(\varphi )]\times \mathcal{H}[1-\frac{|\mathit{\theta }|}{L}sin(\varphi )].\hfill \end{array}$$(A.10)

As E(θ) is not isotropic, in order to obtain the function E(θ) = E(|θ|), we need to calculate the azimuthal average of Eq. (A.10) over all angles ϕ. While the case θ₁, θ₂ ≥ 0 certainly does not hold for all angles ϕ, we can omit the other cases by making use of the symmetry of the problem.

$$\begin{array}{cc}\hfill E(\theta )& =\frac{4}{2\pi }{\int }_0^{\frac{\pi }{2}}\mathrm{d}\varphi E(\mathit{\theta })=\frac{2}{\pi }\{\begin{array}{cc}{\int }_0^{\frac{\pi }{2}}\mathrm{d}\varphi [1-\frac{|\mathit{\theta }|}{L}cos(\varphi )][1-\frac{|\mathit{\theta }|}{L}sin(\varphi )],\hfill & |\mathit{\theta }|\le L\hfill \\ {\int }_{arccos(L/|\mathit{\theta }|)}^{arcsin(L/|\mathit{\theta }|)}\mathrm{d}\varphi [1-\frac{|\mathit{\theta }|}{L}cos(\varphi )][1-\frac{|\mathit{\theta }|}{L}sin(\varphi )],\hfill & L\le |\mathit{\theta }|\le \sqrt{2}L\hfill \\ 0,\hfill & \sqrt{2}L\le \theta \hfill \end{array}\hfill \\ \hfill & =\{\begin{array}{cc}\frac{1}{L^2\pi }[L^2\pi -(4L-\theta )\theta ],\hfill & \theta \le L\hfill \\ \frac{2}{\pi }[4\sqrt{\frac{{\theta }^2}{L^2}-1}-1-\frac{{\theta }^2}{2L^2}-arccos\left(\frac{L}{\theta }\right)+arcsin\left(\frac{L}{\theta }\right)],\hfill & L\le \theta \le \sqrt{2}L\hfill \\ 0,\hfill & \sqrt{2}L\le \theta \hfill \end{array}.\hfill \end{array}$$(A.11)

A.3 Calculation of the shear correlation functions

Given a set of galaxies, we calculate the shear correlation function ${\xi }_{+}^{ij}$ via

$$\begin{array}{c}\hfill {\xi }_{+}^{\mathit{ij}}(\theta )=\frac{{\sum }_{a,b}{w}_a^i{w}_b^j{ϵ}_a^i{ϵ}_b^{j\ast }\mathrm{\Delta }(|{\mathit{\theta }}_a^i-{\mathit{\theta }}_b^i|)}{{\sum }_{a,b}{w}_a^i{w}_b^j\mathrm{\Delta }(|{\mathit{\theta }}_a^i-{\mathit{\theta }}_b^i|)}·\end{array}$$(A.12)

Here, w represents the lensing weight of the galaxy, whereas ϵ is its (complex) ellipticity and θ its position on the sky. We have defined the function Δ as

$$\begin{array}{c}\hfill \mathrm{\Delta }(|{\mathit{\theta }}_a^i-{\mathit{\theta }}_b^i|)=\{\begin{array}{cc}1,\hfill & |{\mathit{\theta }}_a^i-{\mathit{\theta }}_b^j|\in [\theta ,\theta +\mathrm{d}\theta ]\hfill \\ 0,\hfill & \mathrm{else}\hfill \end{array},\end{array}$$(A.13)

where we assume dθ ≪ θ. We define 𝒩 as the number of pointings in the survey and $F_k^i$ as the set of galaxies in pointing k and tomographic redshift bin i. The numerator in Eq. (A.12) then transforms to:

$$\begin{array}{cc}& {\displaystyle \sum _{k,\ell =1}^{\mathcal{N}}\sum _{a\in F_k^i}\sum _{b\in F_{\ell }^j}{w}_a^i{w}_b^j{ϵ}_a^i{ϵ}_b^{j\ast }\mathrm{\Delta }(|{\mathit{\theta }}_a^i-{\mathit{\theta }}_b^i|)=\sum _{k=1}^{\mathcal{N}}\sum _{a\in F_k^i}{w}_a^i\sum _{\ell =1}^{\mathcal{N}}\sum _{b\in F_{\ell }^j}{w}_b^j\mathrm{\Delta }(|{\mathit{\theta }}_a^i-{\mathit{\theta }}_b^i|){ϵ}_a^i{ϵ}_b^{j\ast }}\hfill \\ \hfill =& {\displaystyle \sum _{k=1}^{\mathcal{N}}\sum _{a\in F_k^i}{w}_a^i[\sum _{b\in F_k^j}{w}_b^j\mathrm{\Delta }(|{\mathit{\theta }}_a^i-{\mathit{\theta }}_b^i|){ϵ}_a^i{ϵ}_b^{j\ast }+\sum _{\ell \ne k}\sum _{b\in F_{\ell }^j}{w}_b^j\mathrm{\Delta }(|{\mathit{\theta }}_a^i-{\mathit{\theta }}_b^i|){ϵ}_a^i{ϵ}_b^{j\ast }].}\hfill \end{array}$$(A.14)

When we denote the probability that pointing k is of quantile m by ${\mathcal{P}}_m^k$ and assume that the product ${ϵ}_a^i{ϵ}_b^{j*}$ always equals its expectation value, we can set the numerator as

$$\begin{array}{cc}& {\displaystyle \sum _{k=1}^{\mathcal{N}}\sum _{a\in F_k^i}{w}_a^i\sum _m{\mathcal{P}}_m^k[\stackrel{(\mathrm{A}15.\mathrm{a})}{\overbrace{{\sum }_{b\in F_k^j}{w}_b^j\mathrm{\Delta }(|{\mathit{\theta }}_a^i-{\mathit{\theta }}_b^i|)}}{\xi }_{+,mm}^{\mathit{ij}}(\theta )+\stackrel{(\mathrm{A}15.\mathrm{b})}{\overbrace{{\sum }_{\ell \ne k}{\sum }_{b\in F_{\ell }^j}{w}_b^j\mathrm{\Delta }(|{\mathit{\theta }}_a^i-{\mathit{\theta }}_b^i|)}}\sum _n{\mathcal{P}}_n^{\ell }{\xi }_{+,mn}^{\mathit{ij}}(\theta )].}\hfill \end{array}$$(A.15)