© ESO 2018

1 Introduction

The large scale structure (LSS) of the Universe results from the growth of local density perturbations induced by gravitational collapse within an expanding background (e.g. Peebles 1980). The multi-point correlations and associated spectra that characterise this structure can be measured in real data and combined with theory predictions to constrain cosmological models, including gravity theories and relative amounts of dark matter and dark energy (Astier et al. 2006; Percival et al. 2010; Riess et al. 2009; Hildebrandt et al. 2017; Beutler et al. 2017; Hinton et al. 2017; DES Collaboration 2017). However, likelihood methods of parameter inference also require reliable estimates of the covariances associated with structure observables. Estimates of the covariances can be obtained from the actual data via methods such as jackknife and bootstrap, from simulations and Monte Carlo realisations, and from theoretical predictions that account for the known observational effects (Dodelson & Schneider 2013; Giannantonio et al. 2016; Crocce et al. 2016; O’Connell et al. 2016; Singh et al. 2017; Shirasaki et al. 2017; Blot et al. 2016; Escoffier et al. 2016; Pearson & Samushia 2016; Camacho et al. in prep.). In this work we focus on the latter approach.

Within the so-called halo model paradigm (Cooray & Sheth 2002), the LSS can be characterised by the statistical properties of the Universe building blocks, i.e. dark matter haloes. These haloes are characterised by their abundance, bias, and profiles, all of which can be studied from dark matter N-body simulations and also from high-quality data sets. As galaxy clusters develop within dark matter haloes, they trace the highest density peaks. Their number counts and covariances are very sensitive probes of structure growth and expansion of the Universe (Lima & Hu 2004, 2005, 2007; Schmidt et al. 2009; Aguena & Lima 2016). A number of past and current surveys have detected clusters in multiple wavelengths with hopes to use such detections for cosmological purposes (Miller et al. 2005; Koester et al. 2007; Soares-Santos et al. 2011; Dietrich et al. 2014; Rykoff et al. 2014, 2016; Bleem et al. 2015; Planck Collaboration XXVII 2016; Bayliss et al. 2016).

As we consider scales close to the survey maximum size, the number of modes available decreases significantly, which represents an intrinsic source of uncertainty for the inferred structure properties. The cosmic variance quantifies these uncertainties and may have contributions from scales much larger than the survey itself, in which case we refer to it as super-sample covariance (SSC), whose effects have been studied recently in multiple contexts (e.g. Takada & Hu 2013; Li et al. 2014a,b, 2016; Takahashi et al. 2014; Shirasaki et al. 2017; Hu et al. 2016). SSC has been shown to be particularly important for probe combinations as it couples observables together (Takada & Bridle 2007; Takada & Spergel 2014; Krause & Eifler 2017; Lacasa & Rosenfeld 2016). Within a local patch, large scale modes change the effective average density, which can differ significantly and unpredictably from the true background density, affecting the inferred correlations. In fact, SSC may be the dominant source of errors in Jackknife covariance estimations (Shirasaki et al. 2017).

Finally, the survey geometry or footprint and its selection function, which is characterised by, for example masks and depth maps, also affect the estimation of correlations and covariances (Takahashi et al. 2014). Not only the geometry and selection must be known to good precision, but their properties must be properly propagated into the measured and predicted correlations and covariances. Some predictions can only be directly made under certain approximations, for example for full sky calculations, while for the more realistic case of partial sky coverage, further complicating assumptionsneed to be made. As we attempt to extract maximum information from observations, we may end up with a complicated survey mask containing holes, for example owing to saturated pixels, asteroids, contamination from stars, or simply pixels that do not satisfy a depth criterion. All these effects must be accounted for in a proper cosmological analysis.

In this article, we study the effect of partial sky coverage and arbitrary masks in the estimation of the SSC of cluster counts. We propose a method to estimate this covariance efficiently for an arbitrary mask, and compare this approach to a number of approximated calculations in the literature. We show that our general calculation reduces to the approximated computations in the appropriate limits, but differs from such computations in general. Although beyond the scope of this work, we expect our theoretical estimation can be compared to other methods of estimating the covariance, which also attempt to account for the survey geometry and mask effects such as with simulations. One of the advantages of our method, however, is that it allows for the covariance to be computed as a function of cosmological and nuisance parameters at each step of a Monte Carlo Markov chain, allowing, for example for self-calibration of cluster observable-mass distribution in general cosmological analyses of cluster samples (Lima & Hu 2005). It is also unbiased, contrary to internal covariance estimators (Lacasa & Kunz 2017). Finally, it is numerically much cheaper than running thousands of N-body simulations, the latter needing to be much larger than the survey to capture super-survey modes efficiently.

