© The Authors 2025

1. Introduction

Weak gravitational lensing has emerged as a powerful tool in the field of cosmology, providing valuable insights into the nature of our Universe. In particular, the second-order correlation functions of weak lensing shear have become widely employed for precision cosmological measurements. These correlation functions rely on the matter power spectrum and are especially effective in determining the combined parameter S₈ = σ₈(Ω_m/0.3)^0.5, which represents a combination of the matter density parameter Ω_m and the clustering parameter σ₈. Recent Stage III surveys, including the Hyper Suprime-Cam (HSC) survey (Aihara et al. 2018), the Kilo-Degree Survey (KiDS, Kuijken et al. 2015), and the Dark Energy Survey (DES, Abbott et al. 2016; Becker et al. 2016), have successfully measured S₈ with remarkable precision, thereby establishing cosmic shear as a reliable cosmological tool (Asgari et al. 2021; Amon et al. 2022; Dalal et al. 2023; Li et al. 2023; DES and KiDS Collaboration 2023).

However, intriguing trends have been observed in the results obtained from cosmic shear experiments. These surveys consistently yield lower values of S₈ compared to what would be expected based on measurements of the cosmic microwave background (CMB) using the Planck satellite under the framework of the cosmological standard model. Referred to as the S₈-tension, this discrepancy has a significance of 2 − 3σ, and raises the possibility of issues either in the standard model of cosmology or in the analyses of cosmic shear or the CMB (Di Valentino et al. 2021; Abdalla et al. 2022).

Several effects have been identified as potential sources for the tension observed in S₈ measurements. First, ill-understood astrophysical effects, such as the impact of baryonic matter or intrinsic alignments of galaxies, might bias the measurements from cosmic shear (Semboloni et al. 2011; Troxel & Ishak 2015; Chisari et al. 2019). Second, deviations from the cosmological standard model, such as modified gravity or exotic dark matter, could play a role (Planck Collaboration XIV 2016; Di Valentino et al. 2016; Heimersheim et al. 2020). Third, observational systematics, such as the calibration of photometric redshift estimates and galaxy shape measurements can bias the measured cosmic shear signal (Huterer et al. 2006; Euclid Collaboration 2024a). Finally, the apparent tension could be caused by simplifying assumptions made in the modelling of the cosmic shear signal. Examples of such assumptions include neglecting magnification bias (Unruh et al. 2020; von Wietersheim-Kramsta et al. 2021; Duncan et al. 2022), source obscuration (Hartlap et al. 2011), spatially varying survey depth (Heydenreich et al. 2020), or the intrinsically clustered positions of source galaxies (Yu et al. 2015). Understanding these factors is crucial for comprehending the nature of the S₈-tension, and for assessing the implications for Stage IV surveys such as Euclid (Laureijs et al. 2011; Euclid Collaboration 2025a) or the Vera C. Rubin Legacy Survey of Space and Time (LSST Ivezić et al. 2019). In this study we focus specifically on the clustering of source galaxies with the lensing matter structures, which we refer to as source-lens clustering (SLC).

Conventionally, cosmic shear analyses assume that source galaxies are spatially uniformly distributed within the survey footprint, allowing the cosmic shear field to be sampled randomly across the sky. However, this assumption does not hold in reality. Source galaxies, like all galaxies, trace the underlying density field, meaning their positions correlate with the local matter structures. If the sources are divided into broad redshift bins, some sources in a bin will be at the same redshift as the matter causing a shear signal for other sources in the same bin. Therefore, the lensing signal of the farther sources is correlated to the positions of the closer sources. This correlation contributes to the measured shear correlation functions and, if not taken into account, can, in principle, bias cosmological parameter inference (Bernardeau 1998).

However, the magnitude of this effect is currently unclear. By considering a simple model that incorporates linear galaxy biases and an analytic form of the source redshift distribution, Euclid Collaboration (2024b) predicted that the effect of source-lens clustering would be significant for Stage IV surveys such as Euclid and bias constraints on Ω_m by more than 1σ. However, based on analytical calculations, Krause et al. (2021) expect the effect to be dependent on fourth-order correlations and therefore negligible, at least for Stage III surveys. Yu et al. (2015) quantified the effect with N-body simulations and found a 1–10% effect on the lensing power spectrum. However, they considered only two Gaussian tomographic bins, assumed a linear galaxy bias, and did not include realistic photometric redshift errors in their study. Thus, the magnitude of SLC for galaxies with realistic redshift distributions has not yet been determined. This paper aims to address this gap in knowledge.

A confounding factor in the discussion of SLC is that the non-uniform distribution of source galaxies causes several partially counteracting effects. First, the intrinsic clustering of sources in regions of higher density means that we sample the cosmic shear field predominantly in regions with higher shear signals. This effect, considered for example by Krause et al. (2021), increases the measured cosmic shear signal compared to the theoretical expectation for unclustered sources. Second, the clustering causes a bias on the standard estimator for cosmic shear correlation functions. This estimator is only unbiased for uniformly distributed galaxies. A third effect occurs due to the noise in photometric redshift estimates. Due to the uncertainty of photometric redshifts, galaxies at low true redshifts might be assigned to higher tomographic bins. These galaxies carry less cosmic shear signal and therefore lower the overall signal of a tomographic bin. This decrease in signal is stronger in regions with high foreground matter densities since there are also more galaxies that can be assigned to higher tomographic bins. The correlation between the decrease in signal and the foreground matter distribution causes a decrease in the measured cosmic shear signal. Our goal here is to simultaneously constrain the impact of all these SLC effects on cosmological parameter estimation with cosmic shear.

To achieve this goal, we measured the shear correlation functions for clustered and uniformly distributed source galaxies in cosmological simulations with realistically distributed galaxies and Stage III- and Stage IV-like redshift distributions. We then performed cosmological analyses for the clustered and unclustered cases to assess the effect of source-lens clustering.

The structure of this paper is as follows. In Sect. 2 we review the basics of weak gravitational lensing and cosmic shear and give a theoretical description of the source-lens clustering effect. Section 3 describes the cosmological simulations used and our steps to obtain clustered and unclustered source galaxies. We describe our correlation function measurement and cosmological parameter inference in Sect. 4, and give the resulting parameter estimates in Sect. 5. We conclude with a discussion in Sect. 6.

2. Theoretical background

2.1. Second-order cosmic shear statistics

Cosmological analyses of cosmic shear operate in the weak regime of gravitational lensing (for a review on weak lensing, see, for example, Bartelmann & Schneider 2001). They mostly analyse second-order statistics of the shear γ. The easiest of these statistics, conceptually speaking, are the shear correlation functions ξ₊ and ξ₋, defined by

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

where γ_t and γ_× are the tangential and cross-component of the shear for a galaxy pair with angular separation θ and the brackets denote an ensemble average. Assuming a flat sky, the shear γ is related to the convergence κ by the Kaiser-Squires relation (Kaiser & Squires 1993)

$$\begin{array}{c}\hfill \stackrel{\sim }{\gamma }(\mathit{\ell })={\mathrm{e}}^{2\mathrm{i}{\varphi }_{\ell }}\stackrel{\sim }{\kappa }(\mathit{\ell }),\end{array}$$(2)

where the tilde denotes Fourier transforms and ϕ_ℓ is the polar angle of the wavevector ℓ. Therefore, using the lensing power spectrum C(ℓ), defined with the Dirac delta ‘function’ δ_D, by

$$\begin{array}{c}\hfill {(2\pi )}^2{\delta }_{\mathrm{D}}(\mathit{\ell }+{\mathit{\ell }}^{\prime })C(\ell )=⟨\widehat{\kappa }({\mathit{\ell }}^{\prime })\widehat{\kappa }(\mathit{\ell })⟩,\end{array}$$(3)

the shear correlation functions can be modelled as

$$\begin{array}{c}\hfill {\xi }_{\pm }(\theta )={\displaystyle {\int }_0^{\infty }}\frac{\text{d}\ell \ell }{2\pi }{J}_{0,4}(\ell \theta )C(\ell ),\end{array}$$(4)

where J_i is the i-th order Bessel function. Consequently, the impact of any modelling effect on the shear correlation function can be estimated equivalently on the lensing power spectrum.

The convergence κ is a normalized surface mass density and thus related to the matter density contrast δ(χθ, χ) at angular position θ and comoving distance χ via

$$\begin{array}{c}\hfill \kappa (\mathit{\theta })={\displaystyle {\int }_0^{\infty }}\text{d}\chi \delta (\chi \mathit{\theta },\chi ){\displaystyle {\int }_{\chi }^{\infty }}\text{d}\chi \prime W(\chi ,\chi \prime )p(\chi \prime ;\mathit{\theta }),\end{array}$$(5)

where p(χ; θ) is the probability density of source galaxies at angular position θ and comoving distance χ and W is the lensing efficiency kernel, which is, assuming a flat Universe,

$$\begin{array}{c}\hfill W(\chi ,\chi \prime )=\frac{3{\mathrm{\Omega }}_{\mathrm{m}}{H}_0^2}{2c^2}\frac{\chi }{a(\chi )}\frac{\chi \prime -\chi }{\chi \prime }·\end{array}$$(6)

The probability density p is related to the number density n of source galaxies by

$$\begin{array}{c}\hfill p(\chi ;\mathit{\theta })=\frac{n(\chi \mathit{\theta },\chi )}{\int \text{d}\chi n(\chi \mathit{\theta },\chi )}·\end{array}$$(7)

Usually, one assumes that source galaxies are distributed uniformly across the survey area. Then, p is only a function of the comoving distance χ, and the convergence κ₀(θ) is given by

$$\begin{array}{c}\hfill {\kappa }_0(\mathit{\theta })={\displaystyle {\int }_0^{\infty }}\text{d}\chi \delta (\chi \mathit{\theta },\chi ){\displaystyle {\int }_{\chi }^{\infty }}\text{d}\chi \prime W(\chi ,\chi \prime )p(\chi \prime ).\end{array}$$(8)

Under the extended Limber approximation (Limber 1953; Kaiser 1998; LoVerde & Afshordi 2008), the lensing power spectrum then becomes

$$\begin{array}{c}\hfill C(\ell )={\displaystyle {\int }_0^{\infty }}\text{d}\chi \frac{1}{{\chi }^2}{\left[{\displaystyle {\int }_{\chi }^{\infty }}\text{d}\chi \prime p(\chi \prime )W(\chi ,\chi \prime )\right]}^{2}P(\frac{\ell +1/2}{\chi },\chi ),\end{array}$$(9)

where P(k, χ) is the matter power spectrum. This form of the lensing power spectrum is usually assumed in cosmic shear analyses.

2.2. The source-lens clustering effect

However, the assumption of uniformly distributed source galaxies is not correct. As is the case for any other galaxy population, the sources are tracers of the matter distribution and are therefore spatially clustered. This clustering invalidates Eq. (8), as the source density n_s now depends on angular position θ. Instead, Eq. (5) needs to be used, which leads to different effects outlined in this section.

2.2.1. Estimator bias (EB)

