© The Authors 2023

1. Introduction

Super-sample covariance (SSC) (Hu & Kravtsov 2003; Takada & Hu 2013), also known as sample variance for cluster counts, is an important part of statistical errors or cosmic variance for the current survey of the large-scale structure (LSS) of the Universe (e.g., Hildebrandt et al. 2017), and even more for near future cosmic surveys that aim to deepen our understanding of dark energy and gravity. Specifically, it is expected to be the dominant non-Gaussian contribution to the covariance of weak-lensing statistics (Barreira et al. 2018a), although other covariance terms are also important for galaxy clustering (Lacasa 2018, 2020). It has been shown that for a stage IV survey, the dark energy figure of merit can be degenerate by a factor of two (Gouyou Beauchamps et al. 2022) when the SSC is accounted for in the joint analysis of weak-lensing (WL), photometric galaxy clustering (GC). This joint analysis is also called 3×2pt analysis. However, it is difficult to estimate the SSC because it cannot be calibrated reliably with jackknife or bootstrap methods or simulations with sizes comparable to the survey size (Lacasa & Kunz 2017). Instead, significant progress has recently been made in the analytical modelling of the SSC (e.g., Takada & Spergel 2014; Takahashi et al. 2014; Li et al. 2018; Chan et al. 2018; Lacasa et al. 2018; Wadekar & Scoccimarro 2020) as well as in semi-analytical predictions of the SSC using measurements of the simulated reaction of observables to super-survey modes (e.g., Barreira & Schmidt 2017; Barreira et al. 2018b). An alternative approach that includes the effect of super-survey modes is to account for them with additional nuisance parameters that are marginalised over, that is, to account for them at the likelihood level and not in the covariance (Li et al. 2014; Lacasa & Grain 2019). In this article, we examine the two approaches to the modelling of the impact of super-survey modes (see Sect. 2).

We address four challenges in implementing the SSC in an analysis pipeline: change, speed, reliability, and stability. The first challenge, change, is the difficulty of changing an existing pipeline that may use assumptions that are broken by the SSC (e.g., that the covariance matrix is diagonal) and comparing new results with the previous pipeline version. The second challenge is that the covariance pipeline may be slowed down by the time necessary to predict the SSC and by the inversion of the total covariance matrix (if large), which is no longer diagonal. This challenge is not relevant when the covariance is fixed, but when it is varied with parameters and/or when the robustness of the results is tested with respect to changes in the covariance. The third challenge is about the stability of the results with respect to variations in the theoretical ingredients of the SSC prediction because the SSC literature is quite technical and evolves rapidly. The fourth challenge is a practical one found by the authors: the SSC gives an ill-conditioned covariance matrix so that the inversion becomes unstable when the other covariance term becomes negligible (e.g., when an analysis is pushed in the deep non-linear regime or in a combined probe case with large matrices). In the first section of this article, we propose a framework that solves these four problems. We then develop theoretical and numerical applications of the framework.

In detail, the article is organised as follows. In Sect. 2 we show that SSC can be computed as an update to the precision (i.e. inverse covariance) matrix, present examples for different probes, and show the application in likelihood or Fisher pipelines. In Sect. 3 we use the previous formulae to show the equivalence of three approaches to SSC: at the covariance level, at the likelihood level, and with a quadratic estimator. In Sect. 4 we apply the framework to a mock stage IV photometric galaxy survey, for WL alone or in combination with GC. We conclude in Sect. 5.

2. SSC as an update to the precision matrix

In this section, we show that the inverse covariance matrix can be computed as an update to the noSSC case in a numerically efficient manner. We then show applications to compute the likelihood or the Fisher matrix of an LSS observable.

To do this, we use the Woodbury identity for matrix inversion (Woodbury 1950). This identity states that if A is an invertible (n, n) matrix of known inverse that is updated with a matrix of the form UCV, with U, C, and V the matrices of the respective shapes (n, m), (m, m), and (m, n), and C is invertible, then the inverse of the updated matrix is given by

(1)

More generally, even if C is not invertible, we have

(2)

We use the general form of the Woodbury identity given by Eq. (2) for the case that A and C are symmetric matrices and V = U^T. The original and updated matrices are therefore both symmetric.

It is not immediately obvious that the right-hand side of Eq. (2) is also symmetric, as it should be. However, this becomes clearer when C is invertible, where we can use the equivalent form in Eq. (1), which shows that the right-hand side is symmetric. Then we can note that matrix symmetry is a continuous property, and invertible matrices are dense in matrix space. The property therefore holds by continuity even if C is not invertible.

2.1. General formula

Let 𝒪(i, α) be an abstract LSS probe, where α is the index (multi-index) for the redshift bin (redshift tuple), and i is the index for all other quantities on which the observable depends (e.g., Fourier or angular wavevector, distance or angular separation, mass, size, or peak height bin). This observable is vectorised through an index i_o, with a mapping i_o → (i, α) that is considered implicit hereafter.

We call n_z be the number of indices α and n_o the number of indices i_o, that is, the size of the observable. For instance, if there is the same number n_i of quantities for each redshift bin, then we have n_o = n_i × n_z. This situation can, however, be broken by analysis choices employing a redshift-dependent range for the i quantities, for example a cut k_max(z) or ℓ_max(z).

We use the Sij approximation to SSC (e.g., Hu & Kravtsov 2003 for cluster counts, Lacasa & Grain 2019 for power spectra and other observables). This is a mean-field approximation for the redshift evolution, which is expected valid for observables with compact support in redshift and has been validated for use for cosmological analysis with clusters (Hu & Kravtsov 2003; Aguena & Lima 2018) or galaxies (Lacasa & Grain 2019). In this approximation, the SSC of the probe is given by

(3)

where is the change in matter density in the survey with respect to the global mean. This is also known as the super-survey mode. The S matrix is the covariance matrix of these δ_b in each redshift bin, due to the clustering of matter. ∂_δb𝒪 is the reaction of the probe δ_b, also known as the response. In practice, this response is computed either from a nonlinear theory (e.g., Takada & Hu 2013; Lacasa & Rosenfeld 2016), from an ansatz (Lacasa & Grain 2019), or from simulations (Barreira et al. 2018b).

In the following we often shorten 𝒞_SSC ≡ Cov_SSC, that is,

(4)

where i_o → (i, α) and j_o → (j, β) through the implicit mapping mentioned previously.

