© The Authors 2024

1 Introduction

Over recent decades, we have witnessed a remarkable improvement in the precision of cosmological experiments, and consequently in our grasp of the general properties of the Universe. The Λ cold dark matter (CDM) concordance cosmological model provides an exquisite fit to observational data from both the very early and the very late Universe, but despite its success, the basic components it postulates are poorly understood. Moreover, the nature of the mechanism responsible for the observed accelerated cosmic expansion (Riess et al. 1998; Perlmutter et al. 1999) and that of the component accounting for the vast majority of the matter content, dark matter, are still unknown. Upcoming Stage IV surveys like the Vera C. Rubin Observatory Legacy Survey of Space and Time (LSST, Ivezić et al. 2019), the Nancy Grace Roman Space Telescope (Spergel et al. 2015), and the Euclid mission (Laureijs et al. 2011; Euclid Collaboration 2024) promise to help deepen our understanding of these dark components and the nature of gravity on cosmo-logical scales by providing unprecedented observations of the large-scale structures (LSS) of the Universe.

Because of their high accuracy and precision, these next-generation experiments will require accurate modelling of both the theory and the covariance of the observables under study in order to produce precise and unbiased estimates of the cosmo-logical parameters. Amongst the different theoretical issues to deal with is super-sample covariance (SSC), a form of sample variance arising from the finiteness of the survey area. SSC was first introduced for cluster counts by Hu & Kravtsov (2003), and is sometimes referred to as ‘beat coupling’ (Rimes & Hamilton 2006; Hamilton et al. 2006). In recent years, SSC has received a lot of attention (Takada & Hu 2013; Li et al. 2014; Barreira et al. 2018b; Digman et al. 2019; Bayer et al. 2023; Yao et al. 2024); see also Linke et al. (2024) for an insightful discussion on SSC in real space. Hereafter, Barreira et al. (2018b) is cited as B18.

The effect arises from the coupling between ‘supersurvey’ modes - with wavelength λ larger than the survey typical size $L=V_{\text{s}}^{1/3}$ (where V_s is the volume of the survey) – and short-wavelength (λ < L) modes. This coupling is in turn due to the significant non-linear evolution undergone by low-redshift cosmological probes (contrary to, for example, the cosmic microwave background), which breaks the initial homogeneity of the density field, making its growth position dependent. In Fourier space, this means that modes with different wavenum-ber k = 2π/λ become coupled. The modulation induced by the supersurvey modes is equivalent to a change in the background density of the observed region, which affects and correlates all LSS probes. It is accounted for as an additional, non-diagonal term in the data covariance matrix beyond the Gaussian covari-ance, which is the only term that would exist if the random field under study were Gaussian. Being the most affected by nonlinear dynamics, the smaller scales are heavily impacted by SSC, where the effect is expected to be the dominant source of statistical uncertainty for the two-point statistics of weak-lensing cosmic shear (WL): it has in fact been found to increase conditional uncertainties by up to a factor of about 2 (for a Euclid-like survey, see Barreira et al. 2018a; Gouyou Beauchamps et al. 2022). In the case of photometric galaxy clustering (GCph; again, for a Euclid-like survey), Lacasa & Grain (2019) – hereafter LG19 – found the cumulative signal-to-noise ratio to be decreased by a factor of around 6 at ℓ_max = 2000. These works, however, either do not take into account marginalised uncertainties or the variability of the probe responses, do not include cross-correlations between probes, or do not follow the full specifics (such as modelling of the observables, types of sys-tematics included, binning schemes, sky coverage and so forth) of the Euclid survey detailed below.

There are two aims to the present study. First, we intend to validate the forecast constraints on the cosmological parameters, both including and neglecting the SSC term; these are produced using two independent codes, whose only shared feature is their use of the public Python module PySSC¹,² (LG19) to compute the fundamental elements needed to build the SSC matrix. Second, we investigate the impact of SSC on the marginalised uncertainties and the dark energy figure of merit (FoM), both of which are obtained through a Fisher forecast of the constraining power of Euclid’s photometric observables.

The article is organised as follows: Sect. 2 presents an overview of the SSC and the approximations used to compute it. In Sect. 3 we outline the theoretical model and specifics used to produce the forecasts, while Sect. 4 provides technical details regarding the implementation and validation of the code. In Sect. 5, we then present a study of the impact of SSC on Euclid constraints for different binning schemes and choices of systematic errors and priors. Finally, we present our conclusions in Sect. 6.

2 SSC theory and approximations

2.1 General formalism

Throughout the article, we work with 2D-projected observables, namely the angular Power Spectrum (PS), which in the Limber approximation (Limber 1953; Kaiser 1998) can be expressed as $$C_{ij}^{AB}(\ell )={{\displaystyle \int }}^{\text{}}\text{d}V{W}_i^A(ɀ)W_j^B(ɀ)P_{AB}\left(k_{\ell },ɀ\right),$$(1)

giving the correlation between probes A and B in the redshift bins i and j, as a function of the multipole ℓ; k_ℓ = (ℓ + 1/2)/r(z) is the Limber wavenumber and $W_i^A(z),W_j^B(z)$ are the survey weight functions (WFs), or “kernels”. Here we consider as the element of integration $\text{d}V=r^2(z)\frac{\text{d}r}{\text{d}z}\text{d}z$ which is the comoving volume element per steradian, with r(z) being the comoving distance.

The SSC between two projected observables arises because real observations of the Universe are always limited by a survey window function ℳ(x). Taking ℳ(x) at a given redshift, thus considering only its angular dependence $ℳ(\widehat{n})$³, with $\widehat{n}$ the unit vector on a sphere, we can define the background density contrast as (Lacasa et al. 2018) $${\delta }_{\text{b}}(ɀ)=\frac{1}{{\Omega }_{\text{S}}}{{\displaystyle \int }}^{\text{}} {\text{d}}^2\widehat{n}ℳ(\widehat{\mathit{n}}){\delta }_{\text{m}}[r(ɀ)\widehat{\mathit{n}},ɀ],$$(2)

with $r(z)\widehat{n}=x$. In this equation, ${\delta }_{\text{m}}(x,z)=\left[{\rho }_{\text{m}}(x,z)/{\overline{\rho }}_{\text{m}}(z)-1\right]$ is the matter density contrast, with ρ_m(x, z) the matter density and ${\overline{\rho }}_{\text{m}}(z)$ its spatial average over the whole Universe at redshift z and Ω_S the solid angle observed by the survey.

In other words, δ_b is the spatial average of the density contrast δ_m(x, z) over the survey area: $${\langle {\delta }_{\text{m}}(\mathit{x},ɀ)\rangle }_{\text{universe }}=0,$$(3) $${\langle {\delta }_{\text{m}}(x,ɀ)\rangle }_{\text{survey }}={\delta }_{\text{b}}(ɀ).$$(4)

The covariance of this background density contrast is defined as σ²(z₁, z₂) ≡ 〈δ_b(z₁) δ_b(z₂)〉 and in the full-sky approximation is given by (Lacasa & Rosenfeld 2016) $${\sigma }^2\left({ɀ}_1,{ɀ}_2\right)=\frac{1}{2{\pi }^2}{\displaystyle \int \text{d}}k{k}^2{P}_{\text{mm}}^{\mathrm{lin}}\left(k,{ɀ}_{12}\right){\text{j}}_0\left(k{r}_1\right){\text{j}}_0\left(k{r}_2\right),$$(5)

with $P_{\text{mm}}^{\text{lin}}\left(k,z_{12}\right)\equiv D\left(z_1\right)D\left(z_2\right)P_{\text{mm}}^{\text{lin}}(k,z=0)$ the linear matter cross-spectrum between z₁ and z₂, D(z) the linear growth factor and j₀(kr_i) the first-order spherical Bessel function, and r_i = r(z_i). The use of the linear PS reflects the fact that the SSC is caused by long-wavelength perturbations, which are well described by linear theory. We note that we have absorbed the ${\text{ΩS}}_{-1}^{}$ prefactor of Eq. (2), equal to 4π in full sky, in the dV_i terms, being them the comoving volume element per steradian.

Depending on the portion of the Universe observed, δ_b will be different, and in turn the PS of the considered observables P_AB(k_ℓ, z) (appearing in Eq. (1)) will react to this change in the background density through the probe response ∂P_AB(k_ℓ, z)/∂δ_b.

SSC is then the combination of these two elements, encapsulating the covariance of δ_b and the response of the observables to a change in δ_b; the general expression of the SSC between two projected observables is (Lacasa & Rosenfeld 2016; Schaan et al. 2014; Takada & Hu 2013): $$\begin{array}{l}{\mathrm{Cov}}_{\text{SSC}}\left[C_{ij}^{AB}(\ell ),C_{kl}^{CD}\left({\ell }^{\prime }\right)\right]={\displaystyle \int \text{d}}{V}_1 \text{d}{V}_2{W}_i^A\left({ɀ}_1\right)W_j^B\left({ɀ}_1\right)\\ \times W_k^C\left({ɀ}_2\right)W_l^D\left({ɀ}_2\right)\frac{\partial P_{AB}\left(k_{\ell },{ɀ}_1\right)}{\partial {\delta }_{\text{b}}}\frac{\partial P_{CD}\left(k_{{\ell }^{\prime }},{ɀ}_2\right)}{\partial {\delta }_{\text{b}}}{\sigma }^2\left({ɀ}_1,{ɀ}_2\right).\end{array}$$(6)

We adopt the approximation presented in Lacasa & Grain (2019), which assumes the responses to vary slowly in redshift with respect to σ²(z₁, z₂). We can then approximate the responses with their weighted average over the $W_i^A(z)$ kernels (Gouyou Beauchamps et al. 2022): $$\frac{\partial {\overline{P}}_{AB}\left(k_{\ell },ɀ\right)}{\partial {\delta }_{\text{b}}}=\frac{{\displaystyle \int \text{d}}V{W}_i^A(ɀ)W_j^B(ɀ)\partial P_{AB}\left(k_{\ell },ɀ\right)/\partial {\delta }_{\text{b}}}{{\displaystyle \int \text{d}}V{W}_i^A(ɀ)W_j^B(ɀ)},$$(7)

and pull them out of the integral. The denominator on the right-hand side (r.h.s.) acts as a normalisation term, which we call $I_{ij}^{AB}$. We can further manipulate the above expression by factorising the probe response as $$\frac{\partial P_{AB}\left(k_{\ell },ɀ\right)}{\partial {\delta }_{\text{b}}}=R^{AB}\left(k_{\ell },ɀ\right)P_{AB}\left(k_{\ell },ɀ\right),$$(8)

where R^AB(k_ℓ, z), the “response coefficient”, can be obtained from simulations, as in Wagner et al. (2015a,b); Li et al. (2016); Barreira et al. (2019), or from theory (e.g. via the halo model) as in Takada & Hu (2013); Krause & Eifler (2017); Rizzato et al. (2019). Following LG19, we can introduce the probe response of the angular power spectrum $C_{ij}^{AB}\left(\ell \right)$ in a similar way, using Eq. (1) $$\begin{array}{cc}\frac{\partial C_{ij}^{AB}(\ell )}{\partial {\delta }_{\text{b}}}\hfill & ={{\displaystyle \int }}^{\text{}}\text{d}V{W}_i^A(ɀ)W_j^B(ɀ)\frac{\partial P_{AB}\left(k_{\ell },ɀ\right)}{\partial {\delta }_{\text{b}}}\hfill \\ \hfill & \equiv R_{ij}^{AB}(\ell )C_{ij}^{AB}(\ell ).\hfill \end{array}$$(9)

Substituting Eq. (8) into the r.h.s. of Eq. (7), using Eq. (9) and dividing by the sky fraction observed by the telescope f_sky = Ω_S/4π, we get the expression of the SSC which will be used throughout this work: $$\begin{array}{l}Co{v}_{\text{SSC}}\left[C_{ij}^{AB}(\ell )C_{kl}^{CD}\left({\ell }^{\prime }\right)\right]\simeq f_{sky}^{-1}[R_{ij}^{AB}(\ell )C_{ij}^{AB}(\ell )\hfill \\ \times R_{kl}^{CD}\left({\ell }^{\prime }\right)C_{kl}^{CD}\left({\ell }^{\prime }\right)S_{i,j,k,l}^{A,B,C,D}].\hfill \end{array}$$(10)