As noted by Yu et al. (2015), SLC depends on the estimator used for the shear correlation functions. Different estimators incur different SLC-induced estimator biases. For example, when using an estimator based on pixelized shear maps, the noise in each pixel depends on the number of source galaxies in the pixel. The correlation between this number and the shear signal then leads to a correlation between the signal estimate and the noise (Gatti et al. 2024). Here, though, we consider a catalogue-based estimator for the shear correlation functions which does not require pixelization. This estimator ${\widehat{\xi }}_{\pm }$ is

$$\begin{array}{c}\hfill {\widehat{\xi }}_{\pm }(\theta )=\frac{{\displaystyle {\sum }_i}{\displaystyle {\sum }_j}({\gamma }_{\mathrm{t}}^i{\gamma }_{\mathrm{t}}^j\pm {\gamma }_{\times }^i{\gamma }_{\times }^j)▵(|{\mathit{\vartheta }}_i-{\mathit{\vartheta }}_j|;\theta )}{{\displaystyle {\sum }_i}{\displaystyle {\sum }_j}▵(|{\mathit{\vartheta }}_i-{\mathit{\vartheta }}_j|;\theta )},\end{array}$$(10)

where the sums go over all source galaxies, γ_tⁱ and γ_×ⁱ are the tangential and cross shear of the i-th galaxy and △(ϑ; θ) is one if ϑ lies in the θ-bin and zero otherwise. Using γ_t/×(ϑ, χ) for the tangential (cross) shear of a hypothetical source at angular position ϑ and distance χ, the expectation value of this estimator is

$$\begin{array}{cc}& ⟨{\widehat{\xi }}_{\pm }⟩(|\mathit{\vartheta }-{\mathit{\vartheta }}^{\prime }|)\hfill \\ \hfill & =\langle \{{\displaystyle \int \text{d}}{\chi }_1{\displaystyle \int \text{d}}{\chi }_2{\displaystyle \int {\mathrm{d}}^2}\vartheta p({\chi }_1;\mathit{\vartheta })p({\chi }_2;{\mathit{\vartheta }}^{\prime })\hfill \\ \hfill & \times [{\gamma }_{\mathrm{t}}(\mathit{\vartheta },{\chi }_1){\gamma }_{\mathrm{t}}({\mathit{\vartheta }}^{\prime },{\chi }_2)\pm {\gamma }_{\times }(\mathit{\vartheta },{\chi }_1){\gamma }_{\times }({\mathit{\vartheta }}^{\prime },{\chi }_2)]\}\hfill \\ \hfill & \times {\left[{\displaystyle \int \text{d}}{\chi }_1{\displaystyle \int \text{d}}{\chi }_2{\displaystyle \int {\mathrm{d}}^2}\vartheta p({\chi }_1;\mathit{\vartheta })p({\chi }_2,{\mathit{\vartheta }}^{\prime })\right]}^{-1}\rangle .\hfill \end{array}$$(11)

If we neglect for the moment the correlation between the source density n_s and the shear (which we treat in the next two subsections), we can swap the averaging and the integrations such that the expectation value becomes

$$\begin{array}{cc}& ⟨{\widehat{\xi }}_{\pm }⟩(|\mathit{\vartheta }-\mathit{\vartheta }\prime |)\hfill \\ \hfill & =[{\displaystyle \int \text{d}}{\chi }_1{\displaystyle \int \text{d}}{\chi }_2{\displaystyle \int {\mathrm{d}}^2}\vartheta ⟨p({\chi }_1,\mathit{\vartheta })p({\chi }_2,{\mathit{\vartheta }}^{\prime })⟩\hfill \\ \hfill & \times ⟨{\gamma }_{\mathrm{t}}(\mathit{\vartheta },{\chi }_1){\gamma }_{\mathrm{t}}(\mathit{\vartheta }\prime ,{\chi }_2)\pm {\gamma }_{\times }(\mathit{\vartheta },{\chi }_1){\gamma }_{\times }(\mathit{\vartheta }\prime ,{\chi }_2)⟩]\hfill \\ \hfill & \times {\left[{\displaystyle \int \text{d}}{\chi }_1{\displaystyle \int \text{d}}{\chi }_2{\displaystyle \int {\mathrm{d}}^2}\vartheta ⟨p({\chi }_1,\mathit{\vartheta })p({\chi }_2,\mathit{\vartheta }\prime )⟩\right]}^{-1}.\hfill \end{array}$$(12)

For unclustered sources, p(χ, ϑ) reduces to p(χ) and the expectation value in Eq. (12) becomes exactly ξ_±. This occurs even for clustered sources in the limit of infinitely narrow redshift distributions, that is p(χ, ϑ)∝δ_D(χ − χ′), since then the p-correlations in numerator and denominator cancel. In reality, the shear correlation functions can only be measured for finitely broad p(χ), so in general, ${\widehat{\xi }}_{\pm }$ is biased. We refer to this bias as estimator bias (EB).

The estimator bias is also affected by shear weights. Usually, the measured galaxy ellipticities have an assigned weight that depends on the statistical error or bias of the shape measurement. This weight depends on properties of the considered galaxy, such as magnitude and size, but also on the environment of the galaxy, for example whether it is blended with another object. Consequently, the shear weight is correlated with the source number density, causing further bias to the estimator.

2.2.2. Intrinsic clustering (IC)

For the EB effect we neglected that the shear signal itself is correlated with the source galaxy density. This correlation occurs, because due to source clustering the source redshift distribution is not uniform across the sky but depends on the local matter density contrast. Therefore, the number density n needs to be replaced by

$$\begin{array}{c}\hfill n(\chi )\to n(\chi )[1+{\delta }_{\mathrm{g}}(\chi \mathit{\theta },\chi )],\end{array}$$(13)

where δ_g is the (three-dimensional) galaxy number density contrast. This relation differs from the ‘projected’ SLC ansatz (e.g. Schmidt et al. 2009; Euclid Collaboration 2024b), where δ_g(χ θ, χ) is replaced by a projected density contrast δ_g^p(θ) which no longer depends on χ. However, the projected ansatz only agrees with Eq. (13) for very narrow bins, since then ${\delta }_{\mathrm{g}}^{\mathrm{p}}(\theta )\simeq {\delta }_{\mathrm{g}}(\overline{\chi }\theta ,\overline{\chi })$, where $\overline{\chi }$ is the average distance for galaxies in the redshift bin (Krause et al. 2021).

Under the assumption of a linear galaxy bias b,

$$\begin{array}{c}\hfill n(\chi )\to n(\chi )[1+b\delta (\chi \mathit{\theta },\chi )].\end{array}$$(14)

Applying this change to Eq. (8) means that the true convergence κ is given by

$$\begin{array}{c}\hfill \kappa (\mathit{\theta })={\kappa }_0(\mathit{\theta })+\mathrm{\Delta }\kappa (\mathit{\theta }),\end{array}$$(15)

with

$$\begin{array}{c}\hfill \mathrm{\Delta }\kappa (\mathit{\theta })=b{\displaystyle {\displaystyle {\int }_0^{\infty }}}\text{d}\chi \delta (\chi \mathit{\theta },\chi ){\displaystyle {\int }_{\chi }^{\infty }}\text{d}\chi \prime W(\chi ,\chi \prime )p(\chi \prime )\delta (\chi \prime \mathit{\theta },\chi \prime ).\end{array}$$(16)

The additional term depends on both the density contrast δ(χθ, χ) in front of the source galaxies (i.e. the lens planes), and the density contrast δ(χ′θ, χ′) surrounding the source galaxies (i.e. in the ‘source planes’), which gives the effect its name: source-lens clustering.

As we show in Appendix A, the Δκ induces an additional term ΔC to the lensing power spectrum¹ given by

$$\begin{array}{cc}\hfill \mathrm{\Delta }C(\mathit{\ell })& =b^2{\displaystyle {\int }_0^{\infty }}\text{d}\chi {\displaystyle {\int }_{\chi }^{\infty }}\text{d}\chi \prime \frac{p^2(\chi \prime )}{{\chi \prime }^2}\frac{W^2(\chi ,\chi \prime )}{{\chi }^2}\hfill \\ \hfill & \times {\displaystyle \int \frac{{\mathrm{d}}^{2}L}{{(2\pi )}^2}}P(\frac{L}{\chi },\chi )P(\frac{|\mathit{\ell }-\mathit{L}|}{\chi \prime },\chi \prime ).\hfill \end{array}$$(17)

We note that this expression relies on the assumption of a linear galaxy bias, which might not be fulfilled at small scales. Therefore, instead of simply using Eq. (17) to predict the effect, we measure it in the following using cosmological simulations with realistic galaxy distributions.

2.2.3. Tomographic bin contamination (TC)

Another effect occurs when cosmic shear analyses are performed in tomographic bins divided by photometric redshift estimates. Photometric redshifts of individual galaxies are typically noisy. Therefore, a galaxy associated with a tomographic bin between z₁ and z₂ might have a true redshift well outside of this range. On average, this photometric redshift scatter is taken care of by using a p(χ) for each tomographic bin, which was calibrated with spectroscopic redshift estimates (see e.g. Wright et al. 2020). However, the deviation from the estimated p(χ) along each line-of-sight correlates with the galaxy number density. Along a line-of-sight with an overdensity in the front, there are more galaxies which can potentially be assigned to a higher tomographic bin, thus contaminating this bin. These contaminating galaxies carry a lower lensing signal. Consequently, in denser regions, the cosmic shear signal is suppressed. This effect also impacts weak lensing estimates of galaxy cluster masses (Köhlinger et al. 2015). Due to photometric redshift scatter, cluster member galaxies can be erroneously treated as background source galaxies. They carry no lensing signal due to the galaxy cluster, thus bias the mass estimate low.

All SLC effects are expected to be much smaller than the cosmic shear signal itself. For example, since Eq. (17) depends on the fourth-order moment of the matter distribution, one could assume that this effect is negligible in a realistic weak lensing survey, where shape measurement uncertainties or the distribution of intrinsic galaxy shapes cause significant noise contributions. However, while this argument was made for stage III surveys (Krause et al. 2021), it does not necessarily hold for Euclid. As illustration, we compare in Fig. 1 the shear correlation function ξ₊ for auto-correlations of tomographic bins as measured by the stage III surveys KiDS (from Asgari et al. 2021), DES (from Amon et al. 2022), and HSC (from Hamana et al. 2020) with their reported error bars, to the expectation for Euclid. Due to its greater depth, Euclid will observe more tomographic bins at higher redshift, which carry larger cosmic shear signals. Due to its larger area and higher number density, the noise of the cosmic shear measurements will be significantly reduced, leading to an order of magnitude increase in S/N. This high precision requires us to have at least equally accurate models. The first step to these is understanding how model assumptions impact the cosmological inference.

Fig. 1. Measured and expected weak lensing signals for stage III surveys and Euclid. Upper panel: Shear auto-correlations ξ+ for the stage III surveys KiDS (from Asgari et al. 2021), DES (from Amon et al. 2022), and HSC (from Hamana et al. 2020) with their reported errorbars and the expectation for Euclid. Error bars for Euclid are the square root of the covariance matrix diagonal, calculated as described in Sect. 4.1. The colour indicates the median redshift of the tomographic bins, with darker colours corresponding to higher redshifts. Lower panel: S/N of ξ+, given by dividing the measurement or prediction by the error bars.