We now remark that this SSC matrix can be cast in the following form:

(5)

where S is a square matrix of shape (n_z, n_z) (making it evident that 𝒞_SSC has rank ≤n_z), and U is a rectangular matrix of shape (n_o, n_z) given by

(6)

where the final δ on the right-hand side is a Kronecker symbol.

Therefore, using the Woodbury identity in Eq. (2), the inverse of the total covariance matrix 𝒞_tot = 𝒞_noSSC + 𝒞_SSC is given by

(7)

where we introduced the notation (that we justify below)

(8)

Formula Eq. (7) thus gives the total precision matrix as an update to the noSSC case. It is particularly useful numerically when the SSC dominates the covariance matrix elements: the SSC matrix indeed has a low rank and will thus destabilise the covariance inversion. This was, for instance, noted previously by one of the authors in numerical investigations pushing angular power spectra analyses at high multipoles in the SSC-dominated regime (Lacasa 2020).

2.2. Examples

We now apply the general formula described in the previous subsection to specific class of LSS observables. The case of one-point functions is described in Sect. 2.2.1, the case of two-point functions is set out in Sect. 2.2.2, and the case of other statistics is presented in Sect. 2.2.3.

2.2.1. Counts or one-point functions

The first class of observables includes counts of clusters as a function of (bins of) mass proxy (and redshift), of voids as a function of size, or of WL peaks as a function of peak height. In all cases, the observable can be noted N_obj(i, α): the number of objects as a function of the index i of quantity (mass, size, and height) and redshift index α. For clusters and voids, and except for factors of confusion in bin assignment (e.g., mass errors or photo-z errors), the noSSC (=standard) covariance is Poissonian,

(9)

where the final δ on the right-hand side are Kronecker symbols, that is, the noSSC matrix is diagonal. This diagonality somewhat simplifies Eq. (7),

(10)

which can be written in a more compact form as

(11)

Here, the matrix I(δ_b) has the shape (n_z, n_z), and is given by

(12)

We note that 𝟙 + I(δ_b) S is not symmetric (but S ⋅ (𝟙 + I(δ_b) S)⁻¹ is) and must be inverted numerically. This is not a problem because (i) it is regular (even if S is not invertible, the total is regularised by the presence of the identity matrix), and (ii) it has a small size: the number of redshift bins.

2.2.2. Power spectra or two-point functions

The second class of observable includes the (angular) power spectra or the (angular) correlation function of galaxies, weak-lensing, clusters, cosmic microwave background lensing, HI intensity mapping, and so on as a function of Fourier (harmonic) wave number or (angular) separation and of redshift bin. If we only consider auto-spectra (i.e. correlate only in the same redshift bin), then α is the index of the redshift bin. If we also consider cross-spectra, then α is the index of the redshift pair. Even in the noSSC auto-spectrum case, there is often cross-talk between redshift bins, for example because of photo-z errors, so that the covariance is not diagonal in this direction. However, the noSSC covariance is diagonal in terms of Fourier or harmonic wave number under two conditions: (i) it may be argued that the Gaussian covariance is the dominant term beyond SSC (Barreira et al. 2018a showed this is the case for WL, although Lacasa 2020 reported that other non-Gaussian terms are important for constraints on n_S with GC), and (ii) the wave numbers are binned sufficiently to erase the mode-coupling due to the survey mask (see e.g., Hivon et al. 2002). As for counts, this slightly simplifies Eq. (7), although less so because we still have redshift correlations. For example, for an angular power spectrum C_ℓ(α), the Gaussian covariance is a bloc diagonal multipole by multipole. We call the multipole bloc with the shape (n_z, n_z),

(13)

and denote its inverse with . We thus have

(14)

which can be written in a more compact form as

(15)

where the matrix N_ℓ has the shape (n_z, n_z) and we have

(16)

Here, the matrix I(δ_b) has the shape (n_z, n_z) and is given by

(17)

Equation (15) greatly decreases the numerical cost of inverting 𝒞_tot. For the covariance matrix of shape (n_ℓ × n_z, n_ℓ × n_z), the cost of brute force computation is 𝒪((n_ℓ × n_z)³), while through Eq. (15), we invert only (many) (n_z, n_z) matrices, reducing the cost to . Therefore, we achieve a speedup factor of . For a typical analysis with tens of multipole bins (e.g., we use 40 bins in Sect. 4), this results in an appreciable speedup factor in the hundreds. This speedup means that a larger number of multipole bins can be used instead of binning aggressively to extract as much information as possible from the measurement, for instance, to hunt for localised features.

2.2.3. Other statistics

The formulae can easily be generalised to other statistics. For instance, a combination of angular power spectrum measurements of different probes amounts to generalising the α multi-index to index both probes and redshift tuples. We use this in Sect. 4 for the so-called 3×2pt analysis, which combines WL, GC, and their cross correlation (XC). Higher-order statistics are also possible. For instance, for a bispectrum analysis, i would become a multi-index for triplets of wavenumbers, and α becomes a multi-index for triplets of redshift bins.

In all cases, the formula becomes particularly useful when the noSSC covariance matrix is diagonal in some direction, as it can leverage this to lower the numerical cost of the inversion. For instance, for a combination of power spectra, the Gaussian covariance only correlates with identical wavenumbers, and for bispectra, it only correlates with identical wavenumber triplets.

2.3. Applications

2.3.1. Log likelihood

For a standard Gaussian likelihood with fixed covariance, we have

(18)

where d is the data vector and m is the corresponding model.

Adding SSC can thus be considered as an update to the log-likelihood: lnℒ = lnℒ_noSSC + δlnℒ with

(19)

where . This makes it easy and numerically efficient to update a noSSC pipeline to account for SSC. In the noSSC pipeline, the covariance matrix 𝒞_noSSC should already have been inverted, so no numerical cost is added there. Then we have a matrix operation to downsize the problem into multiplication and inversion of small (n_z, n_z) matrices, which is cheap with respect to all the other operations. Hence the computation of δlnℒ takes a time comparable to that of lnℒ_noSSC and is not a bottleneck of the pipeline.

If the covariance is allowed to vary with model parameters, there is an additional determinant term in the likelihood. The impact of SSC can be computed through the determinant lemma,

(20)

which gives

(21)

We thus keep the property lnℒ = lnℒ_noSSC + δlnℒ, and the additive determinant term is numerically cheap since it is the determinant of a (n_z, n_z) matrix.

