Impact of survey spatial variability on galaxy redshift distributions and the cosmological $3\times2$-point statistics for the Rubin Legacy Survey of Space and Time (LSST)
Qianjun Hang Benjamin Joachimi Eric Charles John Franklin Crenshaw Patricia Larsen Alex I. Malz Sam Schmidt Ziang Yan Tianqing Zhang the LSST Dark Energy Science Collaboration


firstpage–lastpage
Dispelling Four Challenges in Inertial Motion Tracking with One Recurrent Inertial Graph-based Estimator (RING)
    [
    September 9, 2024
===============================================================================================================




§ ABSTRACT 
We investigate the impact of spatial survey non-uniformity on the galaxy redshift distributions for forthcoming data releases of the Rubin Observatory Legacy Survey of Space and Time (LSST). Specifically, we construct a mock photometry dataset degraded by the Rubin OpSim observing conditions, and estimate photometric redshifts of the sample using a template-fitting photo-z estimator, BPZ, and a machine learning method, FlexZBoost. We select the Gold sample, defined as i<25.3 for 10 year LSST data, with an adjusted magnitude cut for each year and divide it into five tomographic redshift bins for the weak lensing lens and source samples. We quantify the change in the number of objects, mean redshift, and width of each tomographic bin as a function of the coadd i-band depth for 1-year (Y1), 3-year (Y3), and 5-year (Y5) data. In particular, Y3 and Y5 have large non-uniformity due to the rolling cadence of LSST, hence provide a worst-case scenario of the impact from non-uniformity. We find that these quantities typically increase with depth, and the variation can be 10-40% at extreme depth values. 

Based on these results and using Y3 as an example, we propagate the variable depth effect to the weak lensing 3×2pt data vector in harmonic space. We find that galaxy clustering is most susceptible to variable depth, causing significant deviations at large scales if not corrected for, due to the depth-dependent number density variations. For galaxy-shear and shear-shear power spectra, we find little impact given the expected LSST Y3 noise.





cosmology: observations – techniques: photometric – large-scale structure of Universe







§ INTRODUCTION



Observational cosmology enters the era of high-precision measurements. For example, weak gravitational lensing, which probes the small distortion of distant galaxy shapes due to the gravity of foreground large-scale structures, is particularly sensitive to the clustering parameter S_8=σ_8√(Ω_ m/0.3). Current weak lensing surveys have measured this parameter to be S_8=0.759^+0.024_-0.021 by the Kilo-Degree Survey <cit.>, S_8=0.759^+0.025_-0.023 by the Dark Energy Survey <cit.>, and S_8=0.760^+0.031_-0.034 (S_8=0.776^+0.032_-0.033) using the shear power spectra (two-point correlation function) by the Hyper Suprime-Cam <cit.>. The constraints are comparable to that measured by <cit.> from the primary cosmic microwave background (CMB), S_8=0.830±0.013, and the recent result from CMB lensing <cit.>, S_8=0.840±0.028, but are interestingly lower by 2-3σ. 
The uncertainties of these measurements are already dominated by systematic errors - without a careful treatment of various systematic effects, the cosmological results can be biased up to a few sigma <cit.>.
The forthcoming Stage IV surveys such as the Rubin Observatory Legacy Survey of Space and Time (LSST) will achieve a combined Figure of Merit ten times as much as the Stage III experiments as mentioned above <cit.>. 
While the high statistical power enables pinning down the nature of such tensions, systematic error needs to be controlled down to sub-percent level to ensure that our results are not biased.

One major systematic uncertainties come from survey non-uniformity. Galaxy samples detected at different survey depth, for example, will have different flux errors and number of faint objects near the detection limit. This could propagate down to systematic errors in redshift distribution and number density fluctuation.
The majority of the LSST footprint will follow the wide-fast-deep (WFD) observing strategy, which means that a large survey region will be covered before building up the survey depth.
At early stages of the survey, fluctuations in observing conditions, such as sky brightness, seeing, and air mass, are expected to be significant across the footprint. These can change the per-visit 5σ limiting magnitude, m_5, leading to depth non-uniformity in the early LSST data <cit.>.


The survey strategy later on could also affect uniformity. LSST will adopt a `rolling cadence', which means that during a fixed period, more frequent revisits will be assigned to a particular area of the sky, whereas the rest of the regions are deprioritized by up to 25% of the baseline observing time. The high- and low-priority regions continue to swap, such that the full footprint is covered with the same exposure time after ten years. 
This can result in different limiting magnitudes across the sky at intermediate stages of rolling.
This strategy greatly advances LSST's potential for time domain science for e.g. denser sampling in light curves. However, it also poses challenges to the analysis of large-scale structure (LSS) probes, which normally prefers a uniform coverage.


Changes in m_5 can change the detected sample of galaxies and its photometric redshifts in two ways. Firstly, a larger m_5 means that fainter, higher redshift galaxies will pass the detection limit. 
This increases the sample size, and could shift the ensemble mean redshift higher. 
These faint galaxies also contain large photometric noise, resulting in larger scatter with respect to the true redshift, hence broadening the redshift distribution. 
Secondly, at fixed magnitude, the signal-to-noise is larger given a larger m_5. This means that, contrary to the previous effect, the scatter in spec-z vs photo-z will be reduced due to the reduced noise.
These effect has been studied previously in a similar context.
The density fluctuation is quantified in <cit.> via 1+δ_ o=(1+δ_ t)(1+δ_ OS), where δ_ o is the observed density contrast, δ_ t is the true density, and δ_ OS is the fluctuation in the observing condition.
The effects on photo-z have been investigated in <cit.> in the context of LSST. They showed that the photo-z quality can change significantly with respect to different observing conditions, although they did not consider tomographic binning. 
<cit.> also quantified the effects for KiDS-1000 data, where the depth varies significantly between different pointings. They showed that by varying the r-band limiting magnitude, a significant amount of high redshift objects can be included in the sample, such that the mean number density can double between the deepest and shallowest pointings, and the average redshift for a tomographic bin can shift by as much as Δ⟨ z⟩∼ 0.2. 

Understanding these effects are important, because weak lensing is particularly sensitive to the mean redshift of the lens and source galaxies. <cit.> demonstrated that this effect is similar to a spatially varying multiplicative bias, and for cosmic shear analysis in configuration space, constraints in the Ω_ m-σ_8 plane can shift up to ∼ 1σ for a KiDS-like survey with the same area as LSST.
<cit.> also derived an analytic expression for anisotropic redshift distributions for galaxy and lensing two-point statistics in Fourier space. They showed that, assuming a spatial variation of scale ℓ_z, the effects are at percent and sub-percent level for the current and forthcoming galaxy surveys, and converge to the uniform case at ℓ≫ℓ_z.


In this paper, we investigate how survey non-uniformity can affect the redshift distribution of tomographic bins for LSST 1-year, 3-year, and 5-year observation (hereafter Y1, Y3, and Y5 respectively).
The LSST Dark Energy Science Collaboration (DESC) Science Requirements Document <cit.> states that the photometric redshifts needs to achieve a precision of ⟨Δ z ⟩ = 0.002(1+z) (0.001(1+z)) for Y1 (Y10) weak lensing analysis, and ⟨Δ z ⟩ = 0.005(1+z) (0.003(1+z)) for Y1 (Y10) large-scale structure analysis. 
Here, using these numbers as a bench mark, we quantify changes in the mean redshift (⟨ z ⟩) and width (σ_z) of tomographic bins, as depth varies.
[Notice that the DESC SRD also provides requirements on the photometric redshift scatter of the full, unbinned sample, σ_Δ z. For weak lensing, this is σ_Δ z=0.006(1+z) (0.003(1+z)) for Y1 (Y10); for large-scale structure analysis, this is σ_Δ z=0.1(1+z) (0.03(1+z)) for Y1 (Y10). Because we do not try to optimize the photometric redshift estimation in this paper, we do not compare our results with the DESC SRD σ_Δ z values.]
We use the up-to-date LSST observing strategy and the simulated 10-year observing conditions for Rubin Observatory <cit.> to quantify the survey non-uniformity, and generate a mock catalogue of true galaxy magnitude in ugrizy, redshift, and ellipticity based on the Roman-Rubin (DiffSky) simulations <cit.>. 
The degradation of photometry and photo-z estimation relies on the public software, Redshift Assessment Infrastructure Layers[<https://github.com/LSSTDESC/RAIL>] <cit.>, which will also be used in the LSST analysis pipeline.
Finally, we propagate these effects to the clustering and weak lensing two-point statistics.




This paper is organized as follows. We describe our simulation datasets in Section <ref> and introduce our methods in Section <ref>. The results are presented in Section <ref>.

We show the variation of the angular power spectra with varying depth effects in Section <ref>.
Finally, we conclude in Section <ref>.







§ SIMULATIONS



This section provides an overview of the simulations used in this work, namely, the Rubin Operation Simulator (OpSim; Section <ref>), which simulates the observing strategy and related properties for Rubin LSST, and the Roman-Rubin simulation (DiffSky; Section <ref>), which provides a truth catalogue complete up to z=3 with realistic galaxy colours.




















 §.§ Rubin Operations Simulator (OpSim)









The Operations Simulator[<https://rubin-sim.lsst.io/>] (OpSim) of the Rubin Observatory is an application that simulates the telescope movements and a complete set of observing conditions across the LSST survey footprint over the 10-year observation period, providing predictions for the LSST performance with respect to various survey strategies. 
We utilise OpSim baseline v3.3, the most recent observing strategy. This strategy involves a rolling cadence that starts after the first year of observation. In subsequent years, parts of the sky will receive more visits than others, enabling higher resolution sampling for time domain science. At the end of the fiducial survey, uniformity will be recovered at the expected 10-year LSST depth. 
The output of OpSim is evaluated by the Metrics Analysis Framework (MAF), a software tool that computes summary statistics (e.g. mean and median of a particular observing condition over a given period) and derived metrics (e.g. coadd 5σ depth) that can be used to assess the performance of the observing strategy, in terms of survey efficiency and various science drivers.

For the purpose of this study, we obtain survey condition maps in HEALPix <cit.> format using the MAF HEALPix slicer with N_ side=128 (corresponding to a pixel size of 755  arcmin^2), using the (RA, Dec) coordinates. We do not choose a higher resolution for the map because we expect that survey conditions vary smoothly on large scales, and this choice of N_ side is enough to capture the variation with the rolling pattern. 
For our purposes, we mainly consider the following quantities in each of the ugrizy filters: extinction-corrected coadd 5σ point source depth (, hereafter m_5^ ex) and the effective full width half maximum seeing (, hereafter θ_ FWHM^ eff) in unit of arcsecond. 
The m_5^ ex is different from the coadd depth, m_5, by the fact that it includes the lost of depth near galactic plane.
The effective seeing, θ_ FWHM^ eff, has a wavelength dependence, with a poorer seeing at bluer filters from Kolmogorov turbulence. 
The MAF also takes into account for increase in PSF size with airmass, X, due to seeing, i.e. θ_ FWHM^ eff∝ X^0.6. However, the MAF does not include the increase in PSF size along the zenith direction with zenith angle, due to differential chromatic refraction.
This quantity is used here to convert point-source depth to that for extended objects. 
We obtain maps of these quantities over the LSST footprint at the end of each full year of observation (e.g. Y3 for nights<1095). 
The coadded depth in each band is computed by assessing the 5σ-depth (in magnitudes) of each visit within each HEALPix pixel, then computing the `stacked' depth.
Maps of θ_ FWHM^ eff contain the median over all visits in a particular HEALPix pixel.

Throughout the paper, we will use Y1, Y3, and Y5 as examples to showcase the impact of spatial variability on photometric redshifts. 
Notice that the choice of Y3 and Y5 are a pessimistic one, because the survey strategy is close to uniformity in Y4 and Y7 where cosmological analysis are expected to be conducted. Hence, this paper provides a worst-case scenario of the severity of the impact from spatial variability.
Also, the Rubin observing strategy is still being decided, and the rolling cadence may move to different times during the survey. 
There are ongoing efforts on recommendations about the observing strategy, and hence the results shown here should be interpreted in light of this particular strategy and years chosen. 
We will focus on the wide-fast-deep (WFD) survey program footprint, and exclude areas with high galactic extinction E(B-V)>0.2 for cosmological studies. Notice that, in practice, additional sky cuts could also be applied (e.g. a depth cut that removes very shallow regions).
Specifically, we will focus on the variation with respect to i-band, the detection band of LSST. 
Figure <ref> shows the spatial variation of the extinction-corrected coadd i-band depth for OpSim baseline v3.3 in Y1, Y3, and Y5. The stripes visible across the footprint in Y3 and Y5 are the characteristics of the rolling cadence. 
The distribution of all OpSim variables are shown in Fig. <ref> for each of the six filters and for selected years of observation. One can see that the coadd depths build up in each band over the years, whereas the distributions of the median effective seeing per visit are relatively unchanged. One can also see a strong skewness in these distributions.





 §.§ Roman-Rubin simulation (DiffSky)



In order to investigate the impact of varying survey conditions on photo-z for LSST, we need a simulated truth catalogue that is complete to beyond the LSST 10-year depth and realistic in colour-redshift space. For this purpose, we use the joint Roman-Rubin simulation v1.1.3. This simulation is an extension of the effort in <cit.>, but with many improvements, including self-consistent, flexible galaxy modelling. The simulation is based on its precursor, CosmoDC2 <cit.>, a synthetic sky catalogue out to z=3 built from the `Outer Rim' N-body cosmological simulation <cit.>. The N-body simulation contains a trillion particles with a box size of (4.225  Gpc)^2.
The galaxies are simulated with Diffsky[<https://github.com/LSSTDESC/lsstdesc-diffsky>], based on two differentiable galaxy models: Diffstar <cit.> and Differentiable Stellar Population Synthesis <cit.>. 
Using Diffstar, one can build a parametric model that links galaxy star formation history with physical parameters in halo mass assembly. Then, with DSPS, one can calculate the SED and photometry of a galaxy as a function of its star formation history, metallicity, dust, and other properties. 
The advantage of this galaxy model is that the distribution in colour-redshift is smooth and more realistic compared to that in CosmoDC2. This is thanks to the separate modeling for different galaxy components, i.e. bulge, disk, and star-forming regions. The SEDs built from these different components with different stellar populations makes the colours more realistic for photo-z estimation. 

We randomly subsample the full simulated catalogue to N=10^6 objects complete to i<26.5 as our truth sample. For each object, we obtain its magnitude in the six LSST bands, true redshift, bulge size s_b, disk sizes s_d, bulge-to-total ratio f_b, and ellipticity e.
We obtain the galaxy semi-major and semi-minor axes, a,b via a=s/√(q) and b= s√(q), where s is the weighted size of the galaxy, s= s_bf_b + s_d(1-f_b), and q is the ratio between the major and minor axes, related to ellipticity via q=(1-e)/(1+e).

One caveat of the current sample is that, at z>1.5, there is an exaggerated bimodal distribution in the g-r colour and redshifts, which is not found in real galaxy data. As a result, the bluest objects in the sample are almost always found at high redshifts. This could be due to the high-redshift SPS models being less well constrained. One direct consequence of this is that, when training a machine learning algorithm to estimate the photo-z, the high-redshift performance may be too optimistic due to this colour-space clustering.








































§ METHODS



This section describes our methodology for generating a mock LSST photometry catalogue for Y1, Y3, and Y5, applying photometric redshift estimation algorithms, and defining metrics to assess the impact of variable depth. Specifically, we describe the degradation process using the LSST error model in Section <ref>, the two photo-z estimators, BPZ and FlexZBoost, in Section <ref>, the tomographic binning strategy in Section <ref>, and the relevant metrics Section <ref>.



 §.§ Degradation of the truth sample





Given a galaxy with true magnitudes m_ t={ugrizy} falling in a HEALPix pixel within the footprint, we `degrade' its magnitude with observing conditions associated with that pixel, and assign a set of `observed' magnitudes m_ o and the associated magnitude error σ_m, o, using the following procedure: (1). Apply galactic extinction. (2). Compute the point-source magnitude error for each object in each filter, using the LSST error model detailed in <cit.>. (3). Compute the correction to obtain the extended-source magnitude errors. (4). Sample from the error and add it to the true magnitudes. Steps (2) - (4) are carried out using the python package [<https://github.com/jfcrenshaw/photerr/tree/main>] <cit.>. We detail each step below.

Firstly, we apply the galactic extinction to each band with the E(B-V) dust map <cit.> via:

    m_ dust=m+[A_λ/E(B-V)]E(B-V),

where for each of the six LSST filters we adopt [A_λ/E(B-V)]={4.81,3.64,2.70,2.06,1.58,1.31}.

Then, we utilize the LSST error model <cit.>  to compute the expected magnitude error, σ_m, per band. 
The magnitude error is related to the noise-to-signal ratio, nsr, via:

    σ_m = 2.5log_10(1+ nsr).

The total nsr consists of two components:

    nsr^2 =  nsr_ sys^2+  nsr_ rand, ext^2,

where nsr_ sys is the systematic error from the instrument read-out and nsr_ rand is the random error arising from observing conditions on the sky, for extended objects. 
Notice that in the high signal-to-noise limit where nsr≪ 1, σ_m∼ nsr, and Eq. <ref> recovers the form in <cit.>. Throughout the paper, we set nsr_ sys≈σ_ sys=0.005, which corresponds to the maximum value allowed from the LSST requirement. 
For point sources, the random component of nsr is given by

    nsr_ rand, pt^2=(0.04-γ)x+γ x^2,

where γ is a parameter that depends on the system throughput. We adopt the default values from <cit.>, γ={0.038, 0.039, 0.039, 0.039, 0.039, 0.039} for ugrizy. x is a parameter that depends on the magnitudes of the object, m, and the corresponding coadd 5σ depth, m_5, in that band:

    log_10 x ≡ 0.4 (m-m_5).

For extended sources, we adopt the expression in <cit.>, where the nsr receives an additional factor related to the ratio between the angular size of the object and that of the PSF:

    nsr_ rand, ext =  nsr_ rand, pt√(A_ ap/A_ psf).

Here, 

    A_ psf=πσ_ psf^2,   σ_ psf = θ_ FWHM^ eff/2.355,


where θ_ FWHM^ eff is the effective FWHM seeing (it is linked to the seeing by θ_ FWHM^ eff = θ_ FWHM X^0.6, where X is the airmass) for a given LSST band. The AP angular size of the object is given by

    A_ ap   =π a_ apb_ ap,
    
        a_ ap    = √(σ_ psf^2+(2.5a)^2), 
    
        b_ ap    = √(σ_ psf^2+(2.5b)^2),

where a, b are the galaxy semi-major and minor axis. 

We make one modification to Eq. <ref>, where we replace the denominator by the mean PSF area, √(⟨ A_ psf⟩), averaged over pixels in the i-band quantiles which we will elaborate shortly. 
In the approximation that nsr_ rand,pt∝ x, the point-source noise is then proportional to θ_ FWHM^ eff (see Eq. <ref>), and so for the extended-source noise, θ_ FWHM^ eff cancels and Eq. <ref> effectively changes the dependence of m_5 on PSF size to that on the extended aperture size. However, in this work, we utilize the median seeing, for which the cancellation may not be exact. 
Naively taking Eq. <ref> could lead to unrealistic cases, where, at fixed depth, nsr_ rand,ext increases with a better seeing. We have tested both scenarios, i.e., using individual A_ psf or the mean ⟨ A_ psf⟩ in Eq. <ref>, and find negligible difference for our main conclusion in the i-band quantiles. However, it does make a significant difference if one were to bin the samples by quantiles of seeing, as investigated in Appendix <ref>.

To obtain the observed magnitudes m_ o, we degrade in flux space, f_ o, by adding a random noise component Δ f drawn from a normal error distribution, Δ f∼𝒩(0, nsr), to the reddened flux f_ dust of the object. Here, nsr is computed by setting m=m_ dust in Eq. <ref>. The flux and magnitude are converted back and forth via

    m_k=-2.5log_10f_k,   k={ dust, o}.

Negative fluxes are set as `non-detection' in that band. The corresponding magnitude error σ_m, o is computed using Eq. <ref> and setting m=m_ o in Eq. <ref>, such that the error de-correlates with the observed magnitude.





To focus on the trend in the depth variation in the detection band, we subdivide pixels in the survey footprint into 10 quantiles in i-band m_5^ ex, where the first quantile (qtl=0) contains the shallowest pixels, and the last quantile (qtl=9) contains the deepest.
Table <ref> shows the mean and standard deviation of each i-band depth quantile. We also show in Table <ref> the mean and standard deviation of all other survey condition maps used in the analysis in each of the i-band depth quantiles.
Within each quantile, we randomly assign each galaxy to a HEALPix pixel in that quantile, with its associated observing conditions { E(B-V), m_5^ ex, θ_ FWHM^ eff} on that pixel for each LSST band, from the OpSim MAF maps. Then, we carry out the above degradation process to our truth sample.
On average, each pixel within each quantile is assigned 121 galaxies.

Notice that there are many other parameters that could affect the photometric errors, e.g. sky background, exposure time, and atmospheric extinction. Following <cit.>, because these quantities only contribute towards m_5, we do not include them otherwise in the degradation, and assume that m_5^ ex completely captures their variation. 
Additionally, we explore the relation between m_5 and these extended quantities using OpSim in Appendix <ref>, and we explore the galaxy redshift distribution dependence with other survey properties in Appendix <ref>.

Finally, we apply an i-band magnitude cut corresponding to the LSST Gold sample selection on the degraded catalogue. For the full 10-year sample this is defined as i<25.3. For data with an observation period of N_ yr years, we adjust the gold cut to i_ lim=25.3 + 2.5log_10(√(N_ yr/10)).
Thus for Y1, Y3, and Y5, we adopt the following gold cuts respectively: i_ lim=24.0, 24.6, 24.9. 
Notice that this is slightly shallower than the definition in the DESC SRD, where the Gold cut is defined as one magnitude shallower than the median coadd m_5. This is due to the fact that OpSim baseline v3.3 has a slightly deeper i-band depth in early years compared to previous expectations. For Y1, the median i-band m_5^ ex is ∼ 25.2, giving a DESC SRD Gold cut to be 0.2 mag deeper than what we adopt here. 
Additionally, for our fiducial sample, we also apply a signal-to-noise cut in i-band: SNR=1/ nsr≥10, although we also look at the case with the full sample. This cut is motivated by the selection of the source sample, where shape measurements typically require a high SNR detection in i-band. In this work, we apply this cut to both the weak lensing and clustering samples.





 §.§ Photo-z estimators





Methods for photometric redshift estimation can be broadly divided into two main categories: template-fitting and machine learning. Template fitting methods assume a set of SED templates for various types of galaxies, and use these to fit the observed magnitudes of the targets. Machine Learning methods, on the other hand, use machine learning algorithms trained on a reference sample, to infer the unknown target redshifts. 
See <cit.> for a review and comparison of the performance of various photo-z estimators in the context of Rubin LSST. In this work, we adopt two algorithms with reasonable performance, a template fitting method, BPZ (Bayesian Photometric Redshifts), and a machine learning method, FlexZBoost. In this work, before applying these redshift estimators, all observed magnitudes are de-reddened, by applying the inverse of Eq. <ref>.



  §.§.§ BPZ (Bayesian Photometric Redshifts)


BPZ <cit.> is a template-based photometric estimation code. Given a set of input templates 𝐭, BPZ computes the joint likelihood P(z,𝐭) for each galaxy with redshift z. A prior P(z,𝐭|m) is included based on the observed magnitude of the galaxy m. For example, the prior restricts bright, elliptical galaxies to lower redshifts. 
For each galaxy, a likelihood P(z,𝐭|c, m) given the galaxy's colour c and magnitude is computed, and by marginalising over the templates, one obtains the per-object redshift probability P(z). 

We use the RAIL interface of the BPZ algorithm, with the list of spectral energy distribution (SED) templates adopted in <cit.>: the CWW+SB4 set introduced by <cit.>, the El, Sbc, Scd & Im from <cit.>, the SB2 & SB3 from <cit.>, and the 25Myr & 15Myr `SSP' from <cit.>. We set the primary observing band set to i-band, and adopt the prior from the original BPZ paper <cit.>, which was used to fit data from the Hubble Deep Field North <cit.>. Notice that these set of SEDs may be different from that in the Roman-Rubin simulation, and the prior distributions may not match exactly. 
The prior mismatch would only affect samples with low signal-to-noise and hence those posteriors are prior-dominated. For the Gold sample considered in this paper, the impact of the prior on the mean difference and scatter of the true and photometric redshifts is expected to be small, although galaxies with broad or bimodal posteriors may end up having a different point estimate (e.g. mode), hence the outlier rate could be slightly higher. We do not include extra SED templates here.

Additionally we compute the odds parameter, defined as 

    odds = ∫_z_ mode-Δ z^z_ mode+Δ zP(z) dz,

where z_ mode is the mode of P(z), and Δ z =ϵ (1+z_ mode) defines an interval around the mode to integrate P(z). The maximum value of odds is 1, which means that the probably density is entirely enclosed within the integration range around the mode, whereas a small odds means that the probability density is diffuse given the range. 
Hence, odds denotes the confidence of the BPZ redshift estimation, and the choice of ϵ essentially sets the criteria. The (1+z_ mode) factor accounts for the fact that larger redshift errors are expected at higher redshifts. We choose ϵ=0.06 as a nominal photo-z scatter, and we use odds as a BPZ `quality control', where a subsample is selected with odds≥0.9, as comparison to the baseline sample. 




  §.§.§ FlexZBoost


FlexZBoost <cit.> is a machine-learning photo-z estimator based on FlexCode <cit.>, a conditional density estimator (CDE) that estimates the conditional probability density p(y|𝐱) for the response or parameters, y, given the features 𝐱.
The algorithm uses basis expansion of univariate y to turn CDE to a series of univariate regression problems. 
Given a set of orthonormal basis functions {ϕ_i(y) }_i, the unknown probability density can be written as an expansion:

    p(y|𝐱)=∑_j β_j(𝐱)ϕ_j(y).

The coefficients β_j(𝐱) can be estimated by a training set (𝐱,y) using regression.
The advantage of FlexCode is the flexibility to apply any regression method towards the CDE.
The main hyper-parameters involved in training is the number of expansion coefficients and those associated with the regression. 
<cit.> found that FZBoost was among the the strongest performing photo-z estimators according to the established performance metrics.

In this paper, we utilise the RAIL interface of the FZBoost algorithm with its default training parameters. 
We construct the training sample by randomly drawing 10% of the degraded objects from each of the deciles, and train each year separately. 
Notice that this training sample is fully representative of the test data, which is not true in practice. Spectroscopic calibration samples typically have a magnitude distribution that is skewed towards the brighter end, and the selection in colour space can be non-trivial depending on the specific dataset used. Although there are methods to mitigate impacts from this incompleteness, such as re-weighting in redshift or colours <cit.>, and, more recently, using training data augmentation from simulations <cit.>, the photo-z performance is not comparable to having a fully representative sample, and one would expect some level of bias and increased scatter depending on the mitigation method adopted. Here, we are interested in whether our results on the non-uniformity impact changes significantly with an alternative photo-z algorithm. We thus leave the more realistic and sophisticated case with training sample imperfection to future work.



  §.§.§ Performance


For both photo-z estimators, we use the mode of the per-object redshift probability, P(z), as the point estimate, z_ phot.
Fig. <ref> shows the scatter in spec-z and photo-z for Y1, Y3, and Y5 with BPZ and FZBoost redshifts, for the shallowest (qtl=0) and the deepest (qtl=9) quantiles in the i-band m_5^ ex respectively. 
The scatter is always larger for the shallower sample in the full sample case (faint purple dots). This is expected following Eqs. <ref> and <ref>, given that the coadd depths in each band are strongly correlated. At fixed magnitude, the larger the m_5, the smaller the photometric error, hence also the smaller the scatter in photo-z. The signal-to-noise cut at SNR≥10 removes some extreme scatter as well as objects from the highest redshifts. This is more obvious for the shallowest sample compared to the deepest, due to the better signal-to-noise for the deepest sample at high redshifts.

There is a significant group of outliers in the BPZ case that are at low redshifts but are estimated to be at z>2, highlighted by the blue contours. By examining individual BPZ posteriors for this group, we find that these objects tend to have very broad or bimodal redshift distributions. This could be a result of the Oxygen line confusion between OII and OIII, and notice that the fraction of this population as well as its location can be influenced by the choice of the BPZ priors. Another possible cause is the spurious bimodal distribution in the colour-redshift space in the Roman-Rubin simulation, as mentioned in Section <ref>.
We see that after applying a strict cut with odds≥0.9, shown by the purple dashed lines enclosing 90% of the sample, the outlier populations are significantly reduced, as expected.
This cut retains 20.4% (27.7%), 25.7% (44.4%), 29.5% (44.0%) of the SNR≥10 sample in qtl=0 (9) for Y1, Y3, and Y5, respectively. We see that this cut further reduces the scatter at z_ phot∼ 1.5. FZBoost in general shows a much better performance, given that the training data is fully representative of the test data.
Table <ref> summarizes these findings for each sample via a few statistics of the distribution of the difference between photo-z and true redshifts: Δ z=(z_ phot-z_ true)/(1+z_ true). Namely, the median bias Median(Δ z), the standard deviation, the normalized Median Absolute Deviation (NMAD) σ_ NMAD=1.48  Median(|Δ z|), and the outlier fraction with outliers defined as |Δ z|>0.15. 





 §.§ Tomographic bins






In weak lensing analysis, the full galaxy catalogue is sub-divided into a `lens' sample and a `source' sample. The lens sample is often limited at lower redshifts, acting as a tracer of the foreground dark matter field which `lenses' the background galaxies. The source sample contains the background galaxies extending to much higher redshifts, whose shapes are measured precisely to construct the shear catalogue. The two samples together allow measurement of the so-called `3×2pt statistics', including galaxy clustering from the lens sample, galaxy-galaxy lensing from the lens galaxies and source shapes, and cosmic shear from the source shapes alone.
Additionally, both the lens and source samples are divided into several tomographic bins, i.e., sub-samples separated with sufficient distinction in redshifts. This further includes evolution information that improves cosmological constraints. 

We adopt the Y1 tomographic bin definitions in the DESC SRD for all of our samples. The lens sample has 5 bins equally spaced in 0.2<z<1.2, with bin width Δ z =0.2, and bin edges defined using z_ phot. 
For source samples, the DESC SRD requires five bins with equal number of galaxies. To do so, we first combine the 10 depth quantiles, and then split the sample into five z_ phot quantiles. 


Notice that in practice, tomographic binning can be determined in different ways, often with the aim of maximizing the signal-to-noise of the two-point measurements. In some cases, clustering algorithm, e.g. random forest, rather than a photo-z estimator, is used to separate samples into broad redshift bins. We refer the interested readers to <cit.> for explorations of optimal tomographic binning strategies for LSST. 
Notice also that, following the DESC SRD, we do not apply additional magnitude cuts for the lens sample. This is done, for example, for the DES Y3 MagLim lens sample, where a selection of i>17.5 and i<4z_ phot+18 is applied <cit.>. These cuts are applied to reduce faint, low-redshift galaxies in the lens sample, such that the photometric redshift calibration is more robust. 
Notice that if the lens samples are selected with a brighter cut, one would expect a different and likely reduced depth variation.
We explore this particular case in Appendix <ref>.

Fig. <ref> shows the normalized true redshift distribution, p(z), of the lens and source tomographic bins for Y3 as an example, split by the BPZ redshifts (with or without odds selection) and the FZBoost redshifts. The dashed lines show the p(z) measured from the shallowest samples, whereas the solid lines show that from the deepest samples. The BPZ case shows more extended tails in each tomographic bin compared to the FZBoost case, and for the source galaxies, a noticeable outlier population at low redshifts in the highest tomographic bin. We see that in most cases, there is a clear difference in p(z) between the shallow and the deep samples: the deep samples seem to shrink the tails, making p(z) more peaky towards the mean redshift (although this is not the case for the odds≥0.9 sample), and their p(z) seems to shift towards higher redshift at the same time. To quantify these changes, we define metrics for the impact of variable depth below.



 §.§ Metrics for impact of variable depth



The first metric is the variation in the number of objects in each sample, N_ gal, as a function of the coadd i-band depth. This is the most direct impact of varying depth: deeper depth leads to more detection of objects within the selection cut. The result is that the galaxy density contrast, δ_g(θ)=[N(θ)-N̅]/N̅, where N(θ) is the per-pixel number count at pixel θ, and N̅ is the mean count over the whole footprint, fluctuates according to the depth variation, leading to a spurious clustering signal in the two-point statistics.
To quantify the relative changes, we measure the average number of objects per tomographic bin across all 10 depth quantiles, N̅_ gal=∑_i N_ gal,iw_i, where i=1,..,10 denotes the depth bin, and w_i∼0.1 is the weight proportional to the number of pixels in that quantile. We quote the change of object number in terms of N_ gal/N̅_ gal.

The second metric quantifies the mean redshift of the tomographic bin as a function of depth:

    ⟨ z ⟩  = ∫ z p(z) dz ,

where p(z) is the true redshift distribution of the galaxy sample in the tomographic bin with normalization ∫ p(z) dz=1. Weak lensing is particularly sensitive to the mean distance to the source sample: the lensing kernel thus differs on patches with different depth. 
Here, we look at the difference between the mean redshift ⟨ z ⟩_i of depth quantile i and that of the full sample,  ⟨ z ⟩_ tot, i.e., Δ⟨ z ⟩≡⟨ z ⟩_i - ⟨ z ⟩_ tot. More specifically, we look at the quantity Δ⟨ z ⟩ /(1+⟨ z ⟩_ tot), where the weighting accounts for the increase in photo-z error towards higher redshifts. This format also allows us to compare with the DESC SRD requirements. 

The third metric quantifies the width of the tomographic bin. This is not a well-defined quantity because the p(z) in many cases deviate strongly from a Gaussian distribution. One could use the variance, or the second moment of the redshift distribution:

    σ_z^2  = ∫ (z-⟨ z ⟩)^2 p(z) dz.

However, this quantity is very sensitive to the tails of the distribution: larger tails of p(z) increases σ_z, even if the bulk of the distribution does not change much.
In our case, the width of the tomographic bin is most relevant for galaxy clustering measurements: the smaller the bin width, the larger the clustering signal. 
Specifically, in the Limber approximation, the galaxy auto-correlation angular power spectrum is given by

    C_ℓ^ gg=∫dχ/χ^2(z)[ H(z)/c p(z)]^2 P_gg(k=ℓ+1/2/χ, z),

where ℓ is the degree of the spherical harmonics, χ is the comoving distance, H(z) is the expansion rate at redshift z, c is the speed of light, k is the 3-dimensional wave vector, and P_gg is the 3-dimensional galaxy power spectrum. 
Assuming that within the tomographic bin, the redshift evolution of galaxy bias is small, and all other functions can be approximated at the mean value at the centre of the bin, the clustering signal is proportional to the integral of the square of the galaxy redshift distribution, p(z). This assumption breaks down if the tomographic bin width is broad, for instances, the combination of all five lens bins.
Hence, we define the following quantity:

    W_z ∫ p^2(z) dz

as the LSS diagnostic metric, which corresponds to changes of the two-point angular power spectrum kernel with respect to changes in p(z).
 This is a useful complement to the second moment, σ_z, because σ_z can be sensitive to the tails of the p(z) distribution caused by a small population of outliers in photo-z; however, the impact of this population could be small for galaxy clustering, which is characterised by W_z. For both of these quantities, we look at the ratio with the overall sample combining all depth quantiles. 

We show all the mean metric quantities in each tomographic bin and each quantile for Y1, Y3, and Y5 in Table <ref> for BPZ and Table <ref> for FZBoost.

Notice that for the p(z)-related quantities, we have used the true redshifts, but in practice, these are not accessible. Rather, unless one uses a Bayesian hierarchical model such as CHIPPR <cit.>, one only has access to the calibrated redshift distribution p_c(z) against some calibration samples via, e.g. a Self-Organizing Map (SOM), which is itself associated with bias and uncertainties that can be impacted by varying depth.
The case we present here thus is idealized, where the calibration produces the perfect true p(z). This allows us to propagate the actual impact of varying depth on p(z) to the 3×2pt data vector, but does not allow us to assess the bias at the level of modelling due to using an `incorrect' p_c(z) that is affected also by the varying depth. 
We leave this more sophisticated case to future work.





§ RESULTS



This section presents our results on the impact of variable depth via three metrics: the number of objects (Section <ref>), mean redshift of the tomographic bin (Section <ref>), and the width of the tomographic bin (Section <ref>).









 §.§ Number of objects



Figure <ref> shows the change in the number of objects, N_ gal, as a function of the i-band extinction-corrected coadd depth, m_5^ ex, compared to the overall mean, for lens and source tomographic bins in Y1, Y3, and Y5.
In general, we find an approximately linear increase of number of objects as the i-band depth increases, with the higher two redshift bins showing the most extreme variation. For the lower redshift bins, the variation can be ∼ 10% compared to the mean value, whereas for bin 5, the variation can be as large as ∼ 40%. The trend does not seem to change much at different observing years. 
This is the result of the i-band gold cut and the high SNR selection. The scatter in magnitudes is larger for the shallower sample, hence given a magnitude cut, the shallower sample will have fewer objects. 
At fixed magnitude, the deeper objects have larger SNR, resulting in more faint galaxies surviving the SNR cut.
Given that the gold cut and SNR at given magnitude evolve with depth in the observation year, we expect the trend to be similar across Y1 to Y5. 
It is interesting to see also that per tomographic bin, the trends for baseline BPZ and FZBoost are similar, despite having quite different features in the photo-z vs spec-z plane. The variation between bins 1 - 4 is slightly larger in the BPZ case.
For the BPZ redshifts, the inclusion of the odds selection increases the variation in object number, especially in the highest redshift bin. The steeper slope might be due to the fact that, objects with larger photometric error from the shallower regions are likely to result in a poorer fit, leading to a smaller odds value. Hence, the odds≥0.9 selection removes more objects from the shallower compared to the baseline case.











 §.§ Mean redshift



Figure <ref> shows the variation in the mean redshift of the tomographic bin, ⟨ z ⟩, as a function of the i-band extinction-corrected coadd depth, m_5^ ex, for lens and source samples in Y1, Y3, and Y5. 
In general, ⟨ z ⟩ increases with the i-band coadd depth. This is expected as more faint, high redshift galaxies that are scattered within the magnitude cut are included in the deeper sample, resulting in an increased high redshift population.
In general, the slope of this relation is similar across tomographic bins for both lens and source samples, with a variation of |Δ z /(1+⟨ z ⟩)|∼0.005 - 0.01. 
This is not true for bin 5 in the source sample, where the variation with depth is noticeably larger. 
This could be explained by this bin containing objects with the highest z_ phot, which are also most susceptible to scatter in the faint end and outliers in the photo-z estimators. This trend becomes more extreme from Y1 to Y5. By reducing outliers with the BPZ odds cut, the variation in source bin 5 is slightly reduced, although still higher than the nominal level.

There are some difference between the BPZ and FZBoost cases: the slope slightly grows from Y1 to Y5 in the BPZ case, whereas it stays consistent in the FZBoost case, but the two cases converge in Y5.
On the same figure, we mark the DESC SRD requirements for photo-z as a dark grey band at Δ z /(1+⟨ z ⟩)=±0.002 and a light grey band at Δ z /(1+⟨ z ⟩)=±0.005. 
The shifts in mean redshift reach the limit of the requirements for Y1, and exceeds the requirement for Y10.












 §.§ Width of the tomographic bin



Figure <ref> shows the change in the tomographic bin width parameters, σ_z and W_z, as defined in Section <ref> for the lens galaxies as a function of the i-band extinction-corrected coadd depth, m_5^ ex, in Y1, Y3, and Y5. 
The width of the tomographic bin can change with depth due to the scatter in the photo-z vs spec-z plane. For example, a deeper sample may have a smaller scatter for the bulk of the sample, but include fainter objects that could result as outliers, resulting a more peaked distribution at the centre with pronounced long tails.

The left two columns of Fig. <ref> show the changes in the second moment, σ_z, for both the BPZ (first column) and FZBoost case (second column).
For BPZ, there is little change in this parameter for Y1 at different depth, but for Y3 and Y5, σ_z increases with depth. Including odds selection reduces the trend, and in some cases reverses it. For FZBoost, the trend is similar to BPZ, but bin 1 shows a particularly large variation by as much as ∼30%. 
This is because σ_z is sensitive to the entire distribution, not just the peak, and outliers at high redshift can significantly impact this parameter. Fig. <ref> shows same p(z) distributions for Y3 in logarithmic scale, where the high redshift outliers are visible. Indeed, one can see an enhanced high-redshift population for bin 1 in the FZBoost case. The odds cut removes most of the outliers, so that σ_z is reflecting the change of the peak width with depth, hence giving the reversed trend. 

The right two columns of Fig. <ref> show the changes in W_z. Given a tomographic bin, a larger W_z means a more peaked redshift distribution, hence a larger clustering signal. 
One can see that W_z is more sensitive to the bulk of the p(z) distribution, as it increases with depth in most bins. We see that the variation in W_z is within 10% from the mean, with the largest variation coming from bins 2, 3, and 4. The highest and lowest tomographic bins, on the other hand, does not change much, despite their σ_z varying significantly with depth. For the BPZ case, adding the additional cut in the odds parameter reduces such trends in general, and the trend in the highest tomographic bin is reversed. 















§ IMPACT ON THE WEAK LENSING 3X2PT MEASUREMENTS




We use the Y3 FZBoost photo-z as an example to showcase the varying depth effects, by propagating the number density and p(z) variation from the previous section into the weak lensing 3×2pt data vector. In Section <ref>, we describe how the mock large-scale structure and weak lensing shear maps are constructed with the inclusion of non-uniformity. In Section <ref>, we show case the measured 3×2pt data vector in both uniform and variable depth case.











 §.§ Mock maps with varying depth






To construct the mock LSST catalogue, we use one of the publicly available Gower Street Simulations <cit.>. This is a suite of 800 N-body cosmological simulations created using PKDGRAV3 <cit.> with various wCDM cosmological parameters. The simulation outputs are saved as 101 lightcones in HEALPix format with N_ side=2048 between 0<z<49. To fill the full sky, the boxes are repeated 8000 times in a 20×20×20 array. For shells z<1.5, though, only three replications are required.
We use the particular simulation with ΛCDM cosmology: w=-1, h=0.70, Ω_m=0.279, Ω_b=0.046, σ_8=0.82, and n_s=0.97. The dark matter density contrast map, δ_m, is computed using particle counts at N_ side=512 (corresponding to a pixel size of 47.2   arcmin^2), and the corresponding lensing convergence map, κ, is produced with Born approximation using [<https://github.com/NiallJeffrey/BornRaytrace>] <cit.>. Finally, the shear map, (γ_1,γ_2) in spherical harmonic space, is produced via

    γ_E,ℓ m = κ_E, ℓ m/ℓ(ℓ + 1)√((ℓ + 2)(ℓ - 1)),

and we transform γ_E,ℓ m as a spin-2 field, γ_ℓ m=γ_E,ℓ m + i γ_B,ℓ m, assuming zero B-mode. For more details see <cit.>.


We construct the lens and source shear maps as follows. In the noiseless case, given a lens (source) redshift distribution, p_i(z), for a tomographic bin i, we construct the lens density (source shear) map by M_i= ∑_j M_j  p_i(z_j)Δ z_j, where j denotes the lightcone shells in the Gower street simulation, M_j denotes the map in this particular shell, and Δ z_j denotes the shell width. 
The noisy maps are generated in the following way.
Lens galaxy counts in tomographic bin i on each pixel θ are drawn from a Poisson distribution. For a shell j, the Poisson mean is μ_j(θ)=n_ gal,j[1+bδ_m,j(θ)], where b is the linear galaxy bias and n_ gal,j=n_ galp_i(z_j)Δ z_j, with n_ gal being the average count per pixel in this tomographic bin. Here, we set b=1 to avoid negative counts in extremely underdens pixels. 
However, notice that in a magnitude-limited survey, the galaxy bias is typically b>1 and evolves with redshift, not to mention the scale-dependence of bias on non-linear scales. 
One approach to sample b>1 is to simply set negative counts to zero. However, this may introduce spurious behaviour in the two-point function of the field.
Given the main purpose here is to propagate the systematic effects due to depth only, we justify our choice by prioritizing the precision of the measured two-point statistics compared to theory inputs.
We assume the ensemble-averaged per-component shape dispersion to be σ_e=⟨√((e_1^2+e_2^2)/2)⟩=0.35, chosen to roughly match that measured in the Stage III lensing surveys <cit.>.
For a tomographic bin i, we first assign source counts in the same way as above, resulting in n̂_ source (θ) galaxies in pixel θ. We then randomly assign shapes drawn from a Gaussian distribution, 𝒩∼(0,σ_e), for each component n̂_ source (θ) times, and we compute the mean shape noise in each pixel. We end up with a shape noise map, which we then add to the true shear map for each tomographic bin.



To imprint the varying depth effects, we divide the footprint into 10 sub-regions containing the pixels in each of the i-band m_5^ ex deciles, and repeat the above procedure with distinct number density and p(z) for both the lens and source galaxies, according to the findings in previous sections. We do not assign depth-varying shape noise, following the finding in <cit.> that the shape noise is only a weak function of depth. 
We also produce the noiseless cases for varying depth. For density contrast, we produce two versions: one with varying p(z) only, and one with additional amplitude modulation δ_m + Δδ, where Δδ+1=N_ gal/N̅_ gal, as shown in Fig. <ref>. The former is to used isolate the effect of varying p(z) only. 


We adopt the cumulative number density of the photometric sample as a function of the i-band limiting magnitude given by the DESC SRD:

    N(<i_ lim)=42.9(1-f_ mask)10^0.359(i_ lim-25)  arcmin^-2,

where f_ mask accounts for the reduction factor for masks due to image defects and bright stars, and f_ mask=0.12 corresponds to a similar level of reduction in HSC Y1 <cit.>. 
Hence, substituting i_ lim=24.6 for LSST Y3, the expected total number density is N(<24.6)=27.1  arcmin^-2.
This is slightly larger but comparable to the HSC Y3 raw number density of N=22.9   arcmin^-2 <cit.> at a similar magnitude cut of i_ lim<24.5 in the cModel magnitude. 
We estimate the total lens galaxy number density for our sample by N_ lens=N(<24.5)f_ LS, where f_ LS=0.90 is the ratio between the total number of lens and source samples (averaged over depth bins) from our degraded Roman-Rubin simulation catalogue, hence N_ lens=24.4   arcmin^-2. 

For each lens tomographic bin, we obtain the following mean number density: 3.93, 6.08, 5.66, 5.71, 3.03   arcmin^-2.
We also explore the case using a MagLim-like lens sample with a much sparser density in Appendix <ref>.
For source sample, it is the effective number density 
n_ eff, rather than the raw number density, that determines the shear signal-to-noise. n_ eff accounts for the down-weighting of low signal-to-noise shape measurements, as defined in e.g. <cit.>. For LSST, n_ eff is estimated for Y1 and Y10 with different scenarios in Table F1 in the DESC SRD. 
In the case adopted for forecasting, where the shapes are measured in i+r and accounting for blending effect, n_ eff is ∼60% of the raw number density for both Y1 and Y10. 
We follow this estimation for Y3, hence adopting n_ eff=16.3   arcmin^-2 for the full source sample, and 3.26   arcmin^-2 for each tomographic bin. This is comparable, but slightly more sparse compared to HSC Y3, where n_ eff=19.9   arcmin^-2 <cit.>. 





Meanwhile, we also generate a uniform sample for comparison, in which the number density and p(z) are given by the mean of the depth quantiles.
We assign uniform weights to lens and source galaxies.





























 §.§ Weak lensing 3x2pt data vector















We use NaMaster <cit.> to measure the 3×2pt data vector in Fourier space: C_ℓ^ gg, C_ℓ^ gγ, and C_ℓ^γγ for the lens and source tomographic bins. NaMaster computes the mixing matrix to account for the masking effects, and produces decoupled band powers. The HEALPix pixel window function correction is also applied when comparing the data with input theory.
We adopt 14 ℓ-bins in range [20,1000] with log spacing. Notice that the maximum ℓ is a conservative choice for C_ℓ^γγ, but is sufficient for our purpose to demonstrate the impact of variable depth on relatively large scales. For galaxy clustering and galaxy-galaxy lensing, we apply an additional scale cut at ℓ_ max=k_ maxχ(⟨ z ⟩)-0.5 following the DESC SRD, where k_ max=0.3h Mpc^-1, and χ(⟨ z ⟩) is the comoving distance at the mean redshift ⟨ z ⟩ of the lens tomographic bin.

We generate theory angular power spectra assuming spatial uniformity with the Core Cosmology Library[<https://github.com/LSSTDESC/CCL>] <cit.>. CCL uses HALOFIT <cit.> non-linear power spectrum and Limber approximation when computing the angular power spectra. 
We compute the Gaussian covariance matrix using NaMaster with theoretical data vectors. The covariance includes mask effects, shot-noise, and shape noise power spectra. It should be noted that this is done assuming uniformity. 

In the varying depth case, the true covariance contains extra variance, due to spatial correlation in the noise with the number count. 


Also, the assumption of a purely Gaussian covariance is not completely true. On very large scales, non-Gaussian mode coupling at scales larger than the survey footprint results in a term called super-sample covariance <cit.>. Here we expect it to be relatively small because of the large sky coverage of LSST. 
On small scales, non-linear structure formation also introduces non-Gaussian terms <cit.>. With the scale cuts adopted in C_ℓ^gg and C_ℓ^gγ we expect that such non-Gaussian contribution to be small.

The galaxy clustering angular power spectra measurements, C_ℓ^ gg, are shown in Fig. <ref>. The tomographic bin number is indicated in the upper right corner as (i,i) for bin i. The measurements for the uniform case are shown as red dots, and that for the varying depth case are shown in purple. The data points are shot-noise-subtracted. 
We see a clear difference between the uniform and the varying depth cases at ℓ<100, and it becomes more significant at higher redshifts. 
The impact at large scales is expected, as the i-band coadd depth varies relatively smoothly and the rolling pattern is imposed at relatively large scales. 

The trend with redshifts is also expected, due to two main reasons. Firstly, the slope d(N_ gal/N̅_ gal)/dm_5 increases slightly with redshift, and is significantly larger for bin 5, as shown in the right middle panel of Fig. <ref>. This means that non-uniformity is most severe in these bins. Secondly, the clustering amplitude increases towards lower redshifts due to structure growth, hence the non-uniformity imprinted in δ_g is less obvious in lower redshift bins. 
In practice, the number density fluctuations are mitigated via the inclusion of the selection weights, w(θ), such that the corrected density field is defined as δ_g(θ) = N(θ) / w(θ) N̅_w, where N̅_w=∑ N(θ) / ∑ w(θ) (see e.g. <cit.>). In addition, these weights will be used to compute the mode coupling matrix and shot noise, such that the varying number density is taken into account in the likelihood analysis.
A more subtle effect is the difference in redshift distribution at different depth. To isolate its impact, we compare the clustering power spectra from the noiseless sample varying p(z) only with that from the noiseless uniform case. The ratio of the measurements are shown as dashed lines in the first panel of Fig. <ref>. We find that once the non-uniformity in number density is removed, the variation in p(z) does not significantly bias the power spectra, and we recover the uniform case at better than 0.5%.

The galaxy-shear and shear-shear power spectra, C_ℓ^ gγ and C_ℓ^γγ, are shown in Figs. <ref> and <ref>, respectively. The source - lens and source - source combinations are indicated on the upper right as (i, j). In both cases, we only show the non-zero E-modes, and we check that the B-modes are consistent with zero. 
For the galaxy-shear case, measurements from combinations i<j are not shown, because we do not include effects such as magnification or intrinsic alignment, hence these measurements are low signal-to-noise or consistent with zero. 
We see that, overall, the impact of variable depth is much smaller compared to galaxy clustering. 
In the galaxy-galaxy shear measurements, only combination (5,5) shows a significant χ^2 in the variable depth case, and the main deviations is at ℓ<100. This could be a joint effect where non-uniformity is largest in the highest redshift bin for both lens and source.



There is negligible difference in the shear-shear measurements for all other combinations given the measurement error.
To look at this further, we take the noiseless case and compute the ratio between measurements from the varying depth sample and the uniform sample. We show some examples along the diagonal, i.e., the (i,i) combinations, in the middle and right panels of Fig. <ref>. The off-diagonal measurements lie mostly within the variation range of the ones shown here. In case of C_ℓ^ gγ, we see that deviations are large at low ℓ when both density and p(z) is non-uniform (shown as solid line); when the density non-uniformity is removed (shown in dashed line), the results are more consistent within 5%. For C_ℓ^γγ, we see that the largest impact is from the highest tomographic bin reaching up to 0.5%. 

These results are consistent with the analytical approach in <cit.>, where, in general, the varying depth effect in the redshift distributions is sub-percent and the weak lensing probes are less susceptible to these variations.
Our results are quite different from <cit.> (hereafter H20) for KiDS cosmic shear analysis in several aspects. H20 found that the largest impact comes from the sub-pointing, small scales, and for a KiDS-like set-up, the difference between the uniform and variable depth cases is 3% - 5% at an angular scale of θ=10  arcmin. Furthermore, the variable depth effect is stronger in lower redshift bins than higher redshift bins. Several differences in the analysis may contribute to these different results. Firstly, the non-uniformity in KiDS is rather different from that considered here: the KiDS footprint consists of many 1  deg^2 pointings, each having distinctive observing conditions due to that each field only received a single visit. This means that survey properties such as depth are weakly correlated at different pointings. One can write down a scale-dependent function, E(θ), to specify the probability of a pair of galaxies falling in the same pointing at each θ, and this essentially gives rise to the scale dependence of the variable depth effect in H20. For LSST, the above assumptions are not true, and E(θ) (if one can write it down) would take a very different form compared with that in KiDS. 
Secondly, due to the single visit, there is a much larger variation in depth, number density, and Δ z in KiDS compared to this work (tomographic bin centre can shift up to Δ z∼ 0.2 in redshift, as shown in Fig.2 of H20). This means that the variable depth effects in KiDS as explored by H20 is significantly larger compared to this work. 
This also explains their redshift dependence, because for KiDS, the average redshift between pointings varies the most in the lowest redshift bins. Lastly, although our ℓ_ max here corresponds to θ∼ 10   arcmin. the results are not directly comparable, as H20 conducted the analysis in real space, i.e. ξ_±(θ). 


To sum up, the largest impact of varying depth comes from galaxy clustering, whereas the impact on weak lensing probes is much smaller. Higher redshift bins are more susceptible due to a higher sensitivity in number density and redshifts with depth. Given the mock LSST Y3 uncertainty, one can clearly detect bias in the power spectrum in galaxy clustering and the galaxy-galaxy shear bin (4,4), while all other combinations do not seem to have detectable impacts. 
Furthermore, once the density non-uniformity is removed, the impact of varying depth is further reduced. There are several ways to mitigate number density variation, such as mode projection <cit.>, template subtraction <cit.>, iterative regression <cit.>, and machine learning methods using neural networks <cit.> and a Self-Organizing Map (SOM) <cit.>. 
See <cit.> for a thorough review. 
Notice that, despite these methods, it is difficult to guarantee a complete removal non-uniformity, and in some cases, clustering signal can also be reduced as a result.
Additional sky cuts to exclude problematic regions can also effectively reduce density variation, at the cost of losing sky coverage.
Finally, for the lens sample, a brighter magnitude cuts can also greatly reduce the variable depth effect (see Appendix <ref> for a MagLim-like lens selection), at the cost of sample sparsity.
Nevertheless, non-uniformity in p(z) only seems to be safely averaged out in the 2-point statics measurements.



  §.§.§ Impact on spectroscopic calibration





Here, we consider another potential source of systematics arising from small spectroscopic calibration fields. Redshift calibration for photometric surveys such as LSST are usually done using small but deep spectroscopic surveys, e.g. C3R2 survey <cit.>. Each field in these surveys has a coverage of a few deg^2. 
Suppose that a calibration field overlaps with a particularly shallow or deep region, the calibration (e.g. a trained SOM) could cause bias to the overall redshift distribution when it is generalised to the whole field. For example, a SOM trained in a shallow region will contain larger noise, which may increase the scatter for the overall sample. The lack of high redshift, fainter objects in the shallow region could also cause bias when the SOM is applied to objects in deeper regions. 

The specific impact will depend on the calibration method and details of the calibration, which is beyond the scope of this paper. 
Here, we qualitatively assess the impact via the difference in the 3×2pt theory vectors computed using the p(z) from a particular quantile and those computed using the mean p(z), as shown in Fig. <ref>. The solid lines show cases from the shallowest quantile, qtl=0, and the dashed lines show cases from the deepest quantile, where qtl=9, highlighting the worst case scenarios. For C_ℓ^gγ and C_ℓ^γγ, only cases where the tracers are in the same bin are shown, but the other lens - source combinations have a comparable variation.
We see that naively taking the p(z) from a quantile and assume it as the p(z) for the full sample can give rise to as much as 10% bias compared to the uniform case.

This effect is reduced by having multiple calibration fields across the LSST footprint. Currently, many of the calibration fields overlaps with the LSST Deept Drilling Field (DDF), which will be much deeper compared to the WFD. Impact of variable depth can then be mitigated via a two-tiered SOM calibration, mapping from the deep to the wide field <cit.>, and synthetic source injection <cit.>, mimicking the degradation of the deep field objects across the LSST footprint, as done in the DES Y3 analysis. 



§ CONCLUSIONS




In this paper, we investigated and quantified the impact of spatial non-uniformity due to survey conditions on redshift distributions in the context of early LSST data. We used the Roman-Rubin simulation as the truth catalogue, and degraded the photometry using the LSST error model implemented in the RAIL package. The degradation utilizes the survey condition maps from the OpSim baseline v3.3 for the 1-year, 3-year, and 5-year LSST data. We run BPZ and FZBoost photometric redshift estimators on the degraded sample and use the photo-z mode to separate the samples into five lens and five source tomographic bins. Finally, we apply the LSST Gold selection and a signal-to-noise cut. Taking the extinction-corrected 5σ coadd depth of the detection band, i-band, as the primary source of non-uniformity, we quantify the impact in terms of three measures: the number of objects, the mean redshift of the tomographic bin, and the tomographic bin width. 
We find that:

    
  * The number of objects increases with the i-band depth in general, and at extreme depth values, the number of objects can vary by a factor of two. The trend is relatively consistent between cases using BPZ and FZBoost, although selecting odds≥ 0.9 for BPZ amplifies the trend. The largest correlation comes from the highest tomographic bin.
    
  * The mean redshift in each bin increases with the i-band depth, with a variation of |Δ z /(1+⟨ z ⟩)|∼0.005 - 0.01. The lens samples show a relatively consistent trend across different tomographic bins, whereas for the source sample, the highest tomographic bin shows the largest variation. This reaches the limit of the requirements of 0.005 for Y1 as listed in the DESC SRD, and exceeds the requirement of 0.003 for Y10. At extreme depth variations, however, deviation in ⟨ z ⟩ could exceed Y1 requirements.
    
  * The width of the lens tomographic bin is measured in terms of σ_z, which is sensitive to the entire redshift distribution, p(z), and W_z, which is sensitive to the peak of p(z), both varying at the level of 10% and slightly increases with year. We find that in general, σ_z increases with the i-band depth due to fainter objects included in the deeper sample. W_z also increases with the i-band depth, due to a more peaked bulk p(z) as a result of higher SNR in deeper samples, although the trend can be reversed in some cases.


As emphasized before, results derived for Y3 and Y5 are with particularly large rolling non-uniformity. Hence, the variations shown should be interpreted as an upper-limit for the early Rubin LSST static science. As shown in Appendix <ref>, if the final LSST lens selection is similar to the DES Y3 MagLim sample with a bright magnitude cut, then the expected variable depth impact will be milder than shown in our baseline cases.

We took the Y3 FZBoost photo-z as an example to propagate the impact of varying depth to the weak lensing 3×2pt measurements. To do this, we used one realization of the Gower Street N-body simulation, and generated lens galaxy maps and source shear maps with spatially varying number density and p(z). We measure the data vector in harmonic space using NaMaster, and also compare them with the theory expectation generated from the Core Cosmology Library. We find that the largest impact is on C_ℓ^ gg with the higher redshift bin measurements significantly biased. C_ℓ^ gγ is less sensitive to varying depth effects, although in the source - lens combination (4,4), there is a visible difference at low ℓ. C_ℓ^γγ shows no significant impact in all source - source combinations from varying depth, given the uncertainties in LSST Y3.
Finally, we also investigate cases where we do not include noise in the lens and source maps. 
The difference between uniform and varying depth cases can be up to a few percent for C_ℓ^ gγ, and less than 0.5% for C_ℓ^γγ. Furthermore, by removing the density non-uniformity, and varying p(z) only with depth, one can reduce the bias in C_ℓ^ gg and C_ℓ^ gγ to sub-percent level. 
Therefore, for early LSST analysis, it is crucial to account for the galaxy density variation, but the impact of varying p(z) seems to be negligible. We leave the investigation of an accurate mitigation strategy of the number density variation to future work.




Our current approach has some caveats. 
Firstly, the fidelity of the colour-redshift relation in the Roman-Rubin simulation at z>1.5 is questionable. As already mentioned, the strong bifurcation of the blue objects at this high redshift may lead to worse (in the case of BPZ) or overly optimisic (in the case of FZBoost) performance when estimating the photo-z.
Secondly, we have adopted an analytic model to obtain the observed magnitudes in each band based on survey conditions. However, in reality, the observed magnitudes and colours also depend on the way they are measured. For example, for extended objects, cModel <cit.> and GAaP <cit.> methods are often applied. Although the photometry will be calibrated, the magnitude error may not be the same for different methods. This could introduce extra scatter in photo-z.
Thirdly, we have only tested on two major photo-z estimators, observing some level of differences in the results. For example, compared to BPZ, the FZBoost samples show more consistency between different tomographic bins regarding to the trend with i-band depth. Therefore, one should take the result as an order of magnitude estimate of the impact, but the specific trends are likely to differ for different photo-z methods.
Moreover, when propagating the effects to the data vector, we have made some simplifications. We considered a galaxy bias of b=1, and did not include systematics such as magnification bias or intrinsic alignments. This choice is to isolate the effect of varying depth on the pure lensing and clustering contribution, but it would be more realistic to include these effects. 
Finally, we have not folded in the effects of blending, i.e. spatially nearby galaxies are detected as one object. This occurs when the surface density is high and the image is crowded, and could be significant for deep photometric surveys such as LSST. 
The level of blending depends on both seeing and depth of the survey, hence, it could correlate with the variable depth effects discussed here.
The impact of blending on photo-z is the inclusion of a small fraction of ill-defined redshifts in the sample, increasing the photo-z scatter. Clustering redshift calibration, which measures galaxy clustering on small scales, can also be affected as these scales are most susceptible to blending.
Moreover, blending can affect shear measurements via e.g. lensing weights, hence introduce impact on galaxy-galaxy lensing and cosmic shear.
As such, <cit.> showed that approximately 12% of the galaxy sample in LSST is unrecognized blends, and can bias S_8 measurement from cosmic shear by 2σ. 

Furthermore, so far our results are based on the p(z) of the true redshifts of the sample. In reality, we do not have access to this, and our theory curve will be based on the calibrated redshift distribution p_c(z), which itself can be impacted by non-uniformity based the calibration method. 
For example, in many weak lensing surveys, a Self-Organizing Map (SOM) is used to calibrate redshifts by training on a photometric sub-sample with spectroscopic counterparts <cit.>. By taking sub-samples from a small calibration field (typically of a few square degrees) located in a particularly shallow region could result in a trained SOM that captures different magnitudes, redshifts, and SNR than that from a deep region, as quanlitatively shown in Section <ref>.
One remedy may come from calibration using clustering redshifts, which takes advantage of galaxy clustering of the target sample with a spectroscopic sample, spliced in thin redshift bins <cit.>. The nonphysical variation with depth will drop out in this method, giving unbiased estimate of p(z).


We have only explored the impact of variable depth on two-point statistics here, but there could be potential impact on statistics beyond two-point. For example, for weak lensing shear, a similar effect in manifestation is source clustering, where the number density of source galaxies n(θ̂, z) is correlated with the measured shear γ(θ̂) for a given direction θ̂ on the sky, because source galaxies are themselves clustered. Impact of source clustering is negligible in two-point statistics for Stage III surveys, but is detected significantly in several higher order statistics in the DES Y3 data <cit.>. Given that the variable depth effect also modulates n(θ̂, z) (hence imprinting a fake `source clustering'), there may be non-negligible impact on higher order statistics with LSST.

We leave these explorations to future work.




§ AUTHOR CONTRIBUTION STATEMENTS















QH: Contributed to the conceptualization, data curation, formal analysis, and writing of the draft. Contributed to the development of RAIL.

BJ: Contributed to the conceptualization, funding acquisition, project administration, and revisions of the text.

EC: Contributed to RAIL software core functionalities used in the analysis. Minor contributions to revisions of text.

JFC: Contributed to the development of RAIL and the  model.

PL: Contributed to the development of the Roman-Rubin Diffsky simulation.

AIM: Contributed to RAIL software core functionalities used in the analysis. Minor contributions to revisions of text.

SS: Contributed to RAIL software used in the analysis, namely the BPZ and FlexZBoost algorithms used in estimation, along with general software development. Minor contributions to revisions of text.

ZY: Contributed to the development of RAIL and the  model; provided reviewing comments for the manuscript.

TZ: Contributed to the development of RAIL, including the LSST error model, and RAIL's core API; provided reviewing comments for the manuscript.




§ ACKNOWLEDGEMENTS


QH and BJ are supported by STFC grant ST/W001721/1 and the UCL Cosmoparticle Initiative.
This paper has undergone internal review in the LSST Dark Energy Science Collaboration. 
The authors thank the internal reviewers, Boris Leistedt and Markus Rau, for their thorough and insightful comments. 
This work also benefited from helpful comments by Pat Burchat, Ofer Lahav, Rachel Mandelbaum, and Andrina Nicola.



The DESC acknowledges ongoing support from the Institut National de 
Physique Nucléaire et de Physique des Particules in France; the 
Science & Technology Facilities Council in the United Kingdom; and the
Department of Energy, the National Science Foundation, and the LSST 
Corporation in the United States.  DESC uses resources of the IN2P3 
Computing Center (CC-IN2P3–Lyon/Villeurbanne - France) funded by the 
Centre National de la Recherche Scientifique; the National Energy 
Research Scientific Computing Center, a DOE Office of Science User 
Facility supported by the Office of Science of the U.S. Department of
Energy under Contract No. DE-AC02-05CH11231; STFC DiRAC HPC Facilities, 
funded by UK BEIS National E-infrastructure capital grants; and the UK 
particle physics grid, supported by the GridPP Collaboration.  This 
work was performed in part under DOE Contract DE-AC02-76SF00515.
JFC acknowledges support from the U.S. Department of Energy, Office of Science, Office of High Energy Physics Cosmic Frontier Research program under Award Number DE-SC0011665.
AIM acknowledges the support of Schmidt Sciences.

We acknowledge the use of arXiv and ADS for referneces, the use of Python libraries and software mentioned in the main text for data analysis and plotting, and Overleaf for the writing of this paper.




§ DATA AVAILABILITY


 The methodology of generating mock LSST photometry with observing conditions is included in the RAIL pipeline[<https://github.com/LSSTDESC/rail_pipelines/tree/main/src/rail/pipelines/examples/survey_nonuniformity>].
 The mock galaxy catalogues with varying depth effect are available upon reasonable request.






mnras
 



















§ COMPARISON OF LSST ERROR MODEL ON DC2





The 5σ depth per visit, m_5, depends on a set of observing conditions in the following way <cit.>: 

    m_5  =     C_m + 0.50 (m_ sky-21) + 2.5log_10(0.7/θ_ eff) + 
        1.25log_10(t_ vis/30)-k_m(X-1),

where C_m is a constant that depend on the overall throughput of the system, m_ sky is the sky brightness in AB mag arcsec^-2, θ_ eff is the seeing in arcsec, t_ vis is the exposure time in seconds,

k is the atmospheric extinction coefficient, and X is the airmass. 
The default values of the parameters in the above equation per band are given in Table 2 in <cit.>. The magnitude error for N-years observation is computed by σ/Nn_ vis, where the mean number of visits per year n_ vis can be derived from Table 1 in <cit.>.

In this Appendix, we compare the LSST error model with the Rubin OpSim output as well as the Data Challenge 2 <cit.> dr6 magnitude error. We perform our tests on the specific OpSim version , and we use the 5-year observing conditions including: coadd 5σ point source depth (), single-visit 5σ point source depth (), sky brightness (), and number of visits ().

We begin by checking Eq. <ref> using OpSim MAF maps over the DC2 footprint. The various survey conditions m_ sky, θ_ eff, and X are taken as the median values over the 5-year period, and other parameters C_m, t_ vis, and k_m are taken as the default values from <cit.>. The results from Eq. <ref> are the 5σ PSF magnitude limit in each band per visit, and we compare it with two quantities: the median 5σ depth map, and the equivalent per-visit depth from the coadded map: m_5=m_5^ coadd-2.5log(√(N_ vis)), where N_ vis is the number of visits at each pixel. The results are shown in Fig. <ref>. We see that in general, m_5 predicted by Eq. <ref> tends to be brighter than that from OpSim, and the difference is larger considering the coadd depth than the median depth. It seems that except for i-band which has a slightly different slope from unity, the difference in all other bands can be fixed by introducing a correction to C_m. For example, for the median m_5 case, the shifts needed are δ C_m={ -0.053, 0.032, -0.063, 0.070, 0.057, 0.027 } for ugrizy, respectively.



We also explicity check whether the dependence of the airmass, seeing, and sky brightness are as expected in Eq. <ref> with the default parameters. This is shown in Fig. <ref>. In all these exercises, we test whether the dependencies of the particular survey condition with m_5 on the ensemble pixels, fixing all other dependence to a constant C which we fit to the ensemble. We see that the airmass and seeing are well captured by Eq. <ref>, although the dependence of m_5 on airmass is weak. The sky brightness relation is less well captured by Eq. <ref> especially for u and g. In general, however, we conclude that in absence of a depth map, one can estimate the unbiased m_5 for Rubin observation conditions using Eq. <ref> with a modification of the C_m parameters for each band.



We then check Eq. <ref> and <ref> with the DC2 DM catalogue, where the magnitude errors are obtained through the detection pipeline, thus supposed to be more realistic. In this case, we directly adopt the coadded depth as m_5. We also compute in the low SNR limit (Eq. <ref>) which allows us to check the fainter magnitudes. For the extended magnitude errors, we compare with the CModel magnitudes in DC2. 
This is shown in Fig. <ref>. We see that there is reasonable agreement for the PSF megnitude errors in most bands, except for the u-band, where the LSST error model predicts larger error compared to that measured in DC2. However, it is also noticeable that the DC2 error seems to be underestimated when comparing the observed magnitude to the truth. It is also noticeable that the LSST error model also predicts consistently larger error at the bright end.
When we add the extended error from the size of the galaxy (Eq. <ref>), we find that the scatter of the magnitude error at fixed magnitude is quite a bit larger than that measured by the cModel in DC2. 
Both the PSF magnitude error and the scatter for the extended error in DC2 can be matched by the LSST error model by simple scaling of the PSF magnitude error by a constant for each band, as well as scaling the galaxy size a_ gal, b_ gal. We emphasise that due to the known issues in the DC2 catalogue, we do not calibrate the LSST error model to DC2 in our analysis. However, it is worth bearing in mind what the differences are, and that one needs to calibrate the model with the real data.


























§ LENS AND SOURCE TOMOGRAPHIC BIN DETAILS





This section includes some supplementary information for the lens and source tomographic bins for the mock photometry sample, as discussed in Section <ref>.

Figure <ref> shows a similar plot as Fig. <ref>, but with the y-axis in logarithmic scale, and extended to z=3. Only tomographic bins 1, 3, and 5 are shown for visual clarity. This scaling enhances the small, high-redshift population for both lens and source galaxies. 








Table <ref> shows the summary statistics on photo-z performance for BPZ and FZBoost at the 10% shallowest i-band coadd depth (qtl=0) and deepest depth (qtl=9) for the 1-year, 3-years, and 5-years mock LSST data. The summary statistics are: median bias, standard deviation (STD), normalized Median Absolute Deviation (NMAD), and outlier fraction.
Tables <ref> and <ref> show the mean values of the various metrics over the depth quantiles, given the gold cut adjusted for each year. The metrics include mean galaxy number N̅_ gal and mean redshift of the tomographic bin ⟨ z ⟩ for both lens and source samples, and additionally the width metrics σ_z and W_z for lens samples. In the BPZ case, we include an additional case where we select objects with odds≥0.9.




§ VARIATION WITH OTHER SURVEY PROPERTIES



In the main analysis, we investigated the trend of galaxy number and redshift distribution as a function of the i-band coadd depth. We considered the i-band depth to be most impactful because it is the detection band, and fluxes in all other bands are measured with forced photometry based on the i-band detection. However, other survey properties can also be important. For example, u-band is important for the quality of photo-z estimation, so extreme variation in the u-band depth could cause additional scatter. The effective seeing could be another important factor, which directly impact the noise-to-signal for extended objects. 
We investigate the variation of galaxy number density and photo-z properties with these other survey properties in this section.




Table <ref> summarizes the mean and standard deviation of the coadd depth in the other five LSST bands and the median effective seeing for Y3 survey properties from Rubin OpSim baseline v3.3. The other years show a similar trend, although Y1 has a larger scatter. We see that there is a strong correlation between the i-band depth and these other survey properties. On average, a deeper i-band quantile also contains deeper coadd depth in all other five bands, as well as a smaller median effective seeing, with more scatter in the latter.  







To check the dependences of other survey properties, we sub-divide each of the i-band deciles into five sub-quantiles of another survey property (such as depth in another band), and check the variation of the metrics, i.e. number of objects N_ gal, mean redshift ⟨ z ⟩, and width of the redshift bin σ_z, with these properties. As a reference, we also compute and compare the variation with sub-quantiles of the i-band depth itself. 
In this section, we show two representative examples for source tomographic bins determined by FZboost photo-z: the sub-quantiles in coadd u-band depth and the i-band seeing, for the fainest, median, and deepest i-band deciles: qtl=0, 5, 9. 
In the results presented here, we over-plot the variation from the i-band depth sub-quantiles (as faint, dashed lines) on top of that from the other survey properties (as solid lines), for visual comparison. That is, one can read off the level of fluctuation from the deepest and shallowest u-band depth sub-bin, for example, and compare it with that from the deepest and shallowest i-band depth sub-bin.
It should be noted, however, that these reference i-band split cases have a different actual x-axis values from those shown in the plots.

The results are shown in Fig. <ref> and <ref> respectively. We see that in general, these trends are consistent with the i-band depth fluctuation for all three metrics: the deeper (smaller) the depth (seeing), the more objects included in the sample, the higher the mean redshift of the tomographic bin, and the larger σ_z. 
Also, qtl=0 has a significantly larger variation compared to qtl=9 in most cases. 
Compared to the trends in the i-band depth sub-bins, we see that the N_ gal variations are always less strong for other properties. This is understood as selections are primarily taken in i-band. The ⟨ z ⟩ variations for the u-band tightly follows the i-band, although the first bin can have slightly larger fluctuations. For seeing, on the other hand, the trend is quite different for qtl=0, where the smallest seeing does not always correspond to a higher mean redshift. This could happen because the seeing is not as well correlated with depth - there are more scatter in the coadd depth and seeing at the faint end. Finally, the variation in σ_z seems to be relatively minor in most cases.

From these exercises, we see that within each i-band decile, the number of objects and p(z) properties can still change significantly with other survey properties such as u-band depth and i-band seeing. Meanwhile, given that these quantities are also quite tightly correlated, we expect that a lot of these variations are also due to the co-variation of the i-band depth. Hence, our main analysis, by splitting into the i-band quantiles, should capture the level of variations of the metrics. However, if one wishes to apply this method in e.g. forward modeling, then co-variation of all bands need to be taken into account.



§ MAGLIM-LIKE LENS SAMPLE









In this section, we explore the impact of variable depth on a lens sample selected with the DES Y3 MagLim cuts <cit.>. Because this sample has a brighter cut, we relax the i-band signal-to-noise limit to SNR≥ 5. The sample is selected with 

    17<i<4 z_ phot+18,

where we use the FZBoost mode redshift as z_ phot. This cut reduces the number of lens sample significantly compared to our fiducial case, resulting in a total sample size of 3.67% of the baseline (Gold cut) lens sample. The true redshift distribution of each tomographic bin is shown in Fig. <ref>, where the dashed lines show those from the shallowest quantile, and the solid lines show those from the deepest. Notice that the distribution is less smooth due to the sparsity of the sample. Overall, thanks to the bright cut, the redshift distribution for each bin has a smaller tail compared to the baseline case, especially for the highest redshift bin. 

Figure <ref> shows the metrics for the variable depth, namely, the galaxy number, mean redshift, and width of the tomographic bin, as a function of the i-band depth. The panels (a) - (d) has the same style as, and should be compared to Fig. <ref> - <ref>. Again, we see a significantly milder, but visible, trend of these metrics with depth, owing to the bright magnitude cut. This shows that the variable depth effect can be greatly reduced, but not completed removed, by introducing a bright cut at the cost of sample size. 

Figure <ref> shows the effect propagated to the galaxy clustering two-point data vector, C_ℓ^gg. We followed the same procedure as in Section <ref>, and set the number density in each bin to be 0.135, 0.117, 0.156, 0.219, 0.267 arcmin^-2 to account for the reduction in the overall number density compared to the fiducial case. The impact of variable depth on C_ℓ^gg is also significantly reduced, especially for (4,4) and (5,5). However, the impact is not negligible still at ℓ<100. 