3. Simulated datasets

We used two sets of simulated data to estimate the impact of source-lens clustering on cosmological analyses. One of these datasets, based on the Marenostrum Institut de Ciències de l’Espai (MICE) simulations, mimics the redshift distribution and tomographic set-up of the KiDS-1000, while the other, based on the Euclid Flagship2 (FS2) simulation, mimics the expected redshift distribution and tomographic set-up of the Euclid wide survey. We describe their creation and properties in the following. Simulated galaxy catalogues with galaxy properties and weak lensing information from both MICE and the FS2 are accessible on CosmoHub² (Carretero et al. 2017; Tallada et al. 2020).

3.1. The Marenostrum Institut de Ciències de l’Espai (MICE) simulations

We created a KiDS-like dataset based on the MICE simulation (Fosalba et al. 2015a), a dark-matter-only N-body simulation. The simulation used a flat ΛCDM cosmology with Ω_m = 0.25, Ω_Λ = 0.75, Ω_b = 0.044, n_s = 0.95, σ₈ = 0.8, h = 0.7. It evolved 4096³ particles with a mass of 2.93 × 10¹⁰ h⁻¹ M_⊙ inside a box with side length 3072 h⁻¹ Mpc from an initial redshift of z = 100 to today. In the dark matter particle distribution, halos were identified using a friends-of-friends halo finder (Crocce et al. 2015). The halos were then populated up to redshift z = 1.4 with galaxies, following a mixture of halo abundance matching and a halo occupation distribution model to match both the spatial and luminosity distribution of observed galaxies (Carretero et al. 2015). We used the second version of this galaxy catalogue, called MICE2.

The MICE2 galaxies have assigned fluxes in many photometric bands. To generate KiDS-like photometry, we used the noiseless simulated photometry from MICE, similar to the approach presented in van den Busch et al. (2020). We chose the filters that are most similar to the KiDS-1000 filter set-up, which consists of the Very Large Telescope (VLT) Survey Telescope (VST) OmegaCAM ugri bands and the Visible and Infrared Survey Telescope for Astronomy (VISTA) VIRCAM ZYJHK_s bands (Kuijken et al. 2019). For the ugriZ bands, we found the Sloan Digital Sky Survey (SDSS) ugriz filters provided in MICE to be most similar to the ones used in KiDS-1000. For the near-infrared imaging in the Y band, the DES y filter is used as a reasonable approximation. The filters of the VISTA JHK_s bands are readily available in MICE. As presented in Fosalba et al. (2015b), we applied a redshift-dependent evolution correction to the MICE magnitudes resulting in a set of magnitudes m^evo. To estimate the observed magnitudes, we approximated the flux noise of KiDS-1000 (Kuijken et al. 2019).

An effective projected galaxy radius size r_eff is estimated taking into account the half-light radius of a galaxy’s bulge and disk and the bulge-to-total emitted flux fraction given in MICE. The estimated projected galaxy size translates into a simulated aperture which mimics the Gaussian Aperture and Point Spread Function (GAaP) apertures used in KiDS (Kuijken 2008). The aperture major and minor axis A_i, x and B_i, x, respectively, are calculated per galaxy i and per photometric filter x through

$$\begin{array}{cc}\hfill A_{i,x}& =min(\sqrt{r_{\mathrm{eff},i}^2+{\sigma }_{\mathrm{PSF},x}^2+r_{\min }^2},2.0\stackrel{″}{.}),\hfill \end{array}$$(18)

$$\begin{array}{cc}\hfill B_{i,x}& =min(\sqrt{{\left({\left[\frac{b}{a}\right]}_i{r}_{\mathrm{eff},i}\right)}^2+{\sigma }_{\mathrm{PSF},x}^2+r_{\min }^2},2.0\stackrel{″}{.}),\hfill \end{array}$$(19)

with the mean seeing σ_PSF, x and the projected axis ratio of the bulge ${\left[\frac{b}{a}\right]}_i$ with the semi-minor and semi-major axis b and a, respectively, which is given for all the MICE sources. The minimum radius r_min for the mocks was set to 0$\stackrel{″}{.}$3, which is the maximum difference in point spread function (PSF) size of the KiDS-1000 data between the different photometric bands. The maximum aperture axis is set to be 2$\stackrel{″}{.}$0, as for the KiDS-1000 data.

Considering the limiting magnitudes m_lim, x of the KiDS observations, we calculated the flux error

$$\begin{array}{c}\hfill \mathrm{\Delta }{f}_{i,x}={10}^{-0.4(m_{\lim ,x}-48.6)}\sqrt{\frac{A_{i,x}{B}_{i,x}}{{\sigma }_{\mathrm{PSF},x}^2}},\end{array}$$(20)

where we included the mock apertures πA_i, xB_i, x to increase the noise for larger apertures compared to the PSF πσ_PSF, x². Given the estimated noise and the evolution-corrected model fluxes f^evo, a flux realization was computed from $f_{i,x}^{\mathrm{obs}}\sim \mathcal{N}(f_{i,x}^{\mathrm{evo}},\mathrm{\Delta }{f}_{i,x})$. Converting to magnitudes $m_{i,x}^{\mathrm{obs}}$ and their errors Δm_i, x, photometric redshift can be computed with the template fitting photo-z algorithm BPZ (Benitez 2000). We excluded all objects with a signal-to-noise ratio of less than 1 in the r band since we label these galaxies as undetected, and therefore they would have no shape measurement. Finally, we performed a k Nearest Neighbours (kNN) matching between the mock catalogue and the KiDS data based on the 9 magnitudes ugrizYJHK_s to assign a mock shape weight recalweight from the real galaxies to the mock ones. This is then used to select the objects with non-zero shape weights, which were used as the source sample in KiDS.

The MICE simulations were designed for application to gravitational lensing surveys, and thus the lensing observables shear and convergence were computed for the simulation as described in Fosalba et al. (2015b) using the approach by Fosalba et al. (2008) under the Born approximation. From this, one obtains full sky maps of the lensing convergence for sources at 265 different redshifts between 0 and 1.4. These κ-maps were converted to maps of the weak lensing shear γ, using first their decomposition in harmonic space,

$$\begin{array}{cc}\hfill \kappa (\mathit{\vartheta })& ={\displaystyle \sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{\ell }}{\widehat{\kappa }}_{\ell m}{Y}_{\ell m}(\mathit{\vartheta }),\hfill \end{array}$$(21)

$$\begin{array}{cc}\hfill \gamma (\mathit{\vartheta })& ={\displaystyle \sum _{\ell =2}^{\infty }\sum _{m=-\ell }^{\ell }}{\widehat{\gamma }}_{\ell m}{}_2{Y}_{\ell m}(\mathit{\vartheta }),\hfill \end{array}$$(22)

where we neglect B-modes and use the Laplace spherical harmonics Y_ℓm and the spin-s-weighted spherical harmonics _sY_ℓm. Then, we apply the inverse Kaiser–Squires relation (Kaiser & Squires 1993),

$$\begin{array}{c}\hfill {\widehat{\gamma }}_{\ell m}=-\sqrt{\frac{(\ell +2)(\ell -1)}{\ell (\ell +1)}}{\widehat{\kappa }}_{\ell m}.\end{array}$$(23)

Each galaxy at angular position ϑ and redshift z in the MICE2 catalogue was then assigned the shear at position ϑ of the weak lensing map for redshift z. The galaxies were also assigned deflected positions, arising from magnification of the sources by the density field. As we are not interested in magnification effects on the cosmic shear signal, we use throughout the original, unmagnified positions of the galaxies.

In the following, we use both the MICE2 catalogue and the shear maps to create clustered and unclustered source galaxy samples, as detailed in Sect. 3.3. For this, it is sufficient for us to use a rather small area, namely the rectangular region with right ascension α ∈ (40°, 50°) and declination δ ∈ (20°, 50°).

3.2. The Euclid flagship simulation 2

Since the MICE galaxies are limited to z < 1.4, it is not possible to construct a Euclid -like sample from them, as Euclid will observe a significant number of galaxies at higher redshifts. Therefore, we used the Euclid flagship simulations 2 (FS2, Euclid Collaboration 2025b) to create a Euclid-like sample.

The FS2 used 16 000³ particles with a mass of 10⁹ h⁻¹ M_⊙ and evolved them in a simulation box with side length 3600 h⁻¹ Mpc with the PKDGRAV3 code (Potter & Stadel 2016). This larger box size allows for ray tracing and galaxy distribution up to z = 3.0. The simulation assumed a flat ΛCDM cosmology with Ω_m = 0.319, Ω_b = 0.049, Ω_Λ = 0.681 − Ω_γ − Ω_ν, with a photon density parameter Ω_γ = 0.00005509, and a neutrino density parameter Ω_ν = 0.00140343. The reduced Hubble constant is h = 0.67, the scalar spectral index of the initial fluctuations n_s = 0.96, and its amplitude A_s = 2.1 × 10⁻⁹, which corresponds to σ₈ ≃ 0.813.

Dark matter halos were found using the friends-of-friends halo finder ROCKSTAR (Behroozi et al. 2013). These halos were populated with central and satellite galaxies according to a halo occupation distribution. Central galaxies were placed at the halo centres, while satellite galaxies are distributed inside the halo following a triaxial Navarro-Frenk-White profile. The galaxies were assigned luminosities with an abundance-matching approach.

The galaxies were assigned photometric redshifts using the nearest-neighbour photometric redshift (NNPZ) pipeline, which estimates redshift probability distributions for each galaxy (see Euclid Collaboration 2020a, for a comparison of this and other redshift estimation algorithms). We used the mode of these probability distributions as photometric redshift estimate. The pipeline used a reference sample of two million objects simulated at the depth of the Euclid calibration fields. Since the full depth of these calibration fields will only be available towards the end of Euclid’s observing run, the accuracy of the derived photometric redshifts is optimistic for the first Euclid data release. However, it is realistic for the final Euclid data releases, which have the highest statistical power and are therefore the most strongly affected by subtle effects such as SLC.

Similarly to the MICE, full-sky shear and convergence maps were computed for the FS2 with the ray-tracing approach by Fosalba et al. (2008). This approach produced shear and convergence maps at 115 source redshifts between z = 2.3 and today, given as healpix maps with N_side = 8192, which corresponds to an angular resolution of $\simeq 0\stackrel{\prime }{.}43$. While the shear maps were computed for the full-sky, galaxy catalogues are (at time of writing) only available for an octant, stretching for right ascensions between 145° and 235° and declinations between 0° and 90°. Therefore, our analysis is limited to this area.

3.3. Creating clustered and unclustered galaxy samples