2.3.2. Fisher matrix

We consider the observable 𝒪 as a vector. The vector of the model parameters is θ with a length n_θ, and the parameters are indexed with Greek indices μ, ν, …. With these notations, the Fisher matrix elements are given by

(22)

where ∂_μ𝒪 is the partial derivative of the observable with respect to the parameter μ, considered as a vector. If we plug in Eq. (7), this yields

(23)

where is a vector with the size n_z. We can thus consider Y as a (n_z, n_θ) matrix, which allows rewriting Eq. (23) in matrix form,

(24)

In particular, for a single redshift bin (n_z = 1), we obtain

(25)

where Y is a covector and is a scalar that can be interpreted as the Fisher information on δ_b as an additional parameter. This interpretation will become clearer in Sect. 3. We note that Eq. (25) recovers a result from Lacasa & Grain (2019), specifically, their Eq. (15) in matrix form, with the mapping , , and S ↔ S_i, i. We note, however, that the latter result is limited to a single redshift bin (which is admittedly of little practical interest), while the general result in Eq. (23) applies to any observable(s) in any number of bins.

In the multi-bin case, I(δ_b) is a (n_z, n_z) matrix that represents the Fisher information on a vector of additional parameters δ_b(α), indeed we have

(26)

2.3.3. Inverse Fisher matrix

The inverse of the Fisher matrix is required to forecast marginalised errors on the parameters. Fortunately, the form of Eq. (24) means that the Woodbury inversion formula can be used again. After some algebra, this yields

(27)

(28)

with

(29)

This can be rewritten in the form that is used throughout this article:

(30)

with

(31)

(32)

whose physical interpretations will become clearer later.

Equation (30) allows the impact of SSC on the parameter constraints to be forecast as an update to the noSSC case. This has advantages in terms of speed for two reasons: (i) The noSSC pipeline does not need to be rewritten nor rerun, its results can be reused, and (ii) if an assumption of diagonality of the noSSC covariance is used, it can also be used for the SSC Fisher elements: it will speed up the computation of Y, and from then on, only cheap matrix multiplications are necessary (matrix sizes: n_z or n_θ).

3. Equivalence of SSC approaches

In this section, we show that three approaches to SSC are equivalent under the conditions detailed in the respective subsections. The first is the classical approach used in the previous section, where SSC is included as an additional covariance term. In the second approach, SSC is instead accounted for by promoting δ_b to be a nuisance parameter (single bin case; Li et al. 2014) or vector of parameters (one per bin, multi-bin case; Lacasa & Grain 2019) that is marginalised over with a prior 𝒩(0, S). The third approach is new. It estimates δ_b with a quadratic estimator (QE) to undo the SSC effect.

3.1. Equivalence between the covariance and likelihood approaches in the single bin case

In the case of a single redshift bin, the inverse Fisher matrix in the covariance approach is a special case of Eq. (30),

(33)

In the likelihood approach, the covariance of the data is the noSSC covariance, and the parameter vector has been extended to include δ_b : θ′=(θ, δ_b). This gives a Fisher matrix,

(34)

where μ and ν are the indices of the model parameters, running from 1 to n_θ.

On top of this, a prior 𝒩(0, S) is imposed on δ_b. Because the Fisher matrix is the inverse covariance of the parameters and because inverse covariances of independent Gaussians add up, imposing this prior is equivalent to adding a Fisher matrix,

(35)

where the δ on the right-hand side are Kronecker symbols.

When we introduce the notations of Sect. 2.3.2, the total Fisher matrix F′=F(freeδ_b)+F(prior) then takes the form

(36)

Then some algebra that is detailed in Appendix B for transparency shows that the inverse Fisher matrix can be written as

(37)

In particular, for the lower right element, we have

(38)

The upper left block for the parameters θ reads

(39)

It is identical to Eq. (33), which shows that the covariance and likelihood approaches give the same constraints on the model parameters, that is, they are equivalent.

It seems natural that the equivalence should extend to the case of several bins, although we did not attempt this calculation.

3.2. Equivalence between the QE and covariance approaches for the power spectrum in the same multipole range

In this section, we denote estimators with an upper hat. For example, is the angular power spectrum estimator. The expectation value is noted with an upper bar, for instance , and can be predicted with a cosmological code as a function of θ the vector of the model parameters. We assume the theoretical predictions to be unbiased.

A QE can be built to estimate the large-scale modes of the density field using their coupling to better measure small-scale modes (see e.g., Li et al. 2020; Darwish et al. 2021). The natural idea is to use such an estimator to measure δ_b, which is the zero mode of the survey, and then combine this with the traditional power spectrum measurement. However, for this zero mode, the QE collapses to a linear combination of power spectrum measurements in full sky. For instance, in terms of angular observables and using notations from the QE literature, Appendix A shows that in the full sky, the δ_b estimator is

(40)

with

(41)

and

(42)

This is consistent with comparable equations from Li et al. (2014), who worked on 3D observables and found an estimator written in terms of P(k) measurements.

If the QE and the power spectrum are estimated on the same field, and if the multipole range used in the QE is included in the multipole range in which the power spectrum is measured, then the QE measurement is entirely determined by the power spectrum measurement. The joint probability is accordingly

(43)

where is fixed at its value in the fiducial cosmology needed to build the QE. As in Sect. 3.1, we note θ′=(θ, δ_b) the vector of parameters. The likelihood then reduces to

(44)

where the first term is the usual likelihood from the power spectrum alone. In the second term, δ_D is a Dirac distribution. This second term has no dependence on the model parameters (e.g., for the weights g_ℓ, ℓ(0), all cosmological parameters are fixed at the fiducial value used to build the estimator). The Fisher matrix is given in terms of the likelihood by

(45)

where ⟨|θ′⟩ denotes the conditional expectation given the values of θ′.

If we call F″ the Fisher matrix of this QE + power spectrum approach, we obtain

(46)

where the second term vanished because the derivatives with respect to model parameters vanish. Thus we have the same Fisher constraints as the likelihood of the power spectrum, as well as the traditional covariance approach for constraints on θ.

3.3. Physical interpretation and conclusion

In this section, we use the two previous sections to give a physical interpretation to some of the quantities introduced in Sects. 2.3.2 and 2.3.3 : I(δ_b), Σ²(δ_b), δI, and Z.