This article is organised as follows. In Sect. 2.2, we introduce the cluster counts covariance and the formalism for an exact SSC computation, and we then consider the case of full sky and flat sky (small angles) limit, providing some comparisonsfor the covariance kernels between those cases. In Sect. 3 we consider the flat sky limit under various approximations that have been proposed in the literature for computing the SCC and compare the results to the exact computation. In Sect. 4, we propose a method to numerically compute SSC in the case of partial sky coverage with an arbitrary survey mask. In Sect. 5 we present our results, applying the method to the case of a geometry similar to that of the Dark Energy Survey (DES), discussing the mask effectsand recovering the flat sky limit. Finally, in Sect. 6 we present our conclusions and perspectives.

In all numerical computations, we take a cosmology for a flat Λ CDM Universe with parameter values h = 0.67, Ω_b h² = 0.022, Ω_c h² = 0.12, w = −1, n_S = 0.96, σ₈ = 0.83. Cluster counts are computed in two bins of redshift in the range z ∈ [0.4, 0.6] with Δz = 0.1, and four bins of mass in the range log[M∕(h⁻¹M_⊙)] ∈ [14, 16] with Δlog[M∕(h⁻¹M_⊙)] = 0.5. The halo mass function is from a fit to simulations from Tinker et al. (2008), the halo bias is from Tinker et al. (2010), and the linear matter power spectrum is from the transfer function by Eisenstein & Hu (1998). We set the following notational conventions: r(z) is the comoving distance, that is dr = c dz∕H(z) with H(z) the Hubble parameter, and dV = r²dr is the comoving volume element per steradian. We use short cuts such as d X₁₂ = dX₁ dX₂ and P_m (k|z₁₂) = P_m(k|z₁, z₂) = G(z₁)G(z₂)P_m(k) with G(z) the linear growth function. Often the limits of redshift or mass integrals are implicitly those of the redshift or mass bin considered.

2 Covariance of large scale structure observables and super-sample covariance

2.1 Covariance of large scale structure observables

LSS observables have a covariance that can be composed of many different terms (e.g. Lacasa & Rosenfeld 2016; Lacasa 2018). The particular case of interest to this article is that of cluster counts, which can be viewed as a halo 1-point correlation function summed over the survey area.

We address the covariance as given by a decomposition in two specific terms: 1-halo and 2-halo (1)

These termsare denoted respectively as shot-noise variance and sample (co)variance in the literature (e.g. Hu & Kravtsov 2003).

The 1-halo or shot-noise variance term is given by, (2)

where Ω_S is the solid angle covered by the survey. In the following discussion, this term is accounted for when showing the total counts covariance. Its dependence on the survey footprint is a trivial scaling by Ω_S ¹ and does not require any particular formalism to be predicted for a survey with large and/or complicated coverage.

We refer to the 2-halo term as SSC, which is described in more detail in Sect. 2.2.

For observables other than cluster number counts (for instance the galaxy angular power spectrum; see e.g. Lacasa & Rosenfeld 2016; Lacasa 2018) the covariance may contain many more terms, whose dependence on the survey geometry may be far from trivial. In this article, we restrict ourselves to studying the geometry dependence of SSC, as needed for cluster number counts analyses.

2.2 Super-sample covariance

SSC is a source of uncertainties for LSS observables coming from modes of size larger than the survey. The effect of these large scale modes on structure comes from the fact that the effective or local matter density averaged within the survey can be different from the true background density averaged over the whole Universe or an ensemble of Universes written as (3)

where δ_b(z) is a background density perturbation induced by large scale modes. Equivalently, defining the density contrasts the ensemble average is ⟨δ_m(x)⟩ = 0 by definition, however the spatial average over the survey is . (4)

We split the survey window function W(x) into its radial W_r(z) and angular pieces, i.e. . The radial part simply specifies the redshift binning and we do not explicitly state this in our description below anymore because it is effectively included with implicit bounds on redshift integrals. On the other hand, the angular window is the main object whose effect on the covariances we want to consider. In principle may depend on redshift as well. In fact, this is very often the case for surveys with significant depth variations across the sky. For simplicity, we keep the angular window only as a function of the angular vector, although the formalism can be easily generalised (see Appendix A).

If the survey angular window subtends a solid angle Ω_S in the sky, and we denote the position vector , the background perturbation δ_b is given by (5)

As a result, all LSS observables respond to such change of background density, becoming correlated. Explicitly, the SSC term of the cross-covariance between twoobservables and is given by (e.g. Lacasa & Rosenfeld 2016) (6)

where 𝔬_i is the density of observable (per comoving volume dV), ∂𝔬_i ∕∂δ_b is its reaction to the change in background δ_b, and σ²(z₁, z₂) is the covariance of δ_b, i.e. (7)

where is the Fourier transform of the survey window function.

In this article, we are mostly interested in the case of cluster number counts N_cl (i_M, i_z) per steradian within bins of mass i_M and bins of redshift i_z ², (8)

where dn_h∕dM is the halo mass function. The SSC of these counts is traditionally denoted simply as sample (co)variance in the cluster literature (Hu & Kravtsov 2003). It is important to notice though that super-horizon modes affect even full sky surveys and must always be accounted for. The SSC for clusters counts is given by (e.g. Lacasa & Rosenfeld 2016) (9)