Both the MICE and the FS2 contain realistically clustered source galaxies. We divided these galaxies into different tomographic bins based on their photometric redshifts. For the KiDS-like set-up, we used the five bins defined by Hildebrandt et al. (2021) and also used by Asgari et al. (2021), while for the Euclid-like set-up, we used the 13 tomographic bins defined by Euclid Collaboration (2021a), shown in Table 1. We show the distributions of true redshifts for these bins in Fig. 2. Both samples exhibit galaxy biases between 0.5 and 1.5 (increasing with redshift), which is comparable to what is expected in a real lensing survey (see Appendix B)

We created two new galaxy catalogues for each tomographic bin based on the shear maps and the original catalogue. The first of these, which we refer to as the ‘clustered’ catalogue, was used purely for validation purposes. To create it, we took the angular positions θ and redshift z of each galaxy in the original catalogue and assigned it the shear γ at position θ on the shear map closest in redshift to z.

The second catalogue we created, which we refer to as ‘unclustered’, is created similarly to the clustered catalogue but crucially without using the galaxies’ angular positions. Instead, for each galaxy in the catalogue we drew a random angular position from a uniform distribution and assigned the shear at this new position and the galaxies’ redshift to the galaxy. In this way, the galaxies in the unclustered catalogue have the same average number density and redshift distribution as those in the clustered and original catalogue but uniformly distributed positions. The SLC effect is thus present only in the clustered catalogue.

To test the impact of the individual SLC effects (see Sect. 2.2), we created two more catalogues from MICE. The first of these consists of taking the true galaxy positions on the sky and assigning to them the shear at positions rotated 90° away with respect to the survey patch center. In this way, the clustering of the source galaxies is no longer correlated to the shear signal. This eliminates the IC and TC effects but retains the EB since the sources still have a non-zero angular correlation function ω.

The last catalogue we created is designed to include the EB and the IC effects but not the TC. For each tomographic bin, we randomly subselected galaxies from the full MICE catalogue such that their true redshifts follow the redshift distribution of the bin but without using the photometric redshifts. In that way, the selection of galaxies in a tomographic bin is no longer dependent on the photometric redshift scatter, eliminating the TC effect.

4. Cosmological analysis

4.1. COSEBI estimation, modelling and covariance

To match the fiducial cosmological analysis of Asgari et al. (2021), we performed our cosmological inference using the Complete Orthogonal Sets of E/B-Integrals (COSEBIs, Schneider et al. 2010). These consist of the E-modes E_n and the B-modes B_n, where, to first order, only E-modes can be generated by weak lensing. They can be obtained as weighted integrals over the shear correlation functions ξ₊ and ξ₋ with

$$\begin{array}{cc}\hfill E_n& =\frac{1}{2}{\displaystyle {\int }_{{\theta }_{\min }}^{{\theta }_{\max }}}\text{d}\theta \theta [T_{+n}(\theta ){\xi }_{+}(\theta )+T_{-n}(\theta ){\xi }_{-}(\theta )],\hfill \end{array}$$(24)

$$\begin{array}{cc}\hfill B_n& =\frac{1}{2}{\displaystyle {\int }_{{\theta }_{\min }}^{{\theta }_{\max }}}\text{d}\theta \theta [T_{+n}(\theta ){\xi }_{+}(\theta )-T_{-n}(\theta ){\xi }_{-}(\theta )],\hfill \end{array}$$(25)

where T_±n(θ) are filter functions defined such that θ is bounded by θ_min and θ_max (Schneider et al. 2010).

To obtain the COSEBIs, we first measured the shear correlation functions ξ₊ and ξ₋ for each combination of the tomographic bins. We used 4000 radial bins spaced between 0$\stackrel{\prime }{.}$5 and 300′ logarithmically for both the KiDS and Euclid-like case. The measurements are conducted with treecorr (Jarvis et al. 2004). The correlation functions are then converted to the first 5 COSEBIs. We chose this number n_max of COSEBIs for two reasons. First, the information content as a function of n_max saturates between n_max = 5 and n_max = 10 (Asgari et al. 2012). Second, the signal-to-noise ratio diminishes with higher n and even for the Euclid-like case, the signal-to-noise ratio for E₆ is lower than 3.8 for all redshift bins. Thus, we do not expect a significant improvement when including more COSEBIs. Nevertheless, determining the optimal choice of n_max, along with the considered angular scales is an important goal for future work. We modelled E_n using

$$\begin{array}{cc}\hfill E_n& ={\displaystyle {\int }_0^{\infty }}\frac{\text{d}\ell \ell }{2\pi }C(\ell )W_n(\ell ),\hfill \end{array}$$(26)

where the W_n are Hankel transforms of the T_±n and given in Asgari et al. (2012). The lensing power spectrum C(ℓ) was modelled as in Joachimi et al. (2021) from Eq. (9). We used the redshift distributions p_z(z) displayed in Fig. 2, but included free parameters δz for shifts in the p_z(z). This δz are the difference between the mean of the estimated p_z(z) of each tomographic bin and the true redshift distribution. We used the linear power spectrum model from CAMB (Lewis et al. 2000; Howlett et al. 2012).

Fig. 2. Redshift distribution of KiDS-1000-like galaxies in the MICE (top) and Euclid-like galaxies in the FS2 (bottom).

We modelled the non-linear matter power spectrum P(k) using HMCode-2020 (Mead et al. 2021) for the KiDS-like analysis. We fixed the baryon feedback parameters η and A_b to 0 and 3.13, respectively, corresponding to the dark-matter-only case. However, we could no use the same prescription for the Euclid-like case. As shown in Fig. 3, the measured shear correlation functions for the KiDS-like case show a deviation of up to 2% to the HMCode-based prediction for the simulation cosmology for sources at redshifts between 0.9 and 1.2, which is well within the measurement uncertainty of KiDS-1000. For Euclid, though, this deviation is of the order of the measurement uncertainty, indicating that the model is not accurate enough. Instead, for the Euclid-set-up, we used the non-linear matter power spectrum predicted by the EuclidEmulator2 (Euclid Collaboration 2021b), which agrees with the measurements within 1%, which is within the uncertainty.

Fig. 3. Fractional difference between ξ+ modelled with HMCode-2020 (red, dashed) and the EuclidEmulator (black, dash-dotted) to the measurements in the FS2 for sources at 0.9 < z < 1.2. The blue area corresponds to the KiDS-1000 uncertainty and the yellow area to the Euclid uncertainty. The correlation function is binned here in 10 logarithmic bins between 1$\stackrel{\prime }{.}$25 and 83′.

Though the EuclidEmulator has better accuracy than the HMCode model, it is calibrated on a smaller range of cosmological parameters than the HMCode. Therefore, we must assume tight priors for the inference in the Euclid-like set-up. Furthermore, we could not use the EuclidEmulator for the KiDS-like set-up since the required tight prior range would bias the cosmological inference.

We further included the impact of intrinsic alignments on the cosmic shear power spectrum via the non-linear alignment (NLA) model (Bridle & King 2007; Kirk et al. 2012; see Equations 13–16 in Joachimi et al. 2021), which gave us the free parameter A_IA characterizing the strength of intrinsic alignments.

For the cosmological inference, we also required a covariance estimate. We modelled the covariance analytically, using the code described in Reischke et al. (2024), which follows the approach in Joachimi et al. (2021). To summarize, we first calculated the three parts of the covariance of the lensing power spectrum, which are the Gaussian part, depending only on the matter power spectrum, the dispersion of intrinsic galaxy shapes and source galaxy number density; the intra-survey non-Gaussian part, depending on the matter trispectrum at modes smaller than the survey; and the super-sample covariance, depending on modes larger than the survey. The sum of these terms was then converted to a covariance estimate for the COSEBIs by convolution with a suitable kernel function following the approach by Asgari et al. (2021). The analytical covariance does not include the impact of SLC. However, we show in Appendix C that this impact is negligible.

For the covariance estimate, we used the component-wise shape noise σ_ϵ = 0.28, 0.27, 0.28, 0.27, and 0.28 for each bin of the KiDS-like set-up. For the Euclid-like set-up, we assumed σ_ϵ = 0.21 for all bins, which is the same as was assumed in Euclid Collaboration (2020b) and Euclid Collaboration (2024b).

4.2. Cosmological parameter estimation

Our cosmological parameter estimation was conducted using cosmosis (Zuntz et al. 2015) and used the priors in Table 2. For the KiDS-1000-like set-up, we sampled the posterior for five cosmological parameters (Ω_c h², Ω_b h², S₈, h, n_s) and six nuisance parameters (A_IA and the redshift distribution shifts for each tomographic bin). We used the same priors as Asgari et al. (2021, see Table 2), except for the parameter Ω_b h², for which the true value in the MICE, 0.255, is very close to the prior boundary in Asgari et al. (2021), 0.26. We therefore extended the prior on the right-hand side.

For the Euclid-like case, we sampled over the power spectrum normalization A_s instead of S₈ to use the fiducial parameter of the EuclidEmulator. We also included two parameters for the equation of state w of dark energy, namely w₀ and w_a defined by

$$\begin{array}{c}\hfill w(z)=w_0+w_a\frac{z}{1+z},\end{array}$$(27)

and the redshift distribution shifts of eight more tomographic bins. Due to the limited range of the EuclidEmulator, we used the tighter priors in Table 2 for the Euclid-like set-up. We also changed the prior for δz. While for the KiDS-like analysis, we assumed a correlated prior on δz, using the correlation matrix by Hildebrandt et al. (2021), for Euclid, we assumed an uncorrelated Gaussian prior. For the Euclid priors, we chose a width of $0.001(1+\overline{z})$.

For both analyses, we assumed a Gaussian likelihood and used the nested sampler nautilus (Lange 2023). We also performed cosmological parameter estimation where we fixed all nuisance parameters to zero. This tests whether the parameter constraints are biased due to SLC if we strictly assume no intrinsic alignments to occur. From the nested sampler, we obtained estimates for the parameter values that maximize the marginal posterior distributions and their 1σ uncertainties. However, these generally differ from the parameter values that maximize the full (high-dimensional) posterior, known as Maximum A-Posteriori (MAP). We therefore ran an optimization procedure using the Nelder–Mead algorithm (Nelder & Mead 1965) to find the MAP for each analysis, following the same approach as in Asgari et al. (2020).

5. Results

5.1. KiDS-1000-like set-up

We first consider the KiDS-1000-like set-up. In the upper right corner of Fig. 4, we show the E_n measured for the MICE simulations using the original, clustered, and unclustered source galaxy catalogues, together with the KiDS-1000 uncertainties. We note that these uncertainties are not error estimates for the measurements in the mock catalogues, since they were obtained from a different survey area and did not include shape noise. We also note that the B-modes B_n were consistent with zero for all three cases. We consider in the lower left half of Fig. 4 the difference between the clustered and the unclustered E_n, and between the original and unclustered E_n normalized by the KiDS-1000 uncertainty.

