© The Authors 2024

1. Introduction

While the Lambda cold dark matter (ΛCDM) model, is currently the best-in-class framework for parameterising the Universe, several open questions remain. A key component that is yet to be fully explained is the acceleration of the expansion of the Universe and its proposed driver: dark energy. A powerful tool for such studies is cosmic shear, namely, a distortion of the ellipticities that we observe for distant galaxies as a result of weak gravitational lensing from the large-scale structure of the Universe (LSS; see e.g. Albrecht et al. 2006).

To date, the most recent generation of cosmic shear surveys (Hikage et al. 2019; Asgari et al. 2021; Abbott et al. 2022) has been able to achieve high-precision cosmology competitive with cosmic microwave background experiments, for a combination of σ₈ and Ω_m (Planck Collaboration VI 2020). Now, upcoming Stage IV surveys (Albrecht et al. 2006) will probe a greater area and depth than previously possible. For example, telescopes such as Euclid¹ (Laureijs et al. 2011), Nancy Grace Roman² (Akeson et al. 2019), and the Vera C. Rubin Observatory³ (LSST Science Collaboration 2009) will achieve more than an order-of-magnitude increase in precision over existing surveys (Euclid Collaboration 2020, heareafter, EC20). We must therefore ensure that any sources of bias in our theoretical formalism are properly accounted for.

In this work, we consider the common approximations made when modelling the cosmic shear angular power spectrum in aggregate. This is of particular importance when deriving the cosmological parameters from a shear-only analysis, but also in a 3×2pt analyses (that includes shear and photometric galaxy clustering analyses, along with their cross-correlation), where modelling the weak-lensing power spectrum sufficiently is also essential. Throughout the literature, the effects considered in this paper have been studied independently, using varying survey and parameter specifications (see references in Table 1). Here, we evaluate them in a consistent framework and quantify their cumulative impact on cosmology inferred from Euclid’s weak lensing probe. As a first step, we review the literature and pinpoint which terms are potentially significant, as well as those for which the impact on the shear power spectrum has not been evaluated. A comparison of the typical magnitudes of the studied corrections is given in Table 1. To create Table 1 we manually read the published numbers from graphs in the referenced papers for the auto-correlation cosmic shear power spectrum for a redshift bin closest to z = 1 for a in each paper (auto-correlation refers to the inter-redshift bin correlation). These values are only approximate due to the inherent inaccuracy of reading value from graphs and varying assumptions the papers⁴.

From our survey of the literature we identify the following as potentially significant systematic effects requiring full analysis, these are as follows. Reduced shear approximation: the effect of assuming that the measured two-point statistics of reduced shear are equal to those of the shear field (Shapiro 2009; Krause & Hirata 2010); it has been previously shown that making this approximation would not lead to precise cosmological parameter estimation using cosmic shear for Euclid (Deshpande et al. 2020a). Magnification bias: the change in the observed number density of sources, due to galaxies at the flux limit of the survey having their flux increased or decreased due to magnification via lensing (Turner et al. 1984); this effect has also shown to significantly bias cosmological information from Euclid (Deshpande et al. 2020a; Duncan et al. 2022) if not accounted for. Additionally, magnification bias must also be accounted for in probes of galaxy clustering. Its impact on the Euclid galaxy clustering probe is discussed in Euclid Collaboration (2022b). Source-lens clustering: the intrinsic clustering of source galaxies correlated with the density field (Bernardeau 1998; Hamana et al. 2002; Yu et al. 2015). Typically, it is assumed source galaxies are distributed homogeneously across the sky. Source obscuration: A reduction in the observed galaxy distribution due to closely-spaced and blended or overlapping source galaxies (Hartlap et al. 2011). Local Universe effects: a possible bias in our measurements of summary statistics of the LSS due to residing in a region with a higher-than-average density (Reischke et al. 2019; Hall 2020). Flat Universe assumption: the impact of assuming that non-flat geometries are sufficiently well represented by modifying the expression for comoving distance and neglecting the additional change in the lensing kernels used to calculate the shear power spectrum (Taylor et al. 2018b).

We made this determination by first excluding the terms which are fourth-order in the lens potential or higher, as these have consistently been shown to be sub-dominant (Cooray & Hu 2002; Shapiro & Cooray 2006). Among these are time delay-lens coupling and deflection-deflection coupling (Bernardeau et al. 2010), which result from foregoing the small-angle and thin-lens approximations and solving the Sachs equation explicitly. Similarly, fourth-order correction terms resulting from relaxing the Born approximation and accounting for line-of-sight coupling of two foreground lenses are negligible (Shapiro & Cooray 2006). Additionally, fourth-order reduced shear corrections (Krause & Hirata 2010) were also neglected. Re-enforcing the sub-dominance of these terms is the fact that the standard reduced shear correction matches forward models and N-body simulations sufficiently well (Dodelson et al. 2006; Deshpande et al. 2020a). We further excluded fourth-order and higher terms resulting from the contribution of dark energy pressure to the lensing potential (Simpson et al. 2010).

Furthermore, we neglected finite beam corrections (Fleury et al. 2017, 2019), which are manifested as a fractional correction to the lensing power spectra on very small scales. This amounts to approximately −(1/3)(ℓθ)², where ℓ is the angular multipole of the power spectrum and θ the mean angular size of galaxies (this was verified by Breton & Fleury 2021 using ray-tracing simulations); for Euclid, this results in a fractional correction of −(ℓ/1.2 × 10⁶)² to the power spectrum.

We also neglected the effects of spatially-varying survey depth (Heydenreich et al. 2020), which can be accounted for directly in the covariance matrix through forward-modelling (Loureiro et al. 2022); although we note that demonstrating that this is sufficient for Euclid-like surveys requires further investigation. This choice could be relaxed in future work. We also note here that one could mitigate the impact of any systematics (including all of the effects discussed in this paper) by adding a corresponding uncertainty in the covariance, at the expense of the accuracy on cosmological parameter constraints. While this is not an approach we take in this paper, it could be investigated in future work.

Of the remaining effects, we neglected those that are significantly smaller than the reduced shear correction. This serves as a good comparison because this correction has been consistently demonstrated to produce biases close to the significance threshold (Shapiro 2009; Deshpande et al. 2020a). Accordingly, we neglected the impact of the Doppler-shift of galaxies on their two-point statistics. Due to the inhomogeneity of the Universe, galaxies have peculiar velocities that affect the measurement of their redshifts. If this is taken into account, it results in an additional contribution to the reduced shear (Bernardeau et al. 2010). As can be seen in Table 1, this effect can be safely neglected as it is two or three orders of magnitude below the reduced shear correction. Additionally, the cosmological parameter biases resulting from it are two orders of magnitude below the level of concern (Deshpande & Kitching 2021).

Likewise, the effect of unequal-time correlators was neglected, as the resulting correction to the angular power spectrum is more than four orders of magnitude smaller than that for the reduced shear (Kitching & Heavens 2017); as illustrated in Table 1. This correction is a consequence of relaxing the equal-time approximation, which approximates the cross-correlation matter power spectrum evaluated at different times as either the power spectrum at a fixed time, or by a geometric mean.

We also did not need to explicitly evaluate the effects of relaxing the Limber, and flat-sky approximations, as the combined corrections for these are an order-of-magnitude smaller than the reduced shear correction over the majority of the range of scales observed by Euclid (Kitching et al. 2017). This is again demonstrated in Table 1. The Limber approximation considers only correlations in the plane of the sky as contributing to the lensing signal, and projects others onto the plane of the sky by replacing spherical Bessel functions with Delta functions (Limber 1953; Kaiser 1998; LoVerde & Afshordi 2008).

Additionally, the Limber approximation is also employed when computing higher order corrections to the angular power spectrum, such as the reduced shear correction. In this work, we did not relax the use of the Limber approximation, as the cosmological parameter biases from this are safely negligible (Deshpande & Kitching 2020).

Another series of corrections that we deemed “safely negligible” stem from corrections to the theoretical expressions describing light propagation (Cuesta-Lazaro et al. 2018). Among these is the effect of second-order corrections to the effective speed of light. This relaxes the assumption that, as lensing potentials are small, the lensing effect can be studied in an effective Minkowskian spacetime and, accordingly, the effective speed of light need only be computed to the first-order. Similarly, a second correction to the effective speed of light presents itself from the energy-momentum tensor. Typically, this quantity is calculated under the assumption that lenses are moving slowly, so that the kinetic contribution to gravity can be ignored. Addressing this creates another correction to the angular power spectrum.

A further effect is that the observed ellipticity is non-linearly related to the shear, but that a linear approximation is often made. For a discussion of the impact of the non-linear relation on the observed shear distribution see Viola et al. (2014). The impact of this effect on the power spectrum (if we assume a linear relation instead of the non-linear relation) is investigated in Krause & Hirata (2010) (Sect. 3.3), who find that the impact is three orders of magnitude smaller than the shear power spectrum⁵.

The temporal-Born approximation is another correction to the description of light propagation. While the correction for the standard Born approximation accounts for the spatial discrepancy between the true perturbed path of a photon from source to observer compared to the mathematically convenient straight one, this discrepancy also produces a temporal one. The photon on the perturbed path will at times be ahead of the photon on the idealised path, and at other times lag behind. Accordingly, the two would encounter different evolutionary stages of the LSS at different times, necessitating a correction in the two-point statistics. The remaining two corrections that fall under this umbrella are the corrections of the Sachs-Wolfe and integrated Sachs-Wolfe effects. The former describes the redshift of an emitted photon due to the source galaxy’s gravitational potential, while the latter encodes the effect on the photon of interaction with the evolving gravitational potential along its path. All of these light propagation corrections are many orders of magnitude below the reduced shear correction, as can be seen in a representative example in Table 1.

Finally, we deemed the flexion corrections to be negligible without requiring explicit calculation. This additional correction term arises from the fact that, for larger sources, the image distortion consists of both shear and a higher-order component labelled flexion (Schneider & Er 2008). This term should be negligible because its effect on the cosmic shear signal will be dependent on third-order or higher-order brightness moments.

Here, we did not consider biases arising from general modelling of unknown shape measurement systematic effects (i.e. multiplicative and additive biases). Instead, we focused only on these well defined theoretical assumptions. For more details on shape measurement effects, we refer to Kitching et al. (2019, 2020, 2021), Kitching & Deshpande (2022). Additionally, we did not evaluate the additional selection effects of flux cuts and size cuts, as the former of these can be calibrated from deep fields, while size cuts are primarily a concern for ground-based telescopes, rather than space-based ones.

In this work, we also did not examine the impact of neglecting effects that are already well-established as requiring evaluation, namely, photometric redshift uncertainties, intrinsic alignments (IA) modelling, baryonic feedback, and modelling of the non-linear component of the matter power spectrum. Determining the exact specification for these is outside of the scope of this work, and each of those effects requires its own through investigation.

This work is structured as follows. In Sect. 2, we detail the theoretical formalism used. We review the basic, first-order cosmic shear angular power spectrum calculation. The expressions for the six correction terms of interest are also detailed, along with the Fisher matrix formalism used to predict cosmological parameter constraints and biases. In Sect. 3, we discuss the modelling and computational specifics used in this work. Finally, we discuss our results in Sect. 4. We show the cosmological parameter biases that result from neglecting the studied corrections, and discuss their implications for Euclid.