where the response of the halo number density n_h = dN_cl∕dV to a change of background density is given by the number-weighted halo-bias averaged within the bin (e.g. Schmidt et al. 2013) as follows: (10)

where b(M, z) is the (first order) halo bias.

Most of the elements for the computation of the covariance pose no particular numerical problem. A notable exception is σ² (z₁, z₂), as it depends on the survey geometry. Numerically tractable formulae are known only for a few cases with simple geometries. First, for the case of full sky, we have (Lacasa & Rosenfeld 2016) (11)

where r_i = r(z_i). Second, an expression for σ²(z₁, z₂) is also known in the flat sky limit (e.g. Hu & Kravtsov 2003; Lima & Hu 2007; Aguena & Lima 2016) (12)

for a cylindrical window function of radius θ_S delineating a survey solid angle , and the wave-vector k = (k_∥, k_⊥) is split into its components parallel and perpendicular to the line of sight, where .

Figure 1 shows the comparison of σ² resulting from the two expressions above for an angular window of radius θ_S = 5 deg. The full sky covariance was rescaled by a factor , where f_SKY = Ω_S∕4π and for θ_S translated in radians.

We see that both covariances share the behaviour of peaking at z₁ = z₂ and decreasing to zero as |z₁ − z₂| increases³. Regarding the amplitudes, we see that rescaling the full sky covariance by the usual 1∕f_SKY factor underpredicts the covariance by a factor ~ 3.4 in this case⁴. Even when rescaling the covariances to the same peak amplitude by hand, we find that their shapes are broadly similar but differ in details. The full sky covariance is more strongly peaked at the centre, then decreases and crosses zero, reaching a negative minima of height ~ 7% of its peak, beforeslowly asymptoting zero. By contrast, the flat sky covariance is a bit broader in its positive part, but its negative minimum is only ~ 0.8% of its peak, and it asymptotes faster to zero. This clearly indicates that SSC is non-trivially related to sky coverage, and a more general approach is required for its accurate computation.

3 Comparison of super-sample covariance approximations

The full Eq. (9) for SSC can be numerically expensive, as it requires a double redshift integral. Furthermore, for each pair of redshift, a double (k_⊥, k_∥) integral is required (e.g. in the flat sky case). To tame down this burden, several approximations have been devised in the literature. Below we consider two approximations in the flat sky regime.

First, in Eq. (9) we can assume that ∂n_h∕∂δ_b varies slowly with redshift within the bins, such that it can be approximated by its bin-averaged value and taken out of the integral. In fact this may be a good approximation for sufficiently narrow bins. In this case, Eq. (9) takes the form (Hu & Kravtsov 2003) (13)

where (14)

and the normalised number-weighted halo bias within the bin is (15)

This approximation implicitly assumes that b(M, z) varies slowlywith redshift within the bin i_z, compared toσ²(z₁, z₂). For a cylindrical window function of height δr in the flat sky case, the sample variance matrix takes the form (Hu & Kravtsov 2003; Lima & Hu 2007; Aguena & Lima 2016) (16)

where the power spectrum P_m(k|z₁₂) is evaluated at the centre of the respective redshift bins. A nice feature of this approximation is that it removes the need to compute a double redshift integral. Furthermore the double integral over (k_⊥, k_∥) only needs to be computed times (where n_z is the number of redshift bins), instead of all redshift pairs (z₁, z₂) needed to compute the redshift integral otherwise. This speeds up considerably the computation of SSC. We call this approximation the Sij method.

Second, another approximation used (e.g. Krause & Eifler 2017) is that σ² (z₁, z₂) is a Dirac delta function at z₁ = z₂, so that the double redshift integral collapses into a single integral (17)

where, following notation from Krause & Eifler (2017), we have (18)

where P_m(k_⊥, z) = P_m(k = k_⊥, z), and we note that the k_∥ integral has disappeared and that σ_b has units of Mpc/h.

The assumption behind this approximation is that the 3D window function W(x) is much wider in the radial direction than in the transverse direction. Thus super-survey modes have k_∥ ≪ k_⊥, and P_m (k) can be taken as approximately constant within the k_∥ integral of Eq. (12). We thus expect the approximation to fare better for wide redshift bins and small angles. By limiting the computation to equal redshift bins i_z = j_z, and reducing the multiple integrals to a single redshift and a single wave-vector k_⊥ integral, the approximation speeds up the SSC computation considerably. Hereafter, we call this approximation the KE method.