Fig. 4. Measurement for the KiDS-like set-up. The x-axis denotes the COSEBI order n, each panel shows the measurement for one combination of tomographic bins, whose numbers are given in the upper right corner of each panel. Upper right: COSEBI E-modes En measured for the original MICE galaxies (black points), the clustered catalogue (red plusses) and the unclustered catalogue (blue crosses) with the KiDS-1000 uncertainty (grey band). Lower left: Difference between En for the original and clustered sources (black points) and between the clustered and unclustered sources (blue crosses), divided by the KiDS-1000 uncertainty. The original sources have shear and position taken directly from the MICE galaxy catalogue, while the clustered sources have their position taken from the catalogue and the shear assigned from the shear maps. The unclustered sources have uniformly distributed positions and shear assigned from the shear maps.

The COSEBIs for the original and clustered source catalogue are identical. This suggests that our method of assigning shears from the MICE shear map to galaxy positions agrees with the original creation of the MICE catalogue.

The impact of SLC is given by the difference between the clustered and the unclustered E_n normalized by the KiDS-1000 uncertainty. The effect shows a dependence on tomographic bins. While SLC decreases the E_n for the lowest tomographic bin, it increases them at higher bins. This increase is stronger for larger n. The deviation between the E_n for the clustered and unclustered sources is strongest for the lowest tomographic bin (up to 20%). However, since the uncertainty is also the largest there, the ratio of deviation to uncertainty peaks for the combination of the third and fourth bin. For all tomographic bins, the SLC has an effect of less than 0.5σ.

Fig. 5. Parameter constraints for the KiDS-like set-up using all parameters (lower corner) and with fixed nuisance parameters (upper corner), either for a data vector including SLC (red and pink) or without SLC (dark blue and cyan). Constraints for all sampled parameters can be found in Tables 3 and 4.

Using the COSEBIs, we first perform a cosmological inference with all nuisance parameters. In Table 3 we report the maxima of the marginalized posteriors together with their 1σ uncertainty and the MAP values. We also give the shifts in the maxima, normalized by the 1σ uncertainties and the shifts in the MAP due to the inclusion of SLC in the data vector. Figure 5 shows the corresponding posteriors for S₈, σ₈, Ω_m and A_IA in dark blue and red. The posteriors for all varied parameters are given in Appendix E.

Similarly to Euclid Collaboration (2024b), we consider shifts of more than 0.25σ significant, since then, if the marginalized posteriors were Gaussian, the contours with and without shift overlap less than 90%. We find that, when including the nuisance parameters, none of the cosmological parameters show a significant shift due to SLC. The largest shift occurs for σ₈, which is shifted by 0.19σ towards higher values. This likely occurs to compensate for the signal increase in higher tomographic bins due to SLC.

The MAP values are shifted by similar absolute amounts as the maxima of the marginalized posteriors. Furthermore, all shifts have the same sign. This indicates that SLC does not have a strong influence on the degeneracy direction between parameters and simply shifts the multidimensional posterior to new values, since otherwise the marginalized posteriors would be affected differently from the MAP.

While the cosmological parameters are relatively unaffected by SLC, the nuisance parameters A_IA and δz₄ are shifted significantly (0.39σ and 0.33σ, respectively) to higher values. These shifts compensate for two different redshift regimes. A larger A_IA decreases the cosmic shear signal, in particular at low redshifts, which is similar to the decrease in E_n for the lowest tomographic bin due to SLC. A larger δz₄ increases the signal for E_n including the fourth tomographic bin, similar to the increase in signal due to SLC for the higher tomographic bins. Therefore, these two nuisance parameters neutralize the effect of SLC on the cosmological parameters, even though their physical meaning is completely unrelated to SLC.

This raises the question of whether the cosmological analysis would still be insensitive to SLC if we excluded the nuisance parameters. This scenario is particularly interesting because future cosmic shear analyses are expected to use informative priors for effects such as intrinsic alignment, constrained from independent observations (e.g. Johnston et al. 2019; Fortuna et al. 2021). These priors are significantly tighter than the broad flat prior used here. To study whether SLC is still unimportant for the analyses when the nuisance parameters are not varied, we perform the second cosmological inference with A_IA and δz_i fixed to zero. The resulting maxima of the marginalized posteriors and the MAP values are given in Table 4.

The posteriors for S₈, σ₈, and Ω_m are shown in Fig. 5 as pink and cyan contours and the full posteriors are in Appendix E. Neglecting SLC would now lead to shifts of 0.42σ in S₈ and 0.31σ in Ω_m, with S₈ shifted to higher and Ω_m shifted to lower values. These increases in the effect of SLC are partly due to the decreased parameter uncertainties compared to the analysis with nuisance parameters, which increase the significance of the parameter shifts. However, the absolute values of the shifts also increase, for example S₈ shifts by 0.7% for fixed nuisance parameters and 0.4% for variable nuisance parameters. This confirms that SLC would bias the cosmological inference stronger if nuisance parameters were fixed.

Both with and without nuisance parameters, the MICE simulation’s true cosmological parameters are within the cosmological inference’s uncertainties. However, in both cases, the values for the unclustered data vector is closer to the true values. This is unsurprising, as the unclustered data vector corresponds to the model assumption of unclustered source galaxies. The remaining difference between the inferred parameters for the unclustered data vector and the true cosmology is likely due to the fact that the simulation is a single realization of the cosmology and the limited accuracy of the model.

We now investigate the impact of the different SLC effects discussed in Sect. 2.2. For this, we show in Fig. 6 the COSEBIs when different parts of the SLC are included, as well as the fractional difference to the signal for unclustered sources. The largest impact on the signal is due to the estimator bias (EB), which suppresses the signal for all tomographic bins. The IC and TC effects partially counteract each other. While the IC increases the signal, as expected from Eq. (17), the contamination by lower redshift galaxies decreases it. The scale dependence of the SLC effects is investigated using shear correlation functions in Appendix D.

Fig. 6. COSEBIs for the KiDS-like set-up when different SLC effects are included, either the full SLC (red crosses), the estimator bias (EB) and intrinsic clustering (IC; orange triangles), only the estimator bias (EB; green triangles), and no SLC (blue dots). Upper right: COSEBIs. Lower left: Fractional difference to COSEBIs for unclustered sources (i.e. without SLC).

To see how the individual SLC effects affect the cosmological parameter constraints, Fig. 7 shows the corresponding parameter constraints. For S₈ the strongest impact occurs due to EB, which lowers the estimated value. This is in line with the effect on the COSEBIs – lowering the cosmic shear signal is equivalent to lowering S₈. The IC slightly increases S₈ again, as expected since it increases the cosmic shear signal. Interestingly, S₈ increases even more, once the TC effect is included. This seems counterintuitive since TC generally lowers the shear signal. However, simultaneously to increasing S₈, the effect also increases the intrinsic alignment parameter A_IA, which lowers the signal. Since A_IA mainly causes a decrease in the signal at low tomographic bins, where the suppression due to TC is strongest, while S₈ mainly increases the signal at higher redshifts, these two parameters together can match the total effect by the contamination.

Fig. 7. Parameter constraints for the KiDS-like set-up, for a data vector including the full SLC, the estimator bias (EB) and intrinsic clustering (IC), only the estimator bias (EB), and without any SLC.

5.2. Euclid-like set-up

Finally, we consider the Euclid-like set-up. Figure 8 shows the measured COSEBIs for the clustered and unclustered set-up, as well as the differences for a subset of the tomographic bins. The full measurements for all tomographic bins are provided in Fig. E.8. The measurements for the clustered and original catalogue agree perfectly, which confirms that the creation of the two mock catalogues is consistent. The difference between the clustered and unclustered case is smaller for the Euclid-like set-up than for the KiDS-1000-like set-up. This is likely due to the narrower redshift distributions p_z(z), with shorter tails and fewer outliers. Consequently, there is less overlap between source galaxies in a tomographic bin and the matter structures lensing other galaxies in the same tomographic bin, which reduces the impact of SLC. Nevertheless, since the statistical uncertainty for Euclid is much smaller than for Stage III surveys such as KiDS, the difference between the two cases approaches 1σ for the highest tomographic bins, while it remains below 0.5σ for KiDS-1000.

Fig. 8. Measurement for the Euclid-like set-up. Upper right: COSEBI E-modes En measured for the original FS2 galaxies (black points), the clustered catalogue (red plusses) and the unclustered catalogue (blue crosses) with the Euclid uncertainty (grey band) for a subset of the 13 tomographic bins. Lower left: Difference between En for the original and clustered sources (black points) and between the clustered and unclustered sources (blue crosses), divided by Euclid uncertainty. The full measurement for all tomographic bins is in Fig. E.8.

Figure 9 shows the parameter constraints for the most important cosmological parameters (red and blue constraints), while the maxima of the marginalized posteriors and the MAP values for all parameters are given in Table 5. We find lower shifts in parameter estimates due to SLC than for the KiDS-like set-up, which is likely due to the narrower redshift bins. This provides an additional motivation for analysing Euclid with more than ten tomographic bins, as already suggested by Euclid Collaboration (2021a).

Fig. 9. Parameter constraints for the Euclid-like set-up using all parameters (lower corner) and with fixed nuisance parameters (upper corner), either for a data vector including SLC (dark blue and cyan) or without SLC (red and pink). Constraints for all sampled parameters can be found in Tables 5 and 6.

Similarly to the KiDS set-up, the nuisance parameters are more strongly affected than the cosmological ones. Several redshift bin shifts δz vary by more than 0.25σ between the cases with and without SLC. Consequently, these parameters again account for the SLC, while the more physical parameters remain stable.

The MAP values show similar absolute shifts due to SLC as the maxima of the marginalized posteriors, which also show the same sign. The shift in S₈ is slightly larger for the MAP than for the marginalized posterior (−2.5 × 10⁻³ compared to −2.0 × 10⁻³) while the shift in Ω_m is slightly smaller (−1.7 × 10⁻³ compared to −3.0 × 10⁻³). Overall, the similarity of the absolute shifts suggests that, similar to the KiDS-set-up, SLC mainly moves the full posterior to a different parameter set and does not strongly influence the degeneracy directions between parameters.

We investigate the effect of neglecting the nuisance parameters in Fig. 9 and Table 6 as pink and cyan contours. Here, we fix the δz and A_IA parameters to zero. In this set-up the shifts due to the unaccounted SLC increase, with the shift for σ₈ and Ω_c h² grazing the significance threshold of 0.25σ and the shift in w₀ becoming larger than this threshold. This increase in significance is partly due to the strongly decreased parameter uncertainties. However, for example for w₀, neglecting the nuisance parameters also increases the absolute value of the shift, from 0.02 to 0.04. Consequently, for Euclid, it is vital to either use nuisance parameters to account for unmodelled biases or to correctly include SLC in the model.

Finally, we also test the impact of the different SLC effects for the Euclid-like set-up. Figure 10 shows the obtained parameter constraints for S₈, σ₈, Ω_m and A_IA when excluding parts of the SLC effects. The effects qualitatively induce the same changes in the parameter constraints as for the KiDS-like set-up. The EB effect decreases S₈, while the IC and TC effects increase it. The intrinsic alignment parameter remains stable under the EB and IC effects but shifts to higher values for the TC effect. Again, the EB has the strongest impact on the S₈ parameter.

