# Redshift evolution and covariances for joint lensing and clustering studies with DESI Y1

Sihan Yuan<sup>1,2,★</sup> Chris Blake<sup>1,3,★</sup> Alex Krolewski,<sup>4,5</sup> Johannes Lange<sup>6</sup>, Jack Elvin-Poole,<sup>7</sup>

Alexie Leauthaud<sup>8</sup>, Joseph DeRose<sup>9</sup>, Jessica Nicole Aguilar,<sup>9</sup> Steven Ahlen<sup>10</sup>

Gillian Beltz-Mohrmann,<sup>11</sup> David Brooks,<sup>12</sup> Todd Claybaugh,<sup>9</sup> Axel de la Macorra<sup>13</sup>, Peter Doel,<sup>12</sup>

Ni Putu Audita Placida Emas,<sup>3</sup> Simone Ferraro<sup>1,9</sup>, Jaime E. Forero-Romero<sup>1,14</sup>

Cristhian Garcia-Quintero<sup>15</sup>, Enrique Gaztañaga,<sup>16</sup> Satya Gontcho A Gontcho<sup>1,9</sup>

Boryana Hadzhiyska<sup>9</sup>, Sven Heydenreich,<sup>8</sup> Klaus Honscheid,<sup>7</sup> Mustapha Ishak<sup>15</sup>, Shahab Joudaki<sup>17</sup>

Eric Jullo<sup>18</sup>, Theodore Kisner<sup>1,9</sup>, Anthony Kremin<sup>1,9</sup>, Andrew Lambert,<sup>9</sup> Martin Landriau<sup>1,9</sup>

Marc Manera<sup>19</sup>, Aaron Meisner<sup>20</sup>, Ramon Miquel,<sup>19</sup> Jundan Nie<sup>21</sup>

Nathalie Palanque-Delabrouille<sup>1,9</sup>, Claire Poppett,<sup>9</sup> Anna Porredon<sup>22</sup>, Mehdi Rezaie<sup>23</sup>

Ashley J. Ross,<sup>7</sup> Graziano Rossi,<sup>24</sup> Rossana Ruggeri<sup>1,3</sup>, EusebioSanchez<sup>25</sup>, Christoph Saulder<sup>26</sup>

Hee-Jong Seo<sup>27</sup>, Joseph Harry Silber<sup>1,9</sup>, Gregory Tarlń<sup>6</sup>, Mariana Vargas-Magaña<sup>28</sup>

Benjamin Alan Weaver,<sup>20</sup> Enia Xhakaj,<sup>8</sup> Zhimin Zhou<sup>21</sup> and Hu Zou<sup>22</sup>

*Affiliations are listed at the end of the paper*

Accepted 2024 July 20. Received 2024 June 17; in original form 2024 March 1

## ABSTRACT

Galaxy–galaxy lensing (GGL) and clustering measurements from the Dark Energy Spectroscopic Instrument Year 1 (DESI Y1) data set promise to yield unprecedented combined-probe tests of cosmology and the galaxy–halo connection. In such analyses, it is essential to identify and characterize all relevant statistical and systematic errors. We forecast the covariances of DESI Y1 GGL + clustering measurements and the systematic bias due to redshift evolution in the lens samples. Focusing on the projected clustering and GGL correlations, we compute a Gaussian analytical covariance, using a suite of N-body and lognormal simulations to characterize the effect of the survey footprint. Using the DESI one percent survey data, we measure the evolution of galaxy bias parameters for the DESI luminous red galaxy (LRG) and bright galaxy survey (BGS) samples. We find mild evolution in the LRGs in $0.4 < z < 0.8$, subdominant to the expected statistical errors. For BGS, we find less evolution for brighter absolute magnitude cuts, at the cost of reduced sample size. We find that for a redshift bin width $\Delta z = 0.1$, evolution effects on DESI Y1 GGL is negligible across all scales, all fiducial selection cuts, all fiducial redshift bins. Galaxy clustering is more sensitive to evolution due to the bias squared scaling. Nevertheless the redshift evolution effect is insignificant for clustering above the 1-halo scale of $0.1h^{-1}$ Mpc. For studies that wish to reliably access smaller scales, additional treatment of redshift evolution is likely needed. This study serves as a reference for GGL and clustering studies using the DESI Y1 sample.

**Key words:** methods: numerical – methods: statistical – galaxies: haloes – large-scale structure of Universe.

## 1 INTRODUCTION

The large-scale structure (LSS) in our Universe carries a wealth of information, and has become an essential probe in constraining fundamental physics and testing the current cosmological paradigm and galaxy formation. Accessing the information contained in the LSS requires summary statistics that capture the salient features in the observed density field. The galaxy autocorrelation function has long been the workhorse summary statistic, and has produced some of

the earliest LSS constraints on the underlying cosmological model (e.g. Percival et al. 2001; Hamilton & Tegmark 2002; Cole et al. 2005; Tegmark et al. 2006; Padmanabhan et al. 2007). In the last two decades, weak gravitational lensing has emerged as another powerful probe of cosmology and gravity that complements the information content of the galaxy autocorrelation function (e.g. Schrabback et al. 2010; Blake et al. 2016; Mandelbaum 2018; Heymans et al. 2021; Abbott et al. 2022; Leauthaud et al. 2022; Dalal et al. 2023; Lange et al. 2023).

Weak gravitational lensing refers to the correlated gravitational distortion induced in background galaxy shapes by foreground LSS, as their light travels towards us. This effect is a powerful cosmological

<footnote>* E-mail: sihany@stanford.edu (SY); cblake@swin.edu.au (CB)</footnote>
probe because it is sensitive to the geometry along the light path and the potential wells of the intervening structure. In principle, we can extract the cosmological information through autocorrelation of the observed shapes of galaxies, commonly referred to as cosmic shear (see Kilbinger 2015 for a review), or by cross-correlation of the shapes of the galaxies with positions of foreground galaxies, often referred to as galaxy-galaxy lensing (GGL; e.g. Brainerd, Blandford & Smail 1996; Mandelbaum et al. 2013). In this paper, we are primarily concerned with the combination of galaxy-galaxy lensing with galaxy autocorrelation functions, which improves the constraining power on cosmology by breaking degeneracies with galaxy bias (e.g. Coupon et al. 2015; Leauthaud et al. 2017; Yuan, Eisenstein & Leauthaud 2020; Pandey et al. 2022; Dvornik et al. 2023).

The Dark Energy Spectroscopic Instrument (DESI) is a spectroscopic galaxy survey with the primary goal of determining the nature of dark energy through the most precise measurement of the expansion history of the universe ever obtained (Levi et al. 2013; DESI Collaboration 2016; Abarashi et al. 2022). The baseline survey began in 2021 and will obtain spectroscopic measurements of 40 million galaxies and quasars in a 14000 square degrees footprint in five years. This represents an order-of-magnitude improvement both in the volume surveyed and the number of galaxies measured over previous surveys such as BOSS and eBOSS (extended Baryon Oscillation Spectroscopic Survey; Ahumada et al. 2020; Alam et al. 2020). The DESI large-scale structure samples are divided into four target classes: the bright galaxy survey (BGS) sample, the luminous red galaxies (LRG), the emission line galaxies, and the quasistellar objects, with increasingly high-redshift kernels. The LRG and BGS samples at lower redshift provide a clean and excellently calibrated foreground lens sample for galaxy-galaxy lensing studies. With these two samples, DESI will provide unprecedentedly high signal-to-noise measurement of galaxy-galaxy lensing over large overlap footprints with photometric surveys such as Hyper-Suprime Cam (HSC) survey (Aihara et al. 2018; Miyazaki et al. 2018; Mandelbaum et al. 2018a, b), the dark energy survey (DES; Dark Energy Survey Collaboration 2016; Park et al. 2016; Drlica-Wagner et al. 2018), the kilo-degree survey (KiDS-1000; Kuijken et al. 2019; Wright et al. 2020; Giblin et al. 2021; Heymans et al. 2021; Hildebrandt et al. 2021), and in the future, Euclid space telescope (Laureijs et al. 2011; Euclid Collaboration 2020) and Vera Rubin Observatories' Legacy Survey of Space and Time (LSST; Ivezic et al. 2008).

In order to robustly analyse the lensing and clustering measurements to the necessary precision in DESI, we need to understand and characterize all the relevant systematics, and accurately determine the covariance between the measured correlations. A leading systematic error arises when the observed lens galaxy samples evolve over the width of the redshift bin, but the model ignores such evolution and makes predictions at a fixed redshift. This simplification rides on the assumption that the targeted sample is sufficient uniform within a redshift bin, but a significantly evolving galaxy bias can break this assumption. Thus, when selecting a lens sample for cosmological analysis, it is essential to test this assumption and characterize any potential biases. Several studies have leveraged existing data to demonstrate significant redshift evolution in the galaxy bias of LRGs (Zhai et al. 2017; Lange et al. 2019; Pandey et al. 2023; Yuan et al. 2024). Specifically, Pandey et al. (2023) found that ignoring the evolution can lead to approximately $1\sigma$ biases in the cosmological constraints in future Rubin LSST analyses.

In this paper, we model the redshift evolution of the DESI BGS and LRG samples using the DESI one-percent survey clustering

measurements (DESI Collaboration 2023a, b). From our best-fitting models, we predict the systematic biases due to redshift binning and compare them to the statistical uncertainties of DESI Y1. We note several key trade-offs: (1) increasing the number of redshift bins reduces the evolution effect within each bin but results in noisier covariances; (2) a higher absolute magnitude cut results in a more uniform sample but at the cost of reduced sample size; (3) a higher maximum redshift in the lens sample increases the sample size but at the cost of increased intrinsic alignment contamination. Finally, we also address the needs of analyses that rely on pre-determined simulation snapshots for modelling. This paper accompanies Lange et al. (2024), which assesses additional theoretical systematic effects on the DESI GGL measurements.

We also present our determination of the analytical covariance for DESI Y1 analyses of projected galaxy-galaxy lensing and clustering statistics. This error forecast is required for both assessing the importance of systematic errors, such as due to lens evolution, as well as performing the future cosmological analysis. Our fiducial covariance is calculated using an analytical approach, following previous treatments by e.g. Singh et al. (2017), Dvornik et al. (2018), Shirasaki & Takada (2018), and Blake et al. (2020). In this paper we use N-body and lognormal simulations, sampling the same footprint as the data, to test and improve the efficacy of this approach.

This paper is organized as follows: In Section 2, we describe the DESI LRG and BGS samples and the photometric source samples used for this analysis. We present the analytic covariance matrix calculation in Section 3, and in Section 4 we describe the pipeline we use to model the redshift evolution of the lens samples. In Section 5, we predict the lens bin biases of the DESI Y1 LRG and BGS samples compared to the statistical noise for a fiducial set of redshift bins. We summarize our findings in Section 6.

## 2 DATA

The fiducial DESI Y1 galaxy-galaxy lensing analysis will utilize the LRG and BGS samples as lenses, and DES, KiDS, and HSC samples as sources. We briefly describe these samples in this section. We also introduce the relevant summary statistics for Y1 analysis: galaxy clustering and galaxy-galaxy lensing.

To model the redshift evolution of the lens samples, we make use of the DESI one percent survey data. DESI observed its one percent survey as the third and final phase of its science validation (SV) programme in 2021 April and May. Observation fields were chosen to be in 20 non-overlapping 'rosettes' selected to cover major data sets from other surveys, including the cosmic evolution survey (COSMOS), HSC, DES deep field, galaxy and mass assembly, Great Observatories Origins Deep Survey, and anticipated deep fields from future LSST and Euclid observations. Each rosette is observed at least 12 times, resulting in over 99 percent completeness for the LRG and BGS samples. The one percent survey essentially produces a smaller but more complete preview version of the upcoming DESI main sample, and is ideal for calibrating galaxy-halo connection models and building high-fidelity mocks. Full details of this sample are presented in DESI Collaboration (2023a, b).

### 2.1 DESI LRG

LRGs are essential for large-scale structure studies. LRGs are specifically selected for observations due to two main advantages: (1) they are bright galaxies with the prominent 4000 Å break in
<figure><img src="image_1.png" /><figcaption>**Figure 1.** The DESI lens density-redshift distribution as seen in the one-per cent sample. We also overplot the source galaxy distribution in each tomographic bin. The three panels correspond to three source galaxy data sets, where the grey vertical dashed lines represent the default lens bins. The short solid vertical lines in the bottom panel mark the ABACUS SUMMIT and AEMULUS simulation snapshots for the reference.</figcaption></figure>