In the above equation, we define $$S_{i,j;k,l}^{A,B,C,D}\equiv {{\displaystyle \int }}^{\text{}}\text{d}{V}_1\text{d}{V}_2\frac{W_i^A\left({ɀ}_1\right)W_j^B\left({ɀ}_1\right)}{I_{ij}^{AB}}\frac{W_k^C\left({ɀ}_2\right)W_l^D\left({ɀ}_2\right)}{I_{kl}^{CD}}{\sigma }^2\left({ɀ}_1,{ɀ}_2\right).$$(11)

The $S_{i,j;k,l}^{A,B;C,D}$ matrix (referred to as S_ijkl from here on) is the volume average of σ²(z₁, z₂), and is a dimensionless quantity. It is computed through the public Python module PySSC, released alongside the above-mentioned LG19. A description of the way this code has been used, and some comments on the inputs to provide and the outputs it produces, can be found in Sect. 4.

The validity of Eq. (10) has been tested in LG19 in the case of GCph and found to reproduce the Fisher matrix (FM, Tegmark et al. 1997) elements and signal-to-noise ratio from the original expression (Eq. 6):

• within 10% discrepancy up to ℓ ≃ 1000 for $R_{ij}^{AB}\left(k_{\ell },z\right)=$ const;

• within 5% discrepancy up to ℓ ≃ 2000 when using the linear approximation in scale for R^AB(k_ℓ, z) provided in Appendix C of the same work.

The necessity to push the analysis to smaller scales, as well as to investigate the SSC impact not only for GCph but also for WL and their cross-correlation, has motivated a more exhaustive characterisation of the probe response functions, which will be detailed in the next section.

Another approximation used in the literature has been presented in Krause & Eifler (2017): the σ²(z₁, z₂) term is considered as a Dirac delta in z₁ = z₂. This greatly simplifies the computation, because the double redshift integral dV₁dV₂ collapses to a single one. This approximation is used by the other two available public codes which can compute the SSC: PyCCL (Chisari et al. 2019) and CosmoLike (Krause & Eifler 2017). Lacasa et al. (2018) compared this approximation against the one used in this work, finding the former to fare better for wide red-shift bins (as in the case of WL), and the latter for narrow bins (as in the case of GCph).

Lastly, we note that in Eq. (10) we account for the sky coverage of the survey through the full-sky approximation by simply dividing by ƒ_sky; in the case of Euclid we have Ω_S = 14 700 deg² ≃ 4.4776 sr, which corresponds to ƒ_sky ≃ 0.356. The validity of this approximation has been discussed in Gouyou Beauchamps et al. (2022), and found to agree at the percent level on the marginalised parameter constraints with the more rigorous treatment accounting for the exact survey geometry when considering large survey areas. For this test, they considered an area of 15 000 deg² and a survey geometry very close to what Euclid will have, i.e. the full sky with the ecliptic and galactic plane removed. Intuitively, the severity of the SSC decays as $f_{\text{sky}}^{-1}$ because larger survey volumes can accommodate more Fourier modes.

We note that we are considering here the maximum sky coverage that Euclid will reach, i.e. the final data release (DR3). For the first data release (DR1), the sky coverage will be significantly lower and the full-sky approximation will not hold. In that case, the partial-sky recipe proposed in Gouyou Beauchamps et al. (2022) should be considered instead.

2.2 Probe response

As mentioned in the previous section, one of the key ingredients of the SSC is the probe response. To compute this term for the probes of interest, we build upon previous works (Wagner et al. 2015a,b; Li et al. 2016; Barreira & Schmidt 2017, B18), and compute the response coefficient of the matter PS as $$R^{\text{mm}}(k,ɀ)=\frac{\partial \ln P_{\text{mm}}(k,ɀ)}{\partial {\delta }_{\text{b}}}=1-\frac{1}{3}\frac{\partial \ln P_{\text{mm}}(k,ɀ)}{\partial \ln k}+G_1^{\text{mm}}(k,ɀ),$$(12)

where $G_1^{\text{mm}}(k,z)$ is called the growth-only response, and is constant and equal to 26/21 in the linear regime and can be computed in the non-linear regime using separate universe simulations, as done in Wagner et al. (2015b), whose results were used in B18 (and in the present work). The latter uses a power law to extrapolate the values of the response for k > k_max, with k_max being the maximum wavenumber at which the power spectrum is reliably measured from the simulations. Further details on this extrapolation, as well as on the redshift and scale dependence of R^mm, can be found respectively in Sect. 2 and the left panel of Fig. 1 of B18. We note that R^mm is the response coefficient of isotropic large-scale density perturbations; we neglect the contribution from the anisotropic tidal-field perturbations to the total response of the power spectrum (and consequently to the SSC), which has been shown in B18 to be subdominant for WL with respect to the first contribution (about 5% of the total covariance matrix at ℓ ≳ 300). While we do not expect this conclusion to change substantially for GCph, we leave an accurate assessment for future work.