Fig. 10. Parameter constraints for Euclid-like set-up, for a data vector including the full SLC, the estimator bias (EB) and intrinsic clustering (IC), only the estimator bias (EB), and without any SLC.

6. Discussion

In this work we investigated the impact of the clustering of source galaxies with lensing matter structures, known as source-lens clustering (SLC), on cosmological parameter inference with second-order statistics of cosmic shear. We considered realistically clustered mock galaxies in the MICE and FS2 simulations, as well as uniformly distributed sources, and measured the COSEBIs for a KiDS-1000-like and a Euclid-like set-up. We then performed a cosmological inference for both clustered and unclustered sources using nested sampling of the posterior and including nuisance parameters for intrinsic alignments and shifts of the source redshift distributions. In summary, we found that if nuisance parameters are taken into account, SLC has only a minor impact on cosmological inference for Stage III, making it an unlikely candidate to resolve the ‘S₈-tension’.

For KiDS-1000, we found little impact of SLC on cosmological parameters if nuisance parameters are taken into account. The SLC causes a small decrease (within 0.5σ) in the cosmic shear signal for the lowest tomographic bin and a similarly small increase for the higher tomographic bins. These shifts are ‘absorbed’ by changes in the intrinsic alignment parameter A_IA and the redshift distribution biases δz, leading to cosmological parameters, such as S₈ and Ω_m, to be shifted by less than 0.2σ. Thus, it is valid for Stage III surveys to neglect SLC in the modelling as long as nuisance parameters for intrinsic alignments and shifts of the tomographic bins are included. We thus confirm the expectation outlined in Krause et al. (2021) that SLC has only a minor impact on cosmological inference with Stage III surveys and can be neglected in their modelling.

This picture changes if nuisance parameters are fixed and not varied with the analysis. In this case, we see that S₈ is shifted by 0.42σ and Ω_m by 0.31σ for the KiDS-1000-like set-up. Consequently, when neglecting nuisance parameters, SLC should be included in the modelling of the cosmic shear signal. Furthermore, we advise against interpreting parameters such as A_IA as ‘quantifying intrinsic alignment’. As shown here, effects completely unrelated to intrinsic alignments can affect this parameter substantially. A further example of this was found by Li et al. (2021), who found that the intrinsic alignment nuisance parameter is strongly correlated with the redshift uncertainty for low-redshift tomographic bins and fixing the assumed redshift uncertainty to zero could lead to an unphysically high inferred A_IA for blue galaxies. Furthermore, Leonard et al. (2024) found that the redshift scaling of the intrinsic alignment amplitude is strongly degenerate with the mean and variance of the redshift distribution. The interplay between redshift distribution uncertainties and intrinsic alignment constraints is further dependent on the considered intrinsic alignment model and can be affected by the considered priors (Fischbacher et al. 2023).

Different intrinsic alignment models might further complicate the interplay between intrinsic alignment parameters and SLC. Here, we only investigated the NLA model, but more complicated models, such as the Tidal Alignment – Tidal Torquing (TATT) model (Blazek et al. 2019), include more free parameters that can further compensate for the impact of SLC.

The fact that unmodelled systematic effects such as SLC change the meaning of the nuisance parameters can be an issue when deciding on priors for these parameters. Ideally, effects such as intrinsic alignment could be constrained from observations using physical models (e.g. Fortuna et al. 2021). These constraints would provide tight priors on the nuisance parameters. However, this is no longer straightforward if the parameters are meant to also capture other effects such as SLC. Consequently, if observational motivated priors for nuisance parameters are used, effects such as SLC need to be carefully included in the modelling.

The parameter uncertainties we find for the KiDS-1000-like set-up are of the same order as those found in Asgari et al. (2021) except for Ω_bh², whose uncertainty is more than twice the one from Asgari et al. (2021). This is due to the larger prior for this parameter we assume here, as Ω_b cannot be constrained from a KiDS-1000-like survey. Since cosmic shear is almost insensitive to the baryonic matter density Ω_b the constraints are completely prior-dominated. Our larger prior thus also leads to a larger posterior. In general, our confidence intervals are slightly larger than those by Asgari et al. (2021), which is likely because we are using the sampler nautilus to sample the posterior, while they use multinest (Feroz et al. 2019), which can underestimate posteriors (Lemos et al. 2023). Furthermore, we also do not vary the parameter A_b for baryonic feedback, which in the simulation is exactly zero.

For KiDS, we investigated the impact of different parts of the SLC effect. We found that the largest effect is due to the estimator bias, which occurs due to the neglect of the source galaxies angular correlation function and decreases the measured signal. The intrinsic clustering of the sources and the contamination of tomographic bins due to photometric redshift errors partly counteract each other, leaving the estimator bias as the most important contribution.

For the Euclid-like set-up, we see a similar effect as for the KiDS-1000-like set-up. However, the absolute shift in the signal due to SLC is smaller. This is likely due to the narrower p_z(z) with shorter tails. These lead to less scattering between tomographic redshift bins, thus reducing the tomographic bin contamination (TC) effect. The p_z(z) are still broader than the typical galaxy correlation length, which is of the order of 1 Mpc, so the intrinsic clustering remains.

We also find that the impact on the cosmological parameter analysis for Euclid, aside from the nuisance parameters, is small. S₈ is shifted only by 0.14σ and w₀ by 0.12σ if nuisance parameters are taken into account. This indicates that SLC might even be neglected for a Stage IV survey. Again, the picture would change if nuisance parameters were fixed in the analysis. Then, the shift in w₀ becomes significant (0.32σ) and shifts in σ₈ and Ω_m also increase.

This small effect of SLC on cosmological parameters is in apparent conflict with Euclid Collaboration (2024b), who found that neglecting SLC can cause biases of more than 1σ on Ω_m and σ₈. However, their analysis relied on modelling SLC analytically, while we directly measured the effect using realistically clustered sources. The analytic modelling naturally included simplifying assumptions, such as a modelled redshift distribution, linear galaxy bias and idealized photometric redshift uncertainties. Our analysis did not require these and could be seen as a more realistic extension of their analysis.

We note that our uncertainties on the cosmological parameters are larger than those quoted in Euclid Collaboration (2024b). This is likely due to three reasons. First, we perform a full sampling of the posterior instead of a Fisher analysis as in Euclid Collaboration (2024b). Fisher analyses can underestimate a posterior’s width, leading to tighter parameter constraints (Wolz et al. 2012). Second, we include more parameters in our analysis: the intrinsic alignment parameter A_IA and the p_z(z)-shifts for the tomographic bins. Including nuisance parameters reduces the constraining power on the cosmological parameters, leading to larger uncertainties. Third, we use a different observable, namely COSEBIs, instead of the lensing power spectrum. Consequently, our analysis uses other (and differently weighted) length scales, and the constraints are not directly comparable. As shown by Asgari et al. (2021), the COSEBIs used here are most sensitive to ℓ between 10 and 1500, while Euclid Collaboration (2024b) considered ℓ up to 5000. COSEBIs also provide a separation into E- and B-modes, while an analysis using the C(ℓ) directly, implicitly assumes vanishing B-modes, which leads to tighter parameter constraints.

We note that the photometric redshift estimates for our simulated Euclid-like sample were derived assuming the full depth of the Euclid calibration fields. These will not be fully available for the first Euclid data releases, so the photometric redshifts will show a larger scatter, increasing the TC effect. However, the statistical uncertainty for the first data release will also be significantly larger due to the smaller survey area, so a subtle effect such as SLC decreases in importance.

Another confounding factor to the SLC that has not been considered here is the impact of blending. Since galaxies have a finite extent on the sky, galaxies that are close by might not be detected as individual objects, which biases shear estimates (Euclid Collaboration: Congedo et al., in prep.). This bias is correlated with the distribution of sources, as stronger source clustering leads to more blended objects, leading to a further correlation between the estimated shear and the source distribution. Studying the impact of blending requires not only simulated sources with a realistic spatial distribution, provided here by the FS2 and MICE simulations but also a realistic (observed) size distribution obtained from image simulations, as created for example by Jansen et al. (2024). We leave the investigation of the impact of blending on SLC to future work.

Accounting for SLC in the theoretical modelling remains a challenge. The intrinsic clustering (IC) effect can be estimated from first principles, as discussed in Krause et al. (2021) and Appendix A. However, this is actually the smallest contributor to the total SLC. Modelling the estimator bias (EB) and TC effects is more involved, as they require realistic models of the (three-dimensional) source galaxy distribution n(χϑ, χ). In principle, this could be measured alongside the sources’ shears, but for photometric samples, the distance estimate will not be accurate enough. Instead, a better understanding of how the scatter in photometric redshift estimates impacts the number density of selected sources, conditioned on the overall galaxy number density, could help in modelling these effects more thoroughly.

Finally, we note that our analysis was concerned with second-order cosmic shear statistics. Euclid will achieve its tightest constraints on cosmological parameters from second-order statistics by combining the cosmic shear analysis with measurements of galaxy clustering and the cross-correlation of shear and galaxy positions. These constraints will be tighter than for cosmic shear alone (e.g., about a factor 2 in Ω_m for a flat ΛCDM setting, Euclid Collaboration 2020b). Consequently, the impact of SLC becomes more significant.

Furthermore, higher-order statistics, which could further improve the cosmological constraints, might be even stronger impacted by SLC. Recently, Gatti et al. (2024) investigated the impact of SLC on higher-order shear statistics and found that it can significantly bias measurements and cosmological inferences. However, this bias is caused partly by the choice of estimator for the higher-order statistics, which relies on shear or convergence maps. Since the SLC impacts the number of sources in each pixel of the map, it also affects the noise of the estimate. The correlation between the source number and the signal then introduces a correlation of the signal variance with the signal itself and, thus, an additional third-order correlation. The effect needs to be modelled for every statistic involving shear or convergence maps but does not occur for statistics measured directly from galaxy catalogues. Thus, we expect SLC to have a negligible effect on higher-order statistics obtained from correlation functions, such as third-order cosmic shear (Schneider et al. 2005; Heydenreich et al. 2023). We leave the proper investigation of SLC for third-order shear statistics measured from correlation functions to future work.

Acknowledgments