First, Eq. (36) shows that I(δ_b) is the direct Fisher information on δ_b. Furthermore, from the QE normalisation in Eq. (41), we can recognise that N = I(δ_b). Computing the variance of the QE gives

(47)

Hence I(δ_b) is the inverse variance of the maximum likelihood estimator for δ_b, which consolidates its interpretation as the raw information on δ_b when all other model parameters are fixed.

Second, Σ²(δ_b) is defined in Eqs. (32), and (38) (in the single bin case) shows that it is the marginalised Fisher constraint on δ_b. In the single bin case, it can be cast in the form

(48)

which makes it clear that there is a difference with the unmarginalised QE error due to (i) the presence of a prior , which adds information, (ii) the negative term −δI, which removes information. In the limit of no prior, we have S → ∞ and . We therefore conclude that δI quantifies the loss of information on δ_b due to the degeneracy with the other model parameters in θ.

Third, the off-diagonal matrix elements of Eq. (37) show that the covector Z quantifies the degeneracy between δ_b and each model parameter in θ. Furthermore, when we write Eq. (39) on the diagonal, we have

(49)

which shows that Z_μ is the (square root of the) increase in variance on parameter μ as a function of the constraint on δ_b. To reflect this, we could have noted this quantity or perhaps , but this would have made the notations of all previous equations in Sect. 3 much heavier and harder to read. We therefore chose to keep a simple Z.

Finally, as a summary of this section, we have shown that accounting for SSC at the likelihood level is formally equivalent, in terms of parameter constraints, to accounting for it traditionally at the covariance level. Combining a QE and a power spectrum does not bring better parameter constraints than the power spectrum alone as long as both estimators act upon the same field and the same multipole range.

The estimator might be used mainly (i) when some observational effects can be dealt with more efficiently (e.g., survey mask or systematics) in the QE than in the power spectrum analysis, (ii) if it allows adding information from another field (e.g., estimate δ_b with galaxies to improve a WL analysis), from another statistic (e.g., a linear, cubic, or log-normal estimator that uses information from the one-point or three-point function or from the log field), or from another survey (see e.g., Digman et al. 2019) when a full combined probe analysis would be too costly.

4. Application to a stage IV weak-lensing and galaxy clustering analysis

In this section, we apply the previous results to the case of a future stage IV photometric galaxy survey that constrains cosmology through WL alone and in combination with GC and XC. We describe the experimental and modelling setup in Sect. 4.1. In Sect. 4.2 we examine how well the survey can constrain the super-sample modes. In Sect. 4.3 we ask whether any theoretical prior on δ_b can be lifted. Finally, in Sect. 4.4 we study the precision to which the SSC responses should be modelled in order to pass certain precision requirements on the reported cosmological errors.

4.1. Setup

For the cosmological modelling, we worked within the w₀w_a CDM model, and in the Fisher forecast, we varied each of its parameters: Ω_m, Ω_b, w₀, w_a, h, n_S, and σ₈. In practice, we mostly focused on the parameters that are most affected by SSC: w₀, w_a, σ₈, and Ω_m. For the prediction of all cosmological quantities, we used the software cosmosis by Zuntz et al. (2015), and for the prediction of the matter power spectrum P(k) on non-linear scales, we used the prescription by Takahashi et al. (2012) and Bird et al. (2012).

In terms of experimental specifications, we followed all the specifications of Euclid Collaboration (2020) for a 15 000 deg² photometric galaxy survey to analyse WL, GC and XC. In slightly more detail, the galaxy redshift distribution is

(50)

with and a total density of 30 galaxies per square arcmin. Then we defined ten tomographic redshift bins such that they are equi-populated, yielding the boundaries z_i = (0.0010, 0.42, 0.56, 0.68, 0.79, 0.90, 1.02, 1.15, 1.32, 1.58, 2.50). The redshifts were then perturbed by photometric errors with the main dispersion and a population of outliers, following the specific details given in Euclid Collaboration (2020). These galaxies are both the source and lens sample, that is, their shape was used for the WL measurement and their position for the GC measurement.

In terms of probe modelling, we also followed Euclid Collaboration (2020). The GC was modelled with a scale-independent bias with a fiducial value . For cosmological constraints, the galaxy bias value in each tomographic bin is a nuisance parameter that needs to be marginalised over. The WL is modelled classically as a tracer of matter, its power is hampered on small scales by the shape noise with standard deviation σ_ϵ = 0.3. Another issue for WL is the contribution from the intrinsic alignments (IA) of galaxies. To model this effect, we followed Euclid Collaboration (2020): the IA power spectrum and cross-spectrum with matter are given by P_II = [ − A(z)]²P_mm and P_mI = −A(z)P_mm, where

(51)

with 𝒞_IA = 0.0134 being a normalisation constant, D(z) the growth factor, and 𝒜_IA a nuisance parameter that controls the amplitude of the IA contribution. The redshift dependence is given by

(52)

where ⟨L⟩(z)/L_*(z) is the ratio of the mean source luminosity and the characteristic scale of the luminosity function (Hirata et al. 2007; Bridle & King 2007). Additionally, we have two nuisance parameters, η_IA and β_IA, making a total of three nuisance parameters for intrinsic alignments, with fiducial values (Euclid Collaboration 2020) (𝒜_IA, η_IA, β_IA) = (1.72, −0.41, 2.17).

For the statistical observable, we used the angular power spectrum for all probes: the WL auto-spectrum, the GC auto-spectrum, and the XC cross-spectrum. For simplicity, we adopted a common range of multipoles from ℓ_min = 10 to ℓ_max = 3000 for all probes. The covariance of these power spectra is described by a classical Gaussian term (including shape noise for WL and shot-noise for GC), and the super-sample covariance term is accounted for through the developments of Sect. 2. Specifically for the SSC, we computed the S matrices using PySSC¹ (Lacasa & Grain 2019), and for the responses, we followed Gouyou Beauchamps et al. (2022) and took ∂_δbC_ℓ = R × C_ℓ with a probe-independent dimensionless response R = 4. Following Barreira et al. (2018a), other non-Gaussian covariance terms have a negligible impact on cosmological constraints from a WL survey, hence we neglected these terms here.

4.2. Information on δ_b

In this section, we compare the information on δ_b brought by the survey, either raw information with I(δ_b) defined by Eq. (8) or marginalised information with I − δI defined by Eq. (29), to the prior information given by the S matrix for the considered redshift bins.

