# Euclid preparation

## XXXVI. Modelling the weak lensing angular power spectrum

Euclid Collaboration: A. C. Deshpande¹, T. Kitching¹, A. Hall², M. L. Brown³, N. Aghanim⁴, L. Amendola⁵, S. Andreon⁶, N. Auricchio⁷, M. Baldi⁸,⁷,⁹, S. Bardelli⁷, R. Bender¹⁰,¹¹, D. Bonino¹², E. Branchini¹³,¹⁴, M. Brescia¹⁵, J. Brinchmann¹⁶, S. Camera¹⁷,¹⁸,¹², G. P. Candini¹, V. Capobianco¹², C. Carbone¹⁹, V. F. Cardone²⁰,²¹, J. Carretero²²,²³, S. Casas²⁴, F. J. Castander²⁵,²⁶, M. Castellano²⁰, S. Cavuoti²⁷,²⁸, A. Cimatti²⁹, R. Cledassou³⁰,³¹,*, G. Congedo², C. J. Conselice³, L. Conversi³²,³³, L. Corcione¹², F. Courbin³⁴, H. M. Courtois³⁵, M. Cropper¹, A. Da Silva³⁶,³⁷, H. Degaudenzi³⁸, M. Douspis⁴, F. Dubath³⁸, C. A. J. Duncan³⁹,³, X. Dupac³², M. Farina⁴⁰, S. Farrens⁴¹, S. Ferriol⁴², P. Fosalba²⁶,²⁵, M. Frailis⁴³, E. Franceschi⁷, M. Fumana¹⁹, S. Galeotta⁴³, B. Garilli¹⁹, B. Gillis², C. Giocoli⁷,⁹, A. Grazian⁴⁴, F. Grupp¹⁰,¹¹, S. V. H. Haugan⁴⁵, H. Hoekstra⁴⁶, W. Holmes⁴⁷, A. Hornstrup⁴⁸,⁴⁹, P. Hudelot⁵⁰, K. Jahnke⁵¹, E. Keihänen⁵², S. Kermiche⁵³, M. Kilbinger⁴¹, M. Kunz⁵⁴, H. Kurki-Suonio⁵⁵,⁵⁶, S. Ligori¹², P. B. Lilje⁴⁵, V. Lindholm⁵⁵,⁵⁶, I. Lloro⁵⁷, E. Maiorano⁷, O. Mansutti⁴³, O. Marggraf⁵⁸, K. Markovic⁴⁷, N. Martinet⁵⁹, F. Marulli⁸,⁷,⁹, R. Massey⁶⁰, S. Mei⁶¹, Y. Mellier⁶²,⁵⁰, M. Meneghetti⁷,⁹, G. Meylan⁵³, L. Moscardini⁸,⁷,⁹, S.-M. Niemi⁶³, J. W. Nightingale⁶⁰, T. Nutma⁴⁶,⁶⁴, C. Padilla²², S. Paltani³⁸, F. Pasian⁴³, K. Pedersen⁶⁵, V. Pettorino⁴¹, S. Pires⁶⁶, G. Polenta⁶⁷, J. Pollack⁶⁸,⁶¹, M. Poncet³⁰, L. A. Popa⁶⁹, F. Raison¹⁰, A. Renzi⁷⁰,⁷¹, J. Rhodes⁴⁷, G. Riccio²⁷, E. Romelli⁴³, M. Roncarelli⁷, E. Rossetti⁷², R. Saglia¹¹,¹⁰, D. Sapone⁷³, B. Sartoris¹¹,⁴³, P. Schneider⁵⁸, T. Schrabback⁷⁴,⁵⁸, A. Secroun⁵³, G. Seidel⁵¹, S. Serrano²⁶,⁷⁵, C. Sirignano⁷⁰,⁷¹, G. Sirri⁹, L. Stanco⁷¹, P. Tallada-Crespi⁷⁶,²³, A. N. Taylor², I. Tereno³⁶,⁷⁷, R. Toledo-Moreo⁷⁸, F. Torradeflot⁷⁶,²³, I. Tutusaus⁷⁹, E. A. Valentijn⁶⁴, L. Valenziano⁷,⁹, T. Vassallo⁴³, Y. Wang⁸⁰, J. Weller¹¹,¹⁰, A. Zacchei⁴³,⁸¹, G. Zamorani⁷, J. Zoubian⁵³, E. Zucca⁷, A. Boucaud⁶¹, E. Bozzo³⁸, C. Colodro-Conde⁸², D. Di Ferdinando⁹, G. Fabbian⁸³,⁸⁴, J. Graciá-Carpio¹⁰, N. Mauri²⁹,⁹, V. Scottez⁵⁰,⁸⁵, M. Tenti⁸⁶, Y. Akrami⁸⁷,⁸⁸,⁸⁹,⁹⁰,⁹¹, C. Baccigalupi⁹²,⁸¹,⁴³,⁹³, A. Balaguera-Antolínez⁸²,⁹⁴, M. Ballardini⁹⁵,⁹⁶,⁷, F. Bernardeau⁹⁷,⁶², A. Biviano⁴³,⁸¹, A. Blanchard⁷⁹, A. S. Borlaff⁹⁸, C. Burigana⁹⁵,⁹⁹,⁸⁶, R. Cabanac⁷⁹, A. Cappi⁷,¹⁰⁰, C. S. Carvalho⁷⁷, G. Castignani⁸,⁷, T. Castro⁴³,⁹³,⁸¹, K. C. Chambers¹⁰¹, A. R. Cooray¹⁰², J. Coupon³⁸, S. Davini¹⁰³, S. de la Torre⁵⁹, G. De Lucia⁴³, G. Desprez³⁸,¹⁰⁴, H. Dole⁴, J. A. Escartin¹⁰, S. Escoffier⁵³, I. Ferrero⁴⁵, F. Finelli⁷,⁸⁶, J. Garcia-Bellido⁸⁷, K. George¹⁰⁵, F. Giacomini⁹, G. Gozaliasl⁵⁵, H. Hildebrandt¹⁰⁶, J. J. E. Kajava¹⁰⁷, V. Kansal⁶⁶, C. C. Kirkpatrick⁵², L. Legrand⁵⁴, A. Loureiro²,⁹¹, J. Macias-Perez¹⁰⁸, M. Magliocchetti⁴⁰, G. Mainetti¹⁰⁹, R. Maoli¹¹⁰,²⁰, M. Martinelli²⁰,²¹, C. J. A. P. Martins¹¹¹,¹⁶, S. Matthew², L. Maurin⁴, R. B. Metcalf⁸,⁷, P. Monaco¹¹²,⁴³,⁹³,⁸¹, G. Morgante⁷, S. Nadathur¹¹³, A. A. Nucita¹¹⁴,¹¹⁵,¹¹⁶, L. Patrizii⁹, A. Peel³⁴, M. Pöntinen⁵⁵, V. Popa⁶⁹, C. Porciani⁵⁸, D. Potter¹¹⁷, A. Pourtsidou²,¹¹⁸, P. Reimberg⁵⁰, Z. Sakr¹¹⁹,⁵,⁷⁹, A. G. Sánchez¹⁰, A. Schneider¹¹⁷, E. Sefusatti⁴³,⁹³,⁸¹, M. Sereno⁷,⁹, A. Shulevski⁴⁶,⁶⁴, A. Spurio Mancini¹, J. Steinwagner¹⁰, R. Teyssier¹²⁰, M. Viel⁹²,⁸¹,⁴³,⁹³, I. A. Zinchenko¹¹, and P. Fleury⁹⁷

(Affiliations can be found after the references)

Received 9 February 2023 / Accepted 12 January 2024

## ABSTRACT

This work considers which higher order modeling effects on the cosmic shear angular power spectra must be taken into account for Euclid. We identified the relevant terms and quantified their individual and cumulative impact on the cosmological parameter inferences from Euclid. We computed the values of these higher order effects using analytic expressions and calculated the impact on cosmological parameter estimations using the Fisher matrix formalism. We reviewed 24 effects and determined the ones that potentially need to be accounted for, namely: the reduced shear approximation, magnification bias, source-lens clustering, source obscuration, local Universe effects, and the flat Universe assumption. After computing these effects explicitly and calculating their cosmological parameter biases, using a maximum multipole of $l = 5000$, we find that the magnification bias, source-lens clustering, source obscuration, and local Universe terms individually produce significant (>0.25σ) cosmological biases in one or more parameters; accordingly, these effects must be accounted for and warrant further investigation. In total, we find biases in $\Omega_m$, $\Omega_b$, $h$, and $\sigma_8$ of 0.73σ, 0.28σ, 0.25σ, and -0.79σ, respectively, for the flat $\Lambda$CDM. For the $w_0w_a$ CDM case, we found biases in $\Omega_m$, $\Omega_b$, $h$, $n_s$, $\sigma_8$, and $w_a$ of 1.49σ, 0.35σ, -1.36σ, 1.31σ, -0.84σ, and -0.35σ, respectively. These are increased relative to the $\Lambda$CDM due to additional degeneracies as a function of redshift and scale.

**Key words.** gravitation – gravitational lensing: weak – cosmological parameters – dark energy

<footnote>* Deceased.</footnote>
# 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 $\sigma_8$ and $\Omega_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<fnref n="1" /> (Laureijs et al. 2011), Nancy Grace Roman<fnref n="2" /> (Akeson et al. 2019), and the Vera C. Rubin Observatory<fnref n="3" /> (LSST Science Collaboration 2009) will achieve more than an order-of-magnitude increase in precision over existing surveys (Euclid Collaboration 2020, hereafter, 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<fnref n="4" />.

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)(\ell\theta)^2$, where $\ell$ is the angular multipole of the power spectrum and $\theta$ 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 $-(\ell/1.2 \times 10^6)^2$ 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

<footnote n="1">¹ https://www.euclid-ec.org/</footnote>

<footnote n="2">² https://roman.gsfc.nasa.gov/</footnote>

<footnote n="3">³ https://www.lsst.org/</footnote>

<footnote n="4">⁴ Every paper made various slightly differing assumptions regarding survey area, depth, and number density; in all cases, the details can be found in the references. In the case that only correction function analyses were available, not power spectrum, we performed a Hankel transform over the correlation functions of the graphs over the angular range available.</footnote>
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 signif-
icantly smaller than the reduced shear correction. This serves
as a good comparison because this correction has been consis-
tently 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 addi-
tional 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 cor-
rection. Additionally, the cosmological parameter biases result-
ing 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 illus-
trated 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 com-
bined 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 negligi-
ble” stem from corrections to the theoretical expressions describ-
ing 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 poten-
tials 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 quan-
tity 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 rela-
tion 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 investi-
gated in Krause & Hirata (2010) (Sect. 3.3), who find that the
impact is three orders of magnitude smaller than the shear power
spectrum<fnref n="5" />.

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 dis-
crepancy between the true perturbed path of a photon from
source to observer compared to the mathematically conve-
nient 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 evo-
lutionary stages of the LSS at different times, necessitating a
correction in the two-point statistics. The remaining two cor-
rections 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 gravita-
tional 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 negligi-
ble without requiring explicit calculation. This additional cor-
rection term arises from the fact that, for larger sources, the
image distortion consists of both shear and a higher-order com-
ponent 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 mod-
elling 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. Deter-
mining 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 cos-
mic shear angular power spectrum calculation. The expressions