We thank the anonymous referee for their helpful suggestions. The Euclid Consortium acknowledges the European Space Agency and a number of agencies and institutes that have supported the development of Euclid, in particular the Agenzia Spaziale Italiana, the Austrian Forschungs förderungsgesellschaft funded through BMK, the Belgian Science Policy, the Canadian Euclid Consortium, the Deutsches Zentrum für Luft- und Raumfahrt, the DTU Space and the Niels Bohr Institute in Denmark, the French Centre National d’Etudes Spatiales, the Fundação para a Ciência e a Tecnologia, the Hungarian Academy of Sciences, the Ministerio de Ciencia, Innovación y Universidades, the National Aeronautics and Space Administration, the National Astronomical Observatory of Japan, the Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the Research Council of Finland, the Romanian Space Agency, the State Secretariat for Education, Research, and Innovation (SERI) at the Swiss Space Office (SSO), and the United Kingdom Space Agency. A complete and detailed list is available on the Euclid web site (www.euclid-ec.org). LL is supported by the Austrian Science Fund (FWF) [ESP 357-N]. The Innsbruck authors acknowledge support from the Austrian Research Promotion Agency (FFG) and the Federal Ministry of the Republic of Austria for Climate Action, Environment, Mobility, Innovation and Technology (BMK) via grants 899537 and 900565. SU acknowledges support from the Max Planck Society and the Alexander von Humboldt Foundation in the framework of the Max Planck-Humboldt Research Award endowed by the Federal Ministry of Education and Research. H. Hildebrandt is supported by a DFG Heisenberg grant (Hi 1495/5-1), the DFG Collaborative Research Center SFB1491, as well as an ERC Consolidator Grant (No. 770935). MA is supported by the UK Science and Technology Facilities Council (STFC) under grant numbers ST/Y002652/1 and the Royal Society under grant numbers RGSR2222268 and ICAR1231094. This work has made use of CosmoHub. CosmoHub has been developed by the Port d’Informació Científica (PIC), maintained through a collaboration of the Institut de Física d’Altes Energies (IFAE) and the Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT) and the Institute of Space Sciences (CSIC & IEEC), and was partially funded by the “Plan Estatal de Investigación Científica y Técnica y de Innovación” program of the Spanish government. Based on data products from observations made with ESO Telescopes at the La Silla Paranal Observatory under programme IDs 177.A-3016, 177.A-3017, 177.A-3018, 179.A-2004, 298.A-5015.

References

Appendix A: Additional term to the convergence power spectrum due to IC source-lens clustering

This appendix calculates the expected additional term ΔC to the lensing power spectrum C(ℓ) due to the IC part of source-lens clustering. As discussed in Sect. 2.2, the effect leads to an additional term Δκ to the lensing convergence. Consequently, the convergence correlation function ξ_κ(|θ₁ − θ₂|) = ⟨κ(θ₁) κ(θ₂)⟩ contains additional terms and is

$$\begin{array}{c}\hfill ⟨\kappa ({\mathit{\theta }}_1)\kappa ({\mathit{\theta }}_2)⟩=⟨{\kappa }_0({\mathit{\theta }}_1){\kappa }_0({\mathit{\theta }}_2)⟩+\underset{\mathcal{A}({\mathit{\theta }}_1,{\mathit{\theta }}_2)}{\underbrace{⟨{\kappa }_0({\mathit{\theta }}_1)\mathrm{\Delta }\kappa ({\mathit{\theta }}_2)}}+⟨\mathrm{\Delta }\kappa ({\mathit{\theta }}_1){\kappa }_0({\mathit{\theta }}_2)⟩+\underset{\mathcal{B}({\mathit{\theta }}_1,{\mathit{\theta }}_2)}{\underbrace{⟨\mathrm{\Delta }\kappa ({\mathit{\theta }}_1)\mathrm{\Delta }\kappa ({\mathit{\theta }}_2)⟩}}.\end{array}$$(A.1)

We are first considering 𝒜, which is

$$\begin{array}{c}\hfill \mathcal{A}({\mathit{\theta }}_1,{\mathit{\theta }}_2)=b{\displaystyle {\int }_0^{\infty }}\text{d}{\chi }_1{\displaystyle {\int }_0^{\infty }}\text{d}{\chi }_2{\displaystyle {\int }_{{\chi }_1}^{\infty }}\text{d}{\chi }_1^{\prime }{\displaystyle {\int }_{{\chi }_2}^{\infty }}\text{d}{\chi }_2^{\prime }n({\chi }_1^{\prime })n({\chi }_2^{\prime })W({\chi }_1,{\chi }_1^{\prime })W({\chi }_2,{\chi }_2^{\prime })⟨\delta ({\chi }_1{\mathit{\theta }}_1,{\chi }_1)\delta ({\chi }_2{\mathit{\theta }}_2,{\chi }_2)\delta ({\chi }_2^{\prime }{\mathit{\theta }}_2,{\chi }_2^{\prime })⟩.\end{array}$$(A.2)

This term vanishes under the Limber approximation. This can be seen qualitatively by considering the three-point correlation function. The correlation function will only give a significant contribution if the densities are evaluated at similar distances. Therefore, χ₁ ≃ χ₂ ≃ χ′₂. However, for χ₂ = χ′₂, the lensing efficiency W(χ₂, χ′₂) vanishes. Consequently, 𝒜 becomes zero.

To see this more clearly, we are Fourier transforming 𝒜 to $\stackrel{\sim }{\mathcal{A}}$, which is

$$\begin{array}{cc}\hfill \stackrel{\sim }{\mathcal{A}}({\mathit{\ell }}_1,{\mathit{\ell }}_2)& ={\displaystyle \int {\mathrm{d}}^2}{\theta }_1{\displaystyle \int {\mathrm{d}}^2}{\theta }_2{\mathrm{e}}^{-\mathrm{i}{\mathit{\theta }}_1·{\mathit{\ell }}_1-\mathrm{i}{\mathit{\theta }}_2·{\mathit{\ell }}_2}\mathcal{A}({\mathit{\theta }}_1,{\mathit{\theta }}_2)\hfill \\ \hfill & =b{\displaystyle {\int }_0^{\infty }}\text{d}{\chi }_1{\displaystyle {\int }_0^{\infty }}\text{d}{\chi }_2{\displaystyle {\int }_{{\chi }_1}^{\infty }}\text{d}{\chi }_1^{\prime }{\displaystyle {\int }_{{\chi }_2}^{\infty }}\text{d}{\chi }_2^{\prime }p({\chi }_1^{\prime })p({\chi }_2^{\prime })W({\chi }_1,{\chi }_1^{\prime })W({\chi }_2,{\chi }_2^{\prime })\hfill \\ \hfill & \times \int \frac{\text{d}{k}_{1z}}{2\pi }\int \frac{\text{d}{k}_{2z}}{2\pi }{\mathrm{e}}^{\mathrm{i}{k}_{1z}({\chi }_1-{\chi }_2^{\prime })}{\mathrm{e}}^{\mathrm{i}{k}_{2z}({\chi }_2-{\chi }_2^{\prime })}{(2\pi )}^2{\delta }_{\mathrm{D}}(\frac{{\mathit{\ell }}_1}{{\chi }_1}+\frac{{\mathit{\ell }}_2}{{\chi }_2}+\frac{{\mathit{\ell }}_2}{{\chi }_2^{\prime }})\hfill \\ \hfill & \times B[\sqrt{\frac{|{\mathit{\ell }}_1{|}^2}{{\chi }_1^2}+k_{1z}^2},\sqrt{\frac{|{\mathit{\ell }}_2{|}^2}{{\chi }_2^2}+k_{2z}^2},\sqrt{\frac{|{\mathit{\ell }}_2{|}^2}{{{\chi }_2^{\prime }}^2}+{(k_{1z}+k_{2z})}^2},{\chi }_1,{\chi }_2,{\chi }_2^{\prime }],\hfill \end{array}$$(A.3)

where B is the (three-dimensional) matter bispectrum. Under the Limber approximation, the bispectrum simplifies as

$$\begin{array}{c}\hfill B[\sqrt{\frac{|{\mathit{\ell }}_1{|}^2}{{\chi }_1^2}+k_{1z}^2},\sqrt{\frac{|{\mathit{\ell }}_2{|}^2}{{\chi }_2^2}+k_{2z}^2},\sqrt{\frac{|{\mathit{\ell }}_2{|}^2}{{{\chi }_2^{\prime }}^2}+{(k_{1z}+k_{2z})}^2},{\chi }_1,{\chi }_2,{\chi }_2^{\prime }]\to B(\frac{{\ell }_1}{{\chi }_1},\frac{{\ell }_2}{{\chi }_2},\frac{{\ell }_2}{{\chi }_2^{\prime }},{\chi }_1,{\chi }_2,{\chi }_2^{\prime }).\end{array}$$(A.4)

Then, $\stackrel{\sim }{\mathcal{A}}$ becomes

$$\begin{array}{cc}& \mathcal{A}({\mathit{\ell }}_1,{\mathit{\ell }}_2)\hfill \\ \hfill & ={(2\pi )}^{2}b{\displaystyle {\int }_0^{\infty }}\text{d}{\chi }_1{\displaystyle {\int }_0^{\infty }}\text{d}{\chi }_2{\displaystyle {\int }_{{\chi }_1}^{\infty }}\text{d}{\chi }_1^{\prime }{\displaystyle {\int }_{{\chi }_2}^{\infty }}\text{d}{\chi }_2^{\prime }p({\chi }_1^{\prime })p({\chi }_2^{\prime })W({\chi }_1,{\chi }_1^{\prime })W({\chi }_2,{\chi }_2^{\prime }){\delta }_{\mathrm{D}}({\chi }_1-{\chi }_2^{\prime }){\delta }_{\mathrm{D}}({\chi }_2-{\chi }_2^{\prime }){\delta }_{\mathrm{D}}(\frac{{\mathit{\ell }}_1}{{\chi }_1}+\frac{{\mathit{\ell }}_2}{{\chi }_2}+\frac{{\mathit{\ell }}_2}{{\chi }_2^{\prime }})\hfill \\ \hfill & \times B(\frac{{\ell }_1}{{\chi }_1},\frac{{\ell }_2}{{\chi }_2},\frac{{\ell }_2}{{\chi }_2^{\prime }},{\chi }_1,{\chi }_2,{\chi }_2^{\prime })\hfill \end{array}$$(A.5)

$$\begin{array}{cc}\hfill & ={(2\pi )}^{2}b{\displaystyle {\int }_0^{\infty }}\text{d}{\chi }_2{\displaystyle {\int }_{{\chi }_2}^{\infty }}\text{d}{\chi }_1^{\prime }p({\chi }_1^{\prime })p({\chi }_2)W({\chi }_2,{\chi }_1^{\prime })W({\chi }_2,{\chi }_2){\delta }_{\mathrm{D}}\left(\frac{{\mathit{\ell }}_1+2{\mathit{\ell }}_2}{{\chi }_2}\right)B(\frac{{\ell }_1}{{\chi }_2},\frac{{\ell }_2}{{\chi }_2},\frac{{\ell }_2}{{\chi }_2},{\chi }_2,{\chi }_2,{\chi }_2)\hfill \\ \hfill & =0,\hfill \end{array}$$(A.6)

where we used in the last step that W(χ₂, χ₂) = 0.

Consequently, the only term that can impact the convergence power spectrum is ℬ. This term is given by

$$\begin{array}{c}\hfill \mathcal{B}({\mathit{\theta }}_1,{\mathit{\theta }}_2)=b^2{\displaystyle {\int }_0^{\infty }}\text{d}{\chi }_1{\displaystyle {\int }_0^{\infty }}\text{d}{\chi }_2{\displaystyle {\int }_{{\chi }_1}^{\infty }}\text{d}{\chi }_1^{\prime }{\displaystyle {\int }_{{\chi }_2}^{\infty }}\text{d}{\chi }_2^{\prime }n({\chi }_1^{\prime })n({\chi }_2^{\prime })W({\chi }_1,{\chi }_1^{\prime })W({\chi }_2,{\chi }_2^{\prime })⟨\delta ({\chi }_1{\mathit{\theta }}_1,{\chi }_1)\delta ({\chi }_2{\mathit{\theta }}_2,{\chi }_2)\delta ({\chi }_1^{\prime }{\mathit{\theta }}_1,{\chi }_1^{\prime })\delta ({\chi }_2^{\prime }{\mathit{\theta }}_2,{\chi }_2^{\prime })⟩.\end{array}$$(A.7)