2. Theoretical formalism

Here, we begin by reviewing the standard first-order calculation of the cosmic shear angular power spectrum. Additional contributions to the lensing signal resulting from IAs and shot-noise are then described. We then detail the analytical forms of the six corrections requiring full evaluation: reduced shear, magnification bias, source-lens clustering, source obscuration, local Universe effects, and the flat Universe assumption. Finally, we review the Fisher matrix formalism used to predict cosmological parameter constraints and biases.

2.1. First-order cosmic shear calculation

As a consequence of weak gravitational lensing by the LSS, the observed ellipticity of distant galaxies is distorted. This change is dependent on the reduced shear, g, according to:

$$\begin{array}{c}\hfill g^{\alpha }(\mathit{\theta })=\frac{{\gamma }^{\alpha }(\mathit{\theta })}{1-\kappa (\mathit{\theta })},\end{array}$$(1)

where θ is the position of the galaxy on the sky, γ is the spin-2 shear with index α ∈ {1, 2} that describes the anisotropic stretching that turns circular distributions of light elliptical, and κ is the convergence – responsible for the isotropic change in the size of the image. Since in the weak lensing regime |κ|≪1, it is standard practice to make the reduced shear approximation, whereby

$$\begin{array}{c}\hfill g^{\alpha }(\mathit{\theta })\approx {\gamma }^{\alpha }(\mathit{\theta }).\end{array}$$(2)

Additionally, the convergence is a projection of the density contrast of the Universe, δ, along the line-of-sight over comoving distance, χ, to the comoving distance to the horizon, χ_h. For a particular tomographic redshift bin i, it is mathematically described by:

$$\begin{array}{c}\hfill {\kappa }_i(\mathit{\theta })={\displaystyle {\int }_0^{{\chi }_{\mathrm{h}}}\mathrm{d}\chi \delta [S_K(\chi )\mathit{\theta },\chi ]W_i(\chi ),}\end{array}$$(3)

where S_K is a function that encodes the effect of the curvature of the Universe, K, on comoving distances according to