We implemented all three SSC methods numerically: (1) full computation from Eqs. (9) and (12), (2) Sij approximation, and (3) KE approximation. Figure 2 shows the ratio of the auto-z covariances to the full computation, plotted as a function of ; i.e. the first four points show the covariance ratio for i_M = i_z = 0 (14 < logM < 14.5 and 0.4 < z < 0.5) and varying j_M, the next four points are for i_M = 1, i_z = 0, etc. This order is simply convenient to plotting this multi-variable function. As visible in Fig. 2 both the Sij and KE approximations underpredict the amplitude of the covariance by 30–35% and by ~15%, respectively. The fact that the Sij approximation underpredicts the covariance conforms to our intuition; indeed, by taking values at the centre of the redshift bins, the approximation neglects the fact that the mass function and matter power spectrum both grow with time, which makes their integral larger⁵. We checked that the Sij method and full computation are in much better agreement for smaller redshift bins, and indeed it can be seen analytically that they reduce to each other in the limit Δ z → 0.

In Fig. 3, we show the correlation matrices for cluster counts obtained from each of the three SSC computations. The data points are ordered with increasing mass then increasing redshift, i.e. the number count index is i_M + n_M i_z⁶.

As visible in Fig. 3, for cross-redshift covariances the full computation gives a positive correlation between our two redshift bins, reaching ~ 13% at maximum for Δz = 0.1. On the other hand, the Sij approximation predicts anti-correlation between the bins, reaching at maximum − 16%, while the KE approximation by construction predicts zero correlation. We checked that the Sij method performs better as the redshift bin width decreases. For instanceat Δz = 0.01, the agreement with the full computation is very good.

As a conclusion of this section, we see that the full computation is necessary to faithfully predict the SSC covariance matrix in the case Δz = 0.1. In fact, thisis a representative binning for current photometric galaxy surveys whose present photo-z errors are ofthis order (e.g. Sánchez et al. 2014).

Fig. 1 Comparison of σ2(z1, z2) for z1 = 0.5 in different cases. In blue the flat sky formula Eq. (12) for a survey angular radius θS = 5 deg. In green the full sky formula Eq. (11) rescaled by a factor 1∕fSKY.

Fig. 2 Ratio of the SSC covariances from the Sij approximation and KE approximation to the full computation. Auto-z covariances are ordered as a function of (see text fordetails).

Fig. 3 Comparison of the cluster counts correlation matrix for different SSC computations, for a survey with angular radius θS = 5 deg. The SSC matrix is shown in 2 blocks for the redshift bins, and each block has 4 entries for the logarithmic mass bins. From left to right: full numerical computation from Eq. (9), Sij approximation, and KE approximation.

4 Super sample covariance for partial sky coverage

In this section, we assume for simplicity that the survey mask is independent of redshift, as we in fact already did implicitly in Eq. (5). It is however straightforward to generalise the formalism to a redshift-dependent mask (see Appendix A).

4.1 Formal derivation

We want to compute the covariance of the background mode (19)

Given that the mask has zero value outside the survey, and given the normalisation of spherical harmonics such that , the background mode of the survey is related to the monopole of the masked matter density as (20)

such that we have (21)

where is the angular power spectrum at ℓ = 0 of the masked matter density field. From C_ℓ pseudo-spectrum methods (e.g. Hivon et al. 2002), we know that the masked power spectrum is related to the true spectrum via the so-called coupling matrix, such that (22)

where (hereafter simply ) is the angular spectrum of matter density in an infinitesimal redshift shell. The coupling matrix is given by (23)

where C_ℓ(W) is the angular power spectrum of the angular window . In our case, given the properties of Wigner symbols, this simplifies to (24)

Furthermore we have and thus (25)

where Ω_S = 4π f_SKY and the matter angular power spectrum is given by (e.g. Campagne et al. 2017) (26)

The combined Eqs. (25) and (26) represent the main analytical result of this article.

4.2 Limiting cases and remarks

It is interesting to examine a few limiting cases of Eqs. (25) and (26). In the full sky case, we have C_ℓ (W) = 4π δ_0,ℓ and f_SKY = 1, and therefore thus we have (27)

recovering indeed Eq. (11).

In partial sky, , such that (28)

i.e. the full sky covariance term is the first in the sum contributing to the partial-sky covariance.

If is scale independent, i.e. independently of ℓ (e.g. for P_m(k|z₁₂) = const., see Sect. 5.2), we have (29)

and therefore , i.e. we obtain the usual f_SKY approximation to partial sky covariance. Conversely the reverse is also true: if the f_SKY approximation holds for any mask, then is constant. Given that for viable cosmological models is not constant, this shows that SSC is a non-trivial source of covariance that cannot be treated by classical approximations.

Now we notice that Eq. (25) can be rewritten as (30)

with (31)

such that in full sky . Applying this to the SSC covariance of two observables, for example number counts, we have (32)

where (33)

Therefore the combined Eqs. (30)–(33), along with the window spectrum C_ℓ (W), allow for the computation of the SSC in partial sky for a general window and binning choice.

Numerically, we can pre-compute and tabulate , and then just change C_ℓ(W) as the mask varies. This enables easier studies of the mask effect such as optimisations of survey strategy (Takahashi et al. 2014), forecasts for improvements as a survey area grows, and comparisons between different surveys. Conversely, for a fixed survey with a well-defined mask, we can pre-compute C_ℓ (W) for the specific survey mask, and change as a function of cosmology within likelihood cosmological analyses, taking full account of geometry, mask, selection, and cosmological dependencies in the covariances, and therefore deriving reliable parameter uncertainties.