their spectra, thus allowing for relatively easy target selection and redshift measurements; and (2) they are highly biased tracers of the large-scale structure, thus yielding a higher signal-to-noise (S/N) per-object for clustering measurements compared to typical galaxies. The LRG SV target selection is defined in Zhou et al. (2020). The DESI LRG sample has a target density of 605 square degrees in the redshift range $0.4 < z < 0.8$, significantly higher than previous LRG surveys (BOSS and eBOSS Dawson et al. 2013, 2016), while the sample also extends to $z \sim 1$. Within the DESI one percent survey, the LRG main sample consists of 89 059 galaxies, 43 269 in the northern footprint and 45 790 in the southern footprint. The shaded blue histogram in Fig. 1 shows the mean spatial density of the LRG sample as a function of redshift. We only consider the LRG

sample in the redshift interval $0.4 < z < 0.8$ for fiducial cosmology analyses.

## 2.2 DESI BGS

During DESI bright time, a BGS sample is observed, which has a ‘Bright’ and ‘Faint’ component, as well as Milky Way stars. In this work, we only consider the higher priority BGS Bright sample for simplicity and do not consider the faint sample as it suffers from complicated incompleteness and systematic effects. The BGS SV target selection is defined in Ruiz-Macias et al. (2021) and Hahn et al. (2023). Within the one percent survey, the BGS Bright sample consists of 115 602 galaxies in redshift range $0.1 < z < 0.4$. The
<figure><img src="image_2.png" /><figcaption>**Figure 2.** The distribution of the DESI one-per cent survey BGS Bright sample in *r*-band absolute magnitude and redshift. The thick lines indicate the redshift bins and magnitude cuts considered in this paper. Two bins are outlined with dashed borders to indicate that these bins are highly incomplete.</figcaption></figure>

mean number density as a function of redshift is shown through the orange shaded histograms in Fig. 1, whereas Fig. 2 shows the full distribution of BGS bright galaxies in the absolute magnitude-redshift plane. The thick blue lines show the fiducial selection cuts we adopt for this analysis. We have three *R*-band magnitude cuts at −19, −20, −21, and three redshift bins 0.1–0.2, 0.2–0.3, 0.3–0.4. The dashed lines show selections that are highly incomplete: $M_R < -19$ in 0.2 < *z* < 0.3, and $M_R < -20$ in 0.3 < *z* < 0.4. These cuts are suboptimal as they do not provide meaningful statistical gains over brighter more complete selections in the same redshift bins. Thus, these two selections are likely not relevant for our final analysis setup. We include these bins in our results for completeness.

Given that the BGS is an apparent-magnitude selected sample, a key consideration when attempting to define a homogeneous subset of lens galaxies is the absolute magnitude cut. For the DESI fiducial large-scale BGS clustering analyses, a conservative magnitude cut of $M_R < -21.5$ is sufficient, as it renders the Poisson noise negligible on the large scales relevant to Baryon acoustic oscillation analysis. However, given galaxy–galaxy lensing analyses are noise-limited rather than sample variance-limited, we intend to retain a much larger fraction of the BGS sample to maximize statistical precision. Fig. 3 shows how different absolute magnitude cuts change the mean number density of the BGS sample as a function of redshift. Applying a $M_R < -19$ cut removes approximately 3.3 per cent of the sample; a $M_R < -20$ cut removes approximately 31 per cent of the sample; $M_R < -21$ cut removes approximately 76 per cent of the sample. The absolute magnitudes are defined in eq. (6) in DESI Collaboration (2023b). Fig. 3 shows the mean number density as a function of redshift for these different absolute magnitude cuts. A brighter magnitude cut results in a more uniform $n(z)$, but at the cost of removing a large fraction of the sample. We will analyse the effects of these choices on the resulting statistical and systematic errors in this study.

## 2.3 Photometric sources

In this study, we consider a cross-correlation analysis of these DESI lens samples with three weak lensing data sets with significant

<figure><img src="image_3.png" /><figcaption>**Figure 3.** The DESI BGS mean density as a function of redshift for different magnitude cuts. The legend also states the BGS sample size for different magnitude cuts, within the fiducial redshift range of 0.1 < *z* < 0.4.</figcaption></figure>

overlap with DESI: the DES Year 3 data set (DES-Y3; Gatti et al. 2021), the current KiDS data set (KiDS-1000; Giblin et al. 2021), and the HSC Year 1 data set (HSC-Y1; Mandelbaum et al. 2018a).

The dark energy survey has used the Dark Energy Camera at the Blanco 4-m telescope to image 5000 square degrees in *grizY* filters. The Y3 sample comprises data taken during the first three years of DES operations, covering the full footprint with approximately half the final exposure time, resulting in a catalogue with depth $i \approx 23$. The DES-Y3 shear catalogue contains approximately 100 million galaxies covering an area of 4143 square degrees, with effective source number density $n_{\text{eff}} = 5.6 \text{ arcmin}^{-2}$ and shape noise $\sigma_e = 0.27$. The sample is divided for tomographic analysis into four photometric redshift bins with median redshifts $z = [0.34, 0.52, 0.74, 0.96]$ (Amon et al. 2022), whose redshift distributions are calibrated by self-organizing maps and further constrained by clustering cross-correlations (Buchs et al. 2019; Myles et al. 2021; Gatti et al. 2022). We adopt these redshift distributions in our forecasts described below.

The kilo-degree survey is a European Southern Observatory public survey using a combination of optical imaging in *ugri* bands from the 2.6-m VLT Survey Telescope with depth $r \approx 25$, and overlapping near-infrared imaging in *ZYJHKs* bands from the VISTA-VIKING survey. The KiDS-1000 lensing catalogue is based on the fourth KiDS data release, comprising 1006 square degrees of imaging data. The sample consists of 21 million galaxies with $n_{\text{eff}} = 6.2 \text{ arcmin}^{-2}$ and $\sigma_e = 0.27$ (Giblin et al. 2021), divided into five photometric bins for tomographic analysis by photometric redshift $z_B$ derived using the BPZ method, with bin limits $z_B = [0.1, 0.3, 0.5, 0.7, 0.9, 1.2]$. The redshift distributions of these source samples were determined by the KiDS collaboration using the self-organizing map (SOM) methodology (Hildebrandt et al. 2021).

Finally, Hyper Suprime-Cam is a wide-field imaging camera on the 8.2-m Subaru telescope, which has been used to conduct a 6-yr deep, multiband imaging survey. The Wide layer, which is designed for weak lensing cosmology, aims at covering 1400 square degrees in *grizy* bands with depth $r \approx 26$. The HSC first-year shear catalogue is based on about 90 nights of HSC Wide data covering 137 square degrees with depth $i \approx 24.5$, and contains about nine million galaxies
with $n_{\text{eff}} = 17.6 \text{ arcmin}^{-2}$ (Mandelbaum et al. 2018a). The source sample is split into tomographic bins by photometric redshifts with bin limits $z = [0.3, 0.6, 0.9, 1.2, 1.5]$, with the redshift distributions for each bin estimated using the COSMOS 30-band photo-z catalogue (Ilbert et al. 2009).

## 2.4 Summary statistics

For this analysis we consider two summary statistics: the projected 2-point galaxy correlation function (2PCF), $w_p$, and the surface mass density contrast, $\Delta\Sigma$, measured by galaxy–galaxy lensing. We start by introducing the 3D 2PCF, which can be compressed into a 2D function due to rotational symmetry. We can express this 2D correlation function in terms of transverse separation $r_p$ and line-of-sight (LoS) separation $r_\pi$, $\xi(r_p, r_\pi)$, which we can compute using the Landy & Szalay (1993) estimator:

$$ \xi(r_p, r_\pi) = \frac{\text{DD} - 2\text{DR} + \text{RR}}{\text{RR}}, \qquad (1) $$

where DD, DR, and RR are the normalized numbers of data–data, data–random, and random–random pair counts in each bin of $(r_p, r_\pi)$. The redshift-space $\xi(r_p, r_\pi)$ represents the full information content of the 2PCF.

To avoid the complexities of modelling the small-scale finger-of-god effect, we can further compress $\xi(r_p, r_\pi)$ to the projected galaxy 2PCF, $w_p$, which is the line-of-sight integral of $\xi(r_p, r_\pi)$,

$$ w_p(r_p) = 2 \int_{0}^{r_{\pi, \text{max}}} \xi(r_p, r_\pi) \, d\pi. \qquad (2) $$

Whilst $w_p$ is strictly less informative than $\xi(r_p, r_\pi)$ as it excludes the velocity information encoded in the LoS clustering, it is easier to model as it avoids non-linear velocity effects, and its covariance matrix is more tractable owing to the data compression. It is common for combined-probe analyses to consider $w_p$ in addition to weak lensing measurements in order to break parameter degeneracies, for example associated with galaxy bias.

The galaxy–galaxy lensing observable we use is the mean surface mass density contrast profile $\Delta\Sigma$, defined as:

$$ \Delta\Sigma(r_p) = \bar{\Sigma}(\langle r_p \rangle) - \bar{\Sigma}(r_p), \qquad (3) $$

where $\bar{\Sigma}(r_\perp)$ is the azimuthally averaged and projected surface mass density at radius $r_p$ and $\bar{\Sigma}(\langle r_p \rangle)$ is the mean projected surface mass density within projected separation $r_p$ (Miralda-Escude 1991; Wilson et al. 2001; Leauthaud et al. 2017):