We first consider the WL-only analysis. The first case of interest is to try to constrain δ_b in a redshift bin by only using the information in this bin. The result is shown in the top panel of Fig. 1, where we plot I(δ_b), I − δI, and 1/S in each of the ten redshift bins. The raw information (solid blue) is an order of magnitude larger than the prior (dash-dotted green), which appears to mean that the data can constrain δ_b very well. However, the marginalised information (dashed orange) falls by ∼3 orders of magnitude compared to the unmarginalised information, making it lower than the prior by ∼2 orders of magnitude. This comes from a high cancellation of I(δ_b) and δI due to high degeneracy between δ_b and cosmological parameters. We found that this is, in fact, not peculiar to δ_b; with a single redshift bin, most cosmological parameters are near perfectly degenerate, and the Fisher matrix has very poor conditioning.

Fig. 1. Information on δb: Raw (I), marginalised (I − δI), and prior (1/S). Top: information in each bin of redshift for WL analyses of each bin separately. Bottom: information in each pair of redshift bins for a joint tomographic WL analysis. The 55 pairs are ordered as (0, 0)(0, 1)...(0, 9)(1, 1)..(1, 9)...(8, 8)(8, 9)(9, 9). The light vertical lines indicate the jump of the first redshift bin, i.e., the leftmost line is at (0, 9) before the jump to (1, 1).

This fact motivated us to pursue the single bin case no longer and instead analyse the tomographic case of a joint analysis of all bins and further pairs of bins together. The result is shown in the bottom panel of Fig. 1, where we plot I(δ_b), I − δI, and 1/S in each of the 55 redshift pairs. Explicitly, each of these quantities is now a (55, 55) matrix, and the figure shows the diagonal of these matrices.

In this tomographic case, the raw information is larger than the prior by an even larger amount, ∼2 orders of magnitude, than in the single bin case. The parameter degeneracies are now all lifted, both for cosmological parameters and for nuisance parameters such as δ_b. We attribute this to the fact that the tomography provides information on the history of expansion and growth. Hence, here the marginalised information I − δI sits almost perfectly on top of the unmarginalised information I(δ_b), with only slight loss in the high-redshift pairs.

Next, we studied the 3×2pt analysis, which combines the WL auto-spectra, the GC auto-spectra, and the XC cross-spectra. There are now 210 bin pairs: 55 for WL, 55 for GC, and 100 for XC because the latter is not symmetric. Figure 2 shows I(δ_b), I − δI, and 1/S in each of the 210 redshift pairs. We chose to order the pairs with WL first, then XC, and then GC; the corresponding intervals are indicated by vertical dashed lines in the figure.

Fig. 2. Information on δb in each redshift bin pair in a tomographic 3×2pt analysis (WL + GC + XC). There are 210 pairs: 55 for WL, 100 for XC, and 55 for GC. The pairs are ordered as in Fig. 1. The only difference is that XC has more pairs because there is no symmetry (i, j)↔(j, i).

The marginalised information is oftentimes indistinguishable from the unmarginalised one and is larger than the prior. Deviations from this behaviour are visible as breaks of the lines in the XC and GC pairs. We cut the y-axis of the plot to zoom in on the interesting region. Otherwise, the blue and orange lines fall extremely low, making the plot unreadable. For these pairs, the two redshift kernels overlap extremely little, so that the signal is practically absent. This does not happen for WL because the WL kernels are broad in redshift, integrating everything between the source and the observer so that all kernels overlap strongly. However, the GC kernels are much more localised, and there is, for instance, practically no overlap between the first and last redshift bins. The cross-correlation XC will also vanish when measured between a low or median WL bin and a high GC bin. We note that the prior (green line) also breaks in these pair cases because PySSC sets the S matrix element to 0 when the overlap of kernels is too small for a reliable computation. In practice, however, this is not an issue because in these pairs, the signal is too low to contain any cosmological information. Essentially, Fig. 2 shows that δ_b can be well constrained for all relevant portions of the data.

A theoretical concern that could come to mind, already for the tomographic WL analysis and even more for the 3×2 pt analysis, is the physical significance of these numerous δ_b parameters: one for each bin pair. Through Eq. (6), they change only the observable for that bin pair; physically, they are the survey matter density contrast averaged over redshift with the product of the two kernels of the pair. The concern is that all these pairs probe a similar redshift range, and that therefore the parameters should be strongly correlated. The answer is that this is indeed the case, and in practice, the correlation of these parameters is quantified by the S matrix, which has high off-diagonal correlation coefficients. When the formulae developed in Sect. 2 for the likelihood or Fisher forecast are used, the δ_b parameters do not appear in the computation, and all the formulae are well behaved numerically. In a different context, for a different kind of analysis that explicitly has δ_b nuisance parameters, there could be numerical problems with their number and high correlation. In that case, Appendix C devises an analysis with fewer less strongly correlated δ_b parameters, at the expense of needing more observable predictions.

4.3. Asymptotic errors

In this section, we determine the change in the Fisher constraints when an infinitely thin or wide prior 𝒩(0, S) is assumed. In the single bin case, this corresponds to the respective cases S → 0 and S → ∞.

The case of an infinitely thin prior is simple. From Eq. (25) (single bin case) or Eq. (24) (multi-bin case), the Fisher constraint asymptote to the noSSC one. This happens because the prior fixes δ_b to 0 in each bin, so that there is no SSC effect.

The case of an infinitely wide prior is more interesting. It corresponds to δ_b being completely free, S → ∞. It might be motivated by the observation of Sect. 4.2 that the data appear to constrain δ_b much better than the prior, so that the latter seems not to be needed.

The Fisher matrix was already given in Eq. (34), but for the parameter vector θ′=(θ, δ_b). After marginalisation over δ_b, the inverse Fisher is

(53)

The corresponding marginalised error bars can be compared to the fiducial case in which S is given by the theoretical prediction from PySSC. The ratios of these two cases for cosmological parameters (after marginalisation over all the others, as well as the nuisance parameters) are given in Table 1. This conservative no-prior analysis would increase the cosmological errors by unacceptably large factors between 1.1 (Ω_b, WL-only) and 5 (w_a, WL-only), with an average of 2.8 for WL-only and 2.5 for 3×2 pt.