4.3 Implementation

When estimating the SSC from Eq. (32), one numerical difficulty is the evaluation of the integral in Eq. (31), given that the Bessel functions highly oscillate with slow damping as kr_i →∞.

For the first multipoles, we may express the Bessel functions in terms of sine and cosine. Through trigonometrical identities, the integralswith products of Bessel function can thus be expressed as a sum of Fourier transforms. Derivations and expressions for the first three multipoles are given in Appendix C. For illustrative purposes we give below the expression for ℓ = 0 (34)

where (35)

is a continuous cosine transform that can be efficiently approximated numerically with a discrete fast Fourier transform (FFT). However, as argued in Appendix C, this method becomes too cumbersome at high ℓ and may also become numerically unstable. We could carry it out only for ℓ = 0, 1, 2.

As an alternative, we decided not to evaluate σ²(z₁, z₂), but instead to switch the order of the integrals over k and z, to compute the covariance directly as (36)

where the kernel is given by (37)

As can be seen in Fig. 4, this redshift integral effectively damps out the Bessel oscillations on scales k > k_peak + 2π∕Δr(i_z), where Δ r(i_z) = r(z_max) − r(z_min) is the width of the redshift bin in terms of comoving distance. This damping makes the k integral Eq. (36) much easier to deal with numerically compared to Eq. (26), since the integrand support is now more compact. So numerically we could just compute the later integrals by brute force, and this is the method we use for all numerical results shown hereafter in this article. We note that this method does not compute σ² (z₁, z₂) as an intermediary product, and thus prevents from comparison of this quantity in partial sky with its alter ego in the full sky or flat sky limits, as in Fig. 1.

We remark that alternatively, Eq. (37) can be expressed as a Hankel transform (38)

via the inversion of the r(z) relation and defining a radial bin i_r from i_z. Equation (38) can be efficiently evaluated using Fourier transform methods such as FFTLog (e.g. Hamilton 2000; Camachoet al. in prep.). Although this is not the numerical method we used in this article, we note it can be a useful approach for future high-precision applications.

The results presented in this article consider the case of cluster counts, and the equations can be straightforwardly generalised to other probes of the LSS; see Appendix B.

Fig. 4 as a functionof k, for the representative case iM = iz = 0 (log M = 14–14.5, z = 0.4–0.5). Left: at the lowest multipole ℓ = 0. Centre: at ℓ = 10. Right: at ℓ = 100.

5 Results

First, as a consistency test, we checked against the full sky covariance matrix computed via given in Eq. (11) (e.g. Lacasa & Rosenfeld 2016). We find good agreement, to 0.8% precision on the auto-redshift covariance, and 7% precision on the cross-redshift covariance.

Second, we show in Fig. 5 the general results for in two representative cases: auto-covariance of the same redshift bin and cross-correlation (i.e. cross-covariance normalised by the corresponding auto-covariances) between two redshift bins. Both plots are for the lowest mass bin (i_M = j_M = 0, corresponding to logM = 14–14.5), although the shape of the curves does not change significantly when taking other mass bins (even with i_M ≠ j_M), only their amplitude changes.

For theauto-covariance (left plot), we see that it first rises with ℓ, to a maximum corresponding to the angular scale of the matter-radiation equality (i.e. the peak of P_m (k)) at that redshift, and then decreases monotonically. This behaviour conforms to our expectations, since the scale dependence is that of a projection of the 3D power spectrum P_m(k). This means that is not constantwith ℓ, and thus, as already discussed in Sect. 4, it means that the f_SKY approximation does not hold for SSC.

For the cross-covariance (right plot), the situation is interesting as we see that the covariance is first negative, increases, crosses zero, and reaches a maximum towards ℓ = 25, and then decreases asymptotically to zero. This means that small surveys have a robustly negligible covariance between redshift bins (at least for this bin width Δz = 0.1), and can thus use block-diagonal covariance matrices to speed up their likelihoods estimations. However for surveys with a large sky coverage, the cross-covariance is non-negligible and depends on the survey area and shape, becoming either positive or negative depending on the mask. In those cases, careful estimation of the SSC is thus critical.

Fig. 5 as a functionof ℓ, in the representative case iM = jM = 0 (log M = 14–14.5). Top: same redshift iz = jz = 0 (z = 0.4–0.5). Bottom: cross-redshift iz = 1 × jz = 2 (z = [0.4–0.5] × [0.5–0.6]).

5.1 Application to a realistic survey mask