The four-point correlation function can be decomposed into its connected and unconnected part,

$$\begin{array}{cc}\hfill ⟨\delta ({\chi }_1{\mathit{\theta }}_1,{\chi }_1)\delta ({\chi }_2{\mathit{\theta }}_2,{\chi }_2)\delta ({\chi }_1^{\prime }{\mathit{\theta }}_1,{\chi }_1^{\prime }\delta ({\chi }_2^{\prime }{\mathit{\theta }}_2,{\chi }_2^{\prime })⟩& =⟨{\delta ({\chi }_1{\mathit{\theta }}_1,{\chi }_1)\delta ({\chi }_2{\mathit{\theta }}_2,{\chi }_2)\delta ({\chi }_1^{\prime }{\mathit{\theta }}_1,{\chi }_1^{\prime })\delta ({\chi }_2^{\prime }{\mathit{\theta }}_2,{\chi }_2^{\prime })⟩}_{\mathrm{c}}\hfill \\ \hfill & +⟨\delta ({\chi }_1{\mathit{\theta }}_1,{\chi }_1)\delta ({\chi }_2^{\prime }{\mathit{\theta }}_2,{\chi }_2^{\prime })⟨\delta ({\chi }_1^{\prime }{\mathit{\theta }}_1,{\chi }_1^{\prime })\delta ({\chi }_2{\mathit{\theta }}_2,{\chi }_2)⟩\hfill \\ \hfill & +⟨\delta ({\chi }_1{\mathit{\theta }}_1,{\chi }_1)\delta ({\chi }_1^{\prime }{\mathit{\theta }}_1,{\chi }_1^{\prime })⟨\delta ({\chi }_2{\mathit{\theta }}_2,{\chi }_2)\delta ({\chi }_2^{\prime }{\mathit{\theta }}_2,{\chi }_2^{\prime })⟩\hfill \\ \hfill & +⟨\delta ({\chi }_1{\mathit{\theta }}_1,{\chi }_1)\delta ({\chi }_2{\mathit{\theta }}_2,{\chi }_2)⟨\delta ({\chi }_1^{\prime }{\mathit{\theta }}_1,{\chi }_1^{\prime })\delta ({\chi }_2^{\prime }{\mathit{\theta }}_2,{\chi }_2^{\prime })⟩.\hfill \end{array}$$(A.8)

Using the same arguments as above, only the last summand can contribute to ℬ, since the lensing efficiency suppresses all correlations with χ_a ≃ χ′_a. Consequently,

$$\begin{array}{c}\hfill \mathcal{B}({\mathit{\theta }}_1,{\mathit{\theta }}_2)=b^2{\displaystyle {\int }_0^{\infty }\text{d}{\chi }_1{\int }_0^{\infty }\text{d}{\chi }_2{\int }_{{\chi }_1}^{\infty }}\text{d}{\chi }_1^{\prime }{\displaystyle {\int }_{{\chi }_2}^{\infty }}\text{d}{\chi }_2^{\prime }n({\chi }_1^{\prime })n({\chi }_2^{\prime })W({\chi }_1,{\chi }_1^{\prime })W({\chi }_2,{\chi }_2^{\prime })⟨\delta ({\chi }_1{\mathit{\theta }}_1,{\chi }_1)\delta ({\chi }_2{\mathit{\theta }}_2,{\chi }_2)⟩⟨\delta ({\chi }_1^{\prime }{\mathit{\theta }}_1,{\chi }_1^{\prime })\delta ({\chi }_2^{\prime }{\mathit{\theta }}_2,{\chi }_2^{\prime })⟩.\end{array}$$(A.9)

This is the full impact of source-lens clustering on the convergence correlation function under the assumptions of the Limber approximation.

To find the impact ΔC(ℓ) on the convergence power spectrum C(ℓ), we Fourier-transform ℬ, so

$$\begin{array}{cc}\hfill \mathrm{\Delta }C({\mathit{\ell }}_1){(2\pi )}^2{\delta }_{\mathrm{D}}({\mathit{\ell }}_1+{\mathit{\ell }}_2)& ={\displaystyle \int {\mathrm{d}}^2}{\theta }_1{\displaystyle \int {\mathrm{d}}^2}{\theta }_2{\mathrm{e}}^{-\mathrm{i}{\mathit{\theta }}_1·{\mathit{\ell }}_1-\mathrm{i}{\mathit{\theta }}_2·{\mathit{\ell }}_2}\mathcal{B}({\mathit{\theta }}_1,{\mathit{\theta }}_2)\hfill \end{array}$$(A.10)

$$\begin{array}{cc}\hfill & =b^2{\displaystyle {\int }_0^{\infty }}\text{d}{\chi }_1{\displaystyle {\int }_0^{\infty }}\text{d}{\chi }_2{\displaystyle {\int }_{{\chi }_1}^{\infty }}\text{d}{\chi }_1^{\prime }{\displaystyle {\int }_{{\chi }_2}^{\infty }}\text{d}{\chi }_2^{\prime }n({\chi }_1^{\prime })n({\chi }_2^{\prime })W({\chi }_1,{\chi }_1^{\prime })W({\chi }_2,{\chi }_2^{\prime })\hfill \\ \hfill & \times {\displaystyle \int \frac{{\mathrm{d}}^2{k}_{\perp }}{{(2\pi )}^2}}{\displaystyle \int \frac{{\mathrm{d}}^2{k}_{\perp }^{\prime }}{{(2\pi )}^2}\int \frac{\text{d}{k}_z}{2\pi }\int \frac{\text{d}{k}_z^{\prime }}{2\pi }}P[|({\mathit{k}}_{\perp },k_z)|]P[|({\mathit{k}}_{\perp }^{\prime },k_z^{\prime })|]{\mathrm{e}}^{\mathrm{i}{k}_z({\chi }_1-{\chi }_2)+\mathrm{i}{k}_z^{\prime }({\chi }_1^{\prime }-{\chi }_2^{\prime })}\hfill \\ \hfill & \times {(2\pi )}^4{\delta }_{\mathrm{D}}({\mathit{k}}_{\perp }{\chi }_1+{\mathit{k}}_{\perp }^{\prime }{\chi }_1^{\prime }-{\mathit{\ell }}_1){\delta }_{\mathrm{D}}(-{\mathit{k}}_{\perp }{\chi }_2-{\mathit{k}}_{\perp }^{\prime }{\chi }_2^{\prime }-{\mathit{\ell }}_2).\hfill \end{array}$$(A.11)

Under the Limber approximation, we replace P[|(k_⊥,k_z)|] by P(k_⊥), so we can execute the k_z integrals, leading to

$$\begin{array}{cc}\hfill \mathrm{\Delta }C({\mathit{\ell }}_1){(2\pi )}^2{\delta }_{\mathrm{D}}({\mathit{\ell }}_1+{\mathit{\ell }}_2)& =b^2{\displaystyle {\int }_0^{\infty }}\text{d}{\chi }_1{\displaystyle {\int }_0^{\infty }}\text{d}{\chi }_2{\displaystyle {\int }_{{\chi }_1}^{\infty }}\text{d}{\chi }_1^{\prime }{\displaystyle {\int }_{{\chi }_2}^{\infty }}\text{d}{\chi }_2^{\prime }n({\chi }_1^{\prime })n({\chi }_2^{\prime })W({\chi }_1,{\chi }_1^{\prime })W({\chi }_2,{\chi }_2^{\prime }){\delta }_{\mathrm{D}}({\chi }_1-{\chi }_2){\delta }_{\mathrm{D}}({\chi }_1^{\prime }-{\chi }_2^{\prime })\hfill \end{array}$$(A.12)

$$\begin{array}{cc}\hfill & \times {\displaystyle \int \frac{{\mathrm{d}}^2{k}_{\perp }}{{(2\pi )}^2}}{\displaystyle \int \frac{{\mathrm{d}}^2{k}_{\perp }^{\prime }}{{(2\pi )}^2}}P(k_{\perp })P(k_{\perp }^{\prime }){(2\pi )}^4{\delta }_{\mathrm{D}}({\mathit{k}}_{\perp }{\chi }_1+{\mathit{k}}_{\perp }^{\prime }{\chi }_1^{\prime }-{\mathit{\ell }}_1){\delta }_{\mathrm{D}}(-{\mathit{k}}_{\perp }{\chi }_2-{\mathit{k}}_{\perp }^{\prime }{\chi }_2^{\prime }-{\mathit{\ell }}_2)\hfill \\ \hfill & =b^2{\displaystyle {\int }_0^{\infty }}\text{d}\chi {\displaystyle {\int }_{\chi }^{\infty }}\text{d}{\chi }^{\prime }{n}_{\mathrm{s}}^2({\chi }^{\prime })W^2(\chi ,{\chi }^{\prime })\frac{1}{{\chi }^{\prime 2}}{\displaystyle \int \frac{{\mathrm{d}}^2{k}_{\perp }}{{(2\pi )}^2}}P(k_{\perp })P\left(\frac{|{\mathit{\ell }}_1-{\mathit{k}}_{\perp }\chi |}{{\chi }^{\prime }}\right){(2\pi )}^2{\delta }_{\mathrm{D}}(-{\mathit{\ell }}_1-{\mathit{\ell }}_2)\hfill \end{array}$$(A.13)

$$\begin{array}{cc}\hfill & ={(2\pi )}^2{\delta }_{\mathrm{D}}({\mathit{\ell }}_1+{\mathit{\ell }}_2)b^2{\displaystyle {\int }_0^{\infty }}\text{d}\chi {\displaystyle {\int }_{\chi }^{\infty }}\text{d}{\chi }^{\prime }{\displaystyle \int \frac{{\mathrm{d}}^{2}L}{{(2\pi )}^2}}\frac{n_{\mathrm{s}}^2({\chi }^{\prime })}{{\chi }^{\prime 2}}\frac{W^2(\chi ,{\chi }^{\prime })}{{\chi }^2}P\left(\frac{L}{\chi }\right)P\left(\frac{|{\mathit{\ell }}_1-\mathit{L}|}{\chi }\right),\hfill \end{array}$$(A.14)

where L = k_⊥ χ. Comparing this expression to Equations (47) and (48) of Krause et al. (2021) shows that

$$\begin{array}{c}\hfill \mathrm{\Delta }C(\mathit{\ell })=\mathrm{\Delta }{C}_{\mathit{EE}}(\mathit{\ell })+\mathrm{\Delta }{C}_{\mathit{BB}}(\mathit{\ell }),\end{array}$$(A.15)