At first, this may seem at odds with the result of Sect. 4.2 that the data appear to constrain δ_b much better than the prior. However, this happens because the data mostly constrain each δ_b independently, while, as discussed at the end of Sect. 4.2, the prior has strong off-diagonal correlations, which encode the fact that the δ_b probe similar redshift ranges. In the no-prior case, we therefore assign much freedom to the data, which blows up the constraints.

4.4. Requirements for reaching robust cosmic errors

The results of Sect. 4.3 pose the more general question of how precisely we need to model the super-sample covariance in order for cosmological constraints to be robust. This question is better posed at the covariance level than at the likelihood level used in Sect. 4.3. We focus on the impact on the error bars, although, as noted by Sellentin & Starck (2019) using an incorrect covariance matrix is also equivalent to analysing a biased data set.

To answer this question, we let the SSC responses vary both in amplitude and scale dependence and forecast the variation in cosmological constraints with respect to the case with the fiducial responses.

4.4.1. Weak-lensing-only case

For the WL case, we introduced two parameters, A₀ and A₁, which change the response as

(54)

with ℓ₀ a pivot scale chosen to be the centre of the multipole range so that the fraction varies by 1 between ℓ_min ∼ 0 and ℓ_max = 2ℓ₀.

Based on the developments of Sect. 2.3, we can compute a Fisher forecast without building and inverting a new covariance matrix for each value of (A₀, A₁) we wish to test. Instead, we first computed the Fisher forecast for the Gaussian covariance matrix and then used Eq. (30) to determine the inverse Fisher matrix. We found that the bottleneck is the computation of I(δ_b), whose brute-force computation takes ∼3 s on a laptop with i7 CPU. By comparison, a covariance inversion takes ∼120 s on the same laptop, which would make the forecast painful to run when many values for A₀ and A₁ need to be tested. This speedup was achieved without using the fact that the Gaussian covariance is diagonal in terms of multipoles.

We can then study how the marginalised errors on cosmological parameters change with A₀ and A₁. The results are shown in Fig. 3 for the most affected parameters, which are in fact the usual targets of WL experiments: the dark energy equation of state, σ₈ and Ω_m.

Fig. 3. Forecasted errors on the dark energy equation of state, σ8 and Ωm for a tomographic WL analysis. Top two rows: constraints as a function of A0 (which changes the amplitude of the SSC response). Bottom two rows: constraints as a function of A1 (which changes the scale dependence of the SSC response). Dotted lines show the value for the fiducial covariance, (A0, A1) = (1, 0), and the dashed lines show ±10% around the dotted line.

We let A₀ run between −1 and 3, that is, between predicting the response with a sign mistake and misestimating the response by a factor of 3. We let A₁ run between −2 and 2, which is a scale dependence that is strong enough to cause the response to vanish at the minimum or maximum bound of the multipole range. The plots also contain dashed horizontal lines to indicate when the errors deviate from the fiducial case by more than 10%, in line with the requirements of stage IV missions (e.g., Euclid, Laureijs et al. 2011).

The results show that the branch in which A₀ is negative is symmetric to the positive branch. This happens because A₀ is squared in the covariance equation, so that an overall sign mistake has no impact. As A₀ increases (in absolute value), so do the cosmological errors. The requirement that errors should have a precision of 10% constrains the allowed value of A₀. The strongest constraints come from Ω_m, giving A₀ ∈ [0.47, 2.82]. We note that this is a quite large interval: response predictions can be under- or over-estimated by a factor larger than 2.

For A₁ the curves are no longer symmetric and have a steep slope near the fiducial case A₁ = 0, so that they easily exceed the 10% requirements. Therefore, satisfying the required precision places strong constraints on A₁; the strongest constraint comes from σ₈, giving A₁ ∈ [ − 0.16, 0.15].

In summary, WL constraints are rather robust to a mismodelling of the amplitude of its response. However, the constraints are more sensitive to the scale dependence of the response, which must be modelled at the 15% level.

4.4.2. 3×2 pt case

For the 3×2 pt case, that is, when WL is combined with GC, including XC, we can first use the same parametrisation Eq. (54), with A₀ and A₁ being the same for all three probes. We can then observe the impact of varying these parameters, which change the amplitude and scale dependence of all probes together. The results are shown in Fig. 4 for the parameters that are usually most affected: dark energy, σ₈, and Ω_m.

Fig. 4. Forecasted errors on the dark energy equation of states, σ8, and Ωm for a 3×2pt analysis. Top two rows: constraints as a function of A0, which changes the amplitude of the SSC response. Bottom two rows: constraints as a function of A1, which changes the scale dependence of the SSC response. Dotted lines show the value for the fiducial covariance, and the dashed lines are ±10% around the dotted line.

The forecasts are generally more robust to these mismodelling parameters. For A₀, the stringent constraints are set by w₀, which requires A₀ ∈ [0.53, 1.70]. This is slightly stricter than the WL case, but still quite large. For A₁, the stringent constraints come from Ω_m, which requires A₁ ∈ [ − 0.44, 0.29]. The scale dependence must therefore be modelled at ∼30% precision; this constraint is twice looser than for the WL-only case.

These results hold for a common mismodelling of all responses of all probes, for example in a separate universe approach, if we have a systematic common to all probes in the simulations, or in a theoretical approach, if we were missing a term in our equations. To study the case when the physics of a single probe is modelled incorrectly, we defined a new model for the responses with new parameters ϵ_p for the amplitude of the response of each probe independently,

(55)

where (ϵ_WL, ϵ_GC, ϵ_XC) modulate the response of the power spectrum of WL, GC, and XC respectively. Figure 5 shows the impact of these parameters on the parameters that are usually most affected (dark energy equation of state parameters, σ₈, and Ω_m), ϵ_WL has little impact on the constraints, which easily stay within the requirements for most of the explored range of values. Only Ω_m discards some values. In particular, the constraints on dark energy are very robust throughout all explored ranges. The error bars are also quite robust to ϵ_GC, although to a lesser extent than for ϵ_WL. For instance, constraints on σ₈ stay within the requirements in all explored ranges. Lower and upper values on ϵ_GC are set by different cosmological parameters: to meet the 10% requirement, Ω_m requires ϵ_GC ≥ 0.01 and w_a requires ϵ_GC ≤ 1.71. The combined allowed range is thus ϵ_GC ∈ [0.01, 1.71], which is quite large. By contrast, in the case of ϵ_XC, the errors exceed the requirements much more easily. The lower constraint is ϵ_XC ≥ 0.69 due to w₀, and the upper constraint is ϵ_XC ≤ 1.34 due to Ω_m. The lower constraint could be relaxed by adopting a slightly looser requirement, as the curves drop just below the −10% line, although we do not enter into these considerations. In summary, the amplitude of this response must be modelled at a precision of ∼30%.