In order to illustrate the SSC method in a realistic case, we used two Healpix (Górski et al. 2005) masks visible in Fig. 6. The first mask is binary and is broadly similar to the footprint of DES⁷,⁸, although wewarn that it does not in any way represent the actual DES survey area, and we do not attempt to draw any particular conclusion for that survey. The second mask represents a more pessimistic case, where we considered that 15% of the survey area had to be discarded due to bright stars, satellite trails, or other systematics. In order to simulate this effect, we simply upgraded the mask to high resolution (nside = 4096), poked random holes in it, then degraded the mask back to the original resolution (nside = 1024).

The angular power spectra C_ℓ(W) for these two masks can be seen on Fig. 7.

At low multipoles, the two spectra differ by a constant multiplicative factor 0.85², which is simply the ratio of the respective sky coverages. A second component in the spectrum of the second mask appears on smaller scales, which is a constant shot noise due to the random holes. However we see that this component is very subdominant, hence we can already expect that the only difference between the covariance of the two masks is due to the different sky coverage.

When implementing the sum over multipoles of Eq. (32), we find that we reached 1% convergence already at ℓ_max = 50. This means that need only be computed on a small number of multipoles for current and future surveys with large sky coverage, rendering the partial sky formalism developed here even more computationally efficient. Comparing the SSC covariances of the two masks, we find that they indeed only differ due to the different sky fraction, but when renormalised by f_SKY they are identicalto numerical precision (0.01% in our case). The total correlation matrix (including shot-noise) for the first mask is shown in Fig. 8.

We see that the SSC is an important contribution to the covariance matrix, dominating the error bars at log M ≤ 14.5. Concerning cross-redshifts, the SSC yields an anti-correlation between the adjacent bins reaching up to −10%. The cross-covariance is thus far from negligible for such a large survey.

The formalism presented in Sect. 4 is thus perfectly adapted to the numerical prediction of SSC, even in the case of a complex survey geometry. In fact we see that such computation is indeed necessary to reproduce the non-trivial behaviour of SSC, yielding, for example non-negligible anti-correlation of redshift bins, in the case presented above.

Fig. 6 Masks used in the analysis. Top: simple footprint, assumed observed uniformly. Bottom: same but with simulated 15% rejection of observations due to systematics (see text for details).

Fig. 7 Angular power spectra of the two masks used in the analysis and shown in Fig. 6.

Fig. 8 Total correlation matrix (SSC + shot-noise) of cluster counts for the first mask. The matrix is organised, as in Fig. 3, in two redshift blocks of increasing mass.

5.2 Flat sky limit

One remaining question is that of the link between the partial sky approach developed in Sect. 4 and the flat sky approximation used previously in the literature and shown in Eq. (12). The two approaches cannot be compared or related at the level of the formalism because the (k_⊥, k_∥) splitting does not apply for sufficiently large angles. However, we can compare what covariance the two formalisms predict in some limits.

First, we can compare the equations analytically in the case of a constant spectrum P_m (k|z₁₂) = cst ≡ 1. In the partial sky formalism we have (39)

and (40)

is independent of ℓ. Thus the SSC is given by (41)

where we made the counts N_cl implicit for the sake of clarity.

Now in the flat sky formalism we have (42)

where we changed variables x = k_⊥θ_Sr₁ and recall . Inserting this into Eq. (9), we again obtain Eq. (41). Therefore the two approaches indeed agree in the flat sky limit for a constant spectrum.

Second, we can compare the results numerically for a mask with small enough sky coverage. To do so, we created a polar cap mask of radius 5 deg⁹. The power spectrum of this mask is shown in Fig. 9.

Compared to Fig. 7, we see that power extends to smaller angular scales or higher multipoles. In this case, we found that we had to extend the sum in Eq. (32) to higher multipoles, ℓ_max = 250 to reach 1% level convergence of the covariance prediction. This stays however numerically tractable, through the observation that the ℓ dependence of is very smooth, especially after ℓ_peak = 25. Thus we can sample this dependence only for a small number of logarithmically spaced multipoles and interpolate when computing the sum in Eq. (32).

We found good agreement between the covariance derived from the partial sky formalism and that derived from the flat sky limit Eq. (12). This is visible, for example in Fig. 10 showing both auto-z covariances as a function of , i.e. the same ordering as in Fig. 2.

The partial sky formalism thus successfully recovers the flat sky approximation in the flat sky limit. Numerically, we even found the partial sky formalism to be faster than the full computation of Eq. (12), and only three times slower than the Sij approximation Eq. (13).

Fig. 9 Angular power spectrum of the 5 deg polar cap mask.

Fig. 10 Comparison of the SSC covariances derived from the partial sky formalism and the flat sky formula for a 5 deg circular sky patch. Auto-z covariances are shown as a function of (see text fordetails).

6 Conclusion