$$\begin{array}{c}\hfill S_K(\chi )=\{\begin{array}{cc}{|K|}^{-1/2}{sin(|K|}^{-1/2}\chi )\hfill & K>0(\mathrm{closed}\mathrm{Universe}),\hfill \\ \chi \hfill & K=0(\mathrm{flat}\mathrm{Universe}),\hfill \\ {|K|}^{-1/2}{sinh(|K|}^{-1/2}\chi )\hfill & K<0(\mathrm{open}\mathrm{Universe}).\hfill \end{array}\end{array}$$(4)

We recall that for the quantity δ[S_K(χ)θ, χ] in Eq. (3) the second χ indicates not only that there is an evaluation at a comoving radius, χ, but also at a conformal time, η = η₀ − χ, meaning that all the integration over χ in this the paper are performed down the background light cone.

The W_i(χ) in Eq. (3) is the lensing projection kernel for tomographic bin i. It takes the form

$$\begin{array}{cc}\hfill W_i(\chi )& =\frac{3}{2}{\mathrm{\Omega }}_{\mathrm{m}}\frac{H_0^2}{c^2}\frac{S_K(\chi )}{a(\chi )}{\displaystyle {\int }_{\chi }^{{\chi }_{\mathrm{h}}}\mathrm{d}{\chi }^{\prime }{n}_i({\chi }^{\prime })\frac{S_K({\chi }^{\prime }-\chi )}{S_K({\chi }^{\prime })},}\hfill \end{array}$$(5)

which is dependent on the dimensionless present-day matter density of the Universe, Ω_m, the speed of light in a vacuum, c, the Hubble constant, H₀, the scale factor of the Universe, a(χ), and the probability distribution of galaxies within the redshift bin, i, namely, n_i(χ).

The spin-2 shear is directly related to the convergence in spherical-harmonic space. For a specified lensing mass distribution, assuming the flat-sky and prefactor-unity approximations (Kitching et al. 2017), and under the small-angle limit, this relationship takes the form:

$$\begin{array}{c}\hfill {\stackrel{\sim }{\gamma }}_i^{\alpha }(\mathit{\ell })=T^{\alpha }(\mathit{\ell }){\stackrel{\sim }{\kappa }}_i(\mathit{\ell }),\end{array}$$(6)

where ℓ is the spherical-harmonic conjugate of θ, with magnitude ℓ and angular component ϕ_ℓ. The functions T^α are two trigonometric weighting functions corresponding to each of the shear components. These take the form:

$$\begin{array}{cc}& T^1(\mathit{\ell })=cos(2{\varphi }_{\ell }),\hfill \end{array}$$(7)

$$\begin{array}{cc}& T^2(\mathit{\ell })=sin(2{\varphi }_{\ell }).\hfill \end{array}$$(8)

In the case of an arbitrary shear field, for example, a field constructed from data, two linear combinations of the individual shear components are pertinent. Specifically, these are a divergence-free B-mode and a curl-free E-mode:

$$\begin{array}{cc}& {\stackrel{\sim }{E}}_i(\mathit{\ell })={\displaystyle \sum _{\alpha }{T}^{\alpha }(\mathit{\ell }){\stackrel{\sim }{\gamma }}_i^{\alpha }(\mathit{\ell }),}\hfill \end{array}$$(9)

$$\begin{array}{cc}& {\stackrel{\sim }{B}}_i(\mathit{\ell })={\displaystyle \sum _{\alpha }\sum _{\beta }{\epsilon }^{\alpha \beta }{T}^{\alpha }(\mathit{\ell }){\stackrel{\sim }{\gamma }}_i^{\beta }(\mathit{\ell }).}\hfill \end{array}$$(10)

Here, the summations are over the shear components, and ε^αβ is the Levi-Civita symbol in the 2D case; such that: ε¹¹ = ε²² = 0 and ε¹² = −ε²¹ = 1.

Assuming that higher order systematic effects in the data have been accounted for, the B-mode of Eq. (10) vanishes. For the remaining E-mode, observables of interest are defined in the form of angular auto and cross-correlation power spectra, $C_{\ell ;ij}^{\gamma \gamma }$, such that:

$$\begin{array}{c}\hfill \langle {\stackrel{\sim }{E}}_i(\mathit{\ell }){\stackrel{\sim }{E}}_j({\mathit{\ell }}^{\prime })\rangle ={(2\pi )}^2{\delta }_{\mathrm{D}}^{(2)}(\mathit{\ell }+{\mathit{\ell }}^{\prime })C_{\ell ;ij}^{\gamma \gamma },\end{array}$$(11)

where the angular brackets on the left-hand-side denote the ensemble average, which (under the assumption of ergodicity) becomes a spatial average, and ${\delta }_{\mathrm{D}}^{(2)}$ is the Dirac delta for 2D. Under the extended Limber approximation (LoVerde & Afshordi 2008), where k = (ℓ+1/2)/S_K(χ), the power spectra themselves are further defined as:

$$\begin{array}{c}\hfill C_{\ell ;ij}^{\gamma \gamma }={\displaystyle {\int }_0^{{\chi }_{\mathrm{h}}}\mathrm{d}\chi \frac{W_i(\chi )W_j(\chi )}{S_K^2(\chi )}{P}_{\delta \delta }(k,\chi ),}\end{array}$$(12)

where P_δδ is the 3D matter power spectrum, and k is the magnitude of the spatial momentum vector, k, which also shares the angular component, ϕ_ℓ. Detailed reviews of this standard calculation can be found in Kilbinger (2015), Munshi et al. (2008), Bartelmann & Schneider (2001).

2.2. Intrinsic alignments and shot noise

When the angular power spectra are actually measured from surveys of galaxies, they contain non-lensing signals together with the pertinent cosmic shear power spectra. It is necessary to model each of these components to ensure accurate cosmological inference. A key non-lensing contribution arises from the fact that galaxies forming close to each other are forming in a similar tidal environment. Consequently, they have intrinsically correlated alignments (Joachimi et al. 2015; Kirk et al. 2015; Kiessling et al. 2015).

The observed ellipticity of an individual source, ϵ, can then, to first-order, be written as a combination of its underlying ellipticity in the absence of any cosmic shear or IA, ϵ^s, the cosmic shear, γ = γ¹ + iγ², and the effect of IA, ϵ^I according to:

$$\begin{array}{c}\hfill ϵ={ϵ}^{\mathrm{s}}+\gamma +{ϵ}^{\mathrm{I}}.\end{array}$$(13)

The angular power spectra corresponding to this observed ellipticity, $C_{\ell ;ij}^{ϵϵ}$, are then the sum of contributions resulting from its components,

$$\begin{array}{c}\hfill C_{\ell ;ij}^{ϵϵ}=C_{\ell ;ij}^{\gamma \gamma }+C_{\ell ;ij}^{\gamma \mathrm{I}}+C_{\ell ;ij}^{\mathrm{I}\gamma }+C_{\ell ;ij}^{\mathrm{II}}+N_{\ell ;ij}^{ϵ},\end{array}$$(14)

in which $C_{\ell ;ij}^{\gamma \gamma }$ are the cosmic shear angular power spectra defined in Eq. (12); in all cases the notation denotes z_i ≤ z_j. The $C_{\ell ;ij}^{\gamma \mathrm{I}}$ are the angular power spectra of correlations between foreground shear and background IA, which are only non-zero if photometric redshift estimates result in the scattering of observed redshifts between bins. On the other hand, the $C_{\ell ;ij}^{\mathrm{I}\gamma }$ arise from the correlation between background shear and foreground IA, and the $C_{\ell ;ij}^{\mathrm{II}}$ represent the auto-correlation of the IA; both must be accounted for. To accomplish this, the non-linear alignment (NLA) model (Bridle & King 2007) can be employed. Under this model, these IA spectra take the form:

$$\begin{array}{cc}& C_{\ell ;ij}^{\mathrm{I}\gamma }={\displaystyle {\int }_0^{{\chi }_{\mathrm{h}}}\frac{\mathrm{d}\chi }{S_K^2(\chi )}[W_i(\chi )n_j(\chi )+n_i(\chi )W_j(\chi )]P_{\delta \mathrm{I}}(k,\chi ),}\hfill \end{array}$$(15)

$$\begin{array}{cc}& C_{\ell ;ij}^{\mathrm{II}}={\displaystyle {\int }_0^{{\chi }_{\mathrm{h}}}\frac{\mathrm{d}\chi }{S_K^2(\chi )}{n}_i(\chi )n_j(\chi )P_{\mathrm{II}}(k,\chi ),}\hfill \end{array}$$(16)

which, in a similar manner to the shear power spectra, are projections of 3D IA power spectra, P_δI and P_II. Both of these are related to the matter power spectrum as follows:

$$\begin{array}{cc}& P_{\delta \mathrm{I}}(k,\chi )=[-\frac{{\mathcal{A}}_{\mathrm{IA}}{\mathcal{C}}_{\mathrm{IA}}{\mathrm{\Omega }}_{\mathrm{m}}}{D(\chi )}]P_{\delta \delta }(k,\chi ),\hfill \end{array}$$(17)

$$\begin{array}{cc}& P_{\mathrm{II}}(k,\chi )=[-\frac{{\mathcal{A}}_{\mathrm{IA}}{\mathcal{C}}_{\mathrm{IA}}{\mathrm{\Omega }}_{\mathrm{m}}}{D(\chi )}{]}^2{P}_{\delta \delta }(k,\chi ),\hfill \end{array}$$(18)

where the product of 𝒜_IA and 𝒞_IA is a free parameter typically set by fitting to simulations or data, and D(χ) is the density perturbation growth factor. We note that the NLA model is a limited description of IAs and, accordingly, it has its own associated modelling uncertainties. Extensions of this model have been proposed (see e.g. Fortuna et al. 2021, EC20). However, investigating the modelling of IAs for Euclid in detail is beyond the scope of this work and merits a separate future investigation.

Finally, the $N_{\mathcal{l};ij}^{ϵ}$ in Eq. (14) is the shot noise term arising from the zero-lag autocorrelation of the unlensed, uncorrelated source ellipticity ϵ^s in Eq. (13) (see e.g. Hu 1999 Eq. (4)). For a survey with equi-populated tomographic redshift bins, such as Euclid (EC20), this is expressed by:

$$\begin{array}{c}\hfill N_{\ell ;ij}^{ϵ}=\frac{{\sigma }_{ϵ}^2}{{\overline{n}}_{\mathrm{g}}/N_{\mathrm{bin}}}{\delta }_{\mathit{ij}}^{\mathrm{K}},\end{array}$$(19)

where ${\sigma }_{ϵ}^2$ is the variance of the observed ellipticities in the survey, ${\overline{n}}_{\mathrm{g}}$ is the surface density of galaxies in the survey, N_bin is the survey’s number of tomographic redshift bins, and ${\delta }_{\mathit{ij}}^{\mathrm{K}}$ is the Kronecker delta; here it indicates that the shot noise vanishes for cross-correlation spectra, as the ellipticities of galaxies at differing redshifts should not be correlated.

2.3. Reduced shear approximation

When relaxing the reduced shear approximation completely and explicitly is intractable, this can be aptly modelled by applying a second-order Taylor expansion (Dodelson et al. 2006; Shapiro 2009; Krause & Hirata 2010; Deshpande et al. 2020a) to Eq. (1), which gives:

$$\begin{array}{c}\hfill g^{\alpha }(\mathit{\theta })={\gamma }^{\alpha }(\mathit{\theta })+({\gamma }^{\alpha }\kappa )(\mathit{\theta })+\mathcal{O}({\kappa }^3).\end{array}$$(20)

Computing the angular E-mode power spectra using this expanded expression results in the standard two-point expression of Eq. (11), plus three-point terms. These additional terms, $\delta ⟨{\stackrel{\sim }{E}}_i(\mathit{\ell }){\stackrel{\sim }{E}}_j({\mathit{\ell }}^{\prime })⟩$, are given by:

$$\begin{array}{cc}& \delta ⟨{\stackrel{\sim }{E}}_i(\mathit{\ell }){\stackrel{\sim }{E}}_j({\mathit{\ell }}^{\prime })⟩={\displaystyle \sum _{\alpha }\sum _{\beta }{T}^{\alpha }(\mathit{\ell })T^{\beta }({\mathit{\ell }}^{\prime })⟨{\stackrel{\sim }{({\gamma }^{\alpha }\kappa )}}_i(\mathit{\ell }){\stackrel{\sim }{\gamma }}_j^{\beta }({\mathit{\ell }}^{\mathbf{\prime }})⟩}\hfill \\ \hfill & +T^{\alpha }({\mathit{\ell }}^{\prime })T^{\beta }(\mathit{\ell })⟨{\stackrel{\sim }{({\gamma }^{\alpha }\kappa )}}_j({\mathit{\ell }}^{\prime }){\stackrel{\sim }{\gamma }}_i^{\beta }(\mathit{\ell })⟩\hfill \\ \hfill & ={(2\pi )}^2{\delta }_{\mathrm{D}}^{(2)}(\mathit{\ell }+{\mathit{\ell }}^{\prime })\delta C_{\ell ;ij}^{\mathrm{RS}},\hfill \end{array}$$(21)

where $\delta C_{\ell ;ij}^{\mathrm{RS}}$ is the corresponding correction to $C_{\ell ;ij}^{\gamma \gamma }$, and is given by:

$$\begin{array}{cc}\hfill \delta C_{\ell ;ij}^{\mathrm{RS}}& ={\displaystyle {\int }_0^{\infty }\frac{{\mathrm{d}}^2{\mathit{\ell }}^{\prime }}{{(2\pi )}^2}cos(2{\varphi }_{{\ell }^{\prime }})B_{\mathit{ij}}^{\kappa \kappa \kappa }(\mathit{\ell },{\mathit{\ell }}^{\prime },-\mathit{\ell }-{\mathit{\ell }}^{\prime }),}\hfill \end{array}$$(22)

where we are always free to choose a coordinate system such that ϕ_ℓ = 0, and accordingly the correction only depends on the magnitude, ℓ. It depends on the two-redshift convergence bispectrum, $B_{\mathit{ij}}^{\kappa \kappa \kappa }$, which is the three-point counterpart of the convergence power spectrum. Higher order terms in the Taylor expansion of Eq. (20) would result in corrections dependent on the matter trispectrum, as well as the Wick contraction terms of O($P_{\delta \delta }^2$). Both types of terms have been shown to be sub-dominant (Cooray & Hu 2002; Shapiro & Cooray 2006; Dodelson et al. 2006; Krause & Hirata 2010; Deshpande et al. 2020a). The latter type of term, although it is of the same perturbative order in the power spectrum as the bispectrum, is still of O(W(χ)⁴), and given that typically χW(χ)≪1, it will still be significantly smaller than the correction of Eq. (22).

Additionally, just as the convergence power spectrum is the projection of the matter power spectrum, the convergence bispectrum is analogously the projection of the matter bispectrum, B_δδδ. Under the Limber approximation, it takes the following form:

$$\begin{array}{cc}& B_{\mathit{ij}}^{\kappa \kappa \kappa }({\mathit{\ell }}_1,{\mathit{\ell }}_2,{\mathit{\ell }}_3)=B_{\mathit{iij}}^{\kappa \kappa \kappa }({\mathit{\ell }}_1,{\mathit{\ell }}_2,{\mathit{\ell }}_3)+B_{\mathit{ijj}}^{\kappa \kappa \kappa }({\mathit{\ell }}_1,{\mathit{\ell }}_2,{\mathit{\ell }}_3)\hfill \\ \hfill & ={\displaystyle {\int }_0^{{\chi }_{\mathrm{h}}}\frac{\mathrm{d}\chi }{S_K^4(\chi )}{W}_i(\chi )W_j(\chi )[W_i(\chi )+W_j(\chi )]}\hfill \\ \hfill & \times B_{\delta \delta \delta }({\mathit{k}}_1,{\mathit{k}}_2,{\mathit{k}}_3,\chi ).\hfill \end{array}$$(23)

For a relaxation of the Limber approximation, we refer to Deshpande & Kitching (2020).

It should also be noted that the use of the reduced shear approximation can produce a B-mode signal contribution. However, it has been demonstrated that this is negligible (Schneider et al. 2002).

2.4. Source-lens clustering

Since, in practice, cosmic shear is only measured where galaxies are present, care must be taken to account for biases from any correlations between background source galaxies and the foreground lensing field. Given that, in reality, tomographic bins must be wide enough to include a sufficient number of galaxies so that shape-measurement noise is minimised, there will be overlap between the source and lensing distributions. The situation is further aggravated by broadening of bins due to photometric redshift uncertainties.

As a consequence of this effect, the observed number density of galaxies used in a given estimator which determines the shear angular power spectra from data is correlated with the intrinsic source galaxy overdensity, ${\delta }_i^g$, such that (Bernardeau 1998; Hamana et al. 2002; Schmidt et al. 2009):

$$\begin{array}{c}\hfill n_i^{\mathrm{obs}}(\mathit{\theta },\chi )=n_i(\chi )[1+{\delta }_i^g(\mathit{\theta })].\end{array}$$(24)

Accordingly, the shear used in the theoretical formalism for inference, is similarly replaced with an ‘observed’ shear,

$$\begin{array}{cc}\hfill {\gamma }_{\mathrm{obs};i}^{\alpha }(\mathit{\theta })& ={\gamma }_i^{\alpha }(\mathit{\theta })+{\gamma }_i^{\alpha }(\mathit{\theta }){\delta }_i^g(\mathit{\theta }).\hfill \end{array}$$(25)

This is similar in form to the Taylor expansion of the reduced shear expressed in Eq. (20), resulting in an analogous correction term, $\delta C_{\ell ;ij}^{\mathrm{SLC}}$, to the angular power spectra,

$$\begin{array}{cc}\hfill \delta C_{\ell ;ij}^{\mathrm{SLC}}& ={\displaystyle {\int }_0^{\infty }\frac{{\mathrm{d}}^2{\mathit{\ell }}^{\prime }}{{(2\pi )}^2}cos(2{\varphi }_{{\ell }^{\prime }})B_{\mathit{ij}}^{\kappa {\delta }^g\kappa }(\mathit{\ell },{\mathit{\ell }}^{\prime },-\mathit{\ell }-{\mathit{\ell }}^{\prime }),}\hfill \end{array}$$(26)

where $B_{ij}^{\kappa {\delta }^g\kappa }$ is now the two-redshift convergence-galaxy bispectrum. By adopting a linear galaxy bias model (so that δ_g = b δ) as used in EC20, and noting that δ^g is the 2D projection of δ_g, the convergence-galaxy bispectrum can also be expressed as a projection of the matter bispectrum:

$$\begin{array}{cc}& B_{\mathit{ij}}^{\kappa {\delta }^g\kappa }({\mathit{\ell }}_1,{\mathit{\ell }}_2,{\mathit{\ell }}_3)=B_{\mathit{iij}}^{\kappa {\delta }^g\kappa }({\mathit{\ell }}_1,{\mathit{\ell }}_2,{\mathit{\ell }}_3)+B_{\mathit{ijj}}^{\kappa {\delta }^g\kappa }({\mathit{\ell }}_{\mathbf{1}},{\mathit{\ell }}_2,{\mathit{\ell }}_3)\hfill \\ \hfill & ={\displaystyle {\int }_0^{{\chi }_{\mathrm{h}}}\frac{\mathrm{d}\chi }{S_K^4(\chi )}[b_i{n}_i(\chi )+b_j{n}_j(\chi )]W_i(\chi )W_j(\chi )}\hfill \\ \hfill & \times B_{\delta \delta \delta }({\mathit{k}}_1,{\mathit{k}}_2,{\mathit{k}}_3,\chi ),\hfill \end{array}$$(27)

where b_i and b_j are the galaxy biases for tomographic bins i and j, respectively. While more complex models of the galaxy bias exist, we proceed with the linear bias in this work, in order to mitigate the already significant computational load of these three-point terms. We note that modelling the galaxy bias requires more complexity at smaller scales, where the SLC effect is most relevant. When ultimately computing this term in the Euclid cosmological analysis, the final Euclid galaxy bias model should be used. The linear galaxy bias for each tomographic bin is given by:

$$\begin{array}{c}\hfill b_i=\sqrt{1+\overline{z_i}},\end{array}$$(28)

where $\overline{z_i}$ is tomographic bin i’s central redshift. For a review of galaxy bias models, we refer to Desjacques et al. (2018).

In addition to this contribution to the E-mode angular power spectra, source-lens clustering produces a B-mode signal as well. This term is comparable to the E-mode correction in magnitude, and accordingly, its detection in the absence of other B-mode contributions could allow for a direct correction of the E-mode signal, rather than requiring the computation of Eq. (26). However, typical B-mode signals are dominated by other contributions (Schneider et al. 2002; Yu et al. 2015).

We note that there are two main references to ‘source-lens clustering’ in the literature. One is the correlation between the shear and density in the source plane (see e.g. Krause et al. 2021, Eqs. (47)–(49)); otherwise referred to as an intrinsic lens-source density correlation. In general, this should be expected to be zero, since the shear should not be correlated with the source density distribution.

The other effect is that described in Schmidt et al. (2009), namely, the sampling effect (or source number-density weighting) caused by the use of estimators (Schmidt et al. 2009 Eq. (3)) for the two-point statistics of the shear field that leads to their Eq. (7) (and Eq. (24) of this paper). This is different from the ‘intrinsic lens-source density’ correlation, as it is caused only as a result of a sampling effect in the estimator. In this case, the form of the correction is exactly the same form as the reduced shear, except the weight functions are different (since the shear is modified by delta rather than kappa). However, as noted by Schmidt et al. (2009) ‘for a sufficiently narrow redshift distribution of source galaxies, this source-lens clustering is negligible’. This is because the distribution of sources and lenses do not overlap’.

In Krause et al. (2021) the second effect is referred to as the ‘source clustering ansatz’. It is noted that ‘in the limit of narrow tomographic bins considered by Schmidt et al., both approaches yield the same power spectrum corrections’ and that their equations are ‘partially cancelled out’ (i.e. account for) biases caused by the sampling effect. They also noted that the DES Y3 galaxygalaxy lensing analysis corrects for sampling effects at the measurement level through so-called ’Boost factors’ that cancel the effect of source number density weighting.

In this paper, we only consider the sampling effect but we note that it may be possible to use ‘de-biased’ estimators to mitigate the effect; this requires further work to demonstrate that such approaches are sufficient for Stage IV experiments. However, it should be noted that in the rest of this paper, the ‘SLC’ terms are dependent on the estimator used in the analysis. Thus, what we are showing here is the worst case where no estimator-level correction has been applied.

Finally, we note that computing the source number density weighting term relies on the bispectrum $B_{ijj}^{\kappa {\delta }^g\kappa }(\mathcal{l})$, whose amplitude will need to be determined from simulations and could therefore could deviate from the ansatz used in this paper (Eq. (27)). This is being further investigated in Linke et al. (prep., and priv. comm.) and could result in the raw (not de-biased) source number density weighting also being small.

Therefore, for all cases we show the results with and without the ‘SLC’ terms to reflect the case that all SLC terms (intrinsic source-lens clustering and source number density weighting) are zero, which is likely to be the case.

2.5. Magnification bias

An additional consequence of gravitational lensing is that the density of galaxies observed by a particular survey is no longer representative of the true underlying galaxy density (Turner et al. 1984). In particular, magnification resulting from the convergence modifies the density in two contrasting ways.

One manifestation of the effect is that individual sources are magnified and as a consequence of this, their flux increases. Accordingly, any sources lying just beyond the flux limit of the survey may have their fluxes increased to the point of then being within the flux limit; increasing the observed density. As sources are magnified, the patch of sky around them is also magnified. This causes the second, competing manifestation. Within the magnified patch of sky, the galaxy density is reduced. The total effect, known as magnification bias, is dependent on the slope of the unlensed galaxy luminosity function. This assumes that the magnification μ > 1.

Assuming that, on our scales of interest, fluctuations in the intrinsic galaxy overdensity are small, and taking into account that, for weak lensing, |κ|≪1, the observed galaxy overdensity for a given tomographic bin, ${\delta }_{\text{obs};i}^g$, is given by (Hui et al. 2007; Schmidt et al. 2009):

$$\begin{array}{c}\hfill {\delta }_{\mathrm{obs};i}^g(\mathit{\theta })={\delta }_i^g(\mathit{\theta })+(5s_i-2){\kappa }_i(\mathit{\theta }),\end{array}$$(29)

where ${\delta }_i^g$ is the intrinsic galaxy overdensity in the absence of magnification or any other systematic effects, and s_i is the slope of the luminosity function for a redshift bin, i. This is given by the derivative of the cumulative galaxy number counts with respect to magnitude, m, evaluated at the survey’s limiting magnitude, m_lim, such that

$$\begin{array}{c}\hfill s_i=\frac{\partial {\log }_{10}\mathfrak{n}(\overline{z_i},m)}{\partial m}{|}_{m_{\lim }},\end{array}$$(30)

where $\mathfrak{n}(\overline{z_i},m)$ is the true, underlying distribution of galaxies, evaluated at the tomographic bin’s central redshift, $\overline{z_i}$. Here, we have suppressed an additional dependence on the wavelength band in which the galaxy is observed. This should be considered when determining the slope from observational data.

Accordingly, Eq. (25) gains an extra term:

$$\begin{array}{cc}\hfill {\gamma }_{\mathrm{obs};i}^{\alpha }& ={\gamma }_i^{\alpha }(\mathit{\theta })+{\gamma }_i^{\alpha }(\mathit{\theta }){\delta }_i^g(\mathit{\theta })+(5s_i-2){\gamma }_i^{\alpha }(\mathit{\theta }){\kappa }_i(\mathit{\theta }).\hfill \end{array}$$(31)

This additional term is near-identical to the second term in Eq. (20), but for the prefactor of (5s_i − 2). Accordingly, it also spawns a correction to the angular power spectra. This correction for magnification bias, $\delta C_{\ell ;ij}^{\mathrm{MB}}$, takes a similar form to the reduced shear correction of Eq. (22),

$$\begin{array}{cc}& \delta C_{\ell ;ij}^{\mathrm{MB}}={\displaystyle {\int }_0^{\infty }\frac{{\mathrm{d}}^2{\mathit{\ell }}^{\prime }}{{(2\pi )}^2}cos(2{\varphi }_{{\ell }^{\prime }})[(5s_i-2)B_{\mathit{iij}}^{\kappa \kappa \kappa }(\mathit{\ell },{\mathit{\ell }}^{\prime },-\mathit{\ell }-{\mathit{\ell }}^{\prime })}\hfill \\ \hfill & +(5s_j-2)B_{\mathit{ijj}}^{\kappa \kappa \kappa }(\mathit{\ell },{\mathit{\ell }}^{\prime },-\mathit{\ell }-{\mathit{\ell }}^{\prime })].\hfill \end{array}$$(32)

Given the similarity of the magnification bias correction to the source-lens clustering and reduced shear corrections, it too would produce a contribution to the B-mode signal. While this has not been explicitly evaluated, we would expect this term (as in the case of its reduced shear and source-lens clustering counterparts) to be sub-dominant.

2.6. Source obscuration

There is another systematic effect which can change the observed galaxy number density. Blending of close galaxy pairs can lead to multiple galaxies being discarded or counted as a lower number than they are (Hartlap et al. 2011). The resulting change in the observed number density of galaxies, Δ n^SO(z, θ), can be modelled by:

$$\begin{array}{cc}\hfill \frac{\mathrm{\Delta }{n}^{\mathrm{SO}}(z,\mathit{\theta })}{n(z,\mathit{\theta })}& =-\pi [{(2\vartheta )}^2{n}_{\text{tot}}+\frac{An(z)}{2-\zeta }{(2\vartheta )}^{2-\zeta }],\hfill \end{array}$$(33)

where n_tot is the total number density of galaxies at all redshifts, n(z) is the observed density of galaxies at redshift z ignoring source-lens clustering, we assume a redshift-independent radius, ϑ, for all galaxies as in Hartlap et al. (2011), and A and ζ are the amplitude and power-law index of a power-law model for the two-point galaxy angular correlation function. This expression is obtained by considering the probability that the centroid of another source lies within 2ϑ of a given one by integrating over the probability that another source centroid lies in an annulus of dθ around the centroid of another one.

Instead, we assume that the blending strategy for Euclid will account for blended pairs sufficiently well, such that the only obscuration of concern is substantial overlap; when the centroid of a source is behind another source (i.e. within ϑ rather than 2ϑ). In this case, we are only concerned with the probability of this overlap. We note that in reality this blending has a complex interaction with shape measurement, but evaluating this is out of the scope of this work. Then, assuming that sources are approximately circular, the probability, dp, of a galaxy at redshift z overlapping with one at redshift z′ is:

$$\begin{array}{c}\hfill \mathrm{d}p(z,z^{\prime },\mathit{\theta })=\pi {\vartheta }^{2}n(z^{\prime },\mathit{\theta })\mathrm{d}z.\end{array}$$(34)

Here, it is also assumed that the expected number of galaxies overlapping with a given galaxy is ≪1, such that the probability of at least one galaxy overlapping a given source (resulting in the removal of that source from the sample) is equal to the probability of just one overlap, which by Poisson statistics is then the expected number of overlaps. Accordingly, the total change in the number of sources at z is then given by:

$$\begin{array}{cc}& \mathrm{\Delta }{n}^{\mathrm{SO}}(z,\mathit{\theta })=-\pi {\vartheta }^{2}n(z,\mathit{\theta }){\displaystyle {\int }_0^{\infty }\mathrm{d}{z}^{\prime }n(z^{\prime },\mathit{\theta })}\hfill \\ \hfill & =-\pi {\vartheta }^{2}n(z,\mathit{\theta })[1+{\delta }^g(\mathit{\theta })]{\displaystyle {\int }_0^{\infty }\mathrm{d}{z}^{\prime }n(z^{\prime })}\hfill \\ \hfill & =-\pi {\vartheta }^{2}n(z,\mathit{\theta })[1+{\delta }^g(\mathit{\theta })]n_{\mathrm{tot}}.\hfill \end{array}$$(35)

We neglect the second term on the right-hand side of Eq. (33), as it specifically accounts for the correlated overlap of galaxies at the same redshift within a fixed disk around the source, in addition to the random one already included. Given that the fractional change is calculated by integrating over redshift slices, that this term would only appear for the slice where z′=z, and that the source obscuration term itself is small (see Sect. 4), we expect it to be safely negligible.

Adopting the tomographic redshift binning approach, for a given source can be obscured by another in the same bin, or by sources in lower redshift bins than the one the source belongs in. The fractional change in the number density of galaxies in redshift bin i then becomes:

$$\begin{array}{cc}\hfill \frac{\mathrm{\Delta }{n}_i^{\mathrm{SO}}(z,\mathit{\theta })}{n_i(z,\mathit{\theta })}& =-\pi {\vartheta }^2{\displaystyle \sum _{q=1}^i[1+{\delta }_q^g(\mathit{\theta })]n_{\mathrm{tot};q}}\hfill \\ \hfill & =-\pi {\vartheta }^2{\displaystyle \sum _{q=1}^i[1+b_q{\delta }^g(\mathit{\theta })]n_{\mathrm{tot};q}}\hfill \\ \hfill & =-\pi {\vartheta }^2{n}_{\mathrm{cumul}.;i}-\pi {\vartheta }^2{\delta }^g{f}_{\mathrm{SO};i},\hfill \end{array}$$(36)

where the n_tot;  q is the total surface density of galaxies for redshift bin q (that is the integral over n_q(z) for bin q), n_cumul.;i is the cumulative total surface density of galaxies for all redshift bins up to and including bin i, and

$$\begin{array}{c}\hfill f_{\mathrm{SO};\mathrm{i}}={\displaystyle \sum _{\rho =1}^i{b}_{\rho }{n}_{\mathrm{tot};\rho }.}\end{array}$$(37)

Including the effect of source obscuration in addition to source-lens clustering and magnification bias, the observed number density for a given tomographic redshift bin, i, becomes

$$\begin{array}{cc}& n_i^{\mathrm{obs}}(\mathit{\theta },\chi )=n_i(\chi )[1+(5s_i-2){\kappa }_i(\mathit{\theta })+{\delta }_i^g(\mathit{\theta })]\hfill \\ \hfill & \times [1-\pi {\vartheta }^2{n}_{\mathrm{cumul}.;i}-\pi {\vartheta }^2{\delta }^g(\mathit{\theta })f_{\mathrm{SO};i}].\hfill \end{array}$$(38)

From here, only terms to first-order in the lensing potential are retained in order to suppress fourth-order or higher terms appearing in the two-point statistic. Then, it can be seen that source obscuration adds prefactors to the base angular shear power spectra: the source-lens clustering correction from Eq. (26) and the magnification bias correction of Eq. (32). Accordingly, source obscuration produces three new correction terms:

$$\begin{array}{cc}& \delta C_{\ell ;ij}^{\mathrm{SO}}=({\pi }^2{\vartheta }^4{n}_{\mathrm{cumul}.;i}{n}_{\mathrm{cumul}.;j}-\pi {\vartheta }^2{n}_{\mathrm{cumul}.;i}\hfill \\ \hfill & -\pi {\vartheta }^2{n}_{\mathrm{cumul}.;j})C_{\ell ;ij}^{\gamma \gamma },\hfill \end{array}$$(39)

$$\begin{array}{cc}\hfill & \delta C_{\ell ;ij}^{\mathrm{SO}-\mathrm{SLC}}={\displaystyle {\int }_0^{\infty }}\frac{{\mathrm{d}}^2{\mathit{\ell }}^{\prime }}{{(2\pi )}^2}cos(2{\varphi }_{{\ell }^{\prime }})\hfill \\ \hfill & \times [({\pi }^2{\vartheta }^4{n}_{\mathrm{cumul}.;i}{n}_{\mathrm{cumul}.;j}-\pi {\vartheta }^2{n}_{\mathrm{cumul}.;i}{b}_i\hfill \\ \hfill & -\pi {\vartheta }^2{n}_{\mathrm{cumul}.;j}{b}_i-\pi {\vartheta }^2(1-n_{\mathrm{cumul}.;j})\hfill \\ \hfill & \times f_{\mathrm{SO};i})B_{\mathit{iij}}^{\kappa {\delta }^g\kappa }({\mathit{\ell }}_1,{\mathit{\ell }}_2,{\mathit{\ell }}_3)+i\leftrightarrow j],\hfill \end{array}$$(40)

$$\begin{array}{cc}& \delta C_{\ell ;ij}^{\mathrm{SO}-\mathrm{MB}}=({\pi }^2{\vartheta }^4{n}_{\mathrm{cumul}.;i}{n}_{\mathrm{cumul}.;j}-\pi {\vartheta }^2{n}_{\mathrm{cumul}.;i}\hfill \\ \hfill & -\pi {\vartheta }^2{n}_{\mathrm{cumul}.;j})\delta C_{\ell ;ij}^{\mathrm{MB}},\hfill \end{array}$$(41)

where i ↔ j indicates a repetition of the preceding bispectrum term and its pre-factor, with all instances of the i and j bin indices exchanged.

2.7. Local Universe effects

A further effect to consider is that the observed two-point statistic at our location may be biased due to local over or under-densities. Accordingly, the angular power spectra must be calculated conditioned on the local density (Hall 2020). The local density contrast, δ₀, can be defined as the matter density contrast smoothed by a top-hat kernel of comoving radius R according to

$$\begin{array}{c}\hfill {\delta }_0(R,\chi )\equiv \frac{3}{4\pi R^3}{\displaystyle \int }{\mathrm{d}}^3\mathit{r}\mathrm{\Theta }(R-|\mathit{r}|)\delta (\mathit{r},\chi ),\end{array}$$(42)

where the matter density contrast is now expressed in terms of spatial distance, r, rather than angle on the sky, and Θ is the Heaviside step-function.

Then, the conditional angular power spectra can be obtained using the Edgeworth expansion for conditional distributions. This calculation is mathematically intensive and, accordingly, it is not reproduced here. The full derivation can be found in Hall (2020). Under the Limber approximation, and assuming ℓ ≫ 1 (and note that cosmic shear is only defined for ℓ ≥ 2), this expression consists of two terms, the standard power spectra of Eq. (12) and a correction term, $\delta C_{\ell ;ij}^{\mathrm{LU}}$, which is defined as

$$\begin{array}{cc}& \delta C_{\ell ;ij}^{\mathrm{LU}}=2\frac{{\delta }_0(R,\chi )}{{\sigma }^2(R,\chi )}{\displaystyle {\int }_0^{{\chi }_{\mathrm{h}}}\mathrm{d}\chi \frac{W_i(\chi )W_j(\chi )}{{\chi }^2}}\hfill \\ \hfill & \times \{[\frac{34}{21}{\xi }_R(\chi )-\frac{4}{21}{\psi }_R(\chi )]P_{\delta \delta }(k,\chi )\hfill \\ \hfill & +[\frac{\chi }{\ell }{\xi }_R^{\prime }(\chi )-\frac{\ell }{\chi }{\mathrm{\Omega }}_R(\chi )]\frac{\partial P_{\delta \delta }(k,\chi )}{\partial k}\frac{1}{\chi }\hfill \\ \hfill & -\frac{4}{7}{\psi }_R(\chi )\frac{{\partial }^2{P}_{\delta \delta }(k,\chi )}{\partial k^2}\frac{1}{{\chi }^2}\},\hfill \end{array}$$(43)

where σ² is the variance of the local density contrast, it is assumed the ratio of the local density contrast to its variance is constant with comoving distance, ${\xi }_R,{\psi }_R,{\xi }_R^{\prime }$, and Ω_R are correlation functions defined in Hall (2020). Also, the expression assumes a flat-geometry, which is valid under the current constraints on Ω_K, as lenses are much less far away than the curvature distance. Additionally, we note that this expression is derived using only the tree-level Eulerian perturbation theory expression for the matter bispectrum.

2.8. The flat Universe assumption

Typically when computing cosmic shear angular power spectra, spatially non-flat universes are accounted for through modifying comoving distances based on curvature, as described by Eq. (4). In practice, however, curvature also modifies the projection kernel (Taylor et al. 2018b).

Under the assumption of a spatially flat Universe, the Poisson equation gives the relationship between the comoving Newtonian gravitational potential, ϕ, and the matter density contrast,

$$\begin{array}{c}\hfill {\mathrm{\nabla }}_{\chi }^2\varphi (\mathit{r},\chi )=\frac{3{\mathrm{\Omega }}_{\mathrm{m}}{H}_0^2}{2c^{2}a(t)}\delta (\mathit{r},\chi ),\end{array}$$(44)

where ${\mathrm{\nabla }}_{\chi }^2$ is the Laplacian for a spatially flat Universe. This allows for the shear angular power spectra to be expressed in terms of the matter power spectrum, as in Eq. (12). However, the matter density contrast has rectilinear coordinates, whereas the lensing potential is defined in terms of angular coordinates (r, θ, φ), from the observer’s frame of reference. Relating the two as above requires expressing the potential in spherical Bessel space as:

$$\begin{array}{c}\hfill {\varphi }_{\ell m}(k)=\sqrt{\frac{2}{\pi }}{\displaystyle \int }{\mathrm{d}}^{3}r\varphi (r)j_{\ell }(kr)Y_{\ell m}(\theta ,\phi ),\end{array}$$(45)

where j_ℓ refers to the spherical Bessel functions and Y_ℓm to the spherical harmonics. Owing to the fact that these spherical Bessel functions and spherical harmonics are eigenfunctions of the Laplacian, the following relationship is obtained:

$$\begin{array}{c}\hfill ({\mathrm{\nabla }}_r^2+k^2)j_{\ell }(kr)Y_{\ell m}(\theta ,\phi )=0.\end{array}$$(46)

This allows for the relation of the lensing potential to the matter density contrast in the spherical harmonic space to be expressed, as well as the eventual calculation of Eq. (12) under the Limber approximation. We refer to Kitching et al. (2017) for a full derivation.

However, in the case of a spatially non-flat Universe, the Laplacian in Eqs. (44)–(46) must be replaced by one corresponding to a curved geometry, ${\nabla }_{S_K}^2$. Accordingly, the projection kernel must also be modified; by replacing the spherical Bessel functions in Eq. (46) with hyper-spherical Bessel functions, ${\mathrm{\Phi }}_{\ell }^{\beta }$, so that

$$\begin{array}{c}\hfill ({\mathrm{\nabla }}_{S_K}^2+k^2){\mathrm{\Phi }}_{\ell }^{\beta }(r)Y_{\ell m}(\theta ,\phi )=0,\end{array}$$(47)

where $\beta =\sqrt{(k^2+K)/|K|}$ (Lesgourgues & Tram 2014). Consequently, the shear angular power spectra for a spatially non-flat Universe, $C_{\ell ;ij}^{\gamma \gamma ;\mathrm{NF}}$, under the Limber approximation, under the Limber approximation, is given by modifying Eq. (12) to be:

$$\begin{array}{c}\hfill C_{\ell ;ij}^{\gamma \gamma }={\displaystyle {\int }_0^{{\chi }_{\mathrm{h}}}\mathrm{d}\chi W_{\ell }^{\mathrm{NF}}(\chi ;K)\frac{W_i(\chi )W_j(\chi )}{S_K^2(\chi )}{P}_{\delta \delta }(k,\chi ),}\end{array}$$(48)

where

$$\begin{array}{c}\hfill W_{\ell }^{\mathrm{NF}}(\chi ;K)={[1-\mathrm{sgn}(K)\frac{{\ell }^2}{{(\ell +1/2)}^2/S_K^2(\chi )+K}]}^{-1/2}.\end{array}$$(49)

Here, sgn(K) is the sign of the curvature, K. Alternatively, for consistency with the previously discussed corrections, this can be expressed as a correction term, $\delta C_{\ell ;ij}^{\mathrm{NF}}$, to the spatially flat angular power spectra such that

$$\begin{array}{c}\hfill \delta C_{\ell ;ij}^{\mathrm{NF}}=C_{\ell ;ij}^{\gamma \gamma ;\mathrm{NF}}-C_{\ell ;ij}^{\gamma \gamma }.\end{array}$$(50)

2.9. The Fisher matrix and bias formalism

The constraining power of cosmological surveys, in terms of the uncertainties on inferred cosmological parameters, is often predicted by using the Fisher matrix formalism. It also allows for the quantification of biases in this inference resulting from neglecting systematic effects within the signal itself. Here, we use this technique to predict how biased cosmological parameters inferred from Euclid would be when the previously discussed systematic effects are neglected. We have followed the conventions of Euclid Collaboration (2020) in all calculations and validated the code with respect to those results⁶.

Explicitly, the Fisher matrix is defined as the expected value of the Hessian of the log likelihood (defined for a Gaussian likelihood, and applied to CMB data in Tegmark et al. 2015). For Stage IV weak lensing cosmology, it has been demonstrated that the likelihood can safely be assumed to be Gaussian (Lin et al. 2020); TaylorNG; UphamGauss; HT22. Accordingly, the Fisher matrix for cosmic shear is defined as:

$$\begin{array}{c}\hfill F_{\mu \nu }={\displaystyle \sum _{{\ell }^{\prime }={\ell }_{\min }}^{{\ell }_{\max }}\sum _{\ell ={\ell }_{\min }}^{{\ell }_{\max }}\sum _{ij,mn}\frac{\partial C_{\ell ;ij}^{ϵϵ}}{\partial {\theta }_{\mu }}{\mathrm{Cov}}^{-1}[C_{\ell ;ij}^{ϵϵ},C_{{\ell }^{\prime };mn}^{ϵϵ}]\frac{\partial C_{{\ell }^{\prime };mn}^{ϵϵ}}{\partial {\theta }_{\nu }},}\end{array}$$(51)

where the μ and ν indices denote element in the Fisher matrix associated with cosmological parameters, θ_μ and θ_ν, respectively, ℓ_min is the minimum angular wavenumber of the survey, ℓ_max is the maximum angular wavenumber used, the sums are over the ℓ-blocks of power spectrum bands, and ${\mathrm{Cov}}^{-1}[C_{\ell ;ij}^{ϵϵ},C_{\ell \prime ;mn}^{ϵϵ}]$ is the inverse of the covariance of the angular power spectra signal.

In practice, this covariance term is non-Gaussian (Barreira et al. 2018a; Takada & Hu 2013; Upham et al. 2022), with an additional contribution arising from the super-sample covariance (SSC; Hu & Kravtsov 2003). This SSC term encapsulates the effects on the covariance of density fluctuations with wavelengths larger than the extent of the galaxy survey. Such fluctuations result in the background density of the survey ceasing to be representative of the underlying density of the Universe. The total covariance is then the sum of the Gaussian, Cov_G and SSC, and Cov_SSC terms, as follows:

$$\begin{array}{cc}& \mathrm{Cov}[C_{\ell ;ij}^{ϵϵ},C_{{\ell }^{\prime };mn}^{ϵϵ}]={\mathrm{Cov}}_{\mathrm{G}}[C_{\ell ;ij}^{ϵϵ},C_{{\ell }^{\prime };mn}^{ϵϵ}]\hfill \\ \hfill & +{\mathrm{Cov}}_{\mathrm{SSC}}[C_{\ell ;ij}^{ϵϵ},C_{{\ell }^{\prime };mn}^{ϵϵ}],\hfill \end{array}$$(52)

where the Gaussian component is given by

$$\begin{array}{c}\hfill {\mathrm{Cov}}_{\mathrm{G}}[C_{\ell ;ij}^{ϵϵ},C_{{\ell }^{\prime };mn}^{ϵϵ}]=\frac{C_{\ell ;im}^{ϵϵ}{C}_{{\ell }^{\prime };jn}^{ϵϵ}+C_{\ell ;in}^{ϵϵ}{C}_{{\ell }^{\prime };jm}^{ϵϵ}}{(2\ell +1)f_{\mathrm{sky}}\mathrm{\Delta }\ell }{\delta }_{\ell {\ell }^{\prime }}^{\mathrm{K}},\end{array}$$(53)

where f_sky is the fraction of the sky observed by the galaxy survey, Δℓ is the bandwidth of the ℓ-modes sampled, and δ^K is the Kronecker delta. Other non-Gaussian terms in the covariance can be neglected (see e.g. Barreira et al. 2018b). The SSC component is well approximated by (Lacasa & Grain 2019)

$$\begin{array}{c}\hfill {\mathrm{Cov}}_{\mathrm{SSC}}[C_{\ell ;ij}^{ϵϵ},C_{{\ell }^{\prime };mn}^{ϵϵ}]\approx R_{\ell }{C}_{\ell ;ij}^{ϵϵ}{R}_{{\ell }^{\prime }}{C}_{{\ell }^{\prime };mn}^{ϵϵ}{S}_{\mathit{ijmn}},\end{array}$$(54)

where S_ijmn is the dimensionless volume-averaged covariance of the background matter density contrast, and R_ℓ is the effective relative response of the observed power spectrum. We assume that there is no interrelation between local Universe effects and the SSC, but this is a caveat that should be verified in the future.

The diagonal of the inverse of the Fisher matrix is used to predict the 1σ uncertainties on each of the parameters. Explicitly, the uncertainty, σ_μ, on parameter θ_μ is given by:

$$\begin{array}{c}\hfill {\sigma }_{\mu }=\sqrt{{F_{\mu \mu }}^{-1}}.\end{array}$$(55)

By extending this formalism, biases on inferred the parameters resulting from neglecting systematic effects can also be predicted (Taylor et al. 2007). For a given systematic, δC_ℓ; ij, the bias, 𝔟_μ, on parameter θ_μ, is given by:

$$\begin{array}{c}\hfill {\mathfrak{b}}_{\mu }={\displaystyle \sum _{\nu }{(F^{-1})}_{\mu \nu }{\mathcal{B}}_{\nu },}\end{array}$$(56)

where

$$\begin{array}{c}\hfill {\mathcal{B}}_{\nu }={\displaystyle \sum _{{\ell }^{\prime }={\ell }_{\min }}^{{\ell }_{\max }}\sum _{\ell ={\ell }_{\min }}^{{\ell }_{\max }}\sum _{ij,mn}\delta C_{\ell ;ij}{\mathrm{Cov}}^{-1}[C_{\ell ;ij}^{ϵϵ},C_{{\ell }^{\prime };mn}^{ϵϵ}]\frac{\partial C_{{\ell }^{\prime };mn}}{\partial {\theta }_{\nu }}.}\end{array}$$(57)

We note that we assume a Gaussian likelihood function but with a correlated covariance matrix (which includes non-Gaussian contributions). The extent to which this assumption is robust to relaxing the Gaussian likelihood assumption was explored in Martinelli et al. (2021), who found good agreement between Fisher matrix (Gaussian likelihood) predictions and full Markov chain Monte Carlo (MCMC) predictions. Additionally, Taylor et al. (2019) found a similar result that also allows for the possible derivation of a fully non-Gaussian likelihood function.

3. Methodology

In this section, we review the computational and modelling specifics used within this investigation. We begin by describing the survey specifications adopted. Then, details are given about our choice of fiducial cosmology, modelling of background quantities, and Fisher matrices. Lastly, we describe modelling choices made in the computation of the magnification bias, source obscuration, and local Universe effect corrections.

3.1. Survey specifications

For Euclid, forecasting specifications are specified in EC20; we adopted these here, but we note that our IA model is simpler (see Sect. 2.2), and we vary Ω_K, rather than the dark energy’s critical density. Specifically, we considered the ‘optimistic’ scenario described in that work, as this is the case where the cosmic shear probe is able to meet its precision goals by itself. Under this scenario, the survey is taken to extend up to ℓ-modes of 5000.

Additionally, the intrinsic variance of observed ellipticities is taken to consist of two components; each with a magnitude of 0.21. Correspondingly, the root mean square (RMS) intrinsic ellipticity variance is ${\sigma }_{ϵ}=\sqrt{2}\times 0.21\approx 0.3$⁷. Euclid is also expected to have a survey area such that f_sky = 0.36. The survey’s galaxy surface density is anticipated to be ${\overline{n}}_{\mathrm{g}}=30$ arcmin⁻².

The cosmic shear probe of Euclid is planned to observe sources between redshifts of 0 and 2.5, and utilise 10 equipopulated tomographic redshift bins with the following edges: {0.001, 0.418, 0.560, 0.678, 0.789, 0.900, 1.019, 1.155, 1.324, 1.576, and 2.50}.

Given that Euclid will use photometric redshifts, the model for the source distributions within these tomographic bins must account for photometric redshift uncertainties. Accordingly, for a particular bin, i, the galaxy redshift distribution, n_i(z), is described as follows:

$$\begin{array}{c}\hfill n_i(z)=\frac{{\int }_{z_i^{-}}^{z_i^{+}}\mathrm{d}{z}_{\mathrm{p}}\mathfrak{n}(z)p_{\mathrm{ph}}(z_{\mathrm{p}}|z)}{{\int }_{z_{\min }}^{z_{\max }}\mathrm{d}z{\int }_{z_i^{-}}^{z_i^{+}}\mathrm{d}{z}_{\mathrm{p}}\mathfrak{n}(z)p_{\mathrm{ph}}(z_{\mathrm{p}}|z)},\end{array}$$(58)

where z_p is measured photometric redshift, $z_i^{-}$ and $z_i^{+}$ are the limits of the ith redshift bin, and z_min and z_max are the redshift limits of the survey itself. Additionally, 𝔫(z) is the underlying distribution of galaxies which here we modeled according to the formalism established in Laureijs et al. (2011):

$$\begin{array}{c}\hfill \mathfrak{n}(z)\propto (\frac{z}{z_0}{)}^2\exp [-(\frac{z}{z_0}{)}^{3/2}],\end{array}$$(59)

where $z_0=z_{\mathrm{m}}/\sqrt{2}$, and z_m = 0.9 is the median redshift of the survey. The remaining function in Eq. (58), p_ph(z_p|z), encapsulates the probability that a source measured to have a photometric redshift of z_p actually has a redshift of z. This distribution takes the form (Kitching et al. 2008):

$$\begin{array}{cc}& p_{\mathrm{ph}}(z_{\mathrm{p}}|z)=\frac{1-f_{\mathrm{out}}}{\sqrt{2\pi }{\sigma }_{\mathrm{b}}(1+z)}\exp \{-\frac{1}{2}[\frac{z-c_{\mathrm{b}}{z}_{\mathrm{p}}-z_{\mathrm{b}}}{{\sigma }_{\mathrm{b}}(1+z)}{]}^2\}\hfill \\ \hfill & +\frac{f_{\mathrm{out}}}{\sqrt{2\pi }{\sigma }_{\mathrm{o}}(1+z)}\exp \{-\frac{1}{2}[\frac{z-c_{\mathrm{o}}{z}_{\mathrm{p}}-z_{\mathrm{o}}}{{\sigma }_{\mathrm{o}}(1+z)}{]}^2\}.\hfill \end{array}$$(60)

Here, the distribution is expressed as the sum of two terms – the first is the uncertainty resulting from multiplicative and additive bias in redshift determination for the fraction of sources with a well measured redshift, whilst the second represents the same, but for a fraction of catastrophic outliers in the sample, f_out. The values used for the individual parameters in this parameterisation match the selection of EC20, and are stated in Table 2, which are fixed throughout our analysis.

3.2. Cosmological modelling and Fisher matrices

Throughout this investigation, we consider the ΛCDM cosmological model and its extension: the w₀w_aCDM model, which also allows for varying dark energy pressure and a separately parameterised dark energy equation of state at early times. The ΛCDM model uses seven parameters, which are defined thusly: the present-day total matter density parameter Ω_m, the present-day baryonic matter density parameter Ω_b, the dimensionless curvature parameter Ω_K = −K(c/H₀)², the Hubble parameter h = H₀/100 km s⁻¹ Mpc⁻¹, the spectral index n_s, the RMS value of density fluctuations on 8 h⁻¹ Mpc scales σ₈, and massive neutrinos with a sum of masses ∑m_ν ≠ 0. The w₀w_aCDM model additionally adds in the present-day value of the dark energy equation of state w₀, and the high-redshift value of the dark energy equation of state, w_a. Typically, the present-day densities Ω_i, i ∈ {m, b, K}, are denoted with an additional subscript 0; we omit this here for brevity. Primarily, we are interested the w₀w_aCDM case when discussing corrections in this investigation, as a key goal of Stage IV surveys is exploring models of dark energy. However, when examining the cosmological parameter biases, we also present the ΛCDM case, for completeness.

The specific values used for each of these parameters were also chosen for consistency with EC20 (shown in Table 3). As in EC20, the value of ∑m_ν ≠ 0 was treated as fixed, and we did not calculate uncertainties or biases for it. When computing biases for all corrections except for the non-flat Universe term, we set Ω_K to 0. Only when testing the significance of the additional non-flat Universe correction term was it set to 0.05. In this case, the selected value serves as the upper-limit of 1σ constraint on the parameter from Planck Collaboration VI (2020).

In cases where the non-flat Universe correction were not evaluated, our Fisher matrices matched those of the Euclid forecasting specification (EC20). Thus, they contain: Ω_m, Ω_b, h, n_s, σ₈ and 𝒜_IA for the ΛCDM case and, additionally, w₀ and w_a, for the w₀w_aCDM case. When the correction for spatially curvature needed to be tested, the Fisher matrices also include Ω_K. It was not necessary to include any further nuisance parameters within the matrix, as EC20 showed that the inclusion of various different nuisance parameters (e.g. those modelling the non-linear part of the matter power spectrum) typically altered the forecasted uncertainties on the cosmological parameters by less than 10%. The S_ijmn were calculated using the publicly available PySSC⁸ code (Lacasa & Grain 2019), with an R_ℓ of 3.

To calculate the cosmological background quantities required for the investigation, including the matter power spectrum and growth factor, we used the CAMB⁹ software package (Lewis et al. 2000). Additionally, we utilised the Halofit (Takahashi et al. 2012) implementation of the non-linear part of the power spectrum, and included additional corrections identified by Bird et al. (2012). Where necessary, we additionally employed Astropy¹⁰ (Astropy Collaboration 2013, 2018) to compute the cosmological distances. The NLA model IA parameters were set to 𝒜_IA = 1.72 and 𝒞_IA = 0.0134, again in accordance with EC20. The required partial derivatives required were computed numerically, using the procedure described in EC20. Throughout this work, all quantities were evaluated for 200 ℓ-bands. The limits for these were logarithmically spaced, with an ℓ_min of 10 and an ℓ_max of 5000.

3.3. Modelling the higher order corrections

To model the matter bispectrum required by the reduced shear, magnification bias, source-lens clustering, and source obscuration corrections, we used the BiHalofit model and code¹¹ (Takahashi et al. 2020). This represents the matter bispectrum using one-halo and three-halo terms, which have been determined via a fitting to N-body simulations.

For the magnification bias correction, we used the slope of the luminosity function as calculated from the fitting formula given in Appendix C of Euclid Collaboration (2022b). This is determined from the Euclid Flagship simulation (Potter et al. 2017) and for the limiting magnitude 24.5 of the VIS instrument (AB in the Euclid VIS band; Cropper et al. 2012). Therefore, it provides the most Euclid specific estimate of this quantity to date. However, once the Euclid survey is in-progress, we note that this quantity should be calculated directly from the observed data. We used a single value for the slope for each tomographic redshift bin. This value was calculated at the central redshift of the bin. The slopes for all bins, together with their central redshifts, can be found in Table 4. We note that the magnification bias from Euclid Collaboration (2022b) was obtained for the n(z) from Euclid Collaboration (2021) but we use the Euclid Collaboration (2020)n(z), however the effect of the small changes in the assumed n(z) should be small, which is consistent with the small differences in the results between this paper and Deshpande et al. (2020a).

In order to evaluate the source obscuration terms, we set the total number of observed galaxies to 2 × 10⁹, so that the total number of galaxies per redshift bin was 2 × 10⁸. We also took the mean galaxy radius to be ϑ = 0.32″ (1.55 × 10⁻⁶ rad). This value is the mean half-light radius of galaxies from the Euclid Flagship mock data (Euclid Collaboration 2022a).

To evaluate the impact of the local Universe correction, we used a smoothing scale of 120 h⁻¹ Mpc. This is the primary value used in Hall (2020), as it is just large enough for the local overdensity to be linear, while still possessing full-sky spherical coverage within the 2M++ galaxy redshift catalogue (Lavaux & Hudson 2011) used to measure the local overdensity. Accordingly, δ₀(R = 120 h⁻¹ Mpc) = 0.045 and δ₀(R = 120 h⁻¹ Mpc)/σ(R = 120 h⁻¹ Mpc) = 0.85. Different choices of smoothing scale result in different values of the local overdensity, δ₀, with some choices consistent with zero. As the LU effect scales linearly with the local overdensity, we do not rule out the LU bias as being exactly zero. Furthermore, it might be expected that the LU effect is mostly subsumed within the SSC uncertainty for any choice of smoothing scale, since it is mostly affected by modes that are outside of the survey. Detailed investigation of these points is beyond the scope of the paper, and our intention here is merely to assess how much bias would result from a nominal amplitude for the local overdensity combined with a fiducial implementation of the SSC covariance. In this work, we compared the magnitude of the studied correction terms to the Gaussian sample variance, ΔC_ℓ/C_ℓ. This was calculated according to Kaiser (1992), and took the form

$$\begin{array}{c}\hfill \mathrm{\Delta }{C}_{\ell }/C_{\ell }=\sqrt{2}{\left[f_{\mathrm{sky}}(2\ell +1)\right]}^{-1/2}.\end{array}$$(61)

4. Results and discussion

This section presents and discusses the computed values for the studied corrections. First, we show the magnitudes of the correction relative to the magnitude of the cosmic shear angular power spectra, and compare them to the sample variance. Throughout this section, we include the SLC and SO-SLC terms; however, due to the lingering uncertainty around the magnitude of the source number density weighting term, we also show results without these terms in Appendix A.

In Fig. 1, the magnitudes of the reduced shear, source-lens clustering, magnification bias, source obscuration, local Universe effect, and non-flat Universe corrections (relative to the angular power spectra0 are shown. The combined correction is also shown, and we note that this is the sum of the signed values of the corrections (rather than the absolute values) that are shown here for comparison. These terms are displayed for the auto-correlations of four redshift bins across the survey’s range; specifically, bins 1, 4, 7, and 10. These particular bins are presented for illustrative purposes, and the remaining bins and cross-correlations display consistent trends. It can be seen that all corrections are typically below sample variance both individually and when combined, with the exception of at small physical scales at the lowest and highest redshifts. A noteworthy detail from this figure is that while individual corrections are either consistently well below sample variance (or higher at low redshifts and reduce significantly at high redshifts or vice versa), the total magnitude of the corrections is consistently high. In Fig. 1 we compare changes in the power spectrum to the “cosmic variance” term. However, this is only for visualisation purposes. The overall result is computed using the covariance as described, that results in specific bias values.

Fig. 1. Absolute magnitudes of the reduced shear, source-lens clustering, magnification bias, local Universe, source obscuration, and non-flat Universe corrections to the shear angular power spectra, relative to those angular power spectra, for Euclid. The corrections to the angular power spectra for four redshift bin auto-correlations are shown as representative examples, spanning the redshift range of Euclid. The remaining auto and cross-correlations exhibit the same patterns. The absolute value of the signed sum of the corrections is also shown. These are all compared to the sample variance, calculated according to Eq. (61). Notably, while the magnitudes of individual corrections are either higher at lower redshifts or vice-versa, the magnitude of the sum of the corrections is consistently high. Additionally, the cross-terms between source obscuration, and magnification bias and source-lens clustering are multiple orders of magnitude below other terms and sample variance, suggesting they are negligible. The remainder of the terms are typically of similar magnitudes across redshifts, suggesting they must all be accounted for. We note that these magnitudes are for both the ΛCDM and w0waCDM cases, as the choice of fiducial values for the latter matches the former, and that the non-flat Universe correction here has been computed for a cosmology with ΩK = 0.05, whilst other corrections are when ΩK = 0. The markers for the SO-MB, SO-SLC, SO, Non-flat Universe, and total lines are only used to distinguish those from the other terms, and do not have any other significance. The symbols (points) are only included to allow a reader to distinguish the lines (in particular if printing in gray-scale) and do not indicate the ℓ-modes where a computation was made; all quantities were evaluated for 200 ℓ-bands, logarithmically spaced, with an ℓmin of 10, and an ℓmax of 5000.

Furthermore, another immediately noticeable feature is that the source obscuration cross terms with magnification bias and source-lens clustering are always multiple orders of magnitude smaller than the sample variance; these are typically the other terms as well; suggesting that these cross-terms are negligible.

Despite the fact that these terms are generally below sample variance because they make contributions consistently across ℓ-modes, they can still cause significant biases in inferred cosmological parameters. The bias in an estimated parameter resulting from neglecting a systematic effect is typically considered significant if it exceeds 25% of the 1σ uncertainty on that parameter (Taylor et al. 2007). This is because at that point, the biased and unbiased 1σ confidence contours overlap by less than 90%.

The predicted cosmological parameter biases resulting from neglecting all corrections except for the non-flat Universe term are stated in Tables 5 and 6, for the ΛCDM and w₀w_aCDM cases, respectively. These are also represented visually in Fig. 2. The biases from the non-flat Universe correction (which requires a fiducial cosmology with non-zero curvature) are stated in Table 7, for both choices of cosmology. From these tables we see that the source-lens clustering, magnification bias, source obscuration, and local Universe terms are individually significant in the ΛCDM case, while instead only the source-lens clustering, magnification bias, and source obscuration corrections are individually of concern in the w₀w_aCDM case. This difference is likely due to the presence of the variable dark energy parameters in the w_ow_aCDM scenario reducing sensitivity to scales where the local Universe term is important. Of these, the source-lens clustering term is particularly concerning, as for this term all but two of the parameters have significant biases.

Fig. 2. Stacked bar chart of cosmological parameter biases resulting from the studied higher-order effects, for the flat ΛCDM case (left) of Table 5, and the flat w0waCDM case (right) of Table 6. The non-flat Universe term is not shown here, due to the different cosmology. Biases are presented here as a fraction of the 1σ parameter uncertainty. A bias is non-negligble if its absolute value reaches or exceeds 0.25σ. ‘RS’ denotes the reduced shear correction, ‘SLC’ is the source-lens clustering term, ‘MB’ is the magnification bias correction, ‘SO’ is the two-point source obscuration correction, ‘SO-MB’ and ‘SO-SLC’ are the source obscuration-magnification bias and source-lens clustering cross terms respectively, and ‘LU’ is the local Universe correction The segments with the dashed outlines show the total parameter biases from these corrections for each parameter.

Also shown in these tables are the combined biases when all of the individually significant corrections are taken into account, and when all of these corrections are taken into account. The full total is also shown in Fig. 2. Owing to the fact that some biases are additive, while others are subtractive, the total biases in the ΛCDM scenario are, in fact, less severe than some of the individual ones; in particular, source-lens clustering. However, the totals are still significant and because they do not strongly resemble any one of the biases uniquely, multiple terms must still be computed. In the w₀w_aCDM case, the magnification bias no longer suppresses the source-lens clustering term; in fact it is adding to it instead, which means that the total biases in this case are more severe than the individual ones. This change likely occurs due to the dark energy terms increasing sensitivity to scales where the opposite component of the magnification bias (e.g. decrease in galaxy number density due to dilution rather than increase due to increased flux) is dominant.

At inference time, the computational load can be reduced by noting that only the individually significant terms (source-lens clustering, magnification bias, source obscuration, and local Universe) need to be computed in both cases, because this total does not significantly differ from the full total. Although, given that the reduced shear correction is also obtained at no additional cost when computing the magnification bias correction, we recommend including this as well. The bias in the two-parameter confidence contours resulting from neglecting the combined effect of the significant biases is shown in Figs. 3 and 4 for the ΛCDM and w₀w_aCDM scenarios, respectively. These figures display the contours in the case where corrections have been made and where they have not been made. As with Tables 5 and 6, we see that the cumulative corrections must be accounted for.

Fig. 3. Projected 1σ and 2σ 2-parameter uncertainty contours for Euclid under a ΛCDM cosmology, with and without correcting for the source-lens clustering, magnification bias, source obscuration, and local Universe terms. These are predicted using the Fisher matrix formalism, using the cosmology specified in Table 3, in the case where ΩK = 0 and is kept fixed. The true location of the constraints is denoted by the blue-dashed contours, while the biased locations if the corrections are not made are given by the solid-gold contours. Significant biases are predicted for Ωm, Ωb, h, and σ8, and their values can be found in Table 5

Fig. 4. Projected 1σ and 2σ 2-parameter uncertainty contours for Euclid, with and without correcting for source-lens clustering, magnification bias, and source obscuration. These are predicted using the Fisher matrix formalism, for the w0waCDM cosmology specified in Table 3, in the case when ΩK = 0 and is kept fixed. The true location of the constraints is denoted by the blue-dashed contours, while the biased locations if the corrections are not made are given by the solid-orange contours. Significant biases are predicted for Ωm, Ωb, ns, σ8, w0, and wa, and their values can be found in Table 6.

We note that many of the modelling specifics used here, for example, the slope of the luminosity function or the smoothing scale of local overdensity, will need to be determined directly from the Euclid survey itself for self-consistency when these corrections are computed at inference time. Accordingly, it would not be meaningful to place constraints on them here, due to the variable survey specifics.

Additionally, we note that the true value of the local density contrast and, accordingly, the radius of the smoothing kernel used to calculate it are still open questions. Accordingly, there is a large uncertainty in the local Universe correction which cannot be meaningfully constrained and it is possible that it may even be zero. Accurate measurements of the local density contrast are required for this.

Similarly, the source obscuration correction as computed in this work represents a worst-case scenario where ’every’ galaxy in the foreground of a given redshift slice has an overlap. This represents an upper limit on the bias from this effect; however, in practice, the number of sources with overlap will be a smaller fraction. Accordingly, it is likely that source obscuration will not result in significant biases for true Euclid observations, particularly if a robust mitigation strategy is employed.

However, the dominant source of quantifiable modelling uncertainty comes from the modelling of the matter bispectrum and given that the bispectrum is not currently well constrained by observations, this model is likely to continue to evolve. Accordingly, it is important to constrain the impact of a change in bispectrum model on these terms. To date, three widely used matter bispectrum models have been produced (Scoccimarro & Couchman 2001; Gil-Marín et al. 2012; Takahashi et al. 2020). As each subsequent model has become more complex and improved upon the accuracy of its predecessor, comparing correction magnitudes using each of these models would not realistically constrain the uncertainty from the bispectrum model.

Instead, it is useful to set a threshold around the latest of these models (Takahashi et al. 2020), within which any change in the model must be contained to avoid producing a significant change in the correction terms. We do this by determining the minimum fractional increase or decrease in the matter bispectrum required across all triangle configurations, for each correction individually, to cause a significant change in any of the cosmological parameter biases, in the w₀w_aCDM case. This corresponds to a change of ±0.25σ in any one of the biases. Given that we are placing multiplicative limits on the change in the bispectrum, for each correction the smallest limits are found when considering the cosmological parameter that already has the largest bias. The resulting limits are stated in Table 8, alongside the parameter which would see the corresponding significant change in its bias. From this, we see that the corrections most susceptible to a change in the bispectrum model are the source-lens clustering and magnification bias terms, as neglecting these already creates the most significant biases individually. Additionally, the two source obscuration cross terms are the least sensitive, requiring a change of an order-of-magnitude. This further reinforces the fact that these terms are safely negligible. We stress that these multiplicative limits are not exhaustive cut-offs on when a bispectrum model would cause a significant change, because a model with a sufficiently large change for only a select sub-range of scales or configurations could still cause a significant difference in the parameter biases. We recommend an explicit revaluation with any future updated bispectrum models, should they non-trivially exceed these thresholds frequently.

Another consideration is how the inclusion of the studied effects would affect the size of the cosmological parameter uncertainty constraints themselves. In this work, we do not explicitly calculate this resulting change. However, it has previously been shown that for the bispectrum-dependent terms, even corrections that cause biases of greater than 1σ result in negligible changes to the uncertainty constraints (Shapiro 2009; Deshpande et al. 2020a).

5. Conclusions

In this investigation, we have examined the higher-order corrections to the cosmic shear angular power spectra that must be modelled when performing inference with Euclid. By first reviewing the literature, we identified 24 correction terms and gathered representative values to facilitate comparison. From these, we identified six corrections which were potentially important for Euclid and evaluated them explicitly. These were: the relaxation of the reduced shear approximation, the source-lens clustering correction, the magnification bias correction, the source obscuration correction, the local Universe correction, and the non-flat Universe correction.

After calculating these corrections, we used the Fisher matrix formalism to predict the biases in cosmological parameter biases if each of these terms were to be neglected, in order to identify which ones are necessary to be modelled for Euclid. This was done for two scenarios: a ΛCDM cosmology, and a w₀w_aCDM cosmology. For the first of these scenarios, we found that the source-lens clustering, magnification bias, source obscuration, and local Universe terms were significant, while for the second case we found that the source-lens clustering, magnification bias, and source obscuration corrections each produced significant biases in multiple parameters individually. The source-lens clustering term was noted as being of particular concern as multiple biases approached or exceeded 1σ. However, in the ΛCDM case we found that when the biases are combined, they frequently suppressed each other, leading the total bias to be lower than many of the individual biases. Despite this, the total biases were still significant, and did not strongly represent the exact biases from any correction individually. In the w₀w_aCDM case, the total of the three biases was higher than the individual terms. Accordingly, we recommend that the source-lens clustering, magnification bias, source obscuration, and local Universe corrections are all taken into account when modelling the shear angular power spectra for Euclid. Additionally, given that the reduced shear correction is obtained at no additional cost when computing the magnification bias correction, we recommend that this too be included.

In the case the source number density weighting term (a contribution to source lens clustering) is zero, we find that for the ΛCDM case we find that none of the overall total bias-to-uncertainty ratios exceed 0.25σ, however this only occurs if all corrections are applied or not applied. Hence, we should either include all corrections or none, but including only some could result in biased results. In the w₀w_aCDM case we find significant biases could be caused by neglecting any of terms, but, in particular, the magnification bias.

To provide some constraints on the predictive ability of this work, we quantified how much the matter bispectrum would have to change by in order to illicit a significant change in the biases predicted in this work. We identified that a ∼20 − 40% increase or decrease in the amplitude of the matter bispectrum at all scales and for all triangle configurations would be required in such a case.

We note that this work did not investigate the impact of higher-order corrections on the IA spectra. Typically, even for effects which, when neglected, produce high biases, namely, ∼O(1σ), the corresponding corrections to the IA spectra cause negligible biases (Deshpande et al. 2020a). Accordingly, an explicit evaluation should not be necessary. Furthermore, we did not consider the impact of baryonic feedback on the bispectrum and, therefore, on its dependent corrections. The impact of this remains poorly understood, with inconsistent findings from different simulations (Semboloni et al. 2013; Barreira et al. 2019). Accordingly, this is beyond the scope of this investigation. However, we note that given the magnitude of change required to the matter bispectrum in order to cause a significant change in the cosmological parameter biases, it is unlikely baryonic feedback would significantly alter the predictions of this investigation.

Given that we find it is necessary to include these higher order terms in the modelling of the shear power spectra, the optimal strategy to do so remains an open question. It has repeatedly been shown that computing these corrections for just one cosmology is relatively time consuming (Deshpande et al. 2020a; Duncan et al. 2022), posing a serious challenge to carrying out computations at the inference time. While it may be possible to sufficiently optimise the required evaluation code, alternate strategies may also prove useful. Scale-cutting techniques such as k-cut cosmic shear (Taylor et al. 2018a) have been shown to mitigate the need to make such corrections without significantly compromising the constraining power of Stage IV surveys (Deshpande et al. 2020b). Alternatively, emulation has recently become a popular tool in cosmology for reducing computation times at the moment of inference by replacing analytical models with emulators (see e.g. recent work emulating the matter power spectrum within Spurio Mancini et al. 2022). Emulators could also be developed directly for these correction terms, or intermediate quantities such as the matter or convergence bispectra. Further work on ℓ-mode weighting and/or k-cut cosmic shear approaches are required to optimise the trade-off between accuracy, precision, and computational expense.

Furthermore, while our analysis here is limited to the angular power spectra, significant corrections for this statistic are also likely to be significant for the two-point correlation function. In fact, due to the mode-mixing that occurs when transforming the power spectra to correlation functions, the effect of the discussed approximations is likely to be more severe. This, combined with the sensitivity of the correlation function to higher ℓ-modes and the additional approximations required (e.g. the flat Hankel transform; Kitching et al. 2017), means that if correlation functions were to be used a similar but separate study would be required to demonstrate modelling of the correlation function to higher order corrections.