<footnote n="5">⁵ In fact there are two definitions of ellipticity: third eccentricity and third flattening (see e.g. Viola et al. 2014) that relate the observed ellipticity to the shear in different ways. Krause & Hirata (2010) find that for third flattening the correction to the power spectrum is zero (since the moments of the third flattening are exactly the moments of the reduced shear as shown by Seitz & Schneider 1997), but that for third eccentricity the effect is non-zero.</footnote>
<table><caption>**Table 1.** List of higher-order correction terms to the shear angular power spectrum resulting from relaxing approximations.</caption><thead><tr><th rowspan="2">Correction</th><th rowspan="2">Reference</th><th colspan="3">Percentage of sample variance</th></tr><tr><th>ℓ = 10</th><th>ℓ = 100</th><th>ℓ = 1000</th></tr></thead><tbody><tr><td><strong>Source-lens clustering</strong></td><td>Yu et al. (2015)</td><td>N/A</td><td>6.0</td><td>57.2</td></tr><tr><td><strong>Reduced shear + magnification bias</strong></td><td>Deshpande et al. (2020a)</td><td>0.4</td><td>1.8</td><td>15.3</td></tr><tr><td>Post-Limber reduced shear</td><td>Deshpande & Kitching (2020)</td><td>0.2</td><td>0.6</td><td>1.9</td></tr><tr><td>Non-linear ellipticity-shear relation</td><td>Krause & Hirata (2010)</td><td>0.2</td><td>0.6</td><td>1.9</td></tr><tr><td>Limber + flat-sky</td><td>Kitching et al. (2017)</td><td>9.8</td><td>3.0</td><td>1.0</td></tr><tr><td>Local Universe effects</td><td>Hall (2020)</td><td>7.8</td><td>24.2</td><td>N/A</td></tr><tr><td>Higher-order reduced shear</td><td>Krause & Hirata (2010)</td><td>O(φ<sup>4</sup>)</td><td>O(φ<sup>4</sup>)</td><td>O(φ<sup>4</sup>)</td></tr><tr><td>Time delay-lens coupling</td><td>Bernardeau et al. (2010)</td><td>O(φ<sup>4</sup>)</td><td>O(φ<sup>4</sup>)</td><td>O(φ<sup>4</sup>)</td></tr><tr><td>Deflection-deflection coupling</td><td>Bernardeau et al. (2010)</td><td>O(φ<sup>4</sup>)</td><td>O(φ<sup>4</sup>)</td><td>O(φ<sup>4</sup>)</td></tr><tr><td>Born approximation</td><td>Cooray & Hu (2002)</td><td>O(φ<sup>4</sup>)</td><td>O(φ<sup>4</sup>)</td><td>O(φ<sup>4</sup>)</td></tr><tr><td>Lensing by dark energy pressure</td><td>Simpson et al. (2010)</td><td>O(φ<sup>4</sup>)</td><td>O(φ<sup>4</sup>)</td><td>O(φ<sup>4</sup>)</td></tr><tr><td>Second-order speed-of-light</td><td>Cuesta-Lazaro et al. (2018)</td><td>&lt;0.01</td><td>&lt;0.01</td><td>&lt;0.01</td></tr><tr><td>Temporal-Born approximation</td><td>Cuesta-Lazaro et al. (2018)</td><td>&lt;0.01</td><td>&lt;0.01</td><td>&lt;0.01</td></tr><tr><td>Finite-beam corrections</td><td>Fleury et al. (2019)</td><td>&lt;0.01</td><td>&lt;0.01</td><td>&lt;0.01</td></tr><tr><td>Doppler-shift</td><td>Deshpande & Kitching (2021)</td><td>&lt;0.01</td><td>&lt;0.01</td><td>&lt;0.01</td></tr><tr><td>Unequal-time correlators</td><td>Kitching & Heavens (2017)</td><td>&lt;0.01</td><td>&lt;0.01</td><td>&lt;0.01</td></tr><tr><td>Sachs-Wolfe effect</td><td>Cuesta-Lazaro et al. (2018)</td><td>&lt;0.01</td><td>&lt;0.01</td><td>&lt;0.01</td></tr><tr><td>Integrated Sachs-Wolfe effect</td><td>Cuesta-Lazaro et al. (2018)</td><td>&lt;0.01</td><td>&lt;0.01</td><td>&lt;0.01</td></tr><tr><td>Flexion correction</td><td>Schneider & Er (2008)</td><td>N/A</td><td>N/A</td><td>N/A</td></tr><tr><td><strong>Flat-geometry assumption</strong></td><td>Taylor et al. (2018b)</td><td>-2.0</td><td>-6.0</td><td>-19.1</td></tr><tr><td><strong>Source obscuration</strong></td><td>Hartlap et al. (2011)</td><td>-2.0</td><td>-6.0</td><td>-19.1</td></tr><tr><td><strong>Spatially-varying survey depth</strong></td><td>Heydenreich et al. (2020)</td><td>-5.9</td><td>-18.1</td><td>-57.2</td></tr></tbody></table>

**Notes.** To illustrate their typical sizes and facilitate comparison, the values of these correction terms at redshift $z \sim 1$ are also stated here. These are provided for $ℓ$-modes 10, 100, and 1000, and as a percentage of sample variance i.e. $100\delta C_ℓ/(\text{sample variance})$. Sample variance here is calculated using the definition of Kaiser (1992), see Eq. (61). These values are taken from the available literature for the case of Stage IV cosmic shear experiments. Unavailable values are represented by 'N/A'. Corrections which have a functional form that is fourth-order in lensing potential, and therefore sub-dominant, are denoted by $O(φ^4)$. The references provided refer to the values stated where available, or to the work describing the correction where explicit values are not available. For the case where the absolute value of an effect is more than four orders of magnitude smaller than the sample variance, it is denoted by $<0.01$. The sign denotes whether the effect is to decrease the power (a negative sign) or to increase (no sign). The rows are ordered in decreasing amplitude for $ℓ = 1000$. These numbers were read from graphs published in the referenced papers, due to the inherent inaccuracy of this approach we quote only one decimal place; in the case that only correlation functions were provided these were converted to power spectrum results using a Hankel transform over the quoted angular range. We note that the Table is ordered by the magnitude of the impact of the effects at $ℓ = 1000$, not the absolute value of the impact, therefore those effects that cause a negative change to the power spectra are at the bottom of the Table. The effects in bold are those whose approximations are investigated in this paper. Where values are the same this is due to rounding and limitations of transcribing values from plots in papers; hence, these values should only be taken as indicative of the amplitude of the effects.

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:

$$g^{\alpha}(\theta) = \frac{\gamma^{\alpha}(\theta)}{1 - \kappa(\theta)}, \qquad (1)$$

where $\theta$ is the position of the galaxy on the sky, $\gamma$ is the spin-2 shear with index $\alpha \in \{1, 2\}$ that describes the anisotropic stretching that turns circular distributions of light elliptical, and $\kappa$ is the convergence – responsible for the isotropic change in the size of the image. Since in the weak lensing regime $|\kappa| \ll 1$, it is standard practice to make the reduced shear approximation, whereby

$$g^{\alpha}(\theta) \approx \gamma^{\alpha}(\theta). \qquad (2)$$

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

$$ \kappa_i(\theta) = \int_0^{\chi_h} d\chi \, \delta[S_K(\chi)\theta, \chi] W_i(\chi), \quad (3) $$

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

$$ S_K(\chi) = \begin{cases} |K|^{-1/2} \sin(|K|^{-1/2}\chi) & K > 0 \text{ (closed Universe),} \\ \chi & K = 0 \text{ (flat Universe),} \\ |K|^{-1/2} \sinh(|K|^{-1/2}\chi) & K < 0 \text{ (open Universe).} \end{cases} \quad (4) $$

We recall that for the quantity $\delta[S_K(\chi)\theta, \chi]$ in Eq. (3) the second $\chi$ indicates not only that there is an evaluation at a comoving radius, $\chi$, but also at a conformal time, $\eta = \eta_0 - \chi$, meaning that all the integration over $\chi$ in this paper are performed down the background light cone.

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

$$ W_i(\chi) = \frac{3}{2} \Omega_m \frac{H_0^2}{c^2} \frac{S_K(\chi)}{a(\chi)} \int_{\chi}^{\chi_h} d\chi' n_i(\chi') \frac{S_K(\chi' - \chi)}{S_K(\chi')}, \quad (5) $$

which is dependent on the dimensionless present-day matter density of the Universe, $\Omega_m$, the speed of light in a vacuum, $c$, the Hubble constant, $H_0$, the scale factor of the Universe, $a(\chi)$, and the probability distribution of galaxies within the redshift bin, $i$, namely, $n_i(\chi)$.

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:

$$ \tilde{\gamma}_i^\alpha(\ell) = T^\alpha(\ell) \tilde{\kappa}_i(\ell), \quad (6) $$

where $\ell$ is the spherical-harmonic conjugate of $\theta$, with magnitude $\ell$ and angular component $\phi_\ell$. The functions $T^\alpha$ are two trigonometric weighting functions corresponding to each of the shear components. These take the form:

$$ T^1(\ell) = \cos(2\phi_\ell), \quad (7) $$