SSC, also often called sample variance, is the dominant error for cluster counts at low cluster masses. For instance, Hu & Cohn (2006) have shown that even for a survey radius θ = 2.4 deg, the sample covariance is of the same order of magnitude as shot noise for cluster counts above log M = 14.2, although shot noise dominates for a threshold logM = 14.4. As shot-noise decreases faster with survey area than SSC, careful predictions of SSC become crucial for current and future surveys covering ever larger sky areas. It becomes even more crucial as these surveys are able to probe lower cluster masses through a higher density of galaxy detections. SSC is also crucial to the analysis of galaxy clustering and lensing shear, where it dominates statistical errors on small scales, and for probe combinations as it has been shown to couple probes (Takada & Bridle 2007; Takada & Spergel 2014; Lacasa & Rosenfeld 2016).

In the case of cluster counts covariance, we examined theoretical SSC computation methods in the flat sky limit, comparing two analytical approximations proposed in the literature to a full computation. We found that both approximations underpredict the auto-z covariance by 15% to 30–35%.

We then presented a harmonic expansion method for efficiently and accurately computing SSC for an arbitrary survey window function. We developed the method in the case of cluster counts, but it can be straightforwardly generalised to other probes such as galaxy clustering or lensing shear. Our derived expressions generalise previous full sky and flat sky equations found in the literature, properly reducing to these equations in the corresponding limits. We have cast the final covariance expression from the partial sky formalism in a way that allows easy modification of the survey mask. This is particularly suitable for comparison of surveys, design of survey strategy, and tests of data cuts due to quality selection criteria or systematics.

When applying our partial sky formalism to a mask broadly similar to the DES footprint with f_SKY ~ 10%, we found a ~−15% cross-z covariance, meaning that the observables in the two redshift bins considered are anti-correlated. Hence the covariance matrix cannot be taken as block diagonal, as is the case for the approximation by Krause & Eifler (2017), which is however restricted to the flat sky limit. We also examined the possibility that the survey area is further reduced by pixel removal due to bright stars or systematics, for example. We found that this does not have important effects on the structure of the SSC matrix, only rescaling its amplitude by the effective survey area.

The results presented in this article thus render possible the theoretical computation of LSS covariances that account for selection and mask effect and also vary as a function of model parameters, as is the case in CMB analyses, for example. The latter parameter dependence is important in the case of likelihood analysis of cluster constraints, as it allows for self-calibration of the cluster observable-mass relation (Lima & Hu 2005; Hu & Cohn 2006; Baxter et al. 2016). For general probes, it can also improve cosmological parameter constraints from likelihood inference analyses compared to methods that either neglect these effects or fix the covariance from data or simulations, thereby avoiding the risk of fixing the covariance at a potentially incorrect cosmology.

Acknowledgements

We thank Flavia Sobreira for providing us with a mask of a DES-like survey footprint. We acknowledge the use of the Healpix package by Górski et al. (2005). F. L. acknowledges support by the Swiss National Science Foundation. M. L. is partially supported by FAPESP and CNPq. M. A. is supported by FAPESP.

Appendix A: Redshift-dependent mask

In the case where the survey angular mask W depends on redshift (e.g. when there are significant depth variations in the sky), the definition of the background mode Eq. (5) is simply changed to (A.1)

where now the survey solid angle Ω_S depends on redshift.

Then Eq. (25) for σ²(z₁, z₂) is changed to (A.2)

where C_ℓ[W(z₁), W(z₂)] is the angular cross-spectrum between the two masks at the two redshifts. This equation can then be integrated over (z₁, z₂) through Eq. (6) to yield the SSC covariance.

We note that in this redshift-dependent case, we can no longer permute the k and z integrals as was done in Sect. 4.3, which yielded a numerically efficient method. Recent advances by Campagne et al. (2017) may however render Eq. (A.2) numerically computable through the AngPow software.

Appendix B: Super-sample covariance for other probes

The main results of this article were derived for the covariance of cluster number counts. However, it is straightforward to generalise the equations to other LSS probes, such as galaxy clustering and lensing shear. In those cases, the covariance equations are simpler for the angular power spectrum C_ℓ than for the angular correlation function. For instance, using Limber’s approximation, Lacasa & Rosenfeld (2016) have shown that the SSC of the galaxy angular power spectrum is given by (B.1)

where .

For lensing shear, we would get a similar equation, replacing the galaxy number density by the lensing selection function and the 3D galaxy power spectrum by the matter power spectrum. The equation can thus be generalised to (B.2)

where the α index refers to either lensing or galaxy, and W_α is the corresponding selection function. In that case, it is straightforward to generalise Eq. (36), even including the possibility of the cross-covariance between galaxy and shear, (B.3)

where (B.4)

and the cross-covariance between cluster counts and either galaxy or shear, (B.5)

which gives us all the equations needed to compute the auto and cross-covariances of cluster counts, galaxy angular power spectrum, and lensing shear power spectrum, i.e. the three main cosmological probes of current and future photometric galaxy surveys.

Appendix C: First multipoles

The spherical Bessel functions j_n(x) obey the recurrence relation (C.1)

such that they can be written analytically in terms of sines, cosines and polynomials, given the initial conditions

The first few spherical Bessel functions are given by

We are trying to compute the following integrals: (C.7)