Fig. 5. For a 3×2pt analysis, forecast errors on dark energy, σ8, and Ωm. Top: as a function of ϵWL, which changes the WL response alone. Centre: as a function of ϵXC, which changes the cross-correlation response alone. Bottom: as a function of ϵGC, which changes the GC response alone. Dotted line: Value for the fiducial covariance. Dashed lines: ±10% around the dotted line.

It is expected that variations in the XC response have the most important effect because it has been shown that this cross correlation is of the utmost importance to control systematics and unlock the full power of the 3×2pt analysis (Tutusaus et al. 2020).

In summary, 3×2pt constraints are more robust to SSC mismodelling than WL-only constraints. The most important quantities to model correctly are the amplitude of the XC response and the scale dependence of all responses. Both need to be calibrated to the ∼30% level.

5. Conclusion

Within the framework of the Sij approximation to SSC, we have presented a new formula in which SSC is an update to the precision (inverse covariance) matrix. This update can be computed numerically in a time comparable to that of the noSSC precision matrix and can take advantage of any sparsity in the noSSC covariance. We showed practical examples for two types of cosmological probes: number counts, and power spectrum. We presented the application of this formula to the computation of the likelihood or a Fisher forecast as an update to the noSSC case. The Fisher matrix (and inverse Fisher matrix) were then used to prove the formal equivalence of three approaches to super-sample covariance: (i) classically at the covariance level, (ii) at the likelihood level through dedicated nuisance parameters, and (iii) with a QE to estimate the super-sample modes. Examination of these three approaches allowed us to give a physical interpretation of the various quantities entering the formula. Although these approaches are formally equivalent, we argue that the second approach, at the likelihood level, is conceptually more interesting for two reasons: (i) the δ_b nuisance parameters can be examined to track systematics, which would cause them to deviate significantly from their prior (e.g., a systematic inducing an excess of power in some bins would push δ_b > 0 for these bins), and (ii) it provides a natural framework to determine whether the data prefer low density (δ_b < 0) at low redshift. This potential explanation has been put forward in the past for acceleration or dark energy or the H₀ tension.

Furthermore, we presented a numerical application of this framework to a mock stage IV weak-lensing survey. Namely, we considered a WL-only and a WL + GC + XC (also known as 3×2 pt) case. We showed that although the δ_b parameters are well constrained, we cannot recommend to remove their theoretical prior because it would strongly degrade the cosmological constraints. We also examined the change induced by mismodelling on the cosmological errors, and especially how much mismodelling can be allowed to still satisfy the requirement of a precision of 10% for stage IV surveys. This is of interest as there are always theoretical or simulation uncertainties in the prediction of these responses. We found that large mispredictions of the response amplitude are allowed for the WL-only case, but the scale dependence must be predicted more precisely at the ∼15% level. For the 3×2pt case, large amplitude mispredictions are possible if it is common to all probes, but a misprediction specific to the cross-spectrum signal (XC = GC − WL power spectrum) is more problematic. This amplitude needs to be calibrated to a precision of ∼30%. The scale dependence of the response must also be predicted precisely, at the ∼30% level, although the requirement is less strict than for the WL-only case. In practice, this is the target of convergence that should be aimed for in the comparison of the response models or simulations.

Acknowledgments

We thank Stéphane Ilić for fruitful exchanges and Joachim Harnois-Déraps for proofreading the manuscript. F.L. acknowledges support from the Swiss National Science Foundation through grant number IZSEZ0_207059, and from Progetto di Eccellenza 2022 of the Physics and Astronomy department at the University of Padova “Physics of the Universe”, within the grant “CosmoGraLSS – Cosmology with Gravitational waves and Large Scale Structure”. JC and LL acknowledge support from a SNSF Eccellenza Professorial Fellowship (No. 186879). AMD acknowledges support from Tommala Foundation for research in gravity as well as Boninchi Foundation. SGB was supported by CNES, and focused on the 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” program (AMX-19-IET-008 -IPhU). AG’s work is supported by the McGill Astrophysics Fellowship funded by the Trottier Chair in Astrophysics, as well as the Canadian Institute for Advanced Research (CIFAR) Azrieli Global Scholars program and the Canada 150 Programme. The calculations in this article made use of https://class-code.net (Blas et al. 2011), https://cosmosis.readthedocs.io/en/latest/ by Zuntz et al. (2015), astropy (Astropy Collaboration 2013, 2018), matplotlib (Hunter 2007), scipy (Jones et al. 2001), and numpy (Oliphant 2006).

References

Appendix A: Quadratic estimator for δ_b

We wish to estimate the background density over the full sky, , which is proportional to the monopole of the matter density field, .

We use notations from the seminal article of Okamoto & Hu (2003) with slight variations to avoid problems at L = 0, and consider a single pair of fields. The estimator for the mode (L, M) has the general form

(A.1)

We apply this to L = M = 0, which leads to an important collapse. The Wigner 3J symbol enforces the triangular conditions |ℓ₁ − ℓ₂|≤L = 0 and m₁ + m₂ − M = 0, so that we have ℓ₁ = ℓ₂ and m₂ = −m₁. We obtain

(A.2)

where we absorbed multiplicative factors into g_ℓ, ℓ(0) and used . Then we use the definition of the power spectrum estimator in full sky , absorb the multiplicative factors again, and obtain

(A.3)

The normalisation A₀ has to be chosen such that the estimator is unbiased. In our case, , so

(A.4)

Our estimator therefore is

(A.5)

(A.6)

where the normalisation follows

(A.7)

We are left with optimising the weight g_ℓ, ℓ(0) (noting that a constant multiplicative factor in its definition is unimportant as it would cancel out in the ratio in Eq. A.6). This optimisation is performed by minimising the variance of the estimator,

(A.8)

where 𝒞_{ℓ, ℓ′} is the covariance of . A significant detail here is that all the brackets mean average at δ_b fixed,

(A.9)

Hence 𝒞_{ℓ, ℓ′} is the noSSC part of the covariance.