$$ T^2(\ell) = \sin(2\phi_\ell). \quad (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:

$$ \tilde{E}_i(\ell) = \sum_{\alpha} T^{\alpha}(\ell) \tilde{\gamma}_i^{\alpha}(\ell), \quad (9) $$

$$ \tilde{B}_i(\ell) = \sum_{\alpha} \sum_{\beta} \varepsilon^{\alpha\beta} T^{\alpha}(\ell) \tilde{\gamma}_i^{\beta}(\ell). \quad (10) $$

Here, the summations are over the shear components, and $\varepsilon^{\alpha\beta}$ is the Levi-Civita symbol in the 2D case; such that: $\varepsilon^{11} = \varepsilon^{22} = 0$ and $\varepsilon^{12} = -\varepsilon^{21} = 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:

$$ \langle \tilde{E}_i(\ell) \tilde{E}_j(\ell') \rangle = (2\pi)^2 \delta_D^{(2)}(\ell + \ell') C_{\ell;ij}^{\gamma\gamma}, \quad (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_D^{(2)}$ is the Dirac delta for 2D. Under the extended Limber approximation (LoVerde & Afshordi 2008), where $k = (\ell + 1/2)/S_K(\chi)$, the power spectra themselves are further defined as:

$$ C_{\ell;ij}^{\gamma\gamma} = \int_0^{\chi_h} d\chi \frac{W_i(\chi)W_j(\chi)}{S_K^2(\chi)} P_{\delta\delta}(k, \chi), \quad (12) $$

where $P_{\delta\delta}$ is the 3D matter power spectrum, and $k$ is the magnitude of the spatial momentum vector, $k$, which also shares the angular component, $\phi_\ell$. 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, $\epsilon$, can then, to first-order, be written as a combination of its underlying ellipticity in the absence of any cosmic shear or IA, $\epsilon^s$, the cosmic shear, $\gamma = \gamma^1 + i\gamma^2$, and the effect of IA, $\epsilon^l$ according to:

$$ \epsilon = \epsilon^s + \gamma + \epsilon^l. \quad (13) $$

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

$$ C_{\ell;ij}^{\epsilon\epsilon} = C_{\ell;ij}^{\gamma\gamma} + C_{\ell;ij}^{\gamma l} + C_{\ell;ij}^{l\gamma} + C_{\ell;ij}^{II} + N_{\ell;ij}^{\epsilon}, \quad (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 \le z_j$. The $C_{\ell;ij}^{\gamma l}$ 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}^{l\gamma}$ arise from the correlation between background shear and foreground IA, and the $C_{\ell;ij}^{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:

$$ C_{\ell;ij}^{l\gamma} = \int_0^{\chi_h} \frac{d\chi}{S_K^2(\chi)} [W_i(\chi)n_j(\chi) + n_i(\chi)W_j(\chi)] P_{\delta l}(k, \chi), \quad (15) $$

$$ C_{\ell;ij}^{II} = \int_0^{\chi_h} \frac{d\chi}{S_K^2(\chi)} n_i(\chi)n_j(\chi) P_{II}(k, \chi), \quad (16) $$

which, in a similar manner to the shear power spectra, are projections of 3D IA power spectra, $P_{\delta l}$ and $P_{II}$. Both of these are related to the matter power spectrum as follows:

$$ P_{\delta l}(k, \chi) = \left[ - \frac{\mathcal{A}_{lA} C_{lA} \Omega_m}{D(\chi)} \right] P_{\delta\delta}(k, \chi), \quad (17) $$

$$ P_{II}(k, \chi) = \left[ - \frac{\mathcal{A}_{lA} C_{lA} \Omega_m}{D(\chi)} \right]^2 P_{\delta\delta}(k, \chi), \quad (18) $$
where the product of AIA and CIA 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_{\ell;ij}^{\epsilon}$ in Eq. (14) is the shot noise term arising from the zero-lag autocorrelation of the un lensed, uncorrelated source ellipticity $\epsilon^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:

$$N_{\ell;ij}^{\epsilon} = \frac{\sigma_{\epsilon}^2}{\bar{n}_g/N_{\text{bin}}} \delta_{ij}^{K}, \qquad (19)$$

where $\sigma_{\epsilon}^2$ is the variance of the observed ellipticities in the survey, $\bar{n}_g$ is the surface density of galaxies in the survey, $N_{\text{bin}}$ is the survey's number of tomographic redshift bins, and $\delta_{ij}^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:

$$g^{\alpha}(\theta) = \gamma^{\alpha}(\theta) + (\gamma^{\alpha}\kappa)(\theta) + \mathcal{O}(\kappa^3). \qquad (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\langle\tilde{E}_i(\ell)\tilde{E}_j(\ell')\rangle$, are given by:

$$\begin{aligned} \delta\langle\tilde{E}_i(\ell)\tilde{E}_j(\ell')\rangle &= \sum_{\alpha} \sum_{\beta} T^{\alpha}(\ell)T^{\beta}(\ell')\langle\langle\tilde{\gamma}^{\alpha\kappa}_i(\ell)\tilde{\gamma}^{\beta}_{j}(\ell')\rangle\rangle \\ &\quad + T^{\alpha}(\ell')T^{\beta}(\ell)\langle\langle\tilde{\gamma}^{\alpha\kappa}_j(\ell')\tilde{\gamma}^{\beta}_{i}(\ell)\rangle\rangle \\ &= (2\pi)^2 \delta_D^{(2)}(\ell + \ell') \delta_{\ell;ij}^{RS}, \end{aligned} \qquad (21)$$

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

$$\delta C_{\ell;ij}^{RS} = \int_{0}^{\infty} \frac{d^2\ell'}{(2\pi)^2} \cos(2\phi_{\ell'}) B_{ij}^{KKK}(\ell, \ell', -\ell - \ell'), \qquad (22)$$

where we are always free to choose a coordinate system such that $\phi_{\ell} = 0$, and accordingly the correction only depends on the magnitude, $\ell$. It depends on the two-redshift convergence bispectrum, $B_{ij}^{KKK}$, 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(\chi)^4)$,

and given that typically $\chi W(\chi) \ll 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_{\delta\delta\delta}$. Under the Limber approximation, it takes the following form:

$$\begin{aligned} B_{ij}^{KKK}(\ell_1, \ell_2, \ell_3) &= B_{ii}^{KKK}(\ell_1, \ell_2, \ell_3) + B_{ijj}^{KKK}(\ell_1, \ell_2, \ell_3) \\ &= \int_0^{\chi h} \frac{d\chi}{S_K^4(\chi)} W_i(\chi) W_j(\chi) [W_i(\chi) + W_j(\chi)] \\ &\quad \times B_{\delta\delta\delta}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3, \chi). \end{aligned} \qquad (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):

$$n_i^{\text{obs}}(\theta, \chi) = n_i(\chi) [1 + \delta_i^g(\theta)]. \qquad (24)$$

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

$$\gamma_{\text{obs};i}^{\alpha}(\theta) = \gamma_{i}^{\alpha}(\theta) + \gamma_{i}^{\alpha}(\theta) \delta_{i}^{g}(\theta). \qquad (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}^{\text{SLC}}$, to the angular power spectra,

$$\delta C_{\ell;ij}^{\text{SLC}} = \int_{0}^{\infty} \frac{d^2 \ell'}{(2\pi)^2} \cos(2\phi_{\ell'}) B_{ij}^{\kappa\delta\gamma\kappa}(\ell, \ell', -\ell - \ell'), \qquad (26)$$

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

$$\begin{aligned} B_{ij}^{\kappa\delta\gamma\kappa}(\ell_1, \ell_2, \ell_3) &= B_{ii}^{\kappa\delta\gamma\kappa}(\ell_1, \ell_2, \ell_3) + B_{ijj}^{\kappa\delta\gamma\kappa}(\ell_1, \ell_2, \ell_3) \\ &= \int_0^{\chi h} \frac{d\chi}{S_K^4(\chi)} [b_i n_i(\chi) + b_j n_j(\chi)] W_i(\chi) W_j(\chi) \\ &\quad \times B_{\delta\delta\delta}(\mathbf{k}_1, \mathbf{k}_2, \mathbf{k}_3, \chi), \end{aligned} \qquad (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:

$$b_i = \sqrt{1 + \bar{z}_i}, \qquad (28)$$

where $\bar{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 galaxy-galaxy 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\kappa\kappa}(\ell)$, 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 un lensed galaxy luminosity function. This assumes that the magnification $\mu > 1$.

Assuming that, on our scales of interest, fluctuations in the intrinsic galaxy overdensity are small, and taking into account that, for weak lensing, $|k| \ll 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):

$$\delta_{\text{obs};i}^g(\theta) = \delta_i^g(\theta) + (5s_i - 2)\kappa_i(\theta), \qquad (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_{\text{lim}}$, such that

$$s_i = \left. \frac{\partial \log_{10} \eta(\bar{z}_i, m)}{\partial m} \right|_{m_{\text{lim}}}, \qquad (30)$$

where $\eta(\bar{z}_i, m)$ is the true, underlying distribution of galaxies, evaluated at the tomographic bin's central redshift, $\bar{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:

$$\gamma_{\text{obs};i}^{\alpha} = \gamma_{i}^{\alpha}(\theta) + \gamma_{i}^{\alpha}(\theta) \delta_{i}^{g}(\theta) + (5s_{i} - 2) \gamma_{i}^{\alpha}(\theta) \kappa_{i}(\theta). \qquad (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}^{\text{MB}}$, takes a similar form to the reduced shear correction of Eq. (22),

$$\delta C_{\ell;ij}^{\text{MB}} = \int_{0}^{\infty} \frac{d^2 \ell'}{(2\pi)^2} \cos(2\phi_{\ell'}) [(5s_i - 2)B_{ijj}^{\text{KKK}}(\ell, \ell', -\ell - \ell') \\ + (5s_j - 2)B_{ijj}^{\text{KKK}}(\ell, \ell', -\ell - \ell')] \qquad (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, $\Delta n^{\text{SO}}(z, \theta)$, can be modelled by:

$$ \frac{\Delta n^{\text{SO}}(z, \theta)}{n(z, \theta)} = -\pi \left[ (2\vartheta)^2 n_{\text{tot}} + \frac{A}{2-\zeta} (2\vartheta)^{2-\zeta} \right], \quad (33) $$

where $n_{\text{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, $\vartheta$, for all galaxies as in Hartlap et al. (2011), and $A$ and $\zeta$ 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\vartheta$ of a given one by integrating over the probability that another source centroid lies in an annulus of $d\vartheta$ 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 $\vartheta$ rather than $2\vartheta$). 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:

$$ dp(z, z', \theta) = \pi \vartheta^2 n(z', \theta) dz. \quad (34) $$

Here, it is also assumed that the expected number of galaxies overlapping with a given galaxy is $\ll 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{aligned} \Delta n^{\text{SO}}(z, \theta) &= -\pi \vartheta^2 n(z, \theta) \int_0^\infty dz' n(z', \theta) \\ &= -\pi \vartheta^2 n(z, \theta) [1 + \delta^g(\theta)] \int_0^\infty dz' n(z') \\ &= -\pi \vartheta^2 n(z, \theta) [1 + \delta^g(\theta)] n_{\text{tot}}. \end{aligned} \quad (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{aligned} \frac{\Delta n_i^{\text{SO}}(z, \theta)}{n_i(z, \theta)} &= -\pi \vartheta^2 \sum_{q=1}^{i} [1 + \delta_q^g(\theta)] n_{\text{tot};q} \\ &= -\pi \vartheta^2 \sum_{q=1}^{i} [1 + b_q \delta^g(\theta)] n_{\text{tot};q} \\ &= -\pi \vartheta^2 n_{\text{cumul.};i} - \pi \vartheta^2 \delta^g f_{\text{SO};i}, \end{aligned} \quad (36) $$

where the $n_{\text{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_{\text{cumul.};i}$ is the cumulative total surface density of galaxies for all redshift bins up to and including bin $i$, and

$$ f_{\text{SO};i} = \sum_{\rho=1}^{i} b_{\rho} n_{\text{tot};\rho}. \quad (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

$$ n_i^{\text{obs}}(\theta, \chi) = n_i(\chi) [1 + (5s_i - 2) \kappa_i(\theta) + \delta_i^g(\theta)] \times [1 - \pi \vartheta^2 n_{\text{cumul.};i} - \pi \vartheta^2 \delta^g(\theta) f_{\text{SO};i}]. \quad (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{aligned} \delta C_{\ell; i j}^{\text{SO}} = & (\pi^2 \vartheta^4 n_{\text{cumul.};i} n_{\text{cumul.};j} - \pi \vartheta^2 n_{\text{cumul.};i} \\ & - \pi \vartheta^2 n_{\text{cumul.};j}) C_{\ell; i j}^{\gamma \gamma} \end{aligned} \quad (39) $$

$$ \begin{aligned} \delta C_{\ell; i j}^{\text{SO-SLC}} = & \int_0^\infty \frac{d^2 \ell'}{(2\pi)^2} \cos(2\phi_{\ell'}) \\ & \times \left[ (\pi^2 \vartheta^4 n_{\text{cumul.};i} n_{\text{cumul.};j} - \pi \vartheta^2 n_{\text{cumul.};i} b_i \right. \\ & \left. - \pi \vartheta^2 n_{\text{cumul.};j} b_i - \pi \vartheta^2 (1 - n_{\text{cumul.};j}) \right. \\ & \left. \times f_{\text{SO};i}) B_{i j}^{\kappa \delta \gamma \kappa}(\ell_1, \ell_2, \ell_3) + i \leftrightarrow j \right], \end{aligned} \quad (40) $$

$$ \begin{aligned} \delta C_{\ell; i j}^{\text{SO-MB}} = & (\pi^2 \vartheta^4 n_{\text{cumul.};i} n_{\text{cumul.};j} - \pi \vartheta^2 n_{\text{cumul.};i} \\ & - \pi \vartheta^2 n_{\text{cumul.};j}) \delta C_{\ell; i j}^{\text{MB}}, \end{aligned} \quad (41) $$

where $i \leftrightarrow 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, $\delta_0$, can be defined as the matter density
contrast smoothed by a top-hat kernel of comoving radius R according to

$$ \delta_0(R, \chi) = \frac{3}{4\pi R^3} \int d^3\mathbf{r} \Theta(R - |\mathbf{r}|) \delta(\mathbf{r}, \chi), \quad (42) $$

where the matter density contrast is now expressed in terms of spatial distance, $\mathbf{r}$, rather than angle on the sky, and $\Theta$ 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 $\ell \gg 1$ (and note that cosmic shear is only defined for $\ell \ge 2$), this expression consists of two terms, the standard power spectra of Eq. (12) and a correction term, $\delta C_{\ell;ij}^{\text{LU}}$, which is defined as

$$ \delta C_{\ell;ij}^{\text{LU}} = 2 \frac{\delta_0(R, \chi)}{\sigma^2(R, \chi)} \int_0^{\chi_h} d\chi \frac{W_i(\chi)W_j(\chi)}{\chi^2} \times \left\{ \left[ \frac{34}{21}\xi_R(\chi) - \frac{4}{21}\psi_R(\chi) \right] P_{\delta\delta}(k, \chi) + \left[ \frac{\chi}{\ell}\xi'_R(\chi) - \frac{\ell}{\chi}\Omega_R(\chi) \right] \frac{\partial P_{\delta\delta}(k, \chi)}{\partial k} \frac{1}{\chi} - \frac{4}{7}\psi_R(\chi) \frac{\partial^2 P_{\delta\delta}(k, \chi)}{\partial k^2} \frac{1}{\chi^2} \right\}, \quad (43) $$

where $\sigma^2$ 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$, and $\Omega_R$ are correlation functions defined in Hall (2020). Also, the expression assumes a flat-geometry, which is valid under the current constraints on $\Omega_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, $\phi$, and the matter density contrast,

$$ \nabla_{\chi}^{2} \phi(\mathbf{r}, \chi) = \frac{3\Omega_{\text{m}}H_{0}^{2}}{2c^{2}a(t)} \delta(\mathbf{r}, \chi), \quad (44) $$

where $\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, \theta, \varphi)$, from the observer's frame of reference. Relating the two as above requires expressing the potential in spherical Bessel space as:

$$ \phi_{\ell m}(k) = \sqrt{\frac{2}{\pi}} \int d^3\mathbf{r} \phi(\mathbf{r}) j_{\ell}(kr) Y_{\ell m}(\theta, \varphi), \quad (45) $$

where $j_\ell$ refers to the spherical Bessel functions and $Y_{\ell 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:

$$ (\nabla_r^2 + k^2) j_{\ell}(kr) Y_{\ell m}(\theta, \varphi) = 0. \quad (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, $\Phi_\ell^\beta$, so that

$$ (\nabla_{S_K}^2 + k^2) \Phi_\ell^\beta(r) Y_{\ell m}(\theta, \varphi) = 0, \quad (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;\text{NF}}$, under the Limber approximation, under the Limber approximation, is given by modifying Eq. (12) to be:

$$ C_{\ell;ij}^{\gamma\gamma} = \int_0^{\chi_h} d\chi W_{\ell}^{\text{NF}}(\chi; K) \frac{W_i(\chi)W_j(\chi)}{S_K^2(\chi)} P_{\delta\delta}(k, \chi), \quad (48) $$

where

$$ W_{\ell}^{\text{NF}}(\chi; K) = \left[ 1 - \operatorname{sgn}(K) \frac{\ell^2}{(\ell + 1/2)^2 / S_K^2(\chi) + K} \right]^{-1/2} . \quad (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}^{\text{NF}}$, to the spatially flat angular power spectra such that

$$ \delta C_{\ell;ij}^{\text{NF}} = C_{\ell;ij}^{\gamma\gamma;\text{NF}} - C_{\ell;ij}^{\gamma\gamma}, \quad (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<fnref n="6" />.

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; Taylor et al. 2019; Upham et al. 2021;

<footnote n="6">⁶ Code for this paper is available at: https://github.com/desh1701/Euclid-Estimators-WL-Approximations-KP.</footnote>
Hall & Taylor 2022). Accordingly, the Fisher matrix for cosmic shear is defined as:

$$F_{\mu\nu} = \sum_{\ell'=\ell_{\min}}^{\ell_{\max}} \sum_{\ell=\ell_{\min}}^{\ell_{\max}} \sum_{ij,mn} \frac{\partial C_{\ell;ij}^{\epsilon\epsilon}}{\partial \theta_{\mu}} \mathrm{Cov}^{-1} [C_{\ell;ij}^{\epsilon\epsilon}, C_{\ell';mn}^{\epsilon\epsilon}] \frac{\partial C_{\ell';mn}^{\epsilon\epsilon}}{\partial \theta_{\nu}}, \quad (51)$$

where the $\mu$ and $\nu$ indices denote element in the Fisher matrix associated with cosmological parameters, $\theta_{\mu}$ and $\theta_{\nu}$, respectively, $\ell_{\min}$ is the minimum angular wavenumber of the survey, $\ell_{\max}$ is the maximum angular wavenumber used, the sums are over the $\ell$-blocks of power spectrum bands, and $\mathrm{Cov}^{-1}[C_{\ell;ij}^{\epsilon\epsilon}, C_{\ell';mn}^{\epsilon\epsilon}]$ 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, $\mathrm{Cov}_G$ and SSC, and $\mathrm{Cov}_{SSC}$ terms, as follows:

$$\mathrm{Cov}[C_{\ell;ij}^{\epsilon\epsilon}, C_{\ell';mn}^{\epsilon\epsilon}] = \mathrm{Cov}_G [C_{\ell;ij}^{\epsilon\epsilon}, C_{\ell';mn}^{\epsilon\epsilon}] + \mathrm{Cov}_{SSC} [C_{\ell;ij}^{\epsilon\epsilon}, C_{\ell';mn}^{\epsilon\epsilon}], \quad (52)$$

where the Gaussian component is given by

$$\mathrm{Cov}_G [C_{\ell;ij}^{\epsilon\epsilon}, C_{\ell';mn}^{\epsilon\epsilon}] = \frac{C_{\ell;im}^{\epsilon\epsilon} C_{\ell';jn}^{\epsilon\epsilon} + C_{\ell;in}^{\epsilon\epsilon} C_{\ell';jm}^{\epsilon\epsilon}}{(2\ell + 1)f_{\mathrm{sky}}\Delta\ell} \delta_{\ell\ell'}, \quad (53)$$

where $f_{\mathrm{sky}}$ is the fraction of the sky observed by the galaxy survey, $\Delta\ell$ is the bandwidth of the $\ell$-modes sampled, and $\delta^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)

$$\mathrm{Cov}_{SSC} [C_{\ell;ij}^{\epsilon\epsilon}, C_{\ell';mn}^{\epsilon\epsilon}] \approx R_{\ell} C_{\ell;ij}^{\epsilon\epsilon} R_{\ell'} C_{\ell';mn}^{\epsilon\epsilon} S_{ijmn}, \quad (54)$$

where $S_{ijmn}$ is the dimensionless volume-averaged covariance of the background matter density contrast, and $R_{\ell}$ 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\sigma$ uncertainties on each of the parameters. Explicitly, the uncertainty, $\sigma_{\mu}$, on parameter $\theta_{\mu}$ is given by:

$$\sigma_{\mu} = \sqrt{F_{\mu\mu}^{-1}}. \quad (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, $\delta C_{\ell;ij}$, the bias, $b_{\mu}$, on parameter $\theta_{\mu}$, is given by:

$$b_{\mu} = \sum_{\nu} (F^{-1})_{\mu\nu} \mathcal{B}_{\nu}, \quad (56)$$

where

$$\mathcal{B}_{\nu} = \sum_{\ell'=\ell_{\min}}^{\ell_{\max}} \sum_{\ell=\ell_{\min}}^{\ell_{\max}} \sum_{ij,mn} \delta C_{\ell;ij} \mathrm{Cov}^{-1} [C_{\ell;ij}^{\epsilon\epsilon}, C_{\ell';mn}^{\epsilon\epsilon}] \frac{\partial C_{\ell';mn}^{\epsilon\epsilon}}{\partial \theta_{\nu}}. \quad (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 $\Omega_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 $\ell$-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_{\epsilon} = \sqrt{2} \times 0.21 \approx 0.3^7$. Euclid is also expected to have a survey area such that $f_{\mathrm{sky}} = 0.36$. The survey’s galaxy surface density is anticipated to be $\bar{n}_g = 30 \ \mathrm{arcmin}^{-2}$.

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:

$$n_i(z) = \frac{\int_{z_i^-}^{z_i^+} dz_p n(z) p_{\mathrm{ph}}(z_p|z)}{\int_{z_{\min}}^{z_{\max}} dz \int_{z_i^-}^{z_i^+} dz_p n(z) p_{\mathrm{ph}}(z_p|z)}, \quad (58)$$

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

$$n(z) \propto \left(\frac{z}{z_0}\right)^2 \exp\left[-\left(\frac{z}{z_0}\right)^{3/2}\right], \quad (59)$$

<footnote n="7">7 We use the specification in EC20, but note that Euclid Collaboration (2019) uses a value of 0.26 per component. Since we are looking at biases caused by the differences in the signal the shot noise component does not affect the reported biases, however the relative significance (bias divided by error) will be lower for a larger shot noise term.</footnote>
<table><caption>**Table 2.** Values of model parameters used in defining the uncertainty of photometric redshift estimates through Eq. (60).</caption><thead><tr><td>Parameter</td><td>Fiducial Value</td></tr></thead><tbody><tr><td>c<sub>b</sub></td><td>1.0</td></tr><tr><td>z<sub>b</sub></td><td>0.0</td></tr><tr><td>σ<sub>b</sub></td><td>0.05</td></tr><tr><td>c<sub>o</sub></td><td>1.0</td></tr><tr><td>z<sub>o</sub></td><td>0.1</td></tr><tr><td>σ<sub>o</sub></td><td>0.05</td></tr><tr><td>f<sub>out</sub></td><td>0.1</td></tr></tbody></table>

Notes. Chosen according to EC20.

where $z_0 = z_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):

$$p_{ph}(z_p|z) = \frac{1 - f_{out}}{\sqrt{2\pi}\sigma_b(1+z)} \exp\left\{-\frac{1}{2}\left[\frac{z - c_b z_p - z_b}{\sigma_b(1+z)}\right]^2\right\} + \frac{f_{out}}{\sqrt{2\pi}\sigma_o(1+z)} \exp\left\{-\frac{1}{2}\left[\frac{z - c_o z_p - z_o}{\sigma_o(1+z)}\right]^2\right\}. \quad (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_0w_a$CDM 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 $\Omega_m$, the present-day baryonic matter density parameter $\Omega_b$, the dimensionless curvature parameter $\Omega_K = -K(c/H_0)^2$, the Hubble parameter $h = H_0/100 \text{ km s}^{-1} \text{ Mpc}^{-1}$, the spectral index $n_s$, the RMS value of density fluctuations on $8h^{-1} \text{ Mpc}$ scales $\sigma_8$, and massive neutrinos with a sum of masses $\sum m_v \neq 0$. The $w_0w_a$CDM model additionally adds in the present-day value of the dark energy equation of state $w_0$, and the high-redshift value of the dark energy equation of state, $w_a$. Typically, the present-day densities $\Omega_i$, $i \in \{m, b, K\}$, are denoted with an additional subscript 0; we omit this here for brevity. Primarily, we are interested in the $w_0w_a$CDM 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 $\sum m_v \neq 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 $\Omega_K$ to 0. Only when testing the significance of the additional

<table><caption>**Table 3.** ΛCDM and $w_0w_a$CDM cosmological parameter fiducial values used in this investigation.</caption><thead><tr><td>Cosmological parameter</td><td>Fiducial value</td></tr></thead><tbody><tr><td>Ω<sub>m</sub></td><td>0.32</td></tr><tr><td>Ω<sub>b</sub></td><td>0.05</td></tr><tr><td>h</td><td>0.67</td></tr><tr><td>n<sub>s</sub></td><td>0.96</td></tr><tr><td>σ<sub>8</sub></td><td>0.816</td></tr><tr><td>w<sub>0</sub></td><td>-1</td></tr><tr><td>w<sub>a</sub></td><td>0</td></tr><tr><td>Σ m<sub>v</sub> (eV)</td><td>0.06</td></tr><tr><td>Ω<sub>K</sub></td><td>{0, 0.05}</td></tr></tbody></table>

Notes. These values correspond to EC20. It should be noted that $\sum m_v \neq 0$ is assumed to be fixed, and uncertainties and biases are not calculated for it. Additionally, two possible values are provided for $\Omega_K$, because the non-zero value must be used when evaluating the non-flat Universe correction. This value is selected using the upper-bound of the Planck Collaboration VI (2020) $1\sigma$ uncertainty.

non-flat Universe correction term was set to 0.05. In this case, the selected value serves as the upper-limit of $1\sigma$ 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: $\Omega_m, \Omega_b, h, n_s, \sigma_8$ and $\mathcal{A}_{IA}$ for the ΛCDM case and, additionally, $w_0$ and $w_a$, for the $w_0w_a$CDM case. When the correction for spatially curvature needed to be tested, the Fisher matrices also include $\Omega_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<fnref n="8" /> code (Lacasa & Grain 2019), with an $R_\ell$ of 3.

To calculate the cosmological background quantities required for the investigation, including the matter power spectrum and growth factor, we used the CAMB<fnref n="9" /> 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<fnref n="10" /> (Astropy Collaboration 2013, 2018) to compute the cosmological distances. The NLA model IA parameters were set to $\mathcal{A}_{IA} = 1.72$ and $C_{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 $\ell$-bands. The limits for these were logarithmically spaced, with an $\ell_{min}$ of 10 and an $\ell_{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 obscu-

<footnote n="8">⁸ https://github.com/fabienlacasa/PySSC</footnote>

<footnote n="9">⁹ https://camb.info/</footnote>

<footnote n="10">¹⁰ http://www.astropy.org</footnote>
<table><caption>**Table 4.** Values of the slope of the luminosity function used in computing the magnification bias correction.</caption><thead><tr><th>Bin i</th><th>Central redshift</th><th>Slope s<sub>i</sub></th></tr></thead><tbody><tr><td>1</td><td>0.2095</td><td>0.108</td></tr><tr><td>2</td><td>0.489</td><td>0.180</td></tr><tr><td>3</td><td>0.619</td><td>0.229</td></tr><tr><td>4</td><td>0.7335</td><td>0.279</td></tr><tr><td>5</td><td>0.8445</td><td>0.335</td></tr><tr><td>6</td><td>0.9595</td><td>0.400</td></tr><tr><td>7</td><td>1.087</td><td>0.480</td></tr><tr><td>8</td><td>1.2395</td><td>0.586</td></tr><tr><td>9</td><td>1.45</td><td>0.753</td></tr><tr><td>10</td><td>2.038</td><td>1.335</td></tr></tbody></table>

Notes. These are calculated at the central redshift of each tomographic bin. The limiting magnitude is taken to be 24.5, and the slopes are calculated with a fitting function (Euclid Collaboration 2022b) determined from the Euclid Flagship simulation.

ration corrections, we used the BiHalofit model and code<fnref n="11" /> (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 \times 10^9$, so that the total number of galaxies per redshift bin was $2 \times 10^8$. We also took the mean galaxy radius to be $\theta = 0.32'' (1.55 \times 10^{-6}$ 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^{-1}$ 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, $\delta_0(R = 120 h^{-1}$ Mpc) = 0.045 and $\delta_0(R = 120 h^{-1}$ Mpc)/$\sigma(R = 120 h^{-1}$ Mpc) = 0.85$. Different choices of smoothing scale result in different values of the local overdensity, $\delta_0$, 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, $\Delta C_\ell/C_\ell$. This was calculated according to Kaiser (1992), and took the form

$$ \Delta C_{\ell}/C_{\ell} = \sqrt{2} [f_{\text{sky}}(2\ell + 1)]^{-1/2}. \quad (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 spectra) 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.

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 $\ell$-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\sigma$ uncertainty on that parameter (Taylor et al. 2007). This is because at that point, the biased and unbiased $1\sigma$ 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 $\Lambda$CDM and $w_0w_a$ CDM cases, respectively. These are also represented visually in Fig. 2. The biases from the non-flat Universe correction (which requires

<footnote n="11">¹¹ http://cosmo.phys.hirosaki-u.ac.jp/takahasi/codes_e.htm</footnote>
<figure><img src="image_1.png" /><figcaption>**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 $w_0w_a$CDM 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 $\Omega_K = 0.05$, whilst other corrections are when $\Omega_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 $\ell$-modes where a computation was made; all quantities were evaluated for 200 $\ell$-bands, logarithmically spaced, with an $\ell_{min}$ of 10, and an $\ell_{max}$ of 5000.</figcaption></figure>

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_0w_a$CDM case. This difference is likely due to the presence of the variable dark energy parameters in the $w_0w_a$CDM 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.

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_0w_a$CDM 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
<table><caption>**Table 5.** Uncertainties on, and biases induced from neglecting the various corrections, in the ΛCDM parameters of Table 3 for *Euclid*.</caption><thead><tr><th>Cosm. param.</th><th>Uncertainty (1σ)</th><th>RS Bias/σ</th><th>SLC Bias/σ</th><th>MB Bias/σ</th><th>SO Bias/σ</th><th>SO-SLC Bias/σ</th><th>SO-MB Bias/σ</th><th>LU Bias/σ</th><th>Total All</th><th>Total Sig.</th></tr></thead><tbody><tr><td>Ω<sub>m</sub></td><td>0.0051</td><td>-0.032</td><td>1.14</td><td>-0.75</td><td>0.36</td><td>0.0032</td><td>0.0047</td><td>0.26</td><td>0.76</td><td>0.73</td></tr><tr><td>Ω<sub>b</sub></td><td>0.021</td><td>0.0035</td><td>0.28</td><td>0.080</td><td>-0.044</td><td>6.5 × 10<sup>-4</sup></td><td>-6.6 × 10<sup>-4</sup></td><td>-0.016</td><td>0.26</td><td>0.28</td></tr><tr><td>h</td><td>0.13</td><td>0.0084</td><td>0.15</td><td>0.19</td><td>-0.11</td><td>2.8 × 10<sup>-4</sup></td><td>-0.0013</td><td>-0.042</td><td>0.23</td><td>0.25</td></tr><tr><td>n<sub>s</sub></td><td>0.029</td><td>0.017</td><td>0.0065</td><td>-0.23</td><td>0.18</td><td>2.1 × 10<sup>-4</sup></td><td>0.0016</td><td>0.038</td><td>2.6 × 10<sup>-4</sup></td><td>-0.07</td></tr><tr><td>σ<sub>8</sub></td><td>0.0072</td><td>0.076</td><td>-1.10</td><td>0.78</td><td>-0.54</td><td>-0.0031</td><td>-0.0050</td><td>-0.22</td><td>-0.78</td><td>-0.79</td></tr></tbody></table>

**Notes.** This table considers the case when $\Omega_K = 0$, and accordingly lists the biases resulting from all corrections except the non-flat Universe correction. That can be found in Table 7. All biases are given as a fraction of the $1\sigma$ uncertainty on each parameter. A bias is considered significant if it reaches or exceeds $0.25\sigma$, as at this point its uncorrected and corrected confidence contours overlap by less than 90%. '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 total biases from the sum of all corrections, as well as the total biases from only the individually significant corrections are also given.

<table><caption>**Table 6.** Uncertainties (along with the biases induced from neglecting the various corrections) on the $w_0w_a$ CDM parameters of Table 3 for *Euclid*.</caption><thead><tr><th>Cosm. Param.</th><th>Uncertainty (1σ)</th><th>RS Bias/σ</th><th>SLC Bias/σ</th><th>MB Bias/σ</th><th>SO Bias/σ</th><th>SO-SLC Bias/σ</th><th>SO-MB Bias/σ</th><th>LU Bias/σ</th><th>Total All</th><th>Total Sig.</th></tr></thead><tbody><tr><td>Ω<sub>m</sub></td><td>0.010</td><td>-0.079</td><td>1.25</td><td>0.12</td><td>0.12</td><td>0.0025</td><td>-0.0022</td><td>0.14</td><td>1.55</td><td>1.49</td></tr><tr><td>Ω<sub>b</sub></td><td>0.021</td><td>-0.0013</td><td>0.30</td><td>0.069</td><td>-0.021</td><td>7.2 × 10<sup>-4</sup></td><td>-5.8 × 10<sup>-4</sup></td><td>-0.0032</td><td>0.35</td><td>0.35</td></tr><tr><td>h</td><td>0.13</td><td>0.0057</td><td>0.11</td><td>0.11</td><td>-0.053</td><td>3.2 × 10<sup>-4</sup></td><td>-6.3 × 10<sup>-4</sup></td><td>-0.015</td><td>0.16</td><td>0.17</td></tr><tr><td>n<sub>s</sub></td><td>0.031</td><td>0.051</td><td>-0.19</td><td>-0.19</td><td>0.026</td><td>-3.0 × 10<sup>-4</sup></td><td>0.0015</td><td>-0.045</td><td>-0.35</td><td>-0.35</td></tr><tr><td>σ<sub>8</sub></td><td>0.012</td><td>0.081</td><td>-1.19</td><td>0.017</td><td>-0.19</td><td>-0.0024</td><td>0.0011</td><td>-0.093</td><td>-1.38</td><td>-1.36</td></tr><tr><td>w<sub>0</sub></td><td>0.13</td><td>-0.076</td><td>0.81</td><td>0.55</td><td>-0.039</td><td>0.0011</td><td>-0.0051</td><td>0.032</td><td>1.28</td><td>1.31</td></tr><tr><td>w<sub>a</sub></td><td>0.36</td><td>0.022</td><td>-0.51</td><td>-0.60</td><td>0.27</td><td>-3.3 × 10<sup>-4</sup></td><td>0.0053</td><td>0.097</td><td>-0.72</td><td>-0.84</td></tr></tbody></table>

**Notes.** This table considers the case when $\Omega_K = 0$, and accordingly lists the biases resulting from all corrections except the non-flat Universe correction. That can be found in Table 7. All biases are given as a fraction of the $1\sigma$ uncertainty on each parameter. A bias is considered significant if it reaches or exceeds $0.25\sigma$, as at this point its uncorrected and corrected confidence contours overlap by less than 90%. '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 total biases from the sum of all corrections, as well as the total biases from only the individually significant corrections are also given.

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_0w_a$ CDM 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.

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 cor-

rections 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,
<figure><img src="image_2.png" /><figcaption>**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 $w_0w_a$ CDM 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-negligible 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.</figcaption></figure>

<table><caption>**Table 7.** Uncertainties on the fiducial ΛCDM and $w_0w_a$ CDM cosmologies of Table 3 when $\Omega_K = 0.05$, and biases induced from neglecting the non-flat Universe correction for the *Euclid* cosmic shear probe.</caption><thead><tr><th>Cosmo. Param.</th><th>ΛCDM Uncertainty (1σ)</th><th>ΛCDM Non-flat Universe Bias/σ</th><th>w<sub>0</sub>w<sub>a</sub>CDM Uncertainty (1σ)</th><th>w<sub>0</sub>w<sub>a</sub>CDM Non-flat Universe Bias/σ</th></tr></thead><tbody><tr><td>Ω<sub>m</sub></td><td>0.0050</td><td>0.043</td><td>0.014</td><td>-0.062</td></tr><tr><td>Ω<sub>b</sub></td><td>0.021</td><td>-0.0070</td><td>0.021</td><td>8.1 × 10<sup>-4</sup></td></tr><tr><td>Ω<sub>K</sub></td><td>0.034</td><td>0.043</td><td>0.062</td><td>-0.054</td></tr><tr><td>h</td><td>0.12</td><td>-0.012</td><td>0.12</td><td>0.012</td></tr><tr><td>n<sub>s</sub></td><td>0.029</td><td>-0.010</td><td>0.030</td><td>0.012</td></tr><tr><td>σ<sub>8</sub></td><td>0.028</td><td>0.026</td><td>0.042</td><td>-0.044</td></tr><tr><td>w<sub>0</sub></td><td>N/A</td><td>N/A</td><td>0.14</td><td>-0.090</td></tr><tr><td>w<sub>a</sub></td><td>N/A</td><td>N/A</td><td>0.56</td><td>0.10</td></tr></tbody></table>

**Notes.** The biases are given as a fraction of the 1σ uncertainty on each parameter. A bias is considered significant if it reaches or exceeds 0.25σ, as at this point its uncorrected and corrected confidence contours overlap by less than 90%. The predicted biases from this correction are well below significance for all parameters.

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_0w_a$ CDM 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
<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Projected $1\sigma$ and $2\sigma$ 2-parameter uncertainty contours for Euclid under a $\Lambda$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 $\Omega_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 $\Omega_m$, $\Omega_b$, $h$, and $\sigma_8$, and their values can be found in Table 5</figcaption></figure>

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\sigma$ 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
<figure><img src="image_4.png" /><figcaption>**Fig. 4.** Projected $1\sigma$ and $2\sigma$ 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 $w_0w_a$ CDM cosmology specified in Table 3, in the case when $\Omega_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 $\Omega_m, \Omega_b, n_s, \sigma_8, w_0$, and $w_a$, and their values can be found in Table 6.</figcaption></figure>

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_0w_a$ CDM 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\sigma$. However, in the ΛCDM case we found that when the biases are combined, they frequently
<table><caption>**Table 8.** Increasing and decreasing multiplicative changes required in the matter bispectrum, across all configurations and all ℓ-modes, to significantly affect the biases from each relevant correction term.</caption><thead><tr><th>Correction</th><th>Min. change param</th><th>Bispec. mult. increase</th><th>Bispec. mult. decrease</th></tr></thead><tbody><tr><td>RS</td><td>Ω<sub>m</sub></td><td>-2.16</td><td>-4.16</td></tr><tr><td>SLC</td><td>Ω<sub>m</sub></td><td>1.20</td><td>0.80</td></tr><tr><td>MB</td><td>w<sub>a</sub></td><td>0.58</td><td>1.41</td></tr><tr><td>SO-SLC</td><td>Ω<sub>m</sub></td><td>100.63</td><td>-99.67</td></tr><tr><td>SO-MB</td><td>w<sub>a</sub></td><td>47.00</td><td>-47.00</td></tr></tbody></table>

**Notes.** A significant change is when any one of the biases on the cosmological parameters changes by ±0.25σ. These values encapsulate the fractional change for a fixed cosmology, in this case the $w_0w_a$CDM cosmology of Table 3, with $\Omega_K = 0$. It should be noted that these values are not exhaustive thresholds, as sufficiently significant changes to the model at a particular sub-range of scales may also be problematic.

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_0w_a$CDM 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_0w_a$CDM 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.

**Acknowledgements.** ACD and AH are supported by the Royal Society. TDK acknowledges funding from the EU's Horizon 2020 programme, grant agreement No 776247. 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 Academy of Finland, the Agenzia Spaziale Italiana, the Belgian Science Policy, the Canadian Euclid Consortium, the French Centre National d'Etudes Spatiales, the Deutsches Zentrum für Luft- und Raumfahrt, the Danish Space Research Institute, the Fundação para a Ciência e a Tecnologia, the Ministerio de Ciencia e Innovación, the National Aeronautics and Space Administration, the National Astronomical Observatory of Japan, the Netherlandse Onderzoeksschool Voor Astronomie, the Norwegian Space Agency, 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. CJM acknowledges FCT and POCH/FSE (EC) support through Investigador FCT Contract 2021.01214.CEECIND/CP1658/CT0001. A complete and detailed list is available on the *Euclid* web site (http://www.euclid-ec.org).

## References

* Abbott, T. M. C., Aguena, M., Alarcon, A., et al. 2022, Phys. Rev. D, 105, 023520

* Akeson, R., Armus, L., Bachelet, E., et al. 2019, arXiv e-prints [arXiv:1902.05569]

* Albrecht, A., Bernstein, G., Cahn, R., et al. 2006, arXiv e-prints [arXiv:astro-ph/0609591]

* Asgari, M., Lin, C.-A., Joachimi, B., et al. 2021, A&A, 645, A104

* Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33

* Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123

* Barreira, A., Krause, E., & Schmidt, F. 2018a, JCAP, 2018, 053

* Barreira, A., Krause, E., & Schmidt, F. 2018b, JCAP, 2018, 015

* Barreira, A., Nelson, D., Pillepich, A., et al. 2019, MNRAS, 488, 2079

* Bartelmann, M., & Schneider, P. 2001, Phys. Rep., 340, 291
Bernardeau, F. 1998, A&A, 338, 375
Bernardeau, F., Bonvin, C., & Vernizzi, F. 2010, Phys. Rev. D, 81, 083002
Bird, S., Viel, M., & Haehnelt, M. G. 2012, MNRAS, 420, 2551
Breton, M.-A., & Fleury, P. 2021, A&A, 655, A54
Bridle, S., & King, L. 2007, New J. Phys., 9, 444
Cooray, A., & Hu, W. 2002, ApJ, 574, 19
Croppier, M., Cole, R., James, A., et al. 2012, in Proc. SPIE, 8442, 84420V
Cuesta-Lazaro, C., Quera-Bofarull, A., Reischke, R., & Schäfer, B. M. 2018, MNRAS, 477, 741
Deshpande, A. C., & Kitching, T. D. 2020, Phys. Rev. D, 101, 103531
Deshpande, A. C., & Kitching, T. D. 2021, Phys. Rev. D, 103, 123510
Deshpande, A. C., Kitching, T. D., Cardone, V. F., et al. 2020a, A&A, 636, A95
Deshpande, A. C., Taylor, P. L., & Kitching, T. D. 2020b, Phys. Rev. D, 102, 083535
Desjacques, V., Jeong, D., & Schmidt, F. 2018, Phys. Rep., 733, 1
Dodelson, S., Shapiro, C., & White, M. 2006, Phys. Rev. D, 73, 023009
Duncan, C. A. J., Harnois-Déraps, J., Miller, L., & Langedijk, A. 2022, MNRAS, 515, 1130
Euclid Collaboration (Martinet, N., et al.) 2019, A&A, 627, A59
Euclid Collaboration (Blanchard, A., et al.) 2020, A&A, 642, A191
Euclid Collaboration (Pocino, A., et al.) 2021, A&A, 655, A44
Euclid Collaboration (Bretonnière, H., et al.) 2022a, A&A, 657, A90
Euclid Collaboration (Lepori, F., et al.) 2022b, A&A, 662, A93
Fleury, P., Larena, J., & Uzan, J.-P. 2017, Phys. Rev. Lett., 119, 191101
Fleury, P., Larena, J., & Uzan, J.-P. 2019, Phys. Rev. D, 99, 023525
Fortuna, M. C., Hoekstra, H., Johnston, H., et al. 2021, A&A, 654, A76
Gil-Marín, H., Wagner, C., Fragkoudi, F., Jimenez, R., & Verde, L. 2012, JCAP, 2012, 047
Hall, A. 2020, Phys. Rev. D, 101, 043519
Hall, A., & Taylor, A. 2022, Phys. Rev. D, 105, 123527
Hamana, T., Colombi, S. T., Thion, A., et al. 2002, MNRAS, 330, 365
Hartlap, J., Hilbert, S., Schneider, P., & Hildebrandt, H. 2011, A&A, 528, A51
Heydenreich, S., Schneider, P., Hildebrandt, H., et al. 2020, A&A, 634, A104
Hikage, C., Oguri, M., Hamana, T., et al. 2019, PASJ, 71, 43
Hu, W. 1999, APJ, 522, L21
Hu, W., & Kravtsov, A. V. 2003, ApJ, 584, 702
Hui, L., Gaztañaga, E., & Loverde, M. 2007, Phys. Rev. D, 76, 103502
Joachimi, B., Cacciato, M., Kitching, T., et al. 2015, Space Sci. Rev., 193, 1
Kaiser, N. 1992, ApJ, 388, 272
Kaiser, N. 1998, ApJ, 498, 26
Kiessling, A., Cacciato, M., Joachimi, B., et al. 2015, Space Sci. Rev., 193, 67
Kilbinger, M. 2015, Rep. Prog. Phys., 78, 086901
Kirk, D., Brown, M., Hoekstra, H., et al. 2015, Space Sci. Rev., 193, 139
Kitching, T. D., & Deshpande, A. C. 2022, Open J. Astrophys., 5, 6
Kitching, T. D., & Heavens, A. F. 2017, Phys. Rev. D, 95, 063522
Kitching, T. D., Taylor, A. N., & Heavens, A. F. 2008, MNRAS, 389, 173
Kitching, T. D., Alsing, J., Heavens, A. F., et al. 2017, MNRAS, 469, 2737
Kitching, T. D., Paykari, P., Hoekstra, H., & Cropper, M. 2019, Open J. Astrophys., 2, 5
Kitching, T. D., Deshpande, A. C., & Taylor, P. L. 2020, Open J. Astrophys., 3, 14
Kitching, T., Deshpande, A., & Taylor, P. 2021, Open J. Astrophys., 4, 17
Krause, E., & Hirata, C. 2010, A&A, 523, A28
Krause, E., Fang, X., Pandey, S., et al. 2021, arXiv e-prints [arXiv:2105.13548]
Lacasa, F., & Grain, J. 2019, A&A, 624, A61
Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv e-prints [arXiv:1110.3193]
Lavaux, G., & Hudson, M. J. 2011, MNRAS, 416, 2840
Lesgourgues, J., & Tram, T. 2014, JCAP, 2014, 032
Lewis, A., Challinor, A., & Lasenby, A. 2000, ApJ, 538, 473
Limber, D. N. 1953, ApJ, 117, 134
Lin, C. H., Harnois-Déraps, J., Eifler, T., et al. 2020, MNRAS, 499, 2977
Loureiro, A., Whittaker, L., Spurio Mancini, A., et al. 2022, A&A, 665, A56
LoVerde, M., & Afshordi, N. 2008, Phys. Rev. D, 78, 123506
LSST Science Collaboration (Abell, P. A., et al.) 2009, arXiv e-prints [arXiv:0912.0201]
Martinelli, M., Tutusaus, I., Archidiacono, M., et al. 2021, A&A, 649, A100
Munshi, D., Valagea, P., van Waerbeke, L., & Heavens, A. 2008, Phys. Rep., 462, 67
Planck Collaboration VI. 2020, A&A, 641, A6
Potter, D., Stadel, J., & Teyssier, R. 2017, Comput. Astrophys. Cosmol., 4, 2
Reischke, R., Schäfer, B. M., Bolejko, K., Lewis, G. F., & Lautsch, M. 2019, MNRAS, 486, 5061
Schmidt, F., Rozo, E., Dodelson, S., Hui, L., & Sheldon, E. 2009, ApJ, 702, 593
Schneider, P., & Er, X. 2008, A&A, 485, 363
Schneider, P., van Waerbeke, L., & Mellier, Y. 2002, A&A, 389, 729

Scoccimarro, R., & Couchman, H. 2001, MNRAS, 325, 1312
Seitz, C., & Schneider, P. 1997, A&A, 318, 687
Semboloni, E., Hoekstra, H., & Schaye, J. 2013, MNRAS, 434, 148
Shapiro, C. 2009, ApJ, 696, 775
Shapiro, C., & Cooray, A. 2006, JCAP, 2006, 007
Simpson, F., Peacock, J. A., & Heavens, A. F. 2010, MNRAS, 402, 2009
Spurio Mancini, A., Piras, D., Alsing, J., Joachimi, B., & Hobson, M. P. 2022, MNRAS, 511, 1771
Takada, M., & Hu, W. 2013, Phys. Rev. D, 87, 123504
Takahashi, R., Sato, M., Nishimichi, T., Taruya, A., & Oguri, M. 2012, ApJ, 761, 152
Takahashi, R., Nishimichi, T., Namikawa, T., et al. 2020, ApJ, 895, 113
Taylor, A., Kitching, T., Bacon, D., & Heavens, A. 2007, MNRAS, 374, 1377
Taylor, P. L., Bernardeau, F., & Kitching, T. D. 2018a, Phys. Rev. D, 98, 083514
Taylor, P. L., Kitching, T. D., McEwen, J. D., & Tram, T. 2018b, Phys. Rev. D, 98, 023522
Taylor, P. L., Kitching, T. D., Alsing, J., et al. 2019, Phys. Rev. D, 100, 023519
Tegmark, M., Taylor, A., & Heavens, A. 2015, ApJ, 78, 086901
Turner, E. L., Ostriker, J. P., & Gott, J. R., III 1984, ApJ, 284, 1
Upham, R. E., Brown, M. L., & Whittaker, L. 2021, MNRAS, 503, 1999
Upham, R. E., Brown, M. L., Whittaker, L., et al. 2022, A&A, 660, A114
Viola, M., Kitching, T. D., & Joachimi, B. 2014, MNRAS, 439, 1909
Yu, Y., Zhang, P., Lin, W., & Cui, W. 2015, ApJ, 803, 46

<footnote n="1">¹ Mullard Space Science Laboratory, University College London, Holmbury St Mary, Dorking, Surrey RH5 6NT, UK
e-mail: t.kitching@ucl.ac.uk</footnote>

<footnote n="2">² Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="3">³ Jodrell Bank Centre for Astrophysics, Department of Physics and Astronomy, University of Manchester, Oxford Road, Manchester M13 9PL, UK</footnote>

<footnote n="4">⁴ Université Paris-Saclay, CNRS, Institut d'astrophysique spatiale, 91405 Orsay, France</footnote>

<footnote n="5">⁵ Institut für Theoretische Physik, University of Heidelberg, Philosophenweg 16, 69120 Heidelberg, Germany</footnote>

<footnote n="6">⁶ INAF-Osservatorio Astronomico di Brera, Via Brera 28, 20122 Milano, Italy</footnote>

<footnote n="7">⁷ INAF-Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Via Piero Gobetti 93/3, 40129 Bologna, Italy</footnote>

<footnote n="8">⁸ Dipartimento di Fisica e Astronomia "Augusto Righi" - Alma Mater Studiorum Università di Bologna, via Piero Gobetti 93/2, 40129 Bologna, Italy</footnote>

<footnote n="9">⁹ INFN-Sezione di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy</footnote>

<footnote n="10">¹⁰ Max Planck Institute for Extraterrestrial Physics, Giessenbachstr. 1, 85748 Garching, Germany</footnote>

<footnote n="11">¹¹ Universitäts-Sternwarte München, Fakultät für Physik, Ludwig-Maximilians-Universität München, Scheinerstrasse 1, 81679 München, Germany</footnote>

<footnote n="12">¹² INAF-Osservatorio Astrofisico di Torino, Via Osservatorio 20, 10025 Pino Torinese (TO), Italy</footnote>

<footnote n="13">¹³ Dipartimento di Fisica, Università di Genova, Via Dodecaneso 33, 16146 Genova, Italy</footnote>

<footnote n="14">¹⁴ INFN-Sezione di Roma Tre, Via della Vasca Navale 84, 00146 Roma, Italy</footnote>

<footnote n="15">¹⁵ Department of Physics “E. Pancini”, University Federico II, Via Cinthia 6, 80126 Napoli, Italy</footnote>

<footnote n="16">¹⁶ Instituto de Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Rua das Estrelas, 4150-762 Porto, Portugal</footnote>

<footnote n="17">¹⁷ Dipartimento di Fisica, Università degli Studi di Torino, Via P. Giuria 1, 10125 Torino, Italy</footnote>

<footnote n="18">¹⁸ INFN-Sezione di Torino, Via P. Giuria 1, 10125 Torino, Italy</footnote>

<footnote n="19">¹⁹ INAF-IASF Milano, Via Alfonso Corti 12, 20133 Milano, Italy</footnote>

<footnote n="20">²⁰ INAF-Osservatorio Astronomico di Roma, Via Frascati 33, 00078 Monteporzio Catone, Italy</footnote>

<footnote n="21">²¹ INFN-Sezione di Roma, Piazzale Aldo Moro, 2 - c/o Dipartimento di Fisica, Edificio G. Marconi, 00185 Roma, Italy</footnote>

<footnote n="22">²² Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, 08193 Bellaterra (Barcelona), Spain</footnote>
23 Port d'Informació Científica, Campus UAB, C. Albareda s/n, 08193 Bellaterra (Barcelona), Spain

24 Institute for Theoretical Particle Physics and Cosmology (TTK), RWTH Aachen University, 52056 Aachen, Germany

25 Institut d'Estudis Espacials de Catalunya (IEEC), Carrer Gran Capità 2-4, 08034 Barcelona, Spain

26 Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans, s/n, 08193 Barcelona, Spain

27 INAF-Osservatorio Astronomico di Capodimonte, Via Moiariello 16, 80131 Napoli, Italy

28 INFN section of Naples, Via Cinthia 6, 80126 Napoli, Italy

29 Dipartimento di Fisica e Astronomia "Augusto Righi" - Alma Mater Studiorum Università di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy

30 Centre National d'Etudes Spatiales - Centre spatial de Toulouse, 18 avenue Edouard Belin, 31401 Toulouse Cedex 9, France

31 Institut national de physique nucléaire et de physique des particules, 3 rue Michel-Ange, 75794 Paris Cedex 16, France

32 ESAC/ESA, Camino Bajo del Castillo, s/n., Urb. Villafranca del Castillo, 28692 Villanueva de la Cañada, Madrid, Spain

33 European Space Agency/ESRIN, Largo Galileo Galilei 1, 00044 Frascati, Roma, Italy

34 Institute of Physics, Laboratory of Astrophysics, Ecole Polytechnique Fédérale de Lausanne (EPFL), Observatoire de Sauverny, 1290 Versoix, Switzerland

35 UCB Lyon 1, CNRS/IN2P3, IUF, IP2I Lyon, 4 rue Enrico Fermi, 69622 Villeurbanne, France

36 Departamento de Física, Faculdade de Ciências, Universidade de Lisboa, Edifício C8, Campo Grande, 1749-016 Lisboa, Portugal

37 Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências, Universidade de Lisboa, Campo Grande, 1749-016 Lisboa, Portugal

38 Department of Astronomy, University of Geneva, ch. d'Ecogia 16, 1290 Versoix, Switzerland

39 Department of Physics, Oxford University, Keble Road, Oxford OX1 3RH, UK

40 INAF-Istituto di Astrofisica e Planetologia Spaziali, via del Fosso del Cavaliere, 100, 00100 Roma, Italy

41 Université Paris-Saclay, Université Paris Cité, CEA, CNRS, Astrophysique, Instrumentation et Modélisation Paris-Saclay, 91191 Gif-sur-Yvette, France

42 University of Lyon, Univ Claude Bernard Lyon 1, CNRS/IN2P3, IP2I Lyon, UMR 5822, 69622 Villeurbanne, France

43 INAF-Osservatorio Astronomico di Trieste, Via G. B. Tiepolo 11, 34143 Trieste, Italy

44 INAF-Osservatorio Astronomico di Padova, Via dell'Osservatorio 5, 35122 Padova, Italy

45 Institute of Theoretical Astrophysics, University of Oslo, PO Box 1029 Blindern, 0315 Oslo, Norway

46 Leiden Observatory, Leiden University, Niels Bohrweg 2, 2333 CA Leiden, The Netherlands

47 Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA, 91109, USA

48 Technical University of Denmark, Elektrovej 327, 2800 Kgs. Lyn-gby, Denmark

49 Cosmic Dawn Center (DAWN), Denmark

50 Institut d'Astrophysique de Paris, 98bis Boulevard Arago, 75014, Paris, France

51 Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany

52 Department of Physics and Helsinki Institute of Physics, Gustaf Hällströmin katu 2, 00014 University of Helsinki, Finland

53 Aix-Marseille Université, CNRS/IN2P3, CPPM, Marseille, France

54 Université de Genève, Département de Physique Théorique and Centre for Astroparticle Physics, 24 quai Ernest-Ansermet, 1211 Genève 4, Switzerland

55 Department of Physics, PO Box 64, 00014 University of Helsinki, Finland

56 Helsinki Institute of Physics, Gustaf Hällströmin katu 2, University of Helsinki, Helsinki, Finland

57 NOVA optical infrared instrumentation group at ASTRON, Oude Hoogeveensedijk 4, 7991PD, Dwingeloo, The Netherlands

58 Universität Bonn, Argelander-Institut für Astronomie, Auf dem Hügel 71, 53121 Bonn, Germany

59 Aix-Marseille Université, CNRS, CNES, LAM, Marseille, France

60 Department of Physics, Institute for Computational Cosmology, Durham University, South Road, DH1 3LE, UK

61 Université Paris Cité, CNRS, Astroparticule et Cosmologie, 75013 Paris, France

62 Institut d'Astrophysique de Paris, UMR 7095, CNRS, and Sorbonne Université, 98 bis boulevard Arago, 75014 Paris, France

63 European Space Agency/ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands

64 Kapteyn Astronomical Institute, University of Groningen, PO Box 800, 9700 AV Groningen, The Netherlands

65 Department of Physics and Astronomy, University of Aarhus, Ny Munkegade 120, 8000 Aarhus C, Denmark

66 AIM, CEA, CNRS, Université Paris-Saclay, Université de Paris, 91191 Gif-sur-Yvette, France

67 Space Science Data Center, Italian Space Agency, via del Politecnico snc, 00133 Roma, Italy

68 CEA Saclay, DFR/IRFU, Service d'Astrophysique, Bât. 709, 91191 Gif-sur-Yvette, France

69 Institute of Space Science, Str. Atomistilor, nr. 409 Măgurele, Ilfov, 077125, Romania

70 Dipartimento di Fisica e Astronomia "G. Galilei", Università di Padova, Via Marzolo 8, 35131 Padova, Italy

71 INFN-Padova, Via Marzolo 8, 35131 Padova, Italy

72 Dipartimento di Fisica e Astronomia, Università di Bologna, Via Gobetti 93/2, 40129 Bologna, Italy

73 Departamento de Física, FCFM, Universidad de Chile, Blanco Encalada 2008, Santiago, Chile

74 Institut für Astro- und Teilchenphysik, Universität Innsbruck, Technikerstr. 25/8, 6020 Innsbruck, Austria

75 Institut de Ciencies de l'Espai (IEEC-CSIC), Campus UAB, Carrer de Can Magrans, s/n Cerdanyola del Vallés, 08193 Barcelona, Spain

76 Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Avenida Complutense 40, 28040 Madrid, Spain

77 Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências, Universidade de Lisboa, Tapada da Ajuda, 1349-018 Lisboa, Portugal

78 Universidad Politécnica de Cartagena, Departamento de Elétrónica y Tecnología de Computadoras, Plaza del Hospital 1, 30202 Cartagena, Spain

79 Institut de Recherche en Astrophysique et Planétologie (IRAP), Université de Toulouse, CNRS, UPS, CNES, 14 Av. Edouard Belin, 31400 Toulouse, France

80 Infrared Processing and Analysis Center, California Institute of Technology, Pasadena, CA 91125, USA

81 IFPU, Institute for Fundamental Physics of the Universe, via Beirut 2, 34151 Trieste, Italy

82 Instituto de Astrofísica de Canarias, Calle Vía Láctea s/n, 38204 San Cristóbal de La Laguna, Tenerife, Spain

83 Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, 10010, New York, NY, USA

84 School of Physics and Astronomy, Cardiff University, The Parade, Cardiff, CF24 3AA, UK

85 Junia, EPA department, 41 Bd Vauban, 59800 Lille, France

86 INFN-Bologna, Via Irnerio 46, 40126 Bologna, Italy

87 Instituto de Física Teórica UAM-CSIC, Campus de Cantoblanco, 28049 Madrid, Spain

88 CERCA/ISO, Department of Physics, Case Western Reserve University, 10900 Euclid Avenue, Cleveland, OH 44106, USA

89 Laboratoire de Physique de l'École Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, 75005 Paris, France

90 Observatoire de Paris, Université PSL, Sorbonne Université, LERMA, 750 Paris, France
91 Astrophysics Group, Blackett Laboratory, Imperial College London, London SW7 2AZ, UK

92 SISSA, International School for Advanced Studies, Via Bonomea 265, 34136 Trieste TS, Italy

93 INFN, Sezione di Trieste, Via Valerio 2, 34127 Trieste TS, Italy

94 Departamento de Astrofísica, Universidad de La Laguna, 38206 La Laguna, Tenerife, Spain

95 Dipartimento di Fisica e Scienze della Terra, Università degli Studi di Ferrara, Via Giuseppe Saragat 1, 44122 Ferrara, Italy

96 Istituto Nazionale di Fisica Nucleare, Sezione di Ferrara, Via Giuseppe Saragat 1, 44122 Ferrara, Italy

97 Institut de Physique Théorique, CEA, CNRS, Université Paris-Saclay 91191 Gif-sur-Yvette Cedex, France

98 NASA Ames Research Center, Moffett Field, CA 94035, USA

99 INAF, Istituto di Radioastronomia, Via Piero Gobetti 101, 40129 Bologna, Italy

100 Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Bd de l'Observatoire, CS 34229, 06304 Nice Cedex 4, France

101 Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu, HI 96822, USA

102 Department of Physics & Astronomy, University of California Irvine, Irvine CA 92697, USA

103 INFN-Sezione di Genova, Via Dodecaneso 33, 16146, Genova, Italy

104 Department of Astronomy & Physics and Institute for Computational Astrophysics, Saint Mary's University, 923 Robie Street, Halifax, Nova Scotia, B3H 3C3, Canada

105 University Observatory, Faculty of Physics, Ludwig-Maximilians-Universität, Scheinerstr. 1, 81679 Munich, Germany

106 Ruhr University Bochum, Faculty of Physics and Astronomy, Astronomical Institute (AIRUB), German Centre for Cosmological Lensing (GCCL), 44780 Bochum, Germany

107 Department of Physics and Astronomy, Vesiilinnantie 5, 20014 University of Turku, Finland

108 Univ. Grenoble Alpes, CNRS, Grenoble INP, LPSC-IN2P3, 53, Avenue des Martyrs, 38000 Grenoble, France

109 Centre de Calcul de l'IN2P3/CNRS, 21 avenue Pierre de Coubertin 69627 Villeurbanne Cedex, France

110 Dipartimento di Fisica, Sapienza Università di Roma, Piazzale Aldo Moro 2, 00185 Roma, Italy

111 Centro de Astrofísica da Universidade do Porto, Rua das Estrelas, 4150-762 Porto, Portugal

112 Dipartimento di Fisica – Sezione di Astronomia, Università di Trieste, Via Tiepolo 11, 34131 Trieste, Italy

113 Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth PO1 3FX, UK

114 Department of Mathematics and Physics E. De Giorgi, University of Salento, Via per Arnesano, CP-193, 73100 Lecce, Italy

115 INFN, Sezione di Lecce, Via per Arnesano, CP-193, 73100 Lecce, Italy

116 INAF-Sezione di Lecce, c/o Dipartimento Matematica e Fisica, Via per Arnesano, 73100 Lecce, Italy

117 Institute for Computational Science, University of Zurich, Winterthurerstrasse 190, 8057 Zurich, Switzerland

118 Higgs Centre for Theoretical Physics, School of Physics and Astronomy, The University of Edinburgh, Edinburgh EH9 3FD, UK

119 Université St Joseph; Faculty of Sciences, Beirut, Lebanon

120 Department of Astrophysical Sciences, Peyton Hall, Princeton University, Princeton, NJ 08544, USA
**Appendix A: Excluding the SLC contributions**

In this appendix, we reproduce the results shown in Figures 1 and 2, except with the SLC and SO-SLC terms removed (discussed in Section 2.4). With future studies, when calibrating the source number density weighting term, these terms could be negligible. In this case, 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/not applied. Hence, one should either include all corrections, or none, but including only some could result in biased results. In the $w_0w_a$ CDM case, we find significant biases could be caused by neglecting any of terms – in particular, the magnification bias.

<figure><img src="image_5.png" /><figcaption>Fig. A.1. Similar to Figure 1 except excluding SLC and SO-SLC terms.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. A.2. Similar to Figure 2 except excluding SLC and SO-SLC terms.</figcaption></figure>