$$ \bar{\Sigma}(\le r_p) = \frac{2}{r_p^2} \int_0^{r_p} \bar{\Sigma}(R') R' dR'. \qquad (4) $$

# 3 ANALYTICAL COVARIANCES

## 3.1 Computation of analytical covariance for projected correlations

The assessment of the significance of systematics due to lens evolution, as well as the statistical analysis of the survey measurements, depends on the covariance of the measured correlation functions, (for example) $\Delta\Sigma_i$ and $\Delta\Sigma_j$ for lens–source combinations $i$ and $j$ between two scales $r_{p,k}$ and $r_{p,l}$,

$$ \begin{aligned} \text{Cov} [\Delta\Sigma_i(r_{p,k}), \Delta\Sigma_j(r_{p,l})] \\ = \langle \Delta\Sigma_i(r_{p,k}) \Delta\Sigma_j(r_{p,l}) \rangle - \langle \Delta\Sigma_i(r_{p,k}) \rangle \langle \Delta\Sigma_j(r_{p,l}) \rangle. \end{aligned} \qquad (5) $$

We use analytical methods to determine the covariances between different statistics, redshifts, and scales for the DESI BGS and LRG

samples and the DES, KiDS, and HSC weak lensing data sets. Our covariance computation follows the methods of Blake et al. (2020) (see also Singh et al. 2017; Dvornik et al. 2018; Shirasaki & Takada 2018) and is fully described in Appendix A; we provide a brief summary here.

We compute a Gaussian analytical covariance including the sample variance, noise, and mixed contributions. We assume a fiducial non-linear matter power spectrum, the measured DESI galaxy redshift distributions in each bin, and the weak lensing source survey configurations in tomographic bins defined for DES-Y3 (Amon et al. 2022), KiDS-1000 (Giblin et al. 2021), and HSC-Y1 (Hikage et al. 2019) including the shape noise and effective source number density in each bin. We evaluate the covariance in 15 logarithmic bins of projected separation in the range $0.08 < r_p < 80 h^{-1} \text{ Mpc}$. For the $w_p$ covariance, we set $r_{\pi, \text{max}} = 100 h^{-1} \text{ Mpc}$. We assume the linear bias factors determined for early DESI data by Prada et al. (2023). The covariance is noise-dominated on small scales (less than several $h^{-1} \text{ Mpc}$ for lensing measurements, depending on the lens and source configuration), so the assumption of linear bias does not have a significant impact on these scales. We exclude the non-Gaussian contributions to the covariance, noting that studies such as Joachimi et al. (2021) show that these are negligible for lensing studies, given current data sets. We evaluate the covariances for the forecast overlap areas of DESI-Y1 and these weak lensing surveys: (716, 456, 142) square degrees for DES-Y3, KiDS-1000, and HSC-Y1 overlapping with BGS, and (845, 455, 153) square degrees with LRGs.

## 3.2 Covariance comparisons with N-body simulations and jack-knife errors

We tested the analytical covariance determination for the surface mass density contrast $\Delta\Sigma(r_p)$ using the Buzzard DESI-Lensing simulations (DeRose et al. prep). These simulations build realistic DESI target and weak lensing source populations within the Buzzard simulation suite (DeRose et al. 2019a). The DESI targets are included via a halo occupation distribution prescription, and the sources are populated by a statistical method which includes photometric redshift errors, magnitude distributions, source weights, and multiplicative shear calibration bias. We refer the reader to DeRose et al. (prep) for a full description of these simulations.

We created mock source and lens samples matched to the angular footprint of the overlapping DESI-Y1 and weak lensing data sets. We determined these ‘survey masks’ from the angular completeness maps of the DESI-Y1 redshift catalogues (where the completeness is defined relative to the final DESI target density), intersected with the footprint of each weak lensing data set, using a Healpixn\_side = 1024 pixelization. The window functions in each case are illustrated for the three weak lensing surveys in Fig. 4, where we note that DESI-Y1 BGS and LRG observations have different completeness maps, such that we treat these two redshift ranges separately. For the case of HSC, we have artificially arranged the six survey regions more compactly in Fig. 4. Tiling the irregular geometries of these intersections within a Buzzard quadrant allows (8, 3, 12) realizations to be extracted for the (KiDS, DES, HSC) footprints, producing a total of (64, 24, 96) realizations across the eight separate Buzzard simulations.

For the different combinations with the source samples, we measured the surface mass density contrast $\Delta\Sigma(r_p)$ around the lenses, using the same measurement code as described in Blake et al. (2020). Specifically, we applied the $\Delta\Sigma(r_p)$ estimator assuming a single spectroscopic redshift distribution for each source sample, and therefore we do not utilize the individual source photo-
<figure><img src="image_4.png" /><figcaption>**Figure 4.** The window functions applied when creating the ensemble of Buzzard mocks including survey overlap masks. Window functions are illustrated for the DESI-Y1 BGS (left-hand column) and LRG observations (right-hand column), using the DESI completeness maps in a healpix side = 1024 pixelization, where the colour bar indicates the DESI completeness level. The different rows display the intersection of these DESI completeness maps with the KiDS-1000, DES-Y3, and HSC-Y1 survey footprints. The axes indicate separation in degrees relative to the field centre, where for the case of HSC-Y1, we have artificially arranged the six widely separated survey regions in a more compact format. The areas shown in the panel titles correspond to the total area of intersection between the DESI and lensing survey footprints.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 5.** The 'tangential' and 'cross' components of the surface mass density contrast $\Delta\Sigma(r_p)$ of the five tomographic source samples of the KiDS-1000 mocks, around the three lens samples of the DESI-Y1 BGS mocks, for the simulations including the overlap mask. The different panels display measurements for different combinations of source and lens samples, as indicated by the captions above and to the right of each panel. We plot the mean and standard deviation of the measurements across 64 Buzzard masked regions. The solid lines show the fiducial cosmological models, where a linear bias factor has been jointly fit to each lens redshift slice using the separation range $r_p > 5 h^{-1} \text{ Mpc}$ (indicated by the un-shaded regions). The units of $\Delta\Sigma(r_p)$ are $h M_{\odot} pc^{-2}$, and the y-axes are scaled by a factor of $r_p$ (in $h^{-1} \text{ Mpc}$) for clarity of presentation.</figcaption></figure>

z estimates. We estimated $\Delta\Sigma(r_p)$ in 15 logarithmically spaced projected separation bins in the range $0.08 < r_p < 80 h^{-1} \text{ Mpc}$, and we corrected our results for the multiplicative shear calibration bias introduced in the mocks. We estimated both the 'tangential' and

'cross' components of $\Delta\Sigma(r_p)$. For the tangential (or 'E-mode') component, for each source–lens pair included in the estimate, the source shape is projected perpendicular to the source–lens separation vector, measuring the contribution of weak gravitational lensing. For
<figure><img src="image_6.png" /><figcaption>**Figure 6.** A comparison between error estimates of the tangential component of the surface mass density contrast $\Delta\Sigma(r_p)$ of the five tomographic source samples of KiDS-1000, around the three lens samples of DESI-BGS, for the Buzzard mocks including survey overlap masks. The different panels display error estimates for different combinations of source and lens samples. We compare the error in $\Delta\Sigma(r_p)$ predicted by the analytical covariance (solid black line), with the standard deviation of the measurements across 64 Buzzard regions (the green band, which indicates the error in the inferred standard deviation arising from the limited number of regions), with the jack-knife error (dashed red line).</figcaption></figure>

the cross (or 'B-mode') component, the source shapes are first rotated by 45°, which serves as a useful systematic error test which we will apply to the real data. We refer the reader to Heydenreich et al. (prep) for a full discussion of $\Delta\Sigma$ estimates in the context of the DESI Y1 data.

We illustrate our results using the example of the masked simulations of DESI-Y1 BGS and KiDS-1000, noting that the analyses of the other survey configurations reach similar conclusions. Fig. 5 displays the tangential and cross $\Delta\Sigma(r_p)$ measurements between the five tomographic source samples of the KiDS-1000 mocks, around the three lens samples of the DESI-Y1 BGS mocks ($0.1 < z < 0.2$, $0.2 < z < 0.3$, and $0.3 < z < 0.4$). We plot the mean and standard deviations of the measurements across the 64 realizations, compared to the fiducial cosmological models in each case, scaled by a linear galaxy bias factor jointly fit to the data at each lens redshift. The models are a good description of the measurements, and (in these mock catalogues) we detect no evidence of a 'cross' component.

Fig. 6 displays a comparison between different estimates of the error in these $\Delta\Sigma$ measurements. In this figure we compare the estimate from the analytical covariance (solid black line), with the standard deviation across the 64 mock realizations (represented by the green band, whose width indicates the 'error in the error' due to the limited number of realizations), with the jack-knife error estimated using 100 jack-knife regions containing equal number of lenses (dashed red line). The jack-knife regions are defined by boundaries of constant right ascension and declination and have average areas between 1.4 and 8.5 square degrees depending on the survey. The analytical covariance is determined from the Gaussian contribution, together with the noise correction using the measured number of source–lens pairs. Generally speaking, the different error estimates agree within 10 per cent to 20 per cent, depending on scale and redshift.

Fig. 7 displays the full correlation matrix of the $\Delta\Sigma$ data vector, when the different redshift and separation bins are concatenated together, where the analytical covariance is displayed in the upper left triangle of both panels, and the numerical and jack-knife covariance in the lower right triangles. We see a similar qualitative structure in the different estimates of the covariance matrix, which features significant correlation between measurements for different source samples and the same lenses, at adjacent large separations, and negligible correlation between measurements for different lens samples. Fig. 8 compares estimates of the correlation coefficients for covariance matrix elements sharing the same scale bin, for different lens and source bins offset from the central diagonal.

In Appendix B, we use lognormal mocks to investigate the effects of the survey footprint on the galaxy–galaxy lensing covariance. The window functions of the source and lens samples can alter the sample variance contribution to the Gaussian covariance on large scales (Kilbinger & Schneider 2004; Sato et al. 2011; Shirasaki et al. 2019; Friedrich et al. 2021; Joachimi et al. 2021). We show that these effects are more significant for the cross-component of the galaxy–galaxy lensing, which may be used in a null test for the presence of B-modes. We find that for the DESI-Y1 footprint geometry shown in Fig. 4, the error in the tangential (cross) component increases with reference to the analytical prediction by 10 per cent (30 per cent) at the scales at which 50 per cent of source–lens pairs are ‘lost’ by the boundary. We describe our corrections to the covariance for these effects in Appendix B.

## 4 MODELLING LENS BIN EVOLUTION

In this section, we describe how we create our model for the redshift evolution of the DESI lens samples. Our pipeline first analyses the DESI one percent survey data to constrain the galaxy bias parameters as a function of redshift. Then we produce high-fidelity redshift-
<figure><img src="image_7.png" /><figcaption>**Figure 7.** The correlation matrix corresponding to the covariance of the tangential ΔΣ measurements of the KiDS-1000 and DESI-Y1 BGS mocks, including the survey masks. The measurements use $N_{\text{tom}} = 5$ source tomographic samples, $N_{\text{lens}} = 3$ lens samples and $N_{\text{sep}} = 15$ separation bins, which are concatenated into a data vector of length $N_{\text{lens}} \cdot N_{\text{tom}} \cdot N_{\text{sep}} = 225$, where the innermost loop is over separations, and the outermost loop is over lens bins. The analytical covariance is depicted in the upper-left triangle of both panels. The lower-right triangles display the numerical covariance derived from 64 realizations (left-hand panel) and the average jack-knife covariance of the 64 realizations (right-hand panel). The correlation matrix is divided into sections corresponding to the different lens bins, indicated by the horizontal and vertical dotted lines.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 8.** The correlation coefficient of the off-diagonal covariance of the tangential ΔΣ measurements of the KiDS-1000 and DESI-Y1 BGS mocks, $r(i, j) = C_{ij} / \sqrt{C_{ii} C_{jj}}$. The different panels display the full set of correlation coefficients for matrix elements sharing the same scale bin, for different lens and source bins offset from the central diagonal. We compare the correlation coefficient predicted by the analytical covariance (solid black line), the Buzzard realizations (green band) and the jack-knife covariance (dashed red line).</figcaption></figure>

dependent mocks using the best-fitting galaxy bias model. Finally, we compute the summary statistics and the biases associated with different lens bin choices, whose significance we assess using the analytical covariances.

## 4.1 HOD model

For galaxy bias modelling, we adopt the Halo Occupation Distribution model (HOD), which probabilistically populates dark matter haloes with galaxies according to a set of halo properties. Statistically,

the HOD can be summarized as a probability distribution $P(n_g|X_h)$, where $n_g$ is the number of galaxies of the given halo, and $X_h$ is some set of halo properties.

In the vanilla HOD model, halo mass is assumed to be the only relevant halo property $X_h = M_h$ (Zheng et al. 2005; Zheng, Coil & Zehavi 2007). This vanilla HOD separates the galaxies into central and satellite galaxies, and assumes the central galaxy occupation follows a Bernoulli distribution whereas the satellites follow a Poisson distribution.

For both the LRG and the luminosity-limited BGS samples, the HOD is well approximated by a vanilla model given by Zheng et al. (2007), Zehavi et al. (2011), and Guo et al. (2015):

$$ \bar{n}_{\text{cent}}^{\text{LRG}}(M) = \frac{f_{\text{ic}}}{2} \text{erfc} \left[ \frac{\log_{10}(M_{\text{cut}}/M)}{\sqrt{2\sigma}} \right], \quad (6) $$

$$ \bar{n}_{\text{sat}}^{\text{LRG}}(M) = \left[ \frac{M - \kappa M_{\text{cut}}}{M_1} \right]^{\alpha} \bar{n}_{\text{cent}}^{\text{LRG}}(M), \quad (7) $$

where the five vanilla parameters characterizing the model are $M_{\text{cut}}$, $M_1$, $\sigma$, $\alpha$, $\kappa$. We describe these parameters as follows: $M_{\text{cut}}$ roughly controls the mean halo mass of central galaxies; $M_1$ sets the typical halo mass that hosts one satellite galaxy; $\sigma$ controls the steepness of the transition from 0 to 1 in the number of central galaxies; $\alpha$ is the power-law index on the number of satellite galaxies; and $\kappa M_{\text{cut}}$ gives the minimum halo mass to host a satellite galaxy. We have added a modulation term $\bar{n}_{\text{cent}}^{\text{LRG}}(M)$ to the satellite occupation function to remove satellites from haloes without centrals.

We have also included an incompleteness parameter $f_{\text{ic}}$, which is a downsampling factor controlling the overall number density of the mock galaxies. This parameter is relevant when trying to account for sample incompleteness and match the observed mean density of the galaxies in addition to clustering measurements. While we expect the one percent survey to be fairly complete in terms of fibre coverage, incompleteness can still arise due to target selection. By definition, $0 < f_{\text{ic}} \le 1$.

For this analysis, we use the ABACUSHOD code to find best-fitting HODs and sample HOD posteriors. ABACUSHOD is a highly efficient
HOD implementation that enables a large set of HOD extensions (Yuan et al. 2021). The code is publicly available as a part of the ABACUSUTILS package at https://github.com/abacusorg/abacusutils. Example usage can be found at https://abacusutils.readthedocs.io/en/latest/hod.html.

The HOD code is implemented on top of the ABACUSSUMMIT simulation suite, which is a set of large, high-accuracy cosmological N-body simulations using the ABACUSN-body code (Garrison, Eisenstein & Pinto 2019; Garrison et al. 2021; Maksimova et al. 2021). For this analysis, we use the base simulation box at Planck 2018 cosmology AbacusSummit_base_c000_ph000 (Planck Collaboration 2020). The box contains 6912³ particles within a (2 h⁻¹Gpc)³ volume, which yields a particle mass of 2.1 × 10⁹ h⁻¹M☉.<fnref n="1" />

The dark matter haloes are identified with the COMPASO halo finder, which is a highly efficient on-the-fly group finder specifically designed for the ABACUSSUMMIT simulations (Hadzhiyska et al. 2022). In addition to determining the number of galaxies per halo, the standard HOD model also dictates the position and velocity of the galaxies. In the vanilla model, the position and velocity of the central galaxy are set to be the same as those of the halo centre, specifically the L2 subhalo centre-of-mass for the COMPASO haloes. For the satellite galaxies, they are randomly assigned to halo particles with uniform weights, each satellite inheriting the position and velocity of its host particle.

To extract the redshift evolution, we conduct standard HOD fits on the sample clustering and number density in multiple redshift bins centred around ABACUSSUMMIT output snapshots. Yuan et al. (2024) describes the LRG HOD fits in detail, whereas the BGS HOD fits are described in Section 5.2 in this paper. The best-fitting parameters are then interpolated as a function of redshift to derive an approximate redshift evolution model, which we then implement on simulations using Dark Emulator (see following subsection) to generate clustering and lensing observables that account for redshift evolution.

## 4.2 Dark emulator

Having obtained the best-fitting HOD parameters in each redshift bin, we interpolate the best-fitting parameters as a function of redshift with simple polynomials to build a simple empirical redshift-dependent galaxy bias model. Then, we apply these models to the Dark Emulator dark matter model to generate realistic redshift-dependent mocks. Dark Emulator essentially provides an emulated model of the halo power spectrum as a function of cosmology, redshift, halo mass, and scale (Nishimichi et al. 2019; Kobayashi et al. 2022). It is trained on a large suite of dark matter only simulations evaluated at 100 different cosmologies.

Dark Emulator then computes the galaxy correlation and galaxy-galaxy lensing statistics by convolving the vanilla HOD with the emulated halo power spectrum and Fourier transforming back to configuration space. Using Dark Emulator, we evaluate the $w_p$ and $\Delta\Sigma$ statistics along a dense grid of lens redshifts ($\Delta z = 0.01$), with input redshift-dependent HOD parameters evaluated from the interpolated HOD(z).

To measure evolution biases from assuming a fixed HOD, we first compute the 'real' measurement $\Delta\Sigma_{evol}$ (or $w_{p, evol}$), which we model by integrating the $\Delta\Sigma$ along the dense redshift grid over the full redshift bin, weighted by the $n(z)$ of the observed sample. Then, we

<table><caption>**Table 1.** LRG HOD posterior means as found in table 3 of Yuan et al. (2024). The masses are in units of $h^{-1}M_{\odot}$.</caption><thead><tr><th>Sample</th><th colspan="2">LRG</th></tr></thead><tbody><tr><td>z range</td><td>0.4 &lt; z &lt; 0.6</td><td>0.6 &lt; z &lt; 0.8</td></tr><tr><td>log M<sub>cut</sub></td><td>12.89</td><td>12.79</td></tr><tr><td>log M<sub>1</sub></td><td>14.08</td><td>13.88</td></tr><tr><td>σ</td><td>0.27</td><td>0.21</td></tr><tr><td>α</td><td>1.20</td><td>1.07</td></tr><tr><td>κ</td><td>0.65</td><td>1.4</td></tr></tbody></table>

compute the naive model measurement $\Delta\Sigma_{fixed}$ (or $w_{p, fixed}$), which assumes a fixed HOD at a fixed redshift. For analytic theories, this fixed redshift can be approximated as the mean redshift of the sample in the corresponding redshift bin. Thus, $\Delta\Sigma_{fixed}$ is computed by simply evaluating our redshift-dependent model at the mean redshift. More precisely,

$$ \Delta\Sigma_{fixed} = \Delta\Sigma_{DE}(z_{mean}, \text{HOD}(z_{mean})), \quad (8) $$

$$ \Delta\Sigma_{evol} = \int dz \, n(z) \Delta\Sigma_{DE}(z, \text{HOD}(z)), \quad (9) $$

where $\Delta\Sigma_{DE}$ is the $\Delta\Sigma$ model from Dark Emulator, $n(z)$ is the lens redshift distribution, and $z_{mean}$ is the mean redshift in the corresponding redshift bin given $n(z)$.<fnref n="2" /> We then use comparisons between these $\Delta\Sigma$ evaluations to determine the impact of these assumptions from different lens bin configurations, characterizing the lens bin bias as,

$$ \delta_{z-evol} = \frac{|\Delta\Sigma_{fixed} - \Delta\Sigma_{evol}|}{\Delta\Sigma_{evol}}. \quad (10) $$

We can repeat this exercise for $w_p$ to calculate the evolution effect on the projected clustering. We note that assuming $z_{mean}$ for the fixed redshift is a conservative scenario, and real analyses may choose to further tune the assumed redshift.

# 5 SYSTEMATICS DUE TO LENS EVOLUTION

In this section, we compare biases due to redshift evolution to the expected statistical error of the DESI Y1 LRG and BGS samples in a set of fiducial redshift bins.

## 5.1 LRG

To model the redshift evolution of the LRG sample, we conduct HOD fits in the two fiducial redshift bins 0.4 < z < 0.6 and 0.6 < z < 0.8. We refer the readers to Section 4.1 for the model description and Yuan et al. (2024) for the detailed HOD analysis. We use the posterior means quoted in table 3 of Yuan et al. (2024) for this analysis. We reproduce the relevant numbers in Table 1.

We build a simple model of LRG redshift evolution HOD(z) by linearly interpolating the best-fitting values of log $M_{cut}$ and log $M_1$ as a function of redshift. The parameters $\alpha$, $\sigma$, and $\kappa$ are poorly constrained and are thus held fixed at their posterior mean values in the lower redshift bin. We then feed the HOD(z) to Dark Emulator to produce a redshift-dependent LRG mock from z = 0.4 to z = 0.8 and compute the biases in $w_p$ and $\Delta\Sigma$ due to redshift evolution using

<footnote n="2">²We note that assuming mean sample redshift is fairly naive, and additional differences between model and data could be accounted for by shifting the mean redshift. Nevertheless, this assumption is useful in presenting a worst case scenario.</footnote>

<footnote n="1">¹For more details, see https://abacussummit.readthedocs.io/en/latest/abacussummit.html</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 9.** The ratio of the error induced by LRG redshift evolution $\delta_{z-\text{evol}}$ divided by the statistical error $\sigma_{\text{stat}}$, as a function of projected scale $r_p$. Note that we show the cumulative measure above a minimum projected scale $r_p$. The two panels correspond to the two fiducial LRG lens bins. The orange curve shows $w_p$, whereas the blue curve shows $\Delta\Sigma$. The width of the two curves indicate the spread between different photometric source samples (DES Y3, KiDS-1000, and HSC). It is clear that for LRGs, the redshift evolution in the fiducial redshift bins is insignificant compared to the statistical error.</figcaption></figure>

<table><caption>**Table 2.** BGS HOD best-fits across four redshift bins and three absolute magnitude cuts. The best-fitting $w_p$ predictions are shown in Fig. 10. Note that we use a non-fiducial redshift range (offset from fiducial bins by 0.05) for the HOD analysis to accommodate the fact that our ABACUS SUMMIT simulation snapshots are at $z = 0.1, 0.2, 0.3, 0.4$. We then use these best-fitting values to derive the redshift evolution of the sample in the fiducial DESI redshift bins. We have fixed three HOD parameters $\sigma = 0.1, \alpha = 0.8, \kappa = 1.0$ for stable fits given the limited S/N in the data vector. We have skipped $M_R < -19$ at $z = 0.4$ since it is the same sample as $M_R < -20$ at $z = 0.4$. The masses are in units of $h^{-1}M_{\odot}$. These best-fitting HODs are also visualized in Fig. 11.</caption><thead><tr><th>z range</th><th>0.05 &lt; z &lt; 0.15</th><th>0.15 &lt; z &lt; 0.25</th><th>0.25 &lt; z &lt; 0.35</th><th>0.35 &lt; z &lt; 0.45</th></tr></thead><tbody><tr><td>Sample</td><td colspan="4">BGS M<sub>R</sub> &lt; -19</td></tr><tr><td>log M<sub>cut</sub></td><td>11.59</td><td>11.87</td><td>12.19</td><td>-</td></tr><tr><td>log M<sub>1</sub></td><td>12.29</td><td>12.68</td><td>13.14</td><td>-</td></tr><tr><td>Sample</td><td colspan="4">BGS M<sub>R</sub> &lt; -20</td></tr><tr><td>log M<sub>cut</sub></td><td>11.89</td><td>12.04</td><td>12.20</td><td>12.39</td></tr><tr><td>log M<sub>1</sub></td><td>12.71</td><td>12.85</td><td>13.15</td><td>13.17</td></tr><tr><td>Sample</td><td colspan="4">BGS M<sub>R</sub> &lt; -21</td></tr><tr><td>log M<sub>cut</sub></td><td>12.59</td><td>12.59</td><td>12.92</td><td>13.12</td></tr><tr><td>log M<sub>1</sub></td><td>13.69</td><td>13.66</td><td>13.50</td><td>13.83</td></tr></tbody></table>

equation (10). Because the $n(z)$ of the LRG is fairly flat, the mean redshift of LRGs is fairly close to the middle of the redshift bin. Specifically, we obtain $z_{\text{mean}} \approx 0.509$ for the $0.4 < z < 0.6$ sample, and $z_{\text{mean}} \approx 0.705$ for the $0.6 < z < 0.8$ sample.

We compare the redshift evolution systematic error $\delta_{z-\text{evol}}$ with the expected statistical error from the same redshift bin, calculated as the cumulative noise-to-signal of the summary statistic, integrating from large to small scales:

$$ \sigma_{\text{stat}}(r_p) = \sqrt{\mathbf{x}^T(r > r_p) C^{-1} \mathbf{x}(r > r_p)}, \quad (11) $$

where $\mathbf{x}$ is the target data vector ($w_p$ or $\Delta\Sigma$), and $C$ is the corresponding analytic covariance matrix forecast in Section 3. One can simply interpret this measure as how accurately we can determine the amplitude of the summary statistics as increasingly smaller scales are incorporated.

Fig. 9 shows the ratio of $\delta_{z-\text{evol}}$ to the statistical error $\sigma_{\text{stat}}$ as a function of the minimum projected scale $r_p$ included in the analysis. We draw a dashed line at $\delta_{z-\text{evol}}/\sigma_{\text{stat}} = 1$, above which redshift evolution becomes a significant systematic effect. The two

panels correspond to the two redshift bins, whereas the two colours correspond to the summary statistics $w_p$ and $\Delta\Sigma$. The widths of the band correspond to the spread between different source samples (DES Y3, KiDS-1000, and HSC). It is clear that for LRGs, the systematic error due to redshift evolution is insignificant given DESI Y1 statistics. This is true for both fiducial lens bins and for all source samples. This statement also holds down to scales as small as $0.1 h^{-1} \text{ Mpc}$.

## 5.2 BGS

For the BGS sample, we perform our analysis not just in terms of redshift binning, but also in terms of cuts in absolute magnitude $M_R$. We follow the same procedure as for the LRG sample to model the BGS redshift evolution for different magnitude cuts. The BGS sample has not been analysed with an HOD model in multiple redshift bins in previous papers, so we conduct an additional set of HOD fits of the BGS sample. To construct HOD($z$), we conduct vanilla HOD fits of the BGS Bright sample at ABACUS SUMMIT
<figure><img src="image_10.png" /><figcaption>**Figure 10.** The DESI one percent survey BGS Bright $w_p$ fits across three different absolute magnitude cuts. The solid lines denote the data with jackknife errors. The dashed lines showcase the best-fitting prediction. The orange and green curves have been shifted up by +10 and +20, respectively for visual clarity. Note that we use a non-fiducial redshift range (offset from fiducial bins by 0.05) for the HOD analysis to accommodate the fact that our ABACUS SUMMIT simulation snapshots are at $z = 0.1, 0.2, 0.3, 0.4$. We then use these best-fitting values to derive the redshift evolution of the sample in the fiducial DESI redshift bins.</figcaption></figure>

snapshots $z = 0.1, 0.2, 0.3, 0.4$, spanning our DESI redshift bins. We repeat the fits for three different absolute magnitude cuts, $M_R < -19$, $M_R < -20$, and $M_R < -21$. We only consider the vanilla HOD model for simplicity, and fit the projected 2PCF $w_p$ within $0.1 < r_p < 32 h^{-1} \text{Mpc}$. We assume Gaussian likelihoods and use the diagonal errors computed with 60 jack-knife regions. The HOD best-fits are summarized in Table 2, and the best-fitting $w_p$ predictions are displayed in Fig. 10. Note that we have shifted the $w_p$ along the y-axis for visual clarity. We have also skipped $M_R < -19$ at $z = 0.4$ since that is the same sample as $M_R < -20$ at $z = 0.4$ (Fig. 2). Note that we are only varying the two mass parameters in these HOD fits because they are much better constrained than the other parameters given the noisy measurements.

Fig. 11 visualizes the evolution of the best-fitting HODs for different magnitude-limited BGS samples. We see clear redshift dependency for all different magnitude cuts. The faintest $M_R < -19$

sample displays the strongest redshift evolution. The $M_R < -21$ sample is less sensitive to redshift and has a similar HOD to that of the LRG samples shown in grey. The dependence of the mass parameters (and linear bias) on absolute magnitude cuts is consistent with findings of Prada et al. (2023; see table 4). The satellite occupation parameters are also similar to that of Prada et al. (2023).

We again build an HOD($z$) model by interpolating parameters $\log M_{\text{cut}}$ and $\log M_1$. We compute $\delta_{z-\text{evol}}$ and $\sigma_{\text{stat}}$ for the BGS sample using Dark Emulator, except we repeat these calculations for different magnitude cuts: $M_R < -19$, $M_R < -20$, and $M_R < -21$.

Fig. 12 shows the ratio of $\delta_{z-\text{evol}}$ to $\sigma_{\text{stat}}$ as a function of scale, but we organize the plots differently than for LRGs. The three columns correspond to the three fiducial redshift bins, whereas the two rows correspond to the two summary statistics. We use three different colours to show three different magnitude cuts. The horizontal dashed lines indicate the threshold where the redshift evolution becomes a
<figure><img src="image_11.png" /><figcaption>**Figure 11.** The DESI one-per cent BGS Bright sample HOD best-fit across three different magnitude-limited subsamples (different panels and colours) and four redshift bins (different line styles). The grey solid and grey dashed lines represent the LRG HOD best-fits from Table 1. We see clear redshift evolution in BGS sample. The $M_R < -21$ sample is qualitatively similar to the LRG sample.</figcaption></figure>

statistically significant effect. The widths on the curves again show the spread between different source samples. We have skipped $M_R < -19$ in $0.3 < z < 0.4$ as it is the same sample as $M_R < -20$ at $0.3 < z < 0.4$. We have also used open hatches instead of solid colours for the two bins where the sample is highly incomplete: $M_R < -19$ in

$0.2 < z < 0.3$, and $M_R < -20$ in $0.3 < z < 0.4$. These two choices are not relevant for the final cosmology analysis as they do not represent significant gains over brighter magnitude cuts in the same redshift bins.

We see that the systematics due to redshift evolution are somewhat more significant for BGS than for LRGs. As expected, the significance of redshift evolution decreases with the limiting absolute magnitude. For lensing, redshift evolution remains insignificant across all scales and all magnitude cuts in the fiducial redshift bins. For clustering, redshift evolution is notably more important due to the bias squared scaling. Specifically, redshift evolution can be significant for $w_p$ at the very small scales $r_p < 1 h^{-1}\text{Mpc}$, but the significance drops as we move to brighter samples and higher completeness. The highly incomplete selections shown in open hatches tend to suffer from worse evolution effects. Thus, for analyses that focus on 2-halo scales and larger, redshift evolution should be insignificant given Y1 statistics. For analyses that aim to leverage the 1-halo scales, more care should be given to evolution effect.

A new type of analyses using simulation-based models to derive cosmological constraints have recently gained popularity. Some recent examples include for example Lange et al. (2022, 2023), Yuan et al. (2022), and Zhai et al. (2022). While promising to utilize smaller scales than traditional analytic methods, these analyses are limited to redshift snapshots available in the simulation outputs. We assess the redshift evolution systematics for these types of analyses in Appendix C.

## 6 DISCUSSION AND CONCLUSIONS

DESI Y1 data are expected to deliver high-precision GGL and clustering measurements and strong constraints on cosmological parameters. It is essential to forecast the statistical errors and characterize the impact of different systematics. In this paper, we present the covariance forecasts for DESI Y1 GGL and clustering measurements and assess the significance of lens sample redshift evolution. Our key findings are:

(i) A Gaussian analytical covariance calculation agrees with the galaxy–galaxy lensing error across an ensemble of mock catalogues within around 10 per cent, across the majority of scales and redshifts.

(ii) Survey footprint effects increase the covariance on the largest scales, especially for the ‘cross’ component of galaxy–galaxy lensing. These effects will be important when using B-mode tests on scales comparable to the survey footprint.

(iii) The DESI LRG sample in fiducial bins of $0.4 < z < 0.6$ and $0.6 < z < 0.8$ shows very mild evolution with insignificant effect on the lensing and clustering measurements.

(iv) The evolution in the BGS sample is dependent on the absolute magnitude cut, with brighter cuts resulting in less evolution and vice versa. Galaxy lensing should not be sensitive to redshift evolution across all scales, magnitude cuts, and redshift bins in our fiducial setup. Galaxy clustering is more sensitive, where evolution effects can be significant on 1-halo scales at $r_p < 1 h^{-1}\text{Mpc}$.

This paper should be a reference for making analysis choices with DESI Y1 samples, especially in selecting redshift bins and magnitude cuts. For LRGs, we find the fiducial redshift bins to be more than sufficient in controlling the redshift evolution effects. For BGS, we find galaxy–galaxy lensing to be safe given fiducial choices. Galaxy clustering is also safe if the analysis does not utilize the smallest scales. For small-scale analyses that want to maximize sample size, additional redshift bins might be needed. Alternatively, one can adopt different magnitude cuts for different redshift ranges. Specifically,
<figure><img src="image_12.png" /><figcaption>**Figure 12.** The ratio of the redshift evolution effect in the fiducial redshift bins to the expected statistical error for BGS sample with different magnitude cuts. The three columns correspond to the three fiducial BGS lens bins. The top row shows the results for $w_p$, whereas the bottom row shows the results for $\Delta\Sigma$. The three colours show the three different absolute magnitude cuts, while the width of the curves represents the spread between different photometric samples. We have skipped $M_R < -19$ in $0.3 < z < 0.4$ as it is the same sample as $M_R < -20$ at $0.3 < z < 0.4$. We have also used open hatches instead of solid colours for the two bins that are highly incomplete (dashed lines in Fig. 2). We see that for galaxy lensing, redshift evolution remains insignificant across all scales and magnitude cuts in the fiducial redshift bins. For galaxy clustering, redshift evolution is insignificant at scales greater than $r_p > 1~h^{-1}\text{Mpc}$. The redshift evolution is less important for brighter samples and tend to be worst for highly incomplete samples.</figcaption></figure>

lower redshift bins tend to see less evolution at fixed luminosity and thus can adopt a fainter magnitude cut. Fig. 12 acts as a reference plot for making these magnitude choices.

There are several caveats with this analysis. For simplicity, we do not characterize or report the uncertainties on the HOD fits. This means that the systematic effects we report in Figs 9 and 12 are of their mean expected amplitudes, not their worst-case amplitudes. However, because the best-fitting HODs produce continuous and expected behaviours as a function of redshift, we expect the uncertainties of the HOD fits to be small and that our conclusions are robust. Similarly, the HODs themselves are calibrated in finite redshift bins assuming a model at a fixed redshift. We designed these bins to be centred on the simulation output in the BGS analysis to reduce this bias. Nevertheless, a more rigorous examination of redshift evolution is needed beyond Y1.

Additionally, we have modelled the redshift dependence with an interpolation scheme between a small number of snapshots. While that is sufficient for identifying the overall trends, it potentially introduces spurious evolution effects due to the low number of interpolated points and errors in the best-fit. We suggest a repeat of this analysis using simulated lightcones and properly implemented HOD($z$) models. Finally, we directly applied the best-fitting HODs from ABACUS SUMMIT to Dark Emulator, which relies on a different set of simulations and different halo find-

ers. This leads to a small systematic bias that is ignored in this paper.

## ACKNOWLEDGEMENTS

This work was supported by grant DE-SC0013718 and under DE-AC02-76SF00515 to SLAC National Accelerator Laboratory, and by the Kavli Institute for Particle Astrophysics and Cosmology. This work was performed in part at the Aspen Center for Physics, which was supported by National Science Foundation grant PHY-2210452.

This material is based upon work supported by the U.S. Department of Energy (DOE), Office of Science, Office of High-Energy Physics, under Contract No. DE-AC02-05CH11231 and DE-SC0010107, and by the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility under the same contract. Additional support for DESI was provided by the U.S. National Science Foundation (NSF), Division of Astronomical Sciences under Contract No. AST-0950945 to the NSF's National Optical-Infrared Astronomy Research Laboratory; the Science and Technology Facilities Council of the United Kingdom; the Gordon and Betty Moore Foundation; the Heising-Simons Foundation; the French Alternative Energies and Atomic Energy Commission (CEA); the National Council of Science and Technology of Mexico (CONACYT); the Ministry of Science and Innovation of Spain (MICINN),
and by the DESI Member Institutions: https://www.desi.lbl.gov/col
laborating-institutions. Any opinions, findings, and conclusions or
recommendations expressed in this material are those of the author(s)
and do not necessarily reflect the views of the U.S. National Science
Foundation, the U.S. Department of Energy, or any of the listed
funding agencies.

The authors are honoured to be permitted to conduct scientific
research on Iolkam Du'ag (Kitt Peak), a mountain with particular
significance to the Tohono O'odham Nation.

This work made use of ASTROPY:<fnref n="3" /> a community-developed core
Python package and an ecosystem of tools and resources for astron-
omy (Astropy Collaboration 2013, 2018, 2022).

DATA AVAILABILITY

The simulation data are available at https://abacussummit.rea
dthedocs.io/en/latest/. The ABACUSHOD code package is pub-
licly available as a part of the ABACUSUTILS package at
https://github.com/abacusorg/abacusutils. Example usage can be
found at https://abacusutils.readthedocs.io/en/latest/hod.html. All
mock products will be made available at https://data.desi.lbl.gov.

All data points shown in the published graph will be also available
at https://doi.org/10.5281/zenodo.10724372.

REFERENCES

Abareshi B. et al., 2022, AJ, 164, 207
Abbott T. M. C. et al., 2022, Phys. Rev. D, 105, 023520
Ahumada R. et al., 2020, ApJS, 249, 3
Aihara H. et al., 2018, PASJ, 70, S8
Alam S., Peacock J. A., Kraljic K., Ross A. J., Comparat J., 2020, MNRAS, 497, 581
Amon A. et al., 2022, Phys. Rev. D, 105, 023514
Astropy Collaboration, 2013, A&A, 558, A33
Astropy Collaboration, 2018, AJ, 156, 123
Astropy Collaboration, 2022, ApJ, 935, 167
Blake C. et al., 2016, MNRAS, 456, 2806
Blake C. et al., 2020, A&A, 642, A158
Brainerd T. G., Blandford R. D., Smail I., 1996, ApJ, 466, 623
Buchs R. et al., 2019, MNRAS, 489, 820
Cole S. et al., 2005, MNRAS, 362, 505
Coupon J. et al., 2015, MNRAS, 449, 1352
DESI Collaboration et al., 2016, preprint (arXiv:1611.00036)
DESI Collaboration, 2023a, AJ, 167, 62
DESI Collaboration, 2023b, AJ, 168, 58
Dalal R. et al., 2023, Phys. Rev. D, 108, 123519
Dark Energy Survey Collaboration, 2016, MNRAS, 460, 1270
Dawson K. S. et al., 2013, AJ, 145, 10
Dawson K. S. et al., 2016, AJ, 151, 44
DeRose J. et al., 2019a, preprint (arXiv:1901.02401)
DeRose J. et al., 2019b, ApJ, 875, 69
Drlica-Wagner A. et al., 2018, ApJS, 235, 33
Dvornik A. et al., 2018, MNRAS, 479, 1240
Dvornik A. et al., 2023, A&A, 675, A189
Euclid Collaboration, 2020, A&A, 642, A191
Friedrich O. et al., 2021, MNRAS, 508, 3125
Garrison L. H., Eisenstein D. J., Pinto P. A., 2019, MNRAS, 485, 3370
Garrison L. H., Eisenstein D. J., Ferrer D., Maksimova N. A., Pinto P. A., 2021, MNRAS, 508, 575
Gatti M. et al., 2021, MNRAS, 504, 4312
Gatti M. et al., 2022, MNRAS, 510, 1223
Giblin B. et al., 2021, A&A, 645, A105

Guo H. et al., 2015, MNRAS, 453, 4368
Hadzhiyska B., Eisenstein D., Bose S., Garrison L. H., Maksimova N., 2022, MNRAS, 509, 501

Hahn C. et al., 2023, AJ, 165, 253

Hamilton A. J. S., Tegmark M., 2002, MNRAS, 330, 506

Heymans C. et al., 2021, A&A, 646, A140

Hikage C. et al., 2019, PASJ, 71, 43

Hildebrandt H. et al., 2021, A&A, 647, A124

Ilbert O. et al., 2009, ApJ, 690, 1236

Ivezic Z. et al., 2008, Serb. Astron. J., 176, 1

Joachimi B. et al., 2021, A&A, 646, A129

Kilbinger M., 2015, Rep. Prog. Phys., 78, 086901

Kilbinger M., Schneider P., 2004, A&A, 413, 465

Kobayashi Y., Nishimichi T., Takada M., Miyatake H., 2022, Phys. Rev. D, 105, 083517

Krause E., Eifler T., 2017, MNRAS, 470, 2100

Kuijken K. et al., 2019, A&A, 625, A2

Landy S. D., Szalay A. S., 1993, ApJ, 412, 64

Lange J. U., Yang X., Guo H., Luo W., van den Bosch F. C., 2019, MNRAS, 488, 5771

Lange J. U., Hearin A. P., Leauthaud A., van den Bosch F. C., Guo H., DeRose J., 2022, MNRAS, 509, 1779

Lange J. U., Hearin A. P., Leauthaud A., van den Bosch F. C., Xhakaj E., Guo H., Wechsler R. H., DeRose J., 2023, MNRAS, 520, 5373

Lange J. U. et al., 2024, preprint (arXiv:2404.09397)

Laureijs R. et al., 2011, preprint (arXiv:1110.3193)

Leauthaud A. et al., 2017, MNRAS, 467, 3024

Leauthaud A. et al., 2022, MNRAS, 510, 6150

Levi M. et al., 2013, preprint (arXiv:1308.0847)

Maksimova N. A., Garrison L. H., Eisenstein D. J., Hadzhiyska B., Bose S., Satterthwaite T. P., 2021, MNRAS, 508, 4017

Mandelbaum R., 2018, ARA&A, 56, 393

Mandelbaum R., Slosar A., Baldauf T., Seljak U., Hirata C. M., Nakajima R., Reyes R., Smith R. E., 2013, MNRAS, 432, 1544

Mandelbaum R. et al., 2018a, PASJ, 70, S25

Mandelbaum R. et al., 2018b, MNRAS, 481, 3170

Marian L., Smith R. E., Angulo R. E., 2015, MNRAS, 451, 1418

Miralda-Escude J., 1991, ApJ, 380, 1

Miyazaki S. et al., 2018, PASJ, 70, S1

Myles J. et al., 2021, MNRAS, 505, 4249

Nishimichi T. et al., 2019, ApJ, 884, 29

Padmanabhan N. et al., 2007, MNRAS, 378, 852

Pandey S. et al., 2022, Phys. Rev. D, 106, 043520

Pandey S., Sánchez C., Jain B., The LSST Dark Energy Science Collaboration, 2023, preprint (arXiv:2310.01315)

Park Y. et al., 2016, Phys. Rev. D, 94, 063533

Percival W. J. et al., 2001, MNRAS, 327, 1297

Planck Collaboration, 2020, A&A, 641, A6

Prada F. et al., 2023, preprint (arXiv:2306.06315)

Ruiz-Macias O. et al., 2021, MNRAS, 502, 4328

Sato M., Takada M., Hamana T., Matsubara T., 2011, ApJ, 734, 76

Schrabback T. et al., 2010, A&A, 516, A63

Shirasaki M., Takada M., 2018, MNRAS, 478, 4277

Shirasaki M., Hamana T., Takada M., Takahashi R., Miyatake H., 2019, MNRAS, 486, 52

Singh S., Mandelbaum R., Seljak U., Slosar A., Vazquez Gonzalez J., 2017, MNRAS, 471, 3827

Tegmark M. et al., 2006, Phys. Rev. D, 74, 123507

Wilson G., Kaiser N., Luppino G. A., Cowie L. L., 2001, ApJ, 555, 572

Wright A. H., Hildebrandt H., van den Busch J. L., Heymans C., 2020, A&A, 637, A100

Xavier H. S., Abdalla F. B., Joachimi B., 2016, MNRAS, 459, 3693

Yuan S., Eisenstein D. J., Leauthaud A., 2020, MNRAS, 493, 5551

Yuan S., Garrison L. H., Hadzhiyska B., Bose S., Eisenstein D. J., 2021, MNRAS, 510, 3301

Yuan S., Garrison L. H., Eisenstein D. J., Wechsler R. H., 2022, MNRAS, 515, 871

Yuan S. et al., 2024, MNRAS, 530, 947

<footnote n="3">³http://www.astropy.org</footnote>
Zehavi I. et al., 2011, ApJ, 736, 59

Zhai Z. et al., 2017, ApJ, 848, 76

Zhai Z. et al., 2022, ApJ, 948, 99

Zheng Z. et al., 2005, ApJ, 633, 791

Zheng Z., Coil A. L., Zehavi I., 2007, ApJ, 667, 760

Zhou R. et al., 2020, Res. Notes Am. Astron. Soc., 4, 181

# APPENDIX A: GAUSSIAN COVARIANCE OF PROJECTED CORRELATIONS

In this section, we detail our determination of the analytical covariance for the projected correlation functions, $\Delta\Sigma(r_p)$ and $w_p(r_p)$.

## A1 Relation of projected correlations to power spectra

The projected mass density around lens galaxies at projected separation $r_p$ is determined by the 3D galaxy-mass cross-correlation function $\xi_{gm}(\mathbf{r})$ as,

$$ \Sigma(\mathbf{r}_p) = \bar{\rho}_m \int_{-\infty}^{\infty} d\mathbf{r}_{\pi} \left[ 1 + \xi_{gm}(\mathbf{r}_p, r_{\pi}) \right], \quad (A1) $$

where $r_\pi$ is the line-of-sight separation and $\bar{\rho}_m$ is the mean cosmic matter density. The surface mass density contrast is then found by averaging the projected density over directions in the plane of the sky defined by $\phi$, such that,

$$ \begin{aligned} \Delta\Sigma(r_p) &= \bar{\Sigma}(\langle r_p \rangle) - \Sigma(r_p) \\ &= \frac{2}{r_p^2} \int_0^{r_p} dr'_p r'_p \Sigma(r'_p) - \Sigma(r_p) \\ &= \bar{\rho}_m \int_{-\infty}^{\infty} dr_{\pi} \left[ \frac{2}{r_p^2} \int_0^{r_p} dr'_p r'_p \int_0^{2\pi} \frac{d\phi'}{2\pi} \xi_{gm}(r'_p, r_{\pi}) \right. \\ & \qquad \left. - \int_0^{2\pi} \frac{d\phi}{2\pi} \xi_{gm}(r_p, r_{\pi}) \right], \end{aligned} \quad (A2) $$

after substituting in equation (A1). We can express this in terms of the galaxy-mass cross-power spectrum by using $\xi_{gm}(\mathbf{r}) = \int \frac{d^3\mathbf{k}}{(2\pi)^3} P_{gm}(\mathbf{k}) e^{-i\mathbf{k}\cdot\mathbf{r}}$, obtaining,

$$ \begin{aligned} \Delta\Sigma(r_p) &= \bar{\rho}_m \int \frac{d^3\mathbf{k}}{(2\pi)^3} P_{gm}(\mathbf{k}) \left[ \int_{-\infty}^{\infty} d\mathbf{r}_{\pi} e^{-i\mathbf{k}\cdot\mathbf{r}_{\pi}} \right] \\ & \quad \times \left[ \frac{2}{r_p^2} \int_0^{r_p} dr'_p r'_p \int_0^{2\pi} \frac{d\phi'}{2\pi} e^{-i\mathbf{k}\cdot\mathbf{r}'_p} \right. \\ & \quad \left. - \int_0^{2\pi} \frac{d\phi}{2\pi} e^{-i\mathbf{k}\cdot\mathbf{r}_p} \right] \\ &= \bar{\rho}_m \int \frac{d^3\mathbf{k}}{(2\pi)^3} P_{gm}(\mathbf{k}) \\ & \quad \times \left[ L_{\parallel} \delta_D(k_{\parallel}) \left[ \frac{2}{r_p^2} \int_0^{r_p} dr'_p r'_p J_0(k_{\perp}r'_p) - J_0(k_{\perp}r_p) \right] \right] \\ &= \bar{\rho}_m \int \frac{d^2\mathbf{k}_{\perp}}{(2\pi)^2} P_{gm}(\mathbf{k}_{\perp}) J_2(k_{\perp}r_p), \end{aligned} \quad (A3) $$

where $\mathbf{k}_{\perp}$ is the 2D Fourier wave vector in the plane of the sky, $k_{\parallel}$ is the wavenumber in the line-of-sight direction, $L_{\parallel}$ is the thickness of the redshift slice, and $J_n$ indicates a Bessel function of the first kind. The projected galaxy correlation function at projected separation $r_p$ is defined by,

$$ w_p(\mathbf{r}_p) = \int_{-r_{\pi, \max}}^{r_{\pi, \max}} d\mathbf{r}_{\pi} \xi_{gg}(\mathbf{r}_p, r_{\pi}), \quad (A4) $$

where $r_{\pi, \max}$ is the upper limit used in the measurement, and $\xi_{gg}(\mathbf{r})$ is the 3D galaxy correlation function. We can express this in terms of the galaxy autopower spectrum by using $\xi_{gg}(\mathbf{r}) = \int \frac{d^3\mathbf{k}}{(2\pi)^3} P_{gg}(\mathbf{k}) e^{-i\mathbf{k}\cdot\mathbf{r}}$, obtaining,

$$ w_p(\mathbf{r}_p) = \int \frac{d^3\mathbf{k}}{(2\pi)^3} P_{gg}(\mathbf{k}) e^{-i\mathbf{k}_{\perp}\cdot\mathbf{r}_p} \int_{-r_{\pi, \max}}^{r_{\pi, \max}} d\mathbf{r}_{\pi} e^{-i\mathbf{k}_{\parallel}\cdot\mathbf{r}_{\pi}}. \quad (A5) $$

Averaging the measurement over directions $\phi$ we find,

$$ w_p(r_p) = 2r_{\pi, \max} \int \frac{d^3\mathbf{k}}{(2\pi)^3} P_{gg}(\mathbf{k}) J_0(k_{\parallel}r_p) J_0(k_{\perp}r_{\pi, \max}). \quad (A6) $$

where $J_0$ is a spherical Bessel function. In the Limber approximation we suppose that only tangential modes contribute to the projected clustering, which allows us to simplify equation (A6) to the form,

$$ \begin{aligned} w_p^{\text{Lim}}(r_p) &= \int \frac{d^2 \mathbf{k}_{\perp}}{(2\pi)^2} P_{gg}(\mathbf{k}_{\perp}) J_0(k_{\perp} r_p) \\ & \quad \times 2r_{\pi, \max} \int_{-\infty}^{\infty} \frac{d k_{\parallel}}{2\pi} J_0(k_{\parallel} r_{\pi, \max}) \\ &= \int \frac{d^2 \mathbf{k}_{\perp}}{(2\pi)^2} P_{gg}(\mathbf{k}_{\perp}) J_0(k_{\perp} r_p). \end{aligned} \quad (A7) $$

## A2 Covariance of projected correlations

We may deduce the covariance of $\Delta\Sigma(r_p)$ using equation (A3) (for similar treatments, see Singh et al. 2017; Dvornik et al. 2018; Shirasaki & Takada 2018; Blake et al. 2020),

$$ \begin{aligned} \text{Cov}[\Delta\Sigma(r_p), \Delta\Sigma(r'_p)] &= \bar{\rho}_m^2 \int \frac{d^2\mathbf{k}_{\perp}}{(2\pi)^2} \int \frac{d^2\mathbf{k}'_{\perp}}{(2\pi)^2} \text{Cov}[P_{gm}(\mathbf{k}_{\perp}), P_{gm}(\mathbf{k}'_{\perp})] \\ & \quad \times J_2(k_{\perp}r_p) J_2(k'_{\perp}r'_p) \\ &= \frac{\bar{\rho}_m^2}{A_s} \int \frac{d^2\mathbf{k}_{\perp}}{(2\pi)^2} \text{Cov}[P_{gm}(\mathbf{k}_{\perp}), P_{gm}(\mathbf{k}'_{\perp})] J_2(k_{\perp}r_p) J_2(k'_{\perp}r'_p), \end{aligned} \quad (A8) $$

where $A_s = \Omega_s \chi_{\text{eff}}^2$ is the projected survey area at the effective lens distance $\chi_{\text{eff}}$.

We now use the relations between the spatial and angular power spectra for a narrow redshift slice,

$$ C_{g\kappa}^{\text{Nar}}(\ell) = \bar{\rho}_m \frac{\overline{\sum_c^{-1}(\chi_{\text{eff}})}}{\chi_{\text{eff}}^2} P_{gm} \left( \frac{\ell}{\chi_{\text{eff}}} \right), \quad (A9) $$

where $\ell$ is the angular wavenumber and the general expression for the covariances between the angular power spectra of different fields (labelled A, B, C, D) and samples (labelled i, j, k, l) (Krause & Eifler 2017),

$$ \begin{aligned} \text{Cov}[C_{AB}^{ij}, C_{CD}^{kl}] = & [C_{AC}^{ik} + \delta_{ik}^{K} \delta_{AC}^{K} N_{A}^{i}] [C_{BD}^{j\ell} + \delta_{jl}^{K} \delta_{BD}^{K} N_{B}^{j}] \\ & + [C_{AD}^{il} + \delta_{il}^{K} \delta_{AD}^{K} N_{A}^{i}] [C_{BC}^{jk} + \delta_{jk}^{K} \delta_{BC}^{K} N_{B}^{j}] \end{aligned} \quad (A10) $$
where N indicates the noise autopower spectrum of a field. Applying equations (A9) and (A10),

$$
\begin{align}
\mathrm{Cov} [P_{\mathrm{gm}}(\mathbf{k}_{\perp}), P_{\mathrm{gm}}(\mathbf{k}_{\perp})] & \\
&= \frac{\chi_{\mathrm{eff}}^4}{\bar{\rho}_m^2 \left[ \frac{\Sigma_c^{-1}(\chi_{\mathrm{eff}})}{\bar{\rho}_m^2} \right]^2} \mathrm{Cov} [C_{g\kappa}^{\mathrm{Nar}}(\mathbf{k}_{\perp}\chi_{\mathrm{eff}}), C_{g\kappa}^{\mathrm{Nar}}(\mathbf{k}_{\perp}\chi_{\mathrm{eff}})] \nonumber \\
&= P_{\mathrm{gm}}^2(k_{\perp}) + \frac{\chi_{\mathrm{eff}}^4}{\bar{\rho}_m^2 \left[ \frac{\Sigma_c^{-1}(\chi_{\mathrm{eff}})}{\bar{\rho}_m^2} \right]^2} [C_{\kappa\kappa}(k_{\perp}\chi_{\mathrm{eff}}) + N_{\kappa}] \nonumber \\
&\quad \times [C_{\mathrm{gg}}^{\mathrm{Nar}}(k_{\perp}\chi_{\mathrm{eff}}) + N_{\mathrm{g}}]. \tag{A11}
\end{align}
$$

The expression for the covariance is then,

$$
\begin{align}
\mathrm{Cov} [\Delta\Sigma(r_p), \Delta\Sigma(r'_p)] & \nonumber \\
&= \frac{1}{A_s} \int \frac{dk_{\perp} k_{\perp}}{2\pi} J_2(k_{\perp} r_p) J_2(k_{\perp} r'_p) \nonumber \\
&\quad \times \left\{ \bar{\rho}_m^2 P_{\mathrm{gm}}^2(k_{\perp}) + \frac{\chi_{\mathrm{eff}}^4}{[\Sigma_c^{-1}(\chi_{\mathrm{eff}})]^2} [C_{\kappa\kappa}(k_{\perp}\chi_{\mathrm{eff}}) + N_{\kappa}] \right. \nonumber \\
&\quad \left. \times [C_{\mathrm{gg}}^{\mathrm{Nar}}(k_{\perp}\chi_{\mathrm{eff}}) + N_{\mathrm{g}}] \right\} \tag{A12}
\end{align}
$$

Similarly for the ‘cross’ component,

$$
\begin{align}
\mathrm{Cov} [\Delta\Sigma_{\times}(r_p), \Delta\Sigma_{\times}(r'_p)] & \nonumber \\
&= \frac{1}{A_s} \frac{\chi_{\mathrm{eff}}^4}{[\Sigma_c^{-1}(\chi_{\mathrm{eff}})]^2} \int \frac{dk_{\perp} k_{\perp}}{2\pi} J_2(k_{\perp} r_p) J_2(k_{\perp} r'_p) N_{\kappa} \nonumber \\
&\quad \times [C_{\mathrm{gg}}^{\mathrm{Nar}}(k_{\perp}\chi_{\mathrm{eff}}) + N_{\mathrm{g}}]. \tag{A13}
\end{align}
$$

The covariance of $w_p(r_p)$ follows from equation (A6),

$$
\begin{align}
\mathrm{Cov} [w_p(r_p), w_p(r'_p)] & \nonumber \\
&= 4r_{\pi, \max}^2 \int \frac{d^3 \mathbf{k}}{(2\pi)^3} \int \frac{d^3 \mathbf{k}'}{(2\pi)^3} \nonumber \\
&\quad \times \mathrm{Cov} [P_{\mathrm{gg}}(\mathbf{k}), P_{\mathrm{gg}}(\mathbf{k}')] J_0(k_{\perp} r_p) J_0(k'_{\perp} r'_p) \nonumber \\
&\quad \times j_0(k_{\parallel} r_{\pi, \max}) j_0(k'_{\parallel} r_{\pi, \max}). \tag{A14}
\end{align}
$$

Using $\mathrm{Cov}[P_{\mathrm{gg}}(\mathbf{k}), P_{\mathrm{gg}}(\mathbf{k}')] = 2 \left[P_{\mathrm{gg}}(\mathbf{k}) + \frac{1}{n_g}\right]^2 \tilde{\delta}_D(\mathbf{k} - \mathbf{k}')$, where the factor of 2 accounts for the fact that the modes $-\mathbf{k}$ and $\mathbf{k}$ are not independent, this becomes,

$$
\mathrm{Cov} [w_p(r_p), w_p(r'_p)] = \frac{8r_{\pi, \max}^2}{V_s} \int \frac{d^3 \mathbf{k}}{(2\pi)^3} \left[ P_{\mathrm{gg}}(\mathbf{k}) + \frac{1}{n_g} \right]^2 \times J_0(k_{\perp} r_p) J_0(k_{\perp} r'_p) J_0^2(k_{\parallel} r_{\pi, \max}), \quad (A15)
$$

where $V_s$ is the survey volume. In the Limber approximation where only tangential modes contribute to the covariance,

$$
\begin{align}
\mathrm{Cov}^{\mathrm{Lim}} [w_p(r_p), w_p(r'_p)] & \nonumber \\
&= \frac{8r_{\pi, \max}^2}{V_s} \int \frac{d^2 \mathbf{k}_{\perp}}{(2\pi)^2} \left[ P_{\mathrm{gg}}(k_{\perp}) + \frac{1}{n_g} \right]^2 J_0(k_{\perp} r_p) J_0(k_{\perp} r'_p) \nonumber \\
&\quad \times \int_{-\infty}^{\infty} \frac{dk_{\parallel}}{2\pi} j_0^2(k_{\parallel} r_{\pi, \max}) \nonumber \\
&= \frac{4r_{\pi, \max}}{V_s} \int \frac{dk_{\perp} k_{\perp}}{2\pi} \left[ P_{\mathrm{gg}}(k_{\perp}) + \frac{1}{n_g} \right]^2 J_0(k_{\perp} r_p) J_0(k_{\perp} r'_p) \nonumber \\
&= \frac{4r_{\pi, \max}}{V_s} \int \frac{dk_{\perp} k_{\perp}}{2\pi} \left[ P_{\mathrm{gg}}^2(k_{\perp}) + \frac{2P_{\mathrm{gg}}(k_{\perp})}{n_g} \right] J_0(k_{\perp} R) J_0(k_{\perp} r'_p) \nonumber \\
&\quad + \frac{2r_{\pi, \max}}{\pi r_p n_g^2 V_s} \delta_D(r_p - r'_p). \tag{A16}
\end{align}
$$

Equation (A16) is identical to equation (31) in Marian, Smith & Angulo (2015), neglecting the term '$w_{gg}^{\mathrm{Lim}}(r_p)$' in their expression. This term arises if we use a Poisson rather than a Gaussian model for the galaxy statistics, and can be neglected on large scales.

### A3 Cross-covariances involving projected correlations

The cross-covariance between $\Delta\Sigma(r_p)$ and $w_p(r_p)$ follows from equations (A3) and (A7),

$$
\begin{align}
\mathrm{Cov} [\Delta\Sigma(r_p), w_p^{\mathrm{Lim}}(r'_p)] & \nonumber \\
&= \frac{\bar{\rho}_m}{A_s} \int \frac{d^2 \mathbf{k}_{\perp}}{(2\pi)^2} \mathrm{Cov} [P_{\mathrm{gm}}(\mathbf{k}_{\perp}), P_{\mathrm{gg}}(\mathbf{k}_{\perp})] J_2(k_{\perp} r_p) J_0(k_{\perp} r'_p). \tag{A17}
\end{align}
$$

We now use the relation for narrow redshift slice,

$$
C_{gg}^{\mathrm{Nar}}(\ell) = \frac{1}{\chi_{\mathrm{eff}}^2 L_{\parallel}} P_{\mathrm{gg}} \left( \frac{\ell}{\chi_{\mathrm{eff}}} \right). \quad (A18)
$$

Applying equations (A18), (A9), and (A10),

$$
\begin{align}
\mathrm{Cov} [P_{\mathrm{gm}}(\mathbf{k}_{\perp}), P_{\mathrm{gg}}(\mathbf{k}_{\perp})] & \nonumber \\
&= \frac{\chi_{\mathrm{eff}}^4 L_{\parallel}}{\bar{\rho}_m \Sigma_c^{-1}(\chi_{\mathrm{eff}})} \mathrm{Cov} [C_{g\kappa}^{\mathrm{Nar}}(k_{\perp}\chi_{\mathrm{eff}}), C_{gg}^{\mathrm{Nar}}(k_{\perp}\chi_{\mathrm{eff}})] \nonumber \\
&= \frac{\chi_{\mathrm{eff}}^4 L_{\parallel}}{\bar{\rho}_m \Sigma_c^{-1}(\chi_{\mathrm{eff}})} \times 2 [C_{gg}^{\mathrm{Nar}}(k_{\perp}\chi_{\mathrm{eff}}) + N_g] C_{g\kappa}^{\mathrm{Nar}}(k_{\perp}\chi_{\mathrm{eff}}) \nonumber \\
&= 2 \left[ P_{\mathrm{gg}}(k_{\perp}) + \frac{1}{n_g} \right] P_{\mathrm{gm}}(k_{\perp}). \tag{A19}
\end{align}
$$

Hence, the cross-covariance expression is,

$$
\begin{align}
\mathrm{Cov} [\Delta\Sigma(r_p), w_p^{\mathrm{Lim}}(r'_p)] & \nonumber \\
&= \frac{2\bar{\rho}_m}{A_s} \int \frac{dk_{\perp} k_{\perp}}{2\pi} \left[ P_{\mathrm{gg}}(k_{\perp}) + \frac{1}{n_g} \right] \nonumber \\
&\quad \times P_{\mathrm{gm}}(k_{\perp}) J_2(k_{\perp} r_p) J_0(k_{\perp} r'_p). \tag{A20}
\end{align}
$$

## APPENDIX B: TESTING FOOTPRINT EFFECTS ON THE COVARIANCE USING LOGNORMAL MOCKS

In this section, we use lognormal mock catalogues to calibrate the impact of the survey footprint on the covariance of galaxy–galaxy lensing statistics. For some previous investigations of this issue, we refer the reader to Kilbinger & Schneider (2004), Sato et al. (2011),
<figure><img src="image_13.png" /><figcaption>**Figure B1.** The impact of survey geometry on the error in the tangential and cross-components of the average shear ($\gamma_t$ and $\gamma_\times$). The upper panel shows the fraction of observed source–lens pairs as a function of angular separation $\theta$, compared to the theoretical pair count that would be observed in the absence of survey boundary effects. Results are shown for KiDS, DES, and HSC DESI overlap footprints, where the band indicates the standard deviation across different pairs of lens and source samples. The angular scale corresponding to losing 50 per cent of pairs is indicated by the vertical dashed line. The middle and lower panel display the ratio of the standard deviation of the average shear measurements across the realizations, to the analytical covariance error. The bands again indicate the standard deviation of results across the different lens and source samples.</figcaption></figure>

Shirasaki et al. (2019), Friedrich et al. (2021), and Joachimi et al. (2021). Footprint effects cannot be readily included in analytical treatments of the covariance, even for the Gaussian term. However, numerical investigations allow the significance of the effects to be explored.

Lognormal realizations are a useful resource for such investigations because a large ensemble of simulations can be cheaply generated (we use 1000 realizations for our analysis), permitting the differences between the covariances for different configurations to be accurately quantified. We generated these mocks using the Full-sky Lognormal Astro-fields Simulation Kit (FLASK; Xavier, Abdalla & Joachimi 2016), adopting the same fiducial cosmology as the Buzzard N-body simulations and a lens redshift distribution and bias evolution matching the DESI configuration. We simulated a weak lensing source distribution with a density of 10 arcmin⁻², sampled from a generic source redshift distribution for $z < 1.85$ which we subsampled into the tomographic bins of the three weak lensing surveys. (We could match the complete source density of

<figure><img src="image_14.png" /><figcaption>**Figure B2.** The distribution of $\chi^2$ values obtained when performing a 'B-mode test' (that $\gamma_\times = 0$) using the Buzzard and FLASK simulations including the KiDS, DES, and HSC overlap footprints. The open and solid histograms correspond, respectively to using the original analytical covariance, and the covariance corrected for survey footprint effects using the FLASK lognormal mocks. The vertical line corresponds to the number of degrees of freedom in each case.</figcaption></figure>

KiDS and DES, and we subsampled the HSC source density by a factor of 2.) We also sampled the mocks to match the angular completeness of each of the DESI and weak lensing survey overlap regions illustrated in Fig. 4.

We measured the tangential and cross-components of the galaxy–galaxy lensing signal for each realization, $\gamma_t$ and $\gamma_\times$, and compared the result to the analytical covariance determination. For the purposes of this comparison, the analytical covariances were computed using the effective convergence and galaxy model power spectra output by the FLASK code, at the corresponding angular resolution, which capture the lognormal effects. These power spectra differ from the models used to determine the analytical covariance of the Buzzard simulations, which is why the numerical FLASK covariance cannot be used directly as the data covariance.

The results of this galaxy–galaxy lensing error comparison are displayed in Fig. B1, for the overlap footprints between DESI and KiDS, DES, and HSC. The upper panel shows the observed source–lens pairs as a function of angular separation, as a fraction of the theoretical pair count that would be observed in the absence of survey boundary effects. The angular scale where 50 per cent of pairs are lost is indicated by the vertical dashed line, which occurs at 3.39°, 4.29°, and 1.72° for the KiDS, DES, and HSC overlaps, respectively. These relative scales agree well with the survey geometries displayed in Fig. 4. The remaining two panels of Fig. B1 display the ratio of the standard deviation of the $\gamma_t$ and $\gamma_\times$ measurements across the realizations, to the prediction of the analytical covariance. The width of the bands indicate the variation across the different combinations of lens and source tomographic bins, which are combined in this plot for clarity. The two covariance determinations agree well on small scales, but we find that the analytical calculation underestimates the survey covariance on large scales, with this effect being most significant for $\gamma_\times$. We can interpret this finding by considering the 'mixing' of E-mode and B-mode convergence power in the
contributions to the Gaussian covariance including a survey footprint: the contribution of the (larger) E-mode power to $\gamma_\times$ exceeds that of the (smaller) B-mode power to $\gamma_t$. We find that at the angular scales noted above where 50 per cent of pairs are lost, the $\gamma_t$ error increases by around 10 per cent or less, and the $\gamma_\times$ error increases by around 30 per cent.

In our calculation of the analytical covariances for these survey footprints, we include these effects by multiplying the diagonal errors by this correction factor, whilst maintaining the same cross-correlation between the off-diagonal elements as predicted by the analytical calculation. (The corrections for the off-diagonal covariances are noisy, as the denominators are small.) For the $\Delta\Sigma$ covariances, we cannot compute the correction directly as a function of $R$, since the FLASK simulations are purely angular. Hence in this case we interpolate the correction at the corresponding angular separation of each lens redshift bin, $\theta = R/\chi(z_l)$.

We tested the efficacy of these corrections by applying them to a 'B-mode test' of the Buzzard and FLASK mocks. This test verifies that the $\gamma_\times$ measurements are consistent with zero, using the assumed covariance combining all source and lens bins. We consider analyses adopting both the original analytical covariance, and the corrected covariance described above. The distribution of $\chi^2$ values across the different realizations, corresponding to these cases, is shown in Fig. B2. We assume the maximum fitting scales indicated in Fig. B1, at which 50 per cent of source–lens pairs are excluded by the survey boundary. We find that if the original (uncorrected) covariance is used, the $\chi^2$ values are overestimated compared to the

number of degrees of freedom, owing to the underestimate of the covariance, such that the $\gamma_\times = 0$ hypothesis is erroneously rejected for many realizations. However, following the correction to the covariance based on the FLASK lognormal mocks, an appropriate $\chi^2$ distribution is recovered again.

## APPENDIX C: SIMULATION-BASED ANALYSES

In this section, we characterize the significance of redshift evolution effects for simulation-based analyses. For DESI, two simulation suites are considered for such analyses: ABACUSSUMMIT (Maksimova et al. 2021) and AEMULUS (DeRose et al. 2019b). The available relevant redshift snapshots are displayed at the bottom of Fig. 1 in short vertical lines. For the BGS sample, the ABACUSSUMMIT snapshots are unfortunately lying at the edge of the fiducial redshift bins. Using these snapshots for the fiducial bins is expected to result in large systematic biases. Thus, it is necessary for simulation-based analyses to use custom non-fiducial redshift bins for the BGS sample.

Given that the relevant ABACUSSUMMIT snapshots are at $z = 0.2, 0.3, 0.4$, the simplest redshift binning scheme for ABACUSSUMMIT-based analyses is $0.15 < z < 0.25$, $0.25 < z < 0.35$, and $0.35 < z < 0.45$. Fig. C1 shows the ratio between the redshift evolution bias and the cumulative noise-to-signal of the BGS sample in this alternative binning scheme and where the theory predictions are made at $z = 0.2, 0.3, 0.4$. We have noted the custom offset bins in the title of the panels. We have again highlighted the incomplete selections with open hatches.

<figure><img src="image_15.png" /><figcaption>**Figure C1.** The ratio of the redshift evolution effect in the offset redshift bins to the expected statistical error for BGS sample in a simulation-based analysis. The three columns correspond to the three new BGS lens bins centred on ABACUSSUMMIT snapshots ($z = 0.2, 0.3, 0.4$). The top row shows the results for $w_p$, whereas the bottom row shows the results for $\Delta\Sigma$. The three colours show the three different absolute magnitude cuts, while the width of the curves represents the spread between different photometric samples. The colour scheme is identical to Fig. 12. We see that redshift evolution is insignificant at scales greater than $r_p > 1 h^{-1}\text{Mpc}$.</figcaption></figure>
Comparing to Fig. 12, we see larger evolution effects due to the fact that we are limited to fixed redshifts instead of the mean sample redshift. However, the qualitative trends in terms of scale and magnitude cuts remain the same, as we expect.

For galaxy–galaxy lensing, we see that the bias due to redshift evolution continues to be insignificant relative to statistical noise. For galaxy clustering, the evolution effects become important at scales $r_p \sim 1 h^{-1} \text{Mpc}$ and below, especially at higher redshift. This suggests that simulation-based analyses that intend to utilize $r_p < 1 h^{-1} \text{Mpc}$ should either carefully account for evolution systematics at the smallest scales, or use narrower redshift bins to reduce the effect of redshift evolution, or change the effective redshift.

This paper has been typeset from a TeX/Latex file prepared by the author.

<footnote n="1">¹Kavli Institute for Particle Astrophysics and Cosmology, Stanford University, 452 Lomita Mall, Stanford, CA 94305, USA</footnote>

<footnote n="2">²SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025, USA</footnote>

<footnote n="3">³Centre for Astrophysics and Supercomputing, Swinburne University of Technology, Hawthorn, VIC 3122, Australia</footnote>

<footnote n="4">⁴Perimeter Institute for Theoretical Physics, 31 Caroline St. North, Waterloo, ON N2L 2Y5, Canada</footnote>

<footnote n="5">⁵Waterloo Centre for Astrophysics, University of Waterloo, 200 University Ave W, Waterloo, ON N2L 3G1, Canada</footnote>

<footnote n="6">⁶Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA</footnote>

<footnote n="7">⁷CCAP, The Ohio State University, Columbus, OH 43210, USA</footnote>

<footnote n="8">⁸Department of Astronomy and Astrophysics, University of California, Santa Cruz, 1156 High Street, Santa Cruz, CA 95065, USA</footnote>

<footnote n="9">⁹Physics Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA</footnote>

<footnote n="10">¹⁰Physics Department, Boston University, 590 Commonwealth Avenue, Boston, MA 02215, USA</footnote>

<footnote n="11">¹¹High-Energy Physics Division, Argonne National Laboratory, 9700 S. Cass Avenue, Argonne, IL 60439, USA</footnote>

<footnote n="12">¹²Department of Physics & Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>

<footnote n="13">¹³Departamento de Física, Centro de Investigación y de Estudios Avanzados del IPN, Av. Politécnico 2508, Col. Sn. Pedro Zacatenco, Del. Gustavo A. Madero, 07010 CDMX, México</footnote>

<footnote n="14">¹⁴Observatorio Astronómico, Universidad de los Andes, Cra. 1 No. 18A-10, Edificio H, CP 111711 Bogotá, Colombia</footnote>

<footnote n="15">¹⁵Department of Physics, The University of Texas at Dallas, Richardson, TX 75080, USA</footnote>

<footnote n="16">¹⁶Institut d'Estudis Espacials de Catalunya (IEEC), E-08034 Barcelona, Spain</footnote>

<footnote n="17">¹⁷Institute of Cosmology & Gravitation, University of Portsmouth, Dennis Sciama Building, Portsmouth PO1 3FX, UK</footnote>

<footnote n="18">¹⁸Aix Marseille University, CNRS, CNES, LAM, Marseille, France</footnote>

<footnote n="19">¹⁹Departament de Física, Serra Húnter, Universitat Autònoma de Barcelona, E-08193 Bellaterra (Barcelona), Spain</footnote>

<footnote n="20">²⁰NSF's NOIRLab, 950 N. Cherry Ave., Tucson, AZ 85719, USA</footnote>

<footnote n="21">²¹National Astronomical Observatories, Chinese Academy of Sciences, A20 Datun Rd., Chaoyang District, Beijing 100012, P. R. China</footnote>

<footnote n="22">²²Faculty of Physics and Astronomy, Ruhr University Bochum, Astronomical Institute (AIRUB), German Centre for Cosmological Lensing, D-44780 Bochum, Germany</footnote>

<footnote n="23">²³Department of Physics, Kansas State University, 116 Cardwell Hall, Manhattan, KS 66506, USA</footnote>

<footnote n="24">²⁴Department of Physics and Astronomy, Sejong University, Seoul 143-747, Korea</footnote>

<footnote n="25">²⁵Centro de Investigaciones Energéticas, CIEMAT, Avenida Complutense 40, E-28040 Madrid, Spain</footnote>

<footnote n="26">²⁶Max Planck Institute for Extraterrestrial Physics, 85748 Garching, Bayern, Deutschland</footnote>

<footnote n="27">²⁷Department of Physics & Astronomy, Ohio University, Athens, OH 45701, USA</footnote>

<footnote n="28">²⁸Instituto de Física, Universidad Nacional Autónoma de México, Cd. de México C.P. 04510, México</footnote>