After some algebra, we find that the minimisation yields

(A.10)

We can therefore take

(A.11)

This also gives the normalisation of the estimator,

(A.12)

In summary, the quadratic estimator takes the form

(A.13)

with given by Eq. A.12 and g_ℓ, ℓ(0) given by Eq. A.11

Appendix B: Inverse Fisher matrix in the likelihood approach

In the likelihood approach, the total Fisher matrix F′=F(freeδ_b)+F(prior) takes the form

(B.1)

For readability, we now write I ≡ I(δ_b). In order to obtain the forecast constraints on θ marginalised over δ_b, we now need to invert F′ and extract the upper left (n_θ, n_θ) block. For this, we shuffle terms around so that we can write F′ in the form with

(B.2)

This makes easy to invert as it is block-diagonal. For the second term, we remark that

(B.3)

which puts F′ in a form amenable to using the Woodbury identity one more time. After some algebra, we find

(B.4)

Shuffling terms around, we obtain the cleaner form

(B.5)

Appendix C: Extension of the Sij approximation to probe-independent δ_bs

A LSS observable 𝒪(i, α) can generally be cast in the form of a light-cone integral of some density 𝔬(i, α),

(C.1)

(We note that for correlation functions of order n ≥ 2, this requires the use of Limber’s approximation.)

In this case, the equation for super-sample covariance is given by (e.g. Lacasa & Rosenfeld 2016)

(C.2)

where σ²(z₁, z₂) = Cov[δ_b(z₁),δ_b(z₂)] is the covariance of the background shifts of infinitesimal redshift slices. Starting from this equation, Lacasa & Grain (2019) derived the general Sij approximation Eq. 3 by assuming that the responses vary slowly with redshift compared to σ²(z₁, z₂). Now this approximation can be refined by cutting the integrals into small pieces and performing the approximation in each piece. Concretely, we start by defining a series of functions Φ_n(z) that sum to unity over a redshift range suitable for the integrals for all considered observables,

(C.3)

An example of such Φ_n is series of non-intersecting top-hat functions that cover [z_min, z_max]. Gaussian-like functions are another possibility. They have the advantage of avoiding ringing effects in later integrals.

The LSS observable then takes the form

(C.4)

and its super-sample covariance is

where we can apply the Sij approximation to Cov_SSC[𝒪(i,α)_m,𝒪(j,β)_n]:

(C.5)

(C.6)

with

(C.7)

and

(C.8)

where we assumed that the Φ_n functions have sufficiently narrow support for probe-specific kernels to be considered constant. In practical tests, we found that a width Δz ∼ 0.1 is small enough. The SSC matrix thus has the same form as Eq. 5,

(C.9)

where now U is a dense rectangular matrix with the shape (n_o, n_Φ),

(C.10)

Thus all the machinery defined in the main text of the article can be used here, as well; Eq. 7 for the inverse covariance, Eq. 19 for the likelihood, and Eq. 24 and Eq. 30 for the Fisher and inverse Fisher matrices.

The key difference is that the S_m, n matrix is now independent of the considered physical probe; indeed, it is the covariance matrix of probe-independent δ_b parameters,

(C.11)

In the likelihood approach to SSC, this δ_b, n is the vector of nuisance parameters over which is to be marginalised. The difference with the main text of the article is that this set of parameters is now independent of the physical probe and does not grow when probes are combined. The number of these parameters is defined by the width of the Φ_n functions, and optimisation can be studied to keep this width as large as possible while maintaining the validity of the Sij approximation.

According to the results of Lacasa & Grain (2019), a width Δz = 0.1 is already small enough to maintain the validity of the approximation. For a stage IV galaxy survey out to z_max = 2, 20 parameters or even fewer are therefore enough to analyse all physical probes. This compares favourably to the 55 δ_b parameters of the WL-only case in the main text or to the 210 parameters of the 3x2pt case. Moreover, this number is small enough not to be a limiting factor in an MCMC analysis, especially in the context of a combined probe analysis where the number of other nuisance parameters is already high. The only downside to this approach is that the pipelines need to be modified to output all the ∂_δb𝒪(i, α)_m, an output with a size 20 times the size of the observable prediction. We note, however, that this needs to be done in general only for the fiducial cosmology for which the covariance is to be computed.

All Tables

All Figures

Fig. 1. Information on δb: Raw (I), marginalised (I − δI), and prior (1/S). Top: information in each bin of redshift for WL analyses of each bin separately. Bottom: information in each pair of redshift bins for a joint tomographic WL analysis. The 55 pairs are ordered as (0, 0)(0, 1)...(0, 9)(1, 1)..(1, 9)...(8, 8)(8, 9)(9, 9). The light vertical lines indicate the jump of the first redshift bin, i.e., the leftmost line is at (0, 9) before the jump to (1, 1).

In the text

	Fig. 2. Information on δb in each redshift bin pair in a tomographic 3×2pt analysis (WL + GC + XC). There are 210 pairs: 55 for WL, 100 for XC, and 55 for GC. The pairs are ordered as in Fig. 1. The only difference is that XC has more pairs because there is no symmetry (i, j)↔(j, i).
In the text

Fig. 3. Forecasted errors on the dark energy equation of state, σ8 and Ωm for a tomographic WL analysis. Top two rows: constraints as a function of A0 (which changes the amplitude of the SSC response). Bottom two rows: constraints as a function of A1 (which changes the scale dependence of the SSC response). Dotted lines show the value for the fiducial covariance, (A0, A1) = (1, 0), and the dashed lines show ±10% around the dotted line.

In the text

Fig. 4. Forecasted errors on the dark energy equation of states, σ8, and Ωm for a 3×2pt analysis. Top two rows: constraints as a function of A0, which changes the amplitude of the SSC response. Bottom two rows: constraints as a function of A1, which changes the scale dependence of the SSC response. Dotted lines show the value for the fiducial covariance, and the dashed lines are ±10% around the dotted line.

In the text

	Fig. 5. For a 3×2pt analysis, forecast errors on dark energy, σ8, and Ωm. Top: as a function of ϵWL, which changes the WL response alone. Centre: as a function of ϵXC, which changes the cross-correlation response alone. Bottom: as a function of ϵGC, which changes the GC response alone. Dotted line: Value for the fiducial covariance. Dashed lines: ±10% around the dotted line.
In the text