The probes considered in the present study are WL, GCph and their cross-correlation (XC); assuming general relativity, the corresponding power spectra are given by the following expressions $$P_{AB}(k,ɀ)=\{\begin{array}{ll}{P}_{\text{mm}}(k,ɀ)\hfill & A=B=\text{L}\hfill \\ b_{(1)}(ɀ)P_{\text{mm}}(k,ɀ)\hfill & A=\text{L},B=\text{G}\hfill \\ b_{(1)}^2(ɀ)P_{\text{mm}}(k,ɀ)\hfill & A=B=\text{G},\hfill \end{array}$$(13)

with (L, G) for (shear, position), P_mm(k, z) the non-linear matter PS and b_(1)(z) the linear, scale-independent and deterministic galaxy bias. A comment is in order about the way we model the galaxy-matter and galaxy-galaxy power spectra. We are indeed using a linear bias, but the non-linear recipe for the matter power spectrum P_mm(k, z). This is reminiscent of the hybrid 1-loop perturbation theory (PT) model adopted by, for example, the DES Collaboration in the analysis of the latest data release (Krause et al. 2021; Pandey et al. 2022), but we drop the higher-order bias terms. This simplified model has been chosen in order to be consistent with the IST:F (Euclid Collaboration 2020, from hereon EC20) forecasts, against which we compare our results (in the Gaussian case) to validate them. We are well aware that scale cuts should be performed to avoid biasing the constraints, but we are here more interested in the relative impact of SSC on the constraints than the constraints themselves. Any systematic error due to the approximate modelling should roughly cancel out in the ratio we compute later on. We note also that we choose to include a perfectly Poissonian shot noise term in the covariance matrix, rather than in the signal, as can be seen in Eq. (25). The responses for the different probes can be obtained in terms⁴ of R^mm(k, z) by using the relations between matter and galaxy PS given above $$R^{\text{gg}}(k,ɀ)=\frac{\partial \ln P_{\text{gg}}(k,ɀ)}{\partial {\delta }_{\text{b}}}=R^{\text{mm}}(k,ɀ)+2b_{(1)}^{-1}(ɀ)\left[b_{(2)}(ɀ)-b_{(1)}^2(ɀ)\right],$$(14)

and similarly for R^gm: $$R^{\text{gm}}(k,ɀ)=\frac{\partial \ln P_{\text{gm}}(k,ɀ)}{\partial {\delta }_{\text{b}}}=R^{\text{mm}}(k,ɀ)+b_{(1)}^{-1}(ɀ)\left[b_{(2)}(ɀ)-b_{(1)}^2(ɀ)\right].$$(15)

Having used the definitions of the first- and second-order galaxy bias, that is, b₍₁₎(z) = (∂n_g/∂δ_b)/n_g and $b_{(2)}(z)=\left({\partial }^2{n}_{\text{g}}/\partial {\delta }_{\text{b}}^2\right)/n_{\text{g}}$, with n_g the total angular galaxy number density, in arcmin⁻². In the following, where there is no risk of ambiguity, we drop the subscript in parenthesis when referring to the first-order galaxy bias – that is, b(z) = b₍₁₎(z) – to shorten the notation, and we indicate the value of the first-order galaxy bias in the i-th red-shift bin with b_i(z). More details on the computation of these terms can be found in Sect. 3.6. We note that Eqs. (14)–(15) are obtained by differentiating a PS model for a galaxy density contrast defined with respect to (w.r.t.) the observed galaxy number density, and so they already account for the fact that the latter also “responds” to the large-scale perturbation δ_b. This is also the reason why $R_{ij}^{\text{GG}}(\ell )$ can have negative values: for galaxy clus-tering, the (number) density contrast δ_gal is measured w.r.t. the observed, local number density ${\overline{n}}_{\text{gal}}:{\delta }_{\text{gal}}=n_{\text{gal}}/{\overline{n}}_{\text{gal}}-1$. The latter also responds to a background density perturbation δ_b, and it can indeed happen that ${\overline{n}}_{\text{gal}}$ grows with δ_b faster than n_gal, which leads to δ_gal decreasing with increasing δ_b (which also implies $\partial C_{ij}^{\text{GG}}(\ell )/\partial {\delta }_{\text{b}}<0)$. We also stress the fact that the second-order galaxy bias appearing in the galaxy-galaxy and galaxy-lensing response coefficients is not included in the signal, following EC20. Once computed in this way, the response coefficient can be projected in harmonic space using Eq. (9), and inserted in Eq. (10) to compute the SSC in the LG19 approximation. The projected $R_{ij}^{AB}(\ell )$ functions are shown in Fig. 1.

Fig. 1 Projected response coefficients for the WL and GCph probes and their cross-correlation for the central redshift bin (0.8 ≲ z ≲ 0.9). The shape and amplitude of the functions for different redshift pairs are analogous. For WL, the baryon acoustic oscillation wiggles are smoothed out by the projection, because the kernels are larger than the GCph ones. The different amplitude of the response is one of the main factors governing the severity of SSC.

3 Forecasts specifics

In order to forecast the uncertainties in the measurement of the cosmological parameters, we follow the prescriptions of the Euclid forecast validation study (EC20), with some updates to the most recent results from the EC, which are used for the third Science Performance Verification (SPV) of Euclid before launch. In particular, the update concerns the fiducial value of the linear bias, the redshift distribution n(z) and the multipole binning.

Once again, the observable under study is the angular PS of probe A in redshift bin i and probe B in redshift bin j, given in the Limber approximation by Eq. (1). The P_AB(k_ℓ, z) multi-probe power spectra are given in Eq. (13); in the following, we refer interchangeably to the probes (WL, XC, GCph) and their auto-and cross-spectra (respectively, LL, GL, GG).

3.1 Redshift distribution

First, we assume that the same galaxy population is used to probe both the WE and the GCph PS. We therefore set $$n_i^{\text{L}}(ɀ)=n_i^{\text{G}}(ɀ)=n_i(ɀ),$$(16)

where $n_i^{\text{L}}(z)\text{ and }{n}_i^{\text{G}}(z)$ are respectively the distribution of sources and lenses in the i-th redshift bin. Then, the same equality applies for the total source and lens number density, ${\overline{n}}^{\text{L}}\text{ and }{\overline{n}}^{\text{G}}$.

A more realistic galaxy redshift distribution than the analytical one presented in EC20 can be obtained from simulations. We use the results from Euclid Collaboration (2021), in which the n(z) is constructed from photometric redshift estimates in a 400 deg² patch of the Flagship 1 simulation (Potter et al. 2017), using the training-based directional neighbourhood fitting (DNF) algorithm (De Vicente et al. 2016).

The training set is a random subsample of objects with true (spectroscopic) redshifts known from the Flagship simulation. We choose the fiducial case presented in Euclid Collaboration (2021), which takes into account a drop in completeness of the spectroscopic training sample with increasing magnitude. A cut in magnitude I_E < 24.5, isotropic and equal for all photometric bands, is applied, corresponding to the optimistic Euclid setting. The DNF algorithm then produces a first estimate of the photo-z, z_mean, using as a metric the objects’ closeness in colour and magnitude space to the training samples. A second estimate of the redshift, z_mc, is computed from a Monte Carlo draw from the nearest neighbour in the DNF metric. The final distributions for the different redshift bins, n_i(z), are obtained by assigning the sources to the respective bins using their Z_mean, and then taking the histogram of the z_mc values in each of the bins – following what has been done in real surveys such as the Dark Energy Survey (Crocce et al. 2019; Hoyle et al. 2018).

As a reference setting, we choose to bin the galaxy distribution into 𝒩_b = 10 equipopulated redshift bins, with edges $$\begin{array}{cc}{ɀ}_{\text{edges }}& =\{0.001,0.301,0.471,0.608,0.731,0.851,\\ & 0.980,1.131,1.335,1.667,2.501\}.\end{array}$$(17)

The total galaxy number density is $\overline{n}=28.73{\text{arcmin}}^{-2}$. As a comparison, this was set to 30 arcmin⁻² in EC20. We note that this choice of redshift binning will be discussed and varied in Sect. 5.3.

3.2 Weight functions

We model the radial kernels, or weight functions, for WL and GCph following once again EC20. Adopting the eNLA (extended non-linear alignment) prescription for modelling the intrinsic alignment (IA) contribution, the weight function ${\mathcal{W}}_i^A(z)$ for the lensing part is given by (see e.g. Kitching et al. 2017; Kilbinger et al. 2017; Taylor et al. 2018b) $${\mathcal{W}}_i^{\text{L}}(ɀ)={\mathcal{W}}_i^{\gamma }(ɀ)-\frac{{\mathcal{A}}_{\text{IA}}{𝒞}_{\text{IA}}{\Omega }_{\text{m},0}{ℱ}_{\text{IA}}(ɀ)}{D(ɀ)}{\mathcal{W}}^{\text{IA}}(ɀ),$$(18)

where we define⁵ $${\mathcal{W}}_i^{\gamma }(ɀ)=\frac{3}{2}{\left(\frac{H_0}{c}\right)}^2{\Omega }_{\text{m},0}(1+ɀ)r(ɀ){\displaystyle {\int }_{ɀ}^{{ɀ}_{\max }}\frac{n_i\left({ɀ}^{\prime }\right)}{\overline{n}}}\left[1-\frac{r(ɀ)}{r\left({ɀ}^{\prime }\right)}\right]\text{d}{ɀ}^{\prime },$$(19)

and $${\mathcal{W}}_i^{\text{IA}}(ɀ)=\frac{1}{c}\frac{n_i(ɀ)}{\overline{n}}H(ɀ).$$(20)

Finally, in Eq. (18), 𝒜_IA is the overall IA amplitude, C_IA a constant, ℱ_IA(z) a function modulating the dependence on redshift, and D(z) is the linear growth factor. More details on the IA modelling are given in Sect. 3.5.

The GCph weight function is equal to the IA one, as long as Eq. (16) holds: $${\mathcal{W}}_i^{\text{G}}(ɀ)={\mathcal{W}}_i^{\text{IA}}(ɀ)=\frac{1}{c}\frac{n_i(ɀ)}{\overline{n}}H(ɀ).$$(21)

Figure 2 shows the redshift dependence of Eqs. (18) and (21), for all redshift bins. We note that we choose to include the galaxy bias term b_i(z) in the PS (see Eq. (13)) rather than in the galaxy kernel, as opposed to what has been done in EC20. This is done to compute the galaxy response as described in Sect. 2.2. However, as the galaxy bias is assumed constant in each bin, the question is of no practical relevance when computing the S_ijkl matrix, since the constant bias cancels out.

We note that the above definitions of the lensing and galaxy kernels $\left({\mathcal{W}}_i^A(z),A=\text{L},\text{G}\right)$ differ from the ones used in LG19. This is simply because of a different definition of the $C_{ij}^{AB}(\ell )$ Limber integral, which is performed in dV in LG19 and in dz in EC20. The mapping between the two conventions is simply given by the expression for the volume element: $$\text{d}V=r^2(ɀ)\frac{\text{d}r}{ \text{d}ɀ} \text{d}ɀ=c\frac{r^2(ɀ)}{H(ɀ)}\text{d}ɀ,$$(22)

and $$W_i^A(ɀ)={\mathcal{W}}_i^A(ɀ)/r^2(ɀ),$$(23)

with A = L, G. In Fig. 2 we plot the values of ${\mathcal{W}}_i^A(z)$ to facilitate the comparison with EC20. As outlined in Appendix A, when computing the S_ijkl matrix through PySSC, the user can either pass the kernels in the form used in LG19 or the one used in EC20 – specifying a non-default convention parameter.

Fig. 2 Kernels and galaxy distribution considered in this work. The first two plots show the kernels, or weight for the two photometric probes. The analytic expressions for these are, respectively, Eqs. (18) (left, WL) and (21) (right, GCph). At high redshifts, the IA term dominates over the shear term in the lensing kernels, making them negative. The rightmost plot shows the redshift distribution per redshift bin for the sources (and lenses), as well as their sum, obtained from the Flagship 1 simulation as described in Sect. 3.1.

3.3 Gaussian covariance

The Gaussian part of the covariance is given by the following expression: $$Co{v}_{\text{G}}\left[{\widehat{C}}_{ij}^{AB}(\ell ),{\widehat{C}}_{kl}^{CD}\left({\ell }^{\prime }\right)\right]={\left[(2\ell +1)f_{\text{sky }}\Delta \ell \right]}^{-1}{\delta }_{\ell {\ell }^{\prime }}^{\text{K}}$$ $$\begin{array}{c}\times \{\left[C_{ik}^{AC}(\ell )+N_{ik}^{AC}(\ell )\right]\left[C_{jl}^{BD}\left({\ell }^{\prime }\right)+N_{jl}^{BD}\left({\ell }^{\prime }\right)\right]\\ +\left[C_{il}^{AD}(\ell )+N_{il}^{AD}(\ell )\right]\left[C_{jk}^{BC}\left({\ell }^{\prime }\right)+N_{jk}^{BC}\left({\ell }^{\prime }\right)\right]\},\end{array}$$(24)

where we use a hat to distinguish the estimators from the true spectra. The noise PS $N_{ij}^{AB}(\ell )$ are, for the different probe combinations, $$N_{ij}^{AB}(\ell )=\{\begin{array}{ll}\left({\sigma }_{ϵ}^2/2{\overline{n}}_i^{\text{L}}\right){\delta }_{ij}^{\text{K}}\hfill & A=B=\text{L}(\text{WL})\hfill \\ 0\hfill & A\ne B\hfill \\ \left(1/{\overline{n}}_i^{\text{G}}\right){\delta }_{ij}^{\text{K}}\hfill & A=B=\text{G}(\text{GCph}).\hfill \end{array}$$(25)

In the above equations, ${\delta }_{ij}^{\text{K}}$ is the Kronecker delta and ${\sigma }_{ϵ}^2$ the variance of the total intrinsic ellipticity dispersion of WL sources, where ${\sigma }_{ϵ}=\sqrt{2}{\sigma }_{ϵ}^{(i)}$, with ${\sigma }_{ϵ}^{(i)}$ being the ellipticity dispersion per component of the galaxy ellipse. Some care is needed when defining the shear noise spectrum: the above equation can then also be written as $N_{ij}^{\text{LL}}(\ell )=\left[{\left({\sigma }_{ϵ}^{(i)}\right)}^2/{\overline{n}}_i^{\text{L}}\right]{\delta }_{ij}^{\text{K}}$, that is, using the ellipticity dispersion per component instead of the total one, which is the appropriate choice for harmonic-space analyses (Hu & Jain 2004; Joachimi & Bridle 2010). We note that the average densities used in Eq. (25) are not the total number densities, but rather those in the ¿-th redshift bin. In the case of 𝒩_b, equipop-ulated redshift bins, they can be simply written as ${\overline{n}}_i^A={\overline{n}}^A/{\mathcal{N}}_{\text{b}}$ for both A = (L, G). Finally, we recall that ƒ_sky is the fraction of the total sky area covered by the survey, while ∆ℓ is the width of the multipole bin centred on a given ℓ. From Sect. 3.1 we have that $\overline{n}=28.73{\mathrm{arcmin}}^{-2}$, while we set σ_є = 0.37 (from the value ${\sigma }_{ϵ}^{(i)}=0.26$ reported in Euclid Collaboration 2019) and ƒ_sky = 0.356 (corresponding to Ω_S = 14 700 deg²). We have now all the relevant formulae for the estimate of the Gaussian and the SSC terms of the covariance matrix. To ease the computation of Eq. (24) we have prepared an optimised Python module, Spaceborne_covg⁶, available as a public repository.

In the context of the present work, we do not consider the other non-Gaussian contribution to the total covariance matrix, the so-called connected non-Gaussian (cNG) term. This additional non-Gaussian term has been shown to be subdominant with respect to the Gaussian and SSC terms for WL both in Barreira et al. (2018a) and in Upham et al. (2022). For what concerns galaxy clustering, Wadekar et al. (2020) showed that the cNG term was subdominant, but this was for a spectroscopic sample so (i) they had a much larger contribution from shot-noise-related terms compared to what is considered here for the Euclid photometric sample, and (ii) they considered larger and more linear scales than in the present study. Lacasa (2020) showed that the cNG term in the covariance matrix of GCph only impacts the spectral index n_s and HOD parameters, but there are a few differences between that analysis and the present work, such as the modelling of galaxy bias. Thus it is still unclear whether the cNG term has a strong impact on cosmological constraints obtained with GCph. Quantifying the impact of this term for the 3×2pt analysis with Euclid settings is left for future work.

3.4 Cosmological model and matter power spectrum

We adopt a flat w₀w_aCDM model, that is, we model the dark energy equation of state with a Chevallier-Polarski-Linder (CPL) parametrisation (Chevallier & Polarski 2001; Linder 2005): $$w(ɀ)=w_0+w_aɀ/(1+ɀ).$$(26)

We also include a contribution from massive neutrinos with total mass equal to the minimum allowed by oscillation experiments (Esteban et al. 2020) ∑m_v = 0.06 eV, which we do not vary in the FM analysis. The vector of cosmological parameters is then $${\theta }_{\text{cosmo}}=\left\{{\text{Ω}}_{\text{m},0},{\text{Ω}}_{\text{b},0},w_0,w_a,h,n_{\text{s}},{\sigma }_8\right\},$$(27)

with Ω_m,0 and Ω_b,0 being respectively the reduced density of total and baryonic matter today, h is the dimensionless Hubble parameter defined as H₀ = 100 h km s⁻¹ Mpc⁻¹ where H₀ is the value of the Hubble parameter today, n_s the spectral index of the primordial power spectrum and σ₈ the root mean square of the linear matter density field smoothed with a sphere of radius 8 h⁻¹ Mpc. We follow EC20 for their fiducial values that are $${\theta }_{\text{cosmo }}^{\text{fid}}=\{0.32,0.05,-1.0,0.0,0.67,0.96,0.816\}.$$(28)

This parameter vector then is used as input for the evaluation of the fiducial linear and non-linear matter PS; for the purpose of validating our forecasts against the EC20 results, we use the TakaBird recipe, that is, the HaloFit version updated in Takahashi et al. (2012) with the Bird et al. (2012) correction for massive neutrinos. For the results shown in this paper, however, we update the non-linear model to the more recent HMCode2020 recipe, (Mead et al. 2021), which includes a baryonic correction parameterised by the log₁₀(T_AGN/K) parameter, characterising the feedback from active galactic nuclei (AGN). This is implemented in CAMB⁷ (Lewis et al. 2000) and, at the time of writing, is planned to be included in CLASS⁸ (Blas et al. 2011) as well. Because of this, we add a further free parameter in the analysis, log₁₀(T_AGN/K), with a fiducial value of 7.75.

3.5 Intrinsic alignment model

We use the eNLA model as in EC20, setting 𝒞_IA = 0.0134 and $${ℱ}_{\text{IA}}(ɀ)={(1+ɀ)}^{{\eta }_{\text{IA}}}{\left[\langle L\rangle (ɀ)/L_{\star }(ɀ)\right]}^{{\beta }_{\text{IA}}},$$(29)

where 〈L〉(z)/L_⋆(z) is the redshift-dependent ratio of the mean luminosity over the characteristic luminosity of WL sources as estimated from an average luminosity function (see e.g. Joachimi et al. 2015, and references therein). The IA nuisance parameters vector is $${\theta }_{\text{IA}}=\left\{{\mathcal{A}}_{\text{IA}},{\eta }_{\text{IA}},{\beta }_{\text{IA}}\right\},$$(30)

with fiducial values – following EC20 $${\theta }_{\text{IA}}^{\text{fid}}=\{1.72,-0.41,2.17\}.$$(31)

All of the IA parameters except for 𝒞_IA are varied in the analysis.

3.6 Linear galaxy bias and multiplicative shear bias

Following EC20 we model the galaxy bias as scale-independent. As for the redshift dependence, we move beyond the simple analytical prescription of EC20 and use the fitting function presented in Euclid Collaboration (2021), obtained from direct measurements from the Euclid Flagship galaxy catalogue, based in turn on the Flagship 1 simulation: $$b(ɀ)=\frac{A{ɀ}^B}{1+ɀ}+C,$$(32)

setting (A, B, C) = (0.81,2.80,1.02).

The galaxy bias is modelled to be constant in each bin with the fiducial value obtained by evaluating Eq. (32) at effective values $z_i^{\text{eff}}$ computed as the median of the redshift distribution considering only the part of the distribution at least larger than 10% of its maximum. We choose to use the median instead of the mean since, for equipopulated bins – as can be seen from the rightmost panel of Fig. 2 – the galaxy distribution of the last bins is highly skewed, and the value of the bias computed at the mean is potentially less accurate; this choice does not, on the other hand, affect the galaxy bias values in the first bins sensibly.

The $z_i^{\text{eff}}$ values obtained in this way are $$\begin{array}{l}{ɀ}^{\text{eff}}=\{0.212,0.363,0.447,0.566,0.682\hfill \\ 0.793,0.910,1.068,1.194,1.628\}.\hfill \end{array}$$(33)

We therefore have 𝒩_b additional nuisance parameters: $${\theta }_{\text{gal}.\text{ bias }}=\left\{b_1,b_2,\dots ,b_{{\mathcal{N}}_{\text{b}}}\right\},$$(34)

with fiducial values $$\begin{array}{l}{\theta }_{\text{gal}.\text{ bias }}^{\text{fid }}=\{1.031,1.057,1.081,1.128,1.187,\hfill \\ 1.258,1.348,1.493,1.628,2.227\}.\hfill \end{array}$$(35)

The modelling of galaxy bias just described is the same used in EC20, with different fiducial values.

We can take a further step forward towards the real data analysis by including the multiplicative shear bias parameters, m, defined as the multiplicative coefficient of the linear bias expansion of the shear field γ (see e.g. Cragg et al. 2023): $$\widehat{\mathit{\gamma }}=(1+m)\mathit{\gamma }+c,$$(36)

with $\widehat{\gamma }$ the measured shear field, γ the true one, m the multiplicative and c the additive shear bias parameters (we do not consider the latter in the present analysis, as we assume it will be corrected in the shear data processing pipeline). The multiplicative shear bias can come from astrophysical or instrumental systemat-ics (such as the effect of the point spread function – PSF), which affect the measurement of galaxy shapes. We take the m_i parameters (one for each redshift bin) as constant and with a fiducial value of 0 in all bins. To include this further nuisance parameter, one just has to update the different angular PS as $$\{\begin{array}{l}{C}_{ij}^{\text{LL}}(\ell )\to \left(1+m_i\right)\left(1+m_j\right)C_{ij}^{\text{LL}}(\ell )\hfill \\ C_{ij}^{\text{GL}}(\ell )\to \left(1+m_j\right)C_{ij}^{\text{GL}}(\ell )\hfill \\ C_{ij}^{\text{GG}}(\ell )\to C_{ij}^{\text{GG}}(\ell ),\hfill \end{array}$$(37)

where m_i is the i-th bin multiplicative bias, and the GCph spectrum is unchanged since it does not include any shear term. We then have $${\theta }_{\text{shear bias }}=\left\{m_1,m_2,\dots ,m_{{\mathcal{N}}_{\text{b}}}\right\}\text{, }$$(38)

with fiducial values $${\theta }_{\text{shear bias }}^{\text{fid }}=\{0,0,\dots ,0\}.$$(39)

Finally, we introduce the ∆z_i parameters to allow for uncertainties over the first moments of the photometric redshift distribution (argued to have the largest impact on the final constraints in Reischke 2024). We then have (Troxel et al. 2018; Abbott et al. 2018; Tutusaus et al. 2020): $$n_i(ɀ)\to n_i\left(ɀ-\Delta {ɀ}_i\right),$$(40)

which adds new entries to our nuisance parameter vector: $${\theta }_{\text{photo-𝔃 }}=\left\{\Delta {ɀ}_1,\Delta {ɀ}_2,\dots ,\Delta {ɀ}_{{\mathcal{N}}_{\text{b}}}\right\},$$(41)

with fiducial values: $${\theta }_{\text{photo-𝔃 }}^{\text{fid }}=\{0,0,\dots ,0\}.$$(42)

These nuisance parameters – unless specified otherwise – are varied in the Fisher analysis so that the final parameters vector is $$\theta ={\theta }_{\text{cosmo }}\cup {\theta }_{\text{IA}}\cup {\theta }_{\text{gal}.\text{ bias }}\cup {\theta }_{\text{shear bias }}\cup {\theta }_{\text{photo- }ɀ},$$

and $${\theta }^{\text{fid }}={\theta }_{\text{cosmo }}^{\text{fid }}\cup {\theta }_{\text{IA}}^{\text{fid }}\cup {\theta }_{\text{gal}.\text{ bias }}^{\text{fid }}\cup {\theta }_{\text{shear bias }}^{\text{fid }}\cup {\theta }_{\text{photo-𝔃 }}^{\text{fid }},$$

both composed of 𝒩_p = 7 + 3 + 3𝒩_b = 3𝒩_b + 10 elements.

3.6.1 Higher-order bias

To compute the galaxy–galaxy and galaxy–galaxy lensing probe response terms (Eqs. (14) and (15)) we need the second-order galaxy bias b₍₂₎(z). To do this, we follow Appendix C of LG19, in which this is estimated following the halo model⁹ as (Voivodic & Barreira 2021; Barreira et al. 2021) $$b_{(i)}(ɀ)={\displaystyle \int \text{d}}M{\Phi }_{\text{MF}}(M,ɀ)b_{(i)}^{\text{h}}(M,ɀ)\langle N\mid M\rangle /n_{\text{gal}}(ɀ),$$(43)

with $$n_{\text{gal}}(ɀ)={\displaystyle \int \text{d}}M{\Phi }_{\text{MF}}(M,ɀ)\langle N\mid M\rangle ,$$(44)

the galaxy number density, Φ_MF(M, z) the halo mass function (HMF), $b_{(i)}^{\text{h}}$(M, z) the i-th order halo bias, and 〈N|M〉 the average number of galaxies hosted by a halo of mass M at redshift z (given by the halo occupation distribution, HOD). These are integrated over the mass range log M ∈ [9,16], with the mass expressed in units of solar masses (we don't include h in our units). The expression for the i-th order galaxy bias (Eq. (43)) is the same as Eq. (C.2) of LG19, but here we are neglecting the scale dependence of the bias evaluating it at k = 0 so that u(k | M = 0, z) = 1, u(k | M, z) being the Fourier Transform of the halo profile. Strictly speaking, this gives us the large-scale bias, but it is easy to check that the dependence on k is negligible over the range of interest.

Although Eq. (43) allows the computation of both the first and second-order galaxy bias, we prefer to use the values of b₍₁₎(z) measured from the Flagship simulation for the selected galaxy sample; this is to maintain consistency with the choices presented at the beginning Sect. 3.6. For each redshift bin, we vary (some of) the HOD parameters to fit the measured b₍₁₎(z), thus getting a model for $b_{(1)}^{\text{h}}\left(z\right)$. We then compute $b_{(2)}^{\text{h}}\left(z\right)$ using as an additional ingredient the following relation between the first and second-order halo bias, which approximates the results from separate universe simulations (Lazeyras et al. 2016) within the fitting range $1\lesssim b_{(1)}^{\text{h}}\lesssim 10$: $$\begin{array}{ll}{b}_{(2)}^{\text{h}}(M,ɀ)\hfill & =0.412-2.143b_{(1)}^{\text{h}}(M,ɀ)\hfill \\ & +0.929{\left[b_{(1)}^{\text{h}}(M,ɀ)\right]}^2+0.008{\left[b_{(1)}^{\text{h}}(M,ɀ)\right]}^3.\hfill \end{array}$$(45)

Finally, we plug the $b_{(2)}^{\text{h}}$ values obtained in this way back into Eq. (43) to get the second-order galaxy bias. The details of the HMF and HOD used and of the fitting procedure are given in Appendix B.

3.7 Data vectors and Fisher matrix

Up to now, we have outlined a fully general approach, without making any assumptions about the data. We now need to set data-related quantities.

First, we assume that we will measure $C_{ij}^{\text{AB}}(\ell )$ in ten equally populated redshift bins over the redshift range (0.001,2.5). When integrating Eq. (1) in dz, z_max must be larger than the upper limit of the last redshift bin to account for the broadening of the bin redshift distribution due to photo-z uncertainties. We have found that the $C_{ij}^{\text{AB}}(\ell )$ stop varying for z_max ≥ 4, which is what we take as the upper limit in the integrals over z. This also means that we need to extrapolate the bias beyond the upper limit of the last redshift bin; we then take its value as constant and equal to the one in the last redshift bin, that is, b(z > 2.501) = b₁₀.

Second, we assume the same multipole limits as in EC20, and therefore examine two scenarios, as follows:

• pessimistic: $$\left({\ell }_{\min },{\ell }_{\max }\right)=\{\begin{array}{ll}(10,1500)\hfill & \text{ for WL }\hfill \\ (10,750)\hfill & \text{ for GCph and XC }\hfill \end{array}\text{, }$$

• optimistic: $$\left({\ell }_{\min },{\ell }_{\max }\right)=\{\begin{array}{cc}(10,5000)\hfill & \text{ for WL }\hfill \\ (10,3000)\hfill & \text{ for GCph and XC }\hfill \end{array}.$$

Then, for the multipole binning, instead of dividing these ranges into N_ℓ (logarithmically equispaced) bins in all cases as is done in EC20, we follow the most recent prescriptions of the EC and proceed as follows:

• we fix the centres and edges of 32 bins (as opposed to 30) in the ℓ range [10,5000] following the procedure described in Appendix C. This will be the ℓ configuration of the optimistic WL case.

• The bins for the cases with ℓ_max < 5000, such as WL pessimistic, GCph, or XC, are obtained by cutting the bins of the optimistic WL case with ℓ_centre > ℓ_max. This means that instead of fixing the number of bins and having different bin centres and edges as done in EC20, we fix the bins’ centres and edges and use a different number of bins, resulting in, for example, ${\mathcal{N}}_{\ell }^{\text{WL}}>{\mathcal{N}}_{\ell }^{\text{GCph}}$.

The number of multipole bins is then ${\mathcal{N}}_{\ell }^{\text{WL}}=26$ and ${\mathcal{N}}_{\ell }^{\text{GCph}}={\mathcal{N}}_{\ell }^{\text{XC}}=22$ in the pessimistic case and ${\mathcal{N}}_{\ell }^{\text{WL}}=32$ and ${\mathcal{N}}_{\ell }^{\text{GCph}}={\mathcal{N}}_{\ell }^{\text{XC}}=29$ in the optimistic case. In all these cases, the angular PS are computed at the centre of the ℓ bin, as done in EC20.

We note that, because of the width of the galaxy – and, more importantly, lensing – kernels, a given fixed ℓ_max will not correspond to a unique k_max value. A more accurate approach could be, for example, to use the k-cut method presented in Taylor et al. (2018a), which leverages the BNT (Bernardeau-Nishimichi-Taruya) transform (Bernardeau et al. 2014) to make the lensing kernels separable in z, hence allowing for a cleaner separation of scales. We leave the investigation of this important open issue to dedicated work.

As mentioned, we consider the different probes in isolation, as well as combine them in the '3×2pt' analysis, which includes three 2-point angular correlation functions (in harmonic space): $C_{ij}^{\text{LL}}(\ell ),C_{ij}^{\text{GL}}(\ell )$ and $C_{ij}^{\text{GG}}(\ell )$. The ℓ binning for the 3×2pt case is the same as for the GCph one.

The covariance matrix and the derivatives of the data vector w.r.t. the model parameters are the only elements needed to compute the FM elements. The one-dimensional data vector C is constructed by simply compressing the redshift and multipole indices (and, in the 3×2pt case, the probe indices) into a single one, which we call p (or q). For Gaussian-distributed data with a parameter-independent covariance, the FM is given by: $$F_{\alpha \beta }=\frac{\partial C}{\partial {\theta }_{\alpha }}{\mathrm{Cov}}^{-1}\frac{\partial C}{\partial {\theta }_{\beta }}={\displaystyle \sum _{pq}\frac{\partial C_p}{\partial {\theta }_{\alpha }}}{\mathrm{Cov}}_{pq}^{-1}\frac{\partial C_q}{\partial {\theta }_{\beta }}.$$(46)

We refer the reader to EC20 for details on the convergence and stability of the Fisher matrix and derivatives computations.

We note that the size of the 3×2pt covariance matrix quickly becomes large. For a standard setting with 𝒩_b = 10 redshift bins, there are respectively (55,100, 55) independent redshift bin pairs for (WL, XC, GCph), to be multiplied by the different 𝒩_ℓ. In general, Cov will be a 𝒩_C × 𝒩_C matrix with $$\begin{array}{l}{\mathcal{N}}_C=\left[{\mathcal{N}}_{\text{b}}\left({\mathcal{N}}_{\text{b}}+1\right)/2\right]\left[{\mathcal{N}}_{\ell }^{\text{WL}}+{\mathcal{N}}_{\ell }^{\text{GCph}}\right]+{\mathcal{N}}_{\text{b}}^2{\mathcal{N}}_{\ell }^{\text{XC}}\\ =\left[{\mathcal{N}}_{\text{b}}\left({\mathcal{N}}_{\text{b}}+1\right)+{\mathcal{N}}_{\text{b}}^2\right]{\mathcal{N}}_{\ell }^{3\times 2\text{pt}},\end{array}$$(47)

where the second line is for the 3×2pt case, which has the same number of ℓ bins for all probes, and $${\mathcal{N}}_C=\left[{\mathcal{N}}_{\text{b}}\left({\mathcal{N}}_{\text{b}}+1\right)/2\right]{\mathcal{N}}_{\ell }^{\text{WL}/\text{GCph}},$$(48)

for the WL and GCph cases. As an example, we have ${\mathcal{N}}_C^{3\times 2\text{pt,opt}}=6090$.

Being diagonal in ℓ, most elements of this matrix will be null in the Gaussian case. As shown in Fig. 3, this is no longer true with the inclusion of the SSC contribution, which makes the matrix computation much more resource-intensive. The use of the Numba JIT compiler¹⁰ can dramatically reduce the CPU (for Central Processing Unit) time from about 260 s to about 2.5 s for the Gaussian + SSC 3×2pt covariance matrix (the largest under study) on a normal laptop working in single-core mode.

Given the highly non-diagonal nature of the Gaussian + SSC covariance, we might wonder whether the inversion of this matrix (which is needed to obtain the FM, see Eq. (46)) is stable. To investigate this, we compute the condition number of the covariance, which is defined as the ratio between its largest and smallest eigenvalues and in this case of order 10¹³. This condition number, multiplied by the standard numpy float64 resolution (2.22 × 10⁻¹⁶), gives us the minimum precision that we have on the inversion of the matrix, of about 10⁻³. This means that numerical noise in the matrix inversion can cause, at most, errors of order 10⁻³ on the inverse matrix. Hence, we consider the inversion to be stable for the purpose of this work.

Fig. 3 Correlation matrix in log scale for all the statistics of the 3×2pt data-vector in the G and GS cases. The positive and negative elements are shown in red and blue, respectively. The Gaussian covariance is block diagonal (i.e. it is diagonal in the multipole indices, but not in the redshift ones; the different diagonals appearing in the plot correspond to the different redshift pair indices, for ℓ1 = ℓ2). The overlap in the WL kernels makes the WL block in the Gaussian + SSC covariance matrix much more dense than the GCph one.

4 Forecast code validation

In order to validate the SSC computation with PySSC, we compare the 1σ forecast uncertainties (which correspond to a 68.3% probability, due to the assumptions of the FM analysis) obtained using two different codes independently developed by two groups, which we call A and B. To produce the FM and the elements needed for its computation (the observables, their derivatives and the covariance matrix), group A uses a private¹¹ code fully written in Python and group B uses CosmoSIS¹² (Jennings et al. 2016). As stated in the introduction, the only shared feature of the two pipelines is the use of PySSC (to compute the S_ijkl matrix). For this reason, and because the SSC is not considered in isolation but added to the Gaussian covariance, we compare the forecast results of the two groups both for the Gaussian and Gaussian + SSC cases.

Following EC20, we consider the results to be in agreement if the discrepancy of each group’s results with respect to the median – which in our case equals the mean – is smaller than 10%. This simply means that the A and B pipelines’ outputs are considered validated against each other if $$\left|\frac{{\sigma }_{\alpha }^i}{{\sigma }_{\alpha }^m}-1\right|<0.1\text{\hspace{1em} for \hspace{1em}}i=\text{A}, \text{B};{\sigma }_{\alpha }^m=\frac{{\sigma }_{\alpha }^A+{\sigma }_{\alpha }^B}{2},$$(49)

with ${\sigma }_{\alpha }^A$ the 1σ uncertainty on the parameter α for group A. The above discrepancies are equal and opposite in sign for A and B.

The marginalised uncertainties are extracted from the FM F_αβ, which is the inverse of the covariance matrix C_αβ of the parameters: (F⁻¹)_αβ = C_αβ. The unmarginalised, or conditional, uncertainties are instead given by ${\sigma }_{\alpha }^{\text{unmarg}.}=\sqrt{1/F_{\alpha \alpha }}$. We then have $${\sigma }_{\alpha }={\sigma }_{\alpha }^{\text{marg}.}=\sqrt{{\left(F^{-1}\right)}_{\alpha \alpha }}.$$(50)

The uncertainties found in the FM formalism constitute lower bounds, or optimistic estimates, on the actual parameters’ uncertainties, as stated by the Cramér-Rao inequality.

In the following, we normalise σ_α by the fiducial value of the parameter θ_α, in order to work with relative uncertainties: ${\overline{\sigma }}_{\alpha }^i={\sigma }_{\alpha }^i/{\theta }_{\alpha }^{\text{fid}};{\overline{\sigma }}_{\alpha }^m={\sigma }_{\alpha }^m/{\theta }_{\alpha }^{\text{fid}}$, again with i = A,B. If a given parameter has a fiducial value of 0, such as w_a, we simply take the absolute uncertainty. The different cases under examination are dubbed 'G', or 'Gaussian', and 'GS', or 'Gaussian + SSC’. The computation of the parameters constraints differs between these two cases only by the covariance matrix used in Eq. (46) to compute the FM: $$\mathrm{Cov}=\{\begin{array}{ll}{\mathrm{Cov}}_{\text{G}}\hfill & \text{ Gaussian }\hfill \\ {\mathrm{Cov}}_{\text{GS}}={\mathrm{Cov}}_{\text{G}}+{\mathrm{CovsSC}}_{\text{SSC}}\hfill & \text{ Gaussian }+\text{SSC}.\hfill \end{array}$$(51)

As mentioned above, we repeat the analysis for both Euclid’s photometric probes taken individually, WL and GCph, as well as for the combination of WL, GCph, and their cross-correlation XC, the 3×2pt.

For the reader wanting to validate their own code, we describe the validation process in Appendix A. Here we sketch the results of the code validation: in Fig. 4, we show the percent discrepancy as defined in Eq. (49) for the 3×2pt case. Similar results have been obtained for the GCph and WL cases, both for the optimistic and pessimistic settings specified in Sect. 3.7. The constraints are all found to satisfy the required agreement level (less than 10% discrepancy with respect to the mean). In light of these results, we consider the two forecasting pipelines validated against each other. All the results presented in this paper are the ones produced by group A.

Fig. 4 Percentage discrepancy of the normalised 1σ uncertainties with respect to the mean for the WL probe, both in the G and GS cases (optimistic settings). The index i = A, B indicates the two pipelines, whilst α indexes the cosmological parameter. The desired agreement level is reached in all cases (WL, GCph probes and pessimistic case not shown).

5 SSC impact on forecasts

We investigate here how the inclusion of SSC degrades the constraints with respect to the Gaussian case. To this end, we look in the following at the quantity $$ℛ(\theta )={\sigma }_{\text{GS}}(\theta )/{\sigma }_{\text{G}}(\theta ),$$(52)

where σ_G(θ) and σ_GS(θ) are the usual marginalised uncertainties on the parameter θ computed, as detailed above, with Gaussian or Gaussian + SSC covariance matrix. We run θ over the set of cosmological parameters listed in Eq. (27), that is, θ ∈ {Ω_m,0, Ω_b,0, w₀, w_a, h, n_s, σ₈}.

In addition, we examine the FoM as defined in Albrecht et al. (2006), a useful way to quantify the joint uncertainty on several parameters. In this work, we parameterise the FoM following EC2O: for two given parameters θ₁ and θ₂, we have $${\text{FoM}}_{{\theta }_1{\theta }_2}=\sqrt{\det \left({\tilde{F}}_{{\theta }_1{\theta }_2}\right)}.$$(53)

This quantity is inversely proportional to the area of the 2σ confidence ellipse in the plane spanned by the parameters (θ₁, θ₂). ${\tilde{F}}_{{\theta }_1{\theta }_2}$ is the Fisher submatrix obtained by marginalising over all the parameters but θ₁ and θ₂, and is computed by inverting F_αβ (that is, taking the parameters’ covariance matrix), removing all the rows and columns but the ones corresponding to θ₁ and θ₂ and reinverting the resulting 2 × 2 matrix. In the following, we mainly focus on the joint uncertainty on the dark energy equation of state parameters w₀ and w_a, and unless specified otherwise we use the notation FoM = FoM_w0wa. However, the FoM can help quantify the joint uncertainty on different sets of parameters, such as (Ω_m,0 − S₈), with S₈ = σ₈(Ω_m,0/0.3)^{0. 5} (see e.g. Abbott et al. 2022).

We also use the notation ℛ(FoM) as a shorthand for FoM_GS/FoM_G. We note that, since we expect the uncertainties to be larger for the GS case, we have ℛ(θ) > 1, and the FoM being inversely proportional to the area of the uncertainty ellipse, ℛ(FoM) < 1.

5.1 Reference scenario

Let us start by considering the case with 𝒩_b = 10 equipopulated redshift bins. To isolate the impact of SSC and gain better physical insight, we begin by computing the conditional uncertainties as described in the last section.

Table 1 gives the values of the ℛ ratios for the different parameters in the optimistic scenarios, for the single probes and their combination.

In accordance with previous results in the literature, we find that the WL constraints are dramatically impacted by the inclusion of SSC: as found in Barreira et al. (2018a) (cf. their Fig. 2), all cosmological parameters are affected, with w_a and h being impacted the least and (Ω_m,0, σ₈, log₁₀(T_AGN/K)) the most. This is because the SSC effect is essentially an unknown shift, or perturbation, in the background density within the survey volume, and is hence degenerate with the parameters which more closely relate to the amplitude of the signal. Being this a nonlinear effect, its impact is also tied to the amount of power on small scales, which in turn is influenced by the baryonic boost, parameterised by log₁₀(T_AGN/K) in our model.

The results in Table 1 also show that GCph is not as strongly affected by SSC – with the exception of the log₁₀(T_AGN/K) constraint. This is an expected result (see e.g. Bayer et al. 2023), mainly driven by the fact that the GCph probe response coefficients are lower (in absolute value) than the WL ones, as can be seen in Fig. 1. This is due to the additional terms that account for the response of the galaxy number density n_g (see Eq. (14)), which is itself affected by the super-survey modes. Additionally, as can be seen from Fig. 2, all WL kernels have non-zero values for z → 0, contrary to the GCph ones. In this limit, the effective volume probed by the survey tends to 0, hence making the variance of the background modes σ² tend to infinity. We thus have a larger S_ijkl matrix, which is one of the main factors driving the amplitude of the SSC. We also note that the importance of bary-onic feedback for GCph depends on the galaxy bias model used; it will likely be reduced when using a non-linear bias expansion, needed for an accurate analysis on small scales (Desjacques et al. 2018). Both topics are still active areas of research, and we leave the characterisation of the SSC impact with the inclusion of higher-order bias terms for future work.

The impact for the full 3×2pt case sits in principle in between the two extremes as a consequence of the data vector containing the strongly affected WL probe, and the less affected GCph one. However, this is clearly apparent only in the case of Ω_m,0, since it is the only cosmological parameter for which the WL constraining power is higher than the GCph one; for the other parameters, the trend resembles very closely the one found for GCph, because of its dominant contribution to the 3×2pt precision.

Lastly, the contribution from the XC probe is again an intermediate case, as can be anticipated by looking at its response coefficient in Fig. 1, so the final impact on the FM elements will be intermediate between the WL and GCph cases, as the ℛ(θ) values in Table 1 indeed show.

Having explored the impact of SSC on the conditional uncertainties, we move on to analyse a more realistic scenario, namely letting the parameters in the analysis free to vary, as opposed to fixing them to their fiducial values.

To compute the 1σ uncertainties in this case, we marginalise over all cosmological and nuisance parameters. We add Gaussian priors of standard deviation σ^p = 5 × 10⁻⁴ on the multiplicative shear bias parameters, and of σ^p = σ_z(1 + z) with σ_z = 0.002 (see Euclid Collaboration 2024 and references therein) on the dz_i parameters. We also include a Gaussian prior of σ^p = 0.06 on log₁₀(T_AGN/K), which roughly matches the prior range recommended in Mead et al. (2021). To add these priors in the FM analysis, it is sufficient to add ${\left({\sigma }_{\alpha }^{\text{p}}\right)}^{-2}$ to the appropriate diagonal elements of the G and GS FMs (${\sigma }_{\alpha }^{\text{p}}$ being the value of the prior on parameter α).

The results of this analysis, which we take as our reference scenario, are shown in Table 2 and in Figs. 5–8. For the WL case, marginalisation over all cosmological and nuisance parameters leads to the result that the SSC has now a very minor impact – of a maximum of about 13% in the optimistic case -on the constraints and the FoM: again, this is in line with what found in Barreira et al. (2018a). We also compute the Ω_m,0 − S₈ FoM, after projecting the FM to the new parameter space with S₈ replacing σ₈ (see e.g. Coe 2009). This drops from 336 to 292 in the optimistic case, corresponding to a ratio of 0.87 (0.92 in the pessimistic case), closely mirroring the results found for FoM_w0wa.

These values are actually easily explained: marginalising over cosmological and nuisance parameters, particularly if these are degenerate with the amplitude of the signal, dilutes the SSC effect in a larger error budget; because of this, it is the relative rather than the absolute impact of SSC that decreases. Indeed, marginalising over additional parameters is formally equivalent to having additional covariance. The parameters that mostly change the SSC impact when marginalised over are the cosmo-logical ones, but the multiplicative nuisance parameters m_i and b_i also play a role: adding these to the set of free parameters introduces a degeneracy between these and the overall amplitude of $C_{ij}^{\text{AB}}(\ell )$. Such a degeneracy is a mathematical one present on the whole ℓ range. As a consequence, the constraints on all the parameters and the FoM are degraded in a way that is independent of the presence of SSC. This is shown in Figs. 5 and 6, which respectively exhibit the relative uncertainty $\overline{\sigma }$ and the FoMs in the G and GS cases for each parameter, if we marginalise or not over nuisance parameters. Letting these free to vary, i.e. marginalising over them, tends to increase the uncertainty on cosmological parameters way more than including SSC, and this is even more true when these nuisance parameters are simply multiplicative such as b_i and m_i.

This of course does not mean that varying more parameters improves the constraints. Indeed, the uncertainties on all parameters increase (hence the FoM decreases) with respect to the case of conditional uncertainties introduced above.

The same reasoning applies to the GCph probe, for which the SSC impact drops now to a subpercent level in the pessimistic case and to about 10% at most in the optimistic case (on w₀ and the two FoMs). As mentioned above, for GCph, one of the reasons behind the observed decrease in the (already low) SSC relative impact is the marginalisation over the galaxy bias parameters, which are perfectly degenerate with the amplitude of the signal and over which we impose no prior.

On the other hand, the results for the 3×2pt case show that the SSC still matters. The additional information carried by the GCph and XC data allows the partial breaking of parameter degeneracies, including those with probe-specific systematics such as m_i and b_i, hence making the scale-dependent increase of the uncertainties due to the inclusion of SSC important again. For this reason, in this case, the 3×2pt does not follow the behaviour of the single probes as seen earlier. In particular, the dark energy FoM, whose increase with respect to present surveys is one of the main objectives of the Euclid mission, is highly degraded – by a factor of about 2 in the optimistic case. This is mainly due to the large impact on the dark energy equation of state parameters, showing the importance of accounting for SSC in upcoming LSS analyses. The same conclusion holds for ${\text{FoM}}_{{\text{Ω}}_{\text{m},0}{S}_8}$, with a slightly larger ℛ value driven by the lower impact on σ₈ w.r.t. w₀ and w_a. In Fig. 7 we show the comparison of the 2D contours for all cosmological parameters between G and GS in the case of the 3×2pt analysis, in the optimistic case; the most impacted parameters in this case are Ω_m,0, w₀ and w_a. In addition, this shows that SSC does not seem to strongly affect the correlations between cosmological parameters.

Table 2 also shows the ratios in the 'GCph diag.' case, in which we neglect the GCph cross-reshift elements (i.e. we take the diagonal of the $C_{ij}^{\text{GG}}(\ell )$ matrix for each ℓ value). This choice is sometimes made because of the large sensitivity of such measurements on the photometric redshift calibration, which is less of an issue for WL due to its broad kernel. The values show the robustness of our conclusions even in this case, with the ℛ values being very close to the standard case both for GCph and for the 3×2pt.

We note that these conclusions depend on the scale cuts imposed in the data vector, as can be seen from the lower FoM decrease (and of the SSC impact as a whole, with some exceptions) in the pessimistic case, in line with what found in Lacasa (2020) for GCph. This is a direct consequence of the larger amount of non-linear modes included in the data vectors, which as mentioned previously are more subject to mode coupling and hence contribute more to SSC. For WL, it should be noted that the diagonal elements of the total covariance matrix are always dominated by the Gaussian contribution, because of the presence of the scale-independent shape noise (see Eq. (24) for A = B = L), which largely dominates over the SSC on small scales. This is consistent with the results of Upham et al. (2022) showing that the diagonal elements of the WL total covariance matrix are more and more dominated by the Gaussian term as we move to higher ℓ. This is also the case for GCph, although the predominance of the Gaussian term along the diagonal is less pronounced because of the smaller contribution of shot noise.

As mentioned above, more sophisticated choices of scale cuts, e.g. through the use of the BNT transform, will allow decreasing the high-k contribution to a given ℓ mode, hence mitigating the SSC impact more than the "hard" angular scale cut considered here (and in EC20). We leave the investigation of this point to a forthcoming publication.

To conclude this section, it is also worth looking at the impact of SSC on the astrophysical nuisance parameters. Indeed, although an issue to be marginalised over when looking at cos-mological ones, the IA and the galaxy bias parameters are of astrophysical interest. We show the impact of SSC on the constraints on these quantities in Fig. 8, as well as on the galaxy bias and the multiplicative shear bias parameters.

For IA-related nuisance parameters, the uncertainty increase due to SSC is lower than 5% when considering WL-only, and below 1% for the full case. The multiplicative shear bias parameters are affected by a similar amount (up to around 10% for the last bin) and in a similar way for WL and the 3×2pt. As for the galaxy bias parameters, the impact is modest for GCph but quite significant (between 30 and 40%) for the 3×2pt. These results are analogous to what was found before for the cosmological parameters: in the marginalised 3×2pt case the impact of SSC is more apparent, since many degeneracies are broken thanks to the probes combination and cross-correlation.

In this case, contrary to the IA parameters, the uncertainty on b_i and m_i in each of the ten redshift bins is significantly affected by SSC. This is because both of these nuisance parameters simply act as a multiplicative factor on the power spectrum and are thus highly degenerated with the effect of SSC. Again, this is related to the fact that the first-order effect of SSC is to modulate the overall clustering amplitude because of a shift in the background density δ_b. For the m_i parameters, we remind that a tight (while realistic) prior is imposed, keeping the SSC uncertainty increase, and hence the ℛ ratio, quite low even for the full 3 × 2pt. We note that going beyond the linear approximation for the modelling of the galaxy bias will add more nuisance parameters, thus degrading the overall constraints on cosmological parameters and further reducing the relative impact of SSC.

Fig. 5 Marginalised and conditional optimistic 1σ uncertainties on the cosmological parameters, relative to their corresponding fiducial values (in percent units), in both the G and GS cases for WL, GCph, and the 3×2pt. We highlight the logarithmic scale on the x-axis, made necessary by the large range of values. The values in these plots have been used to compute the ratios in Tables 1 and 2.

Fig. 6 Same as Fig. 5, but for the dark energy and Ωm,0 − S8 FoM.

Fig. 7 Contour plot for the G and GS constraints, considering the full 3×2pt analysis in the optimistic case, in the reference scenario. The two shaded regions of the ellipses represent the 1 and 2σ contours. For clarity, the nuisance parameters are shown separately in Fig. 8 and the log10(TAGN/K) parameter name has been shortened to log(T).

Fig. 8 Percent increase of the marginalised 1σ uncertainty of the nuisance parameters, for all probe choices, in the optimistic case and for the reference scenario.

5.2 Non-flat cosmologies

In the previous section, we investigate the SSC on the cosmological parameters under the assumption of a flat model. Actually, the requirement on the FoM assessed in the Euclid Red Book (Laureijs et al. 2011) refers to the case with the curvature as an additional free parameter to be constrained, that is, the non-flat w₀w_aCDM model. This is why in EC20 are also reported the marginalised uncertainties for the parameter Ω_DE,0, with a fiducial value ${\text{Ω}}_{\text{DE,0}}^{\text{fid}}=1-{\text{Ω}}_{\text{m},0}$ to be consistent with a flat universe. It is then worth wondering what the impact of SSC is in this case too. This is summarised in Table 3, where we now also include the impact on Ω_DE,0.

A comparison with the results in Table 2 is quite hard if we look at the single parameters. Indeed, opening up the parameter space by removing the flatness assumption introduces additional degeneracy among the parameters controlling the background expansion, which are thus less constrained whether SSC is included or not. We can nevertheless note again that, for the marginalised uncertainties, 3×2pt is still the most impacted probe; the difference between pessimistic and optimistic scenarios is now less evident with ℛ(θ) increasing or decreasing depending on the parameter and probe considered.

Once more, the most affected parameters for WL are (Ω_m,0, σ₈), the uncertainties on which are now further degraded by the fact that they correlate with the parameter Ω_DE,0 which is also affected. Although (w₀, w_a) are also degraded by the SSC, a sort of compensation is at work, so that the overall decrease in the FoM is similar to the case with the flatness prior. The motivations that make GCph much less affected still hold when dropping the flatness prior, explaining the corresponding ℛ(θ) values.

We also note a slight increase of ℛ(FoM) in the 3×2pt optimistic case, meaning a smaller degradation of the FoM due to SSC. The dark energy FoM indeed degrades by 52% (45%) in the non-flat case vs. 50% (58%) for the flat case in the pessimistic (optimistic) scenario, while $Fo{M}_{{\text{Ω}}_{\text{m},0}{S}_8}$ degrades by 47% (35%) in the non-flat case vs. 43% (41%) for the flat case in the pessimistic (optimistic) scenario. This can be qualitatively explained by noting that the decrease of both FoM(G) and FoM(GS) is related to a geometrical degeneracy which is the same on all scales, whether or not they are affected by the increase in uncertainty due to the SSC inclusion.

Overall, these results suggest a dependence of the SSC significance on both the number and type of parameters to be constrained. Qualitatively, we can argue that SSC is more or less important depending on whether the additional parameters (with respect to the reference case of a flat model) introduce degeneracies which are or not scale-dependent and how strong is the degeneracy between these parameters and the amplitude of the power spectrum. To give an example, in future works lens magnification effects should be included in the analysis as these were shown to have a significant impact on cosmological constraints (Unruh et al. 2020). From our results, we can anticipate that the inclusion of magnification-related nuisance parameters will further dilute the impact of SSC.

Fig. 9 Ratio between WL and 3×2pt marginalised uncertainties computed by including or neglecting the SSC contribution, as a function of the number of (equally populated) redshift bins, for the pessimistic and optimistic cases.

5.3 Dependence on redshift binning

The results summarised in Tables 1–3 were obtained for a fixed choice of number and type of redshift bins. We investigate here how the results depend on these settings given that we expect both the G and GS constraints to change as we vary the number and type of bins. We consider again the reference scenario introduced in Sect. 5.1, that is, the case of flat models, marginalising over the full set of nuisance parameters while imposing Gaussian priors over some of them. In this section we only consider the WL and 3×2pt cases, since as we have seen SSC has a modest impact on GCph.

Let us first consider changing the number of redshift bins 𝒩_b. We show the scaling of ℛ(θ) as a function of 𝒩_b for the WL and 3×2pt probes, respectively, in Fig. 9 – for both the pessimistic and optimistic assumptions. The most remarkable result is the weak dependence of ℛ(x) w.r.t. 𝒩_b, as can be inferred from the small y range spanned by the curves – see e.g. the bottom right panel, for the FoM. More specifically, the scaling of ℛ(θ) with 𝒩_b depends on the parameter and the probe one is looking at. It is quite hard to explain the observed trends because of the interplay of different contrasting effects. For instance, a larger number of bins implies a smaller number density in each bin, and hence a larger shot noise. As a consequence, the SSC contribution to the total covariance for the diagonal elements will likely be more and more dominated by the Gaussian component because of the larger shot and shape noise terms. However, this effect also depends on the scale so that, should the SSC be the dominant component on the scales to which a parameter is most sensitive, the impact should still be important. On the other hand, a larger number of bins also comes with a larger number of nuisance parameters which, as shown above, leads to a reduction of the SSC impact. Quantifying which actor plays the major role is hard which explains the variety of trends in the different panels.

As a further modification to the reference settings, we can change how the redshift bins are defined. We have up to now considered equipopulated (EP) bins so that the central bins cover a smaller range in z, because of the larger source number density. As an alternative, we divide the full redshift range into 𝒩_b bins with equal redshift support (‘equidistant’, ED), and recompute the FM forecasts with and without SSC. We show the FoM ratio as a function of the number of bins for EP and ED bins considering WL (left) and 3×2pt (right) probes in the optimistic scenario in Fig. 10. We note that finding the exact number and type of redshift bins used to maximise the constraining power of Euclid is outside the scope of this paper; this effort is indeed brought forward in the context of the SPV exercise.

In order to qualitatively explain these results, let us first consider the WL case. Given that the bins are no longer equipopulated, the number density of galaxies will typically be larger in the lower redshift bins than in the higher ones. As a consequence, the larger the number of bins, the higher the shape noise in the higher redshift bins so that the SSC will be subdominant in a larger number of bins, which explains why its impact decreases (i.e. ℛ(FoM) increases) with 𝒩_b. Nevertheless, the impact of SSC will be larger than in the EP case since SSC will dominate in the low redshift bins which are the ones with the largest signal-to-noise ratio. This effect is small, and the difference in ℛ(FoM) between ED and EP is no larger than 3–5%.

When adding GCph and XC into the game, the impact of SSC is determined by a combination of contrasting effects. On the one hand, we can repeat the same qualitative argument made for WL also for GCph and XC¹³ thus pointing at ℛ(FoM) increasing with 𝒩_b. The larger the number of bins, the narrower they are, and the smaller the cross-correlation between them hence the smaller the Gaussian covariance. This in turn increases the number of elements in the data vector whose uncertainty is dominated by the SSC. Should this effect dominate, we would observe a decrease of ℛ(FoM) with 𝒩_b with the opposite trend if it is the variation of the shape and shot noise to matter the most. This qualitative argument allows us then to roughly explain the non-monotonic behaviour of ℛ(FoM) we see in the right panel of Fig. 10.

It is worth remarking, however, that the overall change of ℛ(FoM) for subsequent ED (or EP) bins over the range in 𝒩_b is smaller than ~5% which is also the maximum value of the difference between ℛ(FoM) values for EP and ED bins once 𝒩_b is fixed.

The analysis in this section motivates us to argue that the constraints and FoM degradation due to SSC have a weak dependence on the redshift binning scheme.

Fig. 10 FoM ratio vs the number of EP and ED гedshift bins for WL (left) and 3×2pt (right) in the optimistic scenario. The blue curves in this plot match the ones in the bottom right panel of Fig 9 in the optimistic case, for the corresponding probes.

5.4 Requirements on prior information

The results in the previous paragraph show that the SSC may dramatically impact the constraints on the cosmological parameters. As a consequence, the 3×2pt FoM is reduced by up to ~ 50% with respect to the case when only the Gaussian term is included in the total covariance. This decrease in the FoM should actually not be interpreted as a loss of information due to the addition of the SSC. On the contrary, one can qualitatively say that removing SSC from the error budget is the same as adding information that is not actually there. It is nevertheless interesting to ask which additional information must be added to recover the Gaussian FoM, which is usually taken as a reference for gauging the potential of a survey. This information can come from priors on the nuisance (or cosmological) parameters. In the following section, we investigate the former option by adding Gaussian priors on the galaxy and multiplicative shear bias parameters.

To this end, we consider the realistic case of a flat model plus the galaxy bias and multiplicative shear bias as nuisance parameters. As a simplifying assumption, we assume that all the 𝒩_b bias values b_i are known with the same percentage uncertainty ε_b = σ_b,i/b_{fid, i}, while we put the same prior σ_m on all the m_i parameters (having set the fiducial value m_fid to 0). We then compute the FoM with and without SSC for the 3×2pt probe in the optimistic scenario and investigate how the ratio ℛ(FoM) scales with (ε_b, σ_m) obtaining the results shown in Fig. 11.

Both in the G and GS cases, the FoM shows little to no sensitivity to priors above ~1% and no sign of saturation even for extremely narrow priors (as little as 0.01%) because of the presence of such a large number of nuisance parameters. For σ_m we observe a similar behaviour, with the difference between the curves corresponding to the two larger priors values being visibly smaller than the difference between the curves corresponding to the smaller ones. This is less true for the GS FoM (orange curves), as the presence of SSC in the error budget decreases the relative overall improvement coming from a smaller uncertainty over the m_i parameters.

A prior on the nuisance parameters increases both the Gaussian and Gaussian + SSC FoM so that one could expect their ratio to be independent of the prior itself. This is not exactly the case since the correlation between different multipoles introduced by SSC alters the way the prior changes the FM elements. As a result, we find a non-flat scaling of ℛ(FoM) as can be seen from the right panel of Fig. 11. The behaviour of ℛ(FoM) with ε_b tells us that FoM_GS increases with decreasing ε_b slower than FoM_G when the galaxy bias is known with an uncertainty smaller than the percent level. Another way to interpret it is that the information gained in the FoM saturates faster when SSC is included: better constraints on ε_b do not bring more information as the SSC now dominates the error budget. However, it is worth stressing that, even for a strong prior on the multiplicative shear bias, the FoM ratio can actually be improved significantly only under the (likely unrealistic) assumption of a subpercent prior on the galaxy bias.

The need for such strong priors comes from the attempt to retrieve the same FoM as a Gaussian case. Alternatively, one can also wonder which additional information must be added through priors to retrieve the idealised FoM value obtained in forecasts that neglect the SSC. In other words, we look for the requirements that must be put on the priors (ε_b, σ_m) in order to make FoM_GS/FoM_ref = 1, where FoM_ref = 617 is the FoM computed for a flat reference case without SSC and with no priors on galaxy bias, but a fiducial prior σ_m = 5 × 10⁻⁴ on the shear bias. The answer to this question is shown in Fig. 12 for the optimistic scenario and 10 equipopulated redshift bins. Some numbers help to better understand how priors can indeed supply the additional information to retrieve the FoM one would obtain in an ideal case where SSC is absent. Solving $${\text{FoM}}_{\text{GS}}\left({\epsilon }_b,{\sigma }_m\right)=f{\text{FoM}}_{\text{ref}}$$

with respect to ε_b, we get $${\epsilon }_b=\{\begin{array}{ll}(0.05,0.03,0.02)\%\hfill & \text{ for }{\sigma }_m=0.5\times {10}^{-4}\hfill \\ (0.04,0.02,0.02)\%\hfill & \text{ for }{\sigma }_m=5\times {10}^{-4}\hfill \\ (0.03,0.02,0.01)\%\hfill & \text{ for }{\sigma }_m=10\times {10}^{-4},\hfill \end{array}$$

where the three values refer to f = (0.8,0.9,1.0). These numbers (and the contours in Fig. 12) show that it is possible to compensate for the degradation due to SSC only by adding strong priors on the galaxy bias, which have a much larger impact on the (G and GS) FoM than strong priors on the multiplicative shear bias. However, it is worth noticing that it is actually easier to obtain priors on the multiplicative shear bias provided a sufficient number of realistic image simulations are produced and fed to the shear measurement code to test its performance. It is therefore worth wondering how much the FoM is restored by improving the prior on m for a fixed one on the bias. We find $$\frac{{\text{FoM}}_{\text{GS}}}{{\text{FoM}}_{\text{ref}}}=\{\begin{array}{ll}(1.24,1.13,1.08)\hfill & \text{ for }{\epsilon }_b=0.01\%\hfill \\ (0.68,0.62,0.59)\hfill & \text{ for }{\epsilon }_b=0.1\%\hfill \\ (0.5,0.45,0.43)\hfill & \text{ for }{\epsilon }_b=1\%,\hfill \end{array}$$

with the three values referring to σ_m = (0.5,5.0,10) × 10⁻⁴. As expected, improving the prior on the multiplicative bias with respect to the fiducial one (which, we remind, is included in FoM_ref) does not help a lot in recovering the constraining power. A very tight prior of around 0.1% prior on the galaxy bias can recover a significant amount of the reference FoM (almost 70%) thanks to the additional information compensating for the presence of SSC.

Investigating whether the priors proposed here can be achieved in practice (e.g. through theoretical bias models tailored to galaxy clustering data or N-body hydrodynamic simulations) is outside the aim of this work. We refer the interested reader to for example Barreira et al. (2021), Zennaro et al. (2022), and Ivanov et al. (2024) for some preliminary results.

Fig. 11 3×2pt FoM in the optimistic scenario with and without SSC as a function of the percentage prior εb on the galaxy bias parameters for σm = (5, 50, 500) × 10−4 (solid, dashed, dotted lines). Left: FoM values. Right. ℛ(FoM), defined as the ratio between the blue and red curves.

Fig. 12 FoMGS contours in the (εb, σm) plane for FoMGS/FoMref going from 0.8 to 1.1 in steps of 0.05 (from right to left).

6 Conclusions

Precision cosmology requires precision computation: previously neglected theoretical contributions must therefore now be taken into account. Motivated by this consideration, we computed and studied the impact of SSC on the Euclid photometric survey, exploring how the different probes and their combination are affected by this additional, non-Gaussian term in the covariance matrix. The analysis of the impact of SSC on the spectroscopic survey, which has been shown to be small in Wadekar et al. (2020) for the Baryon Oscillation Spectroscopic Survey (BOSS) data, is left for future work. We employed a FM analysis, producing forecasts of the 1σ marginalised uncertainties on the measurement of the cosmological parameters of the flat and non-flat w₀w_aCDM cosmological models. We validated two different forecast pipelines against the results of EC20, taking as reference survey the one specified therein, and then updated the galaxy bias and the source redshift distributions according to the most recent versions presented in Euclid Collaboration (2021). The SSC was computed relying on the analytical approximations and numerical routines presented in LG19, interfacing the public code PySSC with two distinct forecast pipelines to validate the constraints. As a further step forward, we build upon the work of LG19 by computing the scale and redshift dependence of the response functions of the different probes, starting from the results of Wagner et al. (2015b) and Barreira et al. (2018b).

We quantify the severity of the impact with the ratio σ_GS/σ_G between the marginalised or conditional uncertainties with and without SSC; this is found to vary significantly between different parameters and probes, and between the number and type of free nuisance parameters, in agreement with recent results (Upham et al. 2022; Barreira et al. 2018a; Lacasa 2020).

The conditional uncertainties show WL to be dramatically impacted by SSC, with all cosmological parameter uncertainties increasing by up to 210% (for Ω_m,0) in the optimistic case. The GCph constraints are less sensitive to the addition of SSC, showing a smaller broadening of the uncertainties for all parameters, namely by up to one order of magnitude with respect to WL. The 3×2pt case sits in between these two, while being in general closer to the GCph results because of its larger constraining power.

When considering marginalised constraints, the relative impact of SSC decreases significantly. In this case, the most impacted parameters for the single probes are mainly Ω_m,0, w₀, w_a, and σ₈. Furthermore, the 3×2pt becomes by far the most impacted probe, precisely because of its power in breaking parameter degeneracies. Indeed, in the reference case and the optimistic scenario, (FoM_w0wa, $Fo{M}_{{\text{Ω}}_{\text{m,}0}{S}_8}$) decrease by (58%, 41%), hinting at the necessity to include SSC in the upcoming Euclid analysis.

These results are the consequence of a complicated interplay between three factors. First, SSC originates from the uncertainty in the determination of the background mean density when measuring it over a finite region. This prevents determination of the overall amplitude of the matter power spectrum, which increases the uncertainty on those parameters that concur in setting its amplitude, mainly Ω_m,0 and σ₈. Second, the elements of the SSC matrix depend on the amplitude of the response functions. Third, the impact depends on how large a contribution the signal receives from the low-z region, where the effective volume probed is smaller, making the variance of the background modes larger. The last two factors are both more severe for WL than for GCph, causing the former probe to be more affected than the latter.

Finally, the deviation of a given element of the GS FM from the Gaussian one also depends on its correlations: in other words, the degradation of the constraints on a given parameter can be large if this is strongly correlated with a parameter severely degraded by SSC. Quantifying the impact of SSC on a single parameter is therefore quite hard in general, and must be investigated on a case-by-case basis taking care of the details of the probe and the way it depends on the parameter of interest.

Nuisance parameters to be marginalised over act as a sort of additional contribution to the covariance. As such, the importance of both the Gaussian and SSC contribution to the overall effective covariance becomes less important when the number of nuisance parameters increases. In order to consider cases that most closely mimic future Euclid data, we opened up the parameter space by adding Ω_DE,0 (i.e. removing the flatness prior). We find that, as long as the additional parameters have a scale-independent degeneracy with the most impacted ones, the relative impact of SSC decreases. We stress, however, that this reduction in the SSC impact has undesired consequences; the marginalised uncertainties on the parameters are definitely worsened, but the degradation is roughly the same whether the SSC is included or not, hence making the ratio σ_GS/σ_G closer to unity for all parameters and probes. This result can be taken as a warning against investing too much effort in refining the estimate of the computationally expensive SSC when no approximations are made. For a Euclid-like survey, the main concern would indeed be the number of nuisance parameters, which makes the impact of the SSC itself less relevant.

We furthermore note that, in light of the recent theoretical developments presented in Lacasa et al. (2023), it appears feasible to include the effect of SSC in the form of nuisance parameters, which would be the value of the density background δ_b in each redshift bin. This approach is interesting as it would reduce the complexity of the data covariance matrix and would allow a simpler interpretation of the effect of SSC and how it is correlated to the other cosmological and nuisance parameters.

Variations in the z binning strategy have contrasting effects: a larger number of bins means a larger number of nuisance parameters (either galaxy bias or multiplicative shear bias for each bin), which leads to a loss of constraining power. Moreover, the larger the number of bins, the larger the Gaussian contribution to the covariance, making the shot and shape noise dominate over SSC for diagonal elements. On the downside, a larger number of bins leads to larger data vectors, thus adding information that can partially compensate for the increase in the covariance. The contrasting effects at play conspire in such a way that the degradation of the FoM due to SSC is found to be approximately independent of the number of redshift bins (cfr. Fig. 10).

An interesting development in this sense is to leverage the SSC dependence on the low-z contribution to investigate whether or not its impact could be mitigated by the use of the BNT transform, which transforms redshift bins in such a way as to increase the separation between the WL kernels. This will be investigated in a forthcoming work.

An alternative strategy is to increase the constraining power by adding information through informative priors, hence recovering the FoM when SSC is incorrectly neglected. We investigate this possibility by quantifying the requirements on the prior information needed to recover the Gaussian FoM. Our results show that the main role is played here by the priors on galaxy bias parameters, while the FoM recovery is less sensitive to the prior on the multiplicative shear bias. However, the galaxy bias must be known to subpercent level in order to recover ~70% of the Gaussian FoM. Investigating whether or not this is possible is outside the scope of this paper. We nevertheless note that such remarkable prior information is the same as stating we are able to model the evolution of the bias with redshift. This is actually quite difficult based on the current knowledge of galaxy formation processes. Alternatively, one could look for an empirical fitting formula as a compromise between the need for strong priors on bias and the number of nuisance parameters.

One part of the covariance modelling not investigated in this work is the geometry of the survey footprint. While it is true that, for the large sky coverage considered, the full-sky approximation for SSC has been shown to suffice (Gouyou Beauchamps et al. 2022), a more realistic treatment accounting for the survey geometry should be considered for the Gaussian term. Still, we expect the main conclusions of this study to hold, as the mode coupling caused by the convolution with the survey mask, generating off-diagonal elements also in the Gaussian covariance, will mainly affect large scales where SSC is subdominant. It is however important to note that small holes in the survey mask (e.g. due to the presence of bright stars) generate mode coupling also on small scales, where we have seen the SSC impact to be most prominent. This will be investigated in future works.

Although some more work is needed to improve the robustness of our results, for example by comparing the different approximations presented in the literature, we can conclude that the effect of including the SSC term in the total covariance matrix of Euclid photometric observables is definitely nonnegligible, especially for WL and 3×2pt. However, the degradation of the constraints on cosmological parameters depends on the particular probe and the number and kind of parameters to constrain. The FoM is nevertheless reduced by 52% (45%) for the 3×2pt probe in the pessimistic (optimistic) scenario where all cosmological (including Ω_DE,0) and nuisance (multiplicative shear bias) parameters are left free to vary. Maximising the power of the actual Euclid photometric data by taking into account the presence of SSC is a daunting task, which we will report on in a forthcoming publication.

Acknowledgements

The computational part of the work has been performed using the Python programming language, interfaced with scientific packages like astropy (Astropy Collaboration 2013, 2018) for cosmological calculations, Numba (Lam et al. 2015) for code speedup, NumPy (Harris et al. 2020) for matrix manipulation, SciPy (Virtanen et al. 2020) for numerical integration and Matplotlib (Hunter 2007) for data visualisation. The authors thank the anonymous referees for their helpful comments that improved the quality of the manuscript. DS would like to thank Raphael Kou for the fruitful discussion on the SSC impact on GCph. SGB was supported by CNES, focused on Euclid mission. The project leading to this publication has received funding from Excellence Initiative of Aix-Marseille University -A*MIDEX, a French “Investissements d’Avenir” programme (AMX-19-IET-008-IPhU). SC acknowledges support from the ‘Departments of Excellence 2018-2022’ Grant (L. 232/2016) awarded by the Italian Ministry of University and Research (MUR). IT acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (Grant agreement No. 863929; project title “Testing the law of gravity with novel large-scale structure observables” and acknowledges support from the Spanish Ministry of Science, Innovation and Universities through grant ESP2017-89838, and the H2020 programme of the European Commission through grant 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 Agenzia Spaziale Italiana, the Austrian Forschungsförderungsgesellschaft funded through BMK, the Belgian Science Policy, the Canadian Euclid Consortium, the Deutsches Zentrum für Luft- und Raumfahrt, the DTU Space and the Niels Bohr Institute in Denmark, the French Centre National d’Etudes Spatiales, the Fundação para a Ciência e a Tecnologia, the Hungarian Academy of Sciences, the Ministerio de Ciencia, Innovación y Universidades, the National Aeronautics and Space Administration, the National Astronomical Observatory of Japan, the Netherlandse Onderzoekschool Voor Astronomie, the Norwegian Space Agency, the Research Council of Finland, the Romanian Space Agency, the State Secretariat for Education, Research, and Innovation (SERI) at the Swiss Space Office (SSO), and the United Kingdom Space Agency. A complete and detailed list is available on the Euclid web site (www.euclid-ec.org).

Appendix A Details of the code validation

In the following, we provide an overview of the steps undertaken to compare and validate the codes used in this work, and some of the lessons learnt in the process.

In order to compute and validate the results we adopt the scheme sketched in Fig. A.1, which highlights the dependency of each main element of the forecast computation on the others. In particular, we have that:

1. The lσ constraints are obtained from the FM through Eq. (50), and the FM is built in turn from the (inverse) covariance matrix and the derivatives of the angular PS $C_{ij}^{AB}(\ell )$ as indicated in Eq. (46).

2. The Gaussian covariance depends on the $C_{ij}^{AB}(\ell )$ through Eq. (24) (and the noise PS, Eq. (25)). The SSC also depends on the $C_{ij}^{AB}(\ell )$, with the added contribution of the $R_{ij}^{AB}(\ell )$ terms and the output of the PySSC module, the S _ijkl matrix – following Eq. (10).

3. The $C_{ij}^{AB}(\ell )$ are constructed by convolving the (non-linear) matter PS with the lensing and galaxy weight functions, as in Eq. (1). The S_ijkl matrix also depends on the weight functions (see Eq. (11)), which are in fact the main external input needed by PySSC, and on the linear matter PS through the σ²(z₁,z₂) term (Eq. (5)). It is to be noted, however, that PySSC computes this PS internally, needing only the specification of a dictionary of cosmological parameters with which to call the Boltzmann solver CLASS through the Python wrapper classy. This means that we also have to make sure that the fiducial value of the parameters used to compute the PS of Eq. (1) are the same ones passed to PySSC (this time to compute the linear PS), in order to work with the same cosmology.

While to compute the constraints we follow the scheme from right to left, starting from the basic ingredients to arrive at the final result, the general trend of the validation is the opposite: we begin by comparing the final results, then work our way back whenever we find disagreement.

We then start the comparison from the σ_α. If a discrepancy larger than 10% is found, we check the quantities they depend on, which in this case are the covariance matrices (see Eq. (46)). If these agree, we check the codes directly. If these disagree, we iterate the process by checking the subsequent element in the scheme (in this case the S_ijkl matrix and the $C_{ij}^{AB}(\ell )$, until agreement is found. Essentially, this means that the disagreement in the outputs of the codes at each step can either come from the inputs, or from the codes themselves. Once the cause of the discrepancy is found and fixed, the computation is repeated and the process can start again.

The pipelines under comparison are both written in the Python language. One of them requires as external inputs the weight functions, the angular $C_{ij}^{AB}(\ell )$ and their derivatives with respect to the cosmological parameters; whilst the other produces these through the use of CosmoSIS¹⁴ (Jennings et al. 2016), and hence needs no external inputs but the vectors of fiducial cosmological and nuisance parameters. For the reader wishing to repeat the validation, we list below some of the lessons learnt in the code comparison process.

• PySSC needs as input the WL and GCph kernels of Eqs. (18) and (21), as well as their argument, the redshift values. The code then uses this redshift array to perform the necessary integrals in dV through Simpson’s rule. The user is responsible for sampling the kernels on a sufficiently fine z grid [O(10⁴) values have been found to be sufficient in the present case] to make sure these integrals are performed accurately.

• The latest version of PySSC accepts a convention parameter. This specifies whether the kernels are in the form used in LG19 (convention = 0) or the one prescribed in EC20 (convention = 1). The two differ by a l/r²(z) factor, as shown in Eq. (23). Passing the kernels in the EC20 form without changing the parameter’s value from 0 – the default – to 1 will obviously yield incorrect results.

• The ordering of the S_ijkl matrix’s elements depends on the ordering chosen when passing the input kernels to PySSC – whether $W_i^{\text{L}}(z)$ first and $W_i^{\text{G}}(z)$ second or vice versa. This must be kept in mind when implementing Eq. (10).

• The GCph constraints can show a discrepancy greater the 10% for the dark energy equation of state parameters w₀ and w_a even when the corresponding covariance is found to be in good agreement. This discrepancy is due to GCph being less numerically stable because of the lower constraining power compared to the other probes, and because the bias model considered has a strong degeneracy with σ₈, making the numerical derivatives unstable (see e.g. Casas et al. 2023). Since this is a known issue, not coming from the SSC computation, and the covariance matrices and angular PS show good agreement, we choose to overcome the problem by using, for GCph, one code to compute both sets of parameter constraints (that is, we run one FM evaluation code with as input the covariance matrices from both groups).

Fig. A.1 Some of the most important elements examined in the comparison. The arrows show the ordering followed to produce the parameters constraints, which is opposite to the one followed to validate the code. The derivatives of the PS with respect to the cosmological parameters, entering the final step of the computation, are not shown.

Appendix B High order bias from halo model

As described in Sect. 3.6.1, the higher-order bias b₍₂₎(z) has been estimated using the halo model. In the following, we provide further details on the input quantities, and how we set the relevant parameters.

A key role is played by the halo mass function Φ_MF(M, z), which we model as $${\Phi }_{\text{MF}}(M,z)=\frac{{\overline{\rho }}_{\text{m}}}{M}f(v)\frac{\text{d\hspace{0.17em}}\ln {\sigma }^{-1}}{ \text{d}M},$$(B.1)

with M the halo mass, ${\overline{\rho }}_{\text{m}}$ the mean matter density, v = δ_c/σ(M,z), δ_c = 1.686 the critical overdensity for collapse, and σ(M, z) the variance of linear perturbation smoothed with a top-hat filter of radius $R={\left[3M/\left(4\pi {\overline{\rho }}_{\text{m}}\right)\right]}^{1/3}$. We follow Tinker et al. (2010), setting $$f(v)={\mathcal{N}}_{\text{MF}}\left[1+{\left({\beta }_{\text{MF}}v\right)}^{-2{\varphi }_{\text{MF}}}\right]v^{2{\text{M}}_{\text{MF}}}\exp \left(-{\gamma }_{\text{MF}}{v}^2/2\right),$$(B.2)

where N_MF is a normalisation constant, and the halo mass function fitting parameters β_MF, η_MF,γ_MF and ϕ_MF – not to be confused with Φ_MF(M, z) – scale with redshift as illustrated in Eqs. (9–13) of the above-mentioned paper.

The other quantity needed is the average number of galaxies hosted by a halo of mass M at redshift z. This is given by $$\langle N\mid M\rangle (M)=N_{\text{cen}}(M)\left[1+N_{\text{sat}}(M)\right],$$(B.3)