For j_ℓ(x) of the form A_ℓ(x)sin(x)∕x − B_ℓ(x)cos(x)∕x, this yields (C.8)

where the four integrals are Fourier (sine or cosine) transforms and can be computed numerically with FFTs. Let us note

Then Eq. (C.8) can be rewritten as (C.13)

with A_ℓ,i = A_ℓ(kr_i). For ℓ = 0 we have A₀ = 1 and B₀ = 0, thus (C.14)

For ℓ = 1 we have A₁ = 1∕x and B₁ = 1. Therefore (C.15)

If we define , the above equations yield

For ℓ = 2 we have A₂ = 3∕x² − 1 and B₂ = 3∕x, thus

We note that A_ℓ and B_ℓ follow the same recurrence relation Eq. (C.1) as the spherical Bessel functions j_n. As such we can look at whether there is a recurrence relation that would allow us to get the analytical formula for a general . Our efforts in this direction have shown only partially fruitful and are described below. We have (C.20)

and then

Now we define the following quantities

and decompose , and onto the basis

Then Eq. (C.22) gives a first recurrence relation (C.28)

By computing and , we can see that we have the two other recurrence relations

The system is closed through the initial conditions

The following properties can easily be shown by recurrence:

With a bit more work, we can in principle solve for and ; since ,

This can be put in the form (C.40)

for

Finally, Eq. (C.40) has the following solution: (C.44)

which in principle gives the solution for and . We can insert this into Eq. (C.28) for i = 2 to get a closed recurrence relation for and solve for it. In turn, this can be inserted into Eqs. (C.29) and (C.30) for i = 4 to get a closed recurrence relation for , and so on. So this represents a procedure for solving this set of equations and get analytical expressions for . Similarly, we could derive recurrence relations for , , and . In practice this represents a daunting task, as even the formula for ℓ = 2 (Eq. (C.19)) is already cumbersome.

Moreover, we may expect that these analytical formulae are doomed to present numerical instability at high ℓ. Indeed, the expansion of the spherical Bessel functions in terms of powers of 1∕x as in Eqs. (C.4)–(C.6) is ill-advised at high ℓ, as it leads to delicate cancellations of the numerous terms, in particular for x ≤ ℓ. We thus expect that our analytical formulae for C_ℓ also leads to delicate cancellations that may be numerically unstable. Another way of seeing this is to notice that, at high ℓ, we have from Limber’s approximation that , and is thus a decreasing function of ℓ. However we saw that and thus C_ℓ always contains a term . This term hence needs to be (at least partially) cancelled by higher order terms, and this cancellation needs to be increasingly precise at high ℓ, since P_m (k) is a steep decreasing function of k.

Given these issues, in practice, we implemented this low multipole method only for ℓ = 0, 1, 2, and used these to check our results obtained from the numerical method described in Sect. 4.3.

References

All Figures

	Fig. 1 Comparison of σ2(z1, z2) for z1 = 0.5 in different cases. In blue the flat sky formula Eq. (12) for a survey angular radius θS = 5 deg. In green the full sky formula Eq. (11) rescaled by a factor 1∕fSKY.
In the text

	Fig. 2 Ratio of the SSC covariances from the Sij approximation and KE approximation to the full computation. Auto-z covariances are ordered as a function of (see text fordetails).
In the text

	Fig. 3 Comparison of the cluster counts correlation matrix for different SSC computations, for a survey with angular radius θS = 5 deg. The SSC matrix is shown in 2 blocks for the redshift bins, and each block has 4 entries for the logarithmic mass bins. From left to right: full numerical computation from Eq. (9), Sij approximation, and KE approximation.
In the text

	Fig. 4 as a functionof k, for the representative case iM = iz = 0 (log M = 14–14.5, z = 0.4–0.5). Left: at the lowest multipole ℓ = 0. Centre: at ℓ = 10. Right: at ℓ = 100.
In the text

	Fig. 5 as a functionof ℓ, in the representative case iM = jM = 0 (log M = 14–14.5). Top: same redshift iz = jz = 0 (z = 0.4–0.5). Bottom: cross-redshift iz = 1 × jz = 2 (z = [0.4–0.5] × [0.5–0.6]).
In the text

	Fig. 6 Masks used in the analysis. Top: simple footprint, assumed observed uniformly. Bottom: same but with simulated 15% rejection of observations due to systematics (see text for details).
In the text

	Fig. 7 Angular power spectra of the two masks used in the analysis and shown in Fig. 6.
In the text

	Fig. 8 Total correlation matrix (SSC + shot-noise) of cluster counts for the first mask. The matrix is organised, as in Fig. 3, in two redshift blocks of increasing mass.
In the text

	Fig. 9 Angular power spectrum of the 5 deg polar cap mask.
In the text

	Fig. 10 Comparison of the SSC covariances derived from the partial sky formalism and the flat sky formula for a 5 deg circular sky patch. Auto-z covariances are shown as a function of (see text fordetails).
In the text