Forecasting Super-Sample Covariance in Future Weak Lensing Surveys with SuperSCRAM

TL;DR

This paper introduces a formalism and a public code, SuperSCRAM, to forecast and mitigate super-sample covariance effects in future weak lensing surveys, improving parameter constraint accuracy.

Contribution

It presents a Fisher matrix-based approach and a practical implementation for estimating and reducing super-sample covariance in cosmological surveys.

Findings

Super-sample covariance can significantly inflate error volumes in parameter space.

Mitigation strategies can reduce super-sample covariance effects by more than half.

Application to WFIRST and LSST data demonstrates practical effectiveness.

Abstract

The observable universe contains density perturbations on scales larger than any finite volume survey. Perturbations on scales larger than a survey can measure degrade its power to constrain cosmological parameters. The dependence of survey observables such as the weak lensing power spectrum on these long-wavelength modes results in super-sample covariance. Accurately forecasting parameter constraints for future surveys requires accurately accounting for the super-sample effects. If super-sample covariance is in fact a major component of the survey error budget, it may be necessary to investigate mitigation strategies that constrain the specific realization of the long-wavelength modes. We present a Fisher matrix based formalism for approximating the magnitude of super-sample covariance and the effectiveness of mitigation strategies for realistic survey geometries. We implement our formalism in the public code SuperSCRAM: Super-Sample Covariance Reduction and Mitigation. We illustrate SuperSCRAM with an example application, where the modes contributing to super-sample covariance in the WFIRST weak lensing survey are constrained by the low-redshift galaxy number counts in the wider LSST footprint. We find that super-sample covariance increases the volume of the error ellipsoid in 7D cosmological parameter space by a factor of 4.5 relative to Gaussian statistical errors only, but our simple mitigation strategy more than halves the contamination, to a factor of 2.0.

Tables

Table 1. Table 1 : Results from a well converged run of SuperSCRAM with a sample WFIRST and LSST survey footprints. with k cut = 0.0814 ​ Mpc − 1 subscript 𝑘 cut 0.0814 superscript Mpc 1 k_{\text{cut}}=0.0814\;\text{Mpc}^{-1} as described in Appendix  A , with the product of eigenvalues, eigenvalue in most contaminated direction, and the coefficients describing the most contaminated directions, as described in subsection  2.3 . The convergence of Π ( i ) ​ λ ( i ) subscript Π 𝑖 superscript 𝜆 𝑖 \Pi_{(i)}\lambda^{(i)} with respect to k cut subscript 𝑘 cut k_{\text{cut}} is ≲ 𝒪 ​ ( 1 % ) less-than-or-similar-to absent 𝒪 percent 1 \lesssim\mathcal{O}(1\%) . The mitigation reduces the volume of parameter space contaminated by SSC by approximately a factor of 5. The most contaminated direction does not change drastically after mitigation.

	${\Pi }_{\left(i\right)}{\lambda }^{\left(i\right)}$	${\lambda }^{\text{top}}$	$\widehat{n_s}$	$\widehat{{\Omega }_m{h}^2}$	$\widehat{{\Omega }_b{h}^2}$	$\widehat{{\Omega }_{\text{de}}{h}^2}$	$\widehat{ln\left(A_s\right)}$	$\widehat{w_0}$	$\widehat{w_a}$
$C_{\text{SSC}}$	20.26	9.94	329.3	5051.1	-8425.4	-390.9	269.4	-206.3	-58.6
$C_{\text{mit}}$	4.02	2.93	271.0	4319.2	-7132.2	-346.2	229.2	-176.5	-48.9

Equations

$$F_{ij}^{\text{tot}}=F_{ij}^{\text{prior}}+F_{ij}^{\text{survey}},$$

$$F_{ij}^{\text{survey}}=\frac{\partial O^I}{\partial {\Theta }^i}{F}_{IJ}^{\text{tot}}\frac{\partial O^J}{\partial {\Theta }^j},$$

$$C_{\text{tot}}^{IJ}=C_{\text{g}}^{IJ}+C_{\text{ng}}^{IJ}+C_{\text{SSC}}^{IJ}.$$

$$C_{\text{SSC}}^{IJ}=\frac{\partial O^I}{\partial {\delta }^{\alpha }}(F_{\text{SSC}}^{\alpha \beta }{)}^{-1}\frac{\partial O^J}{\partial {\delta }^{\beta }}.$$

$$F_{\alpha \beta }^{\text{SSC+mit}}=F_{\alpha \beta }^{\text{SSC}}+\sum _n{F}_{\alpha \beta }^n.$$

$${\mathcal{M}}_j^i=(C_{\text{SSC}}^{ik})(C_{\text{g}}^{-1}{)}_{kj},$$

$${\tilde{\mathcal{M}}}_{\hat{i}\hat{j}}=(L^{-1}{)}_{k\hat{i}}{C}_{\text{SSC}}^{kl}(L^{-1}{)}_{l\hat{j}},$$

$${\tilde{C}}_{\text{SSC},\text{mit}}^{(i)(j)}\equiv {\tilde{v}}_k^{(i)}{C}_{\text{mit}}^{kl}{\tilde{v}}_l^{(j)},$$

$$n_1$$

$$=\frac{3}{(r_{\text{max}}^3-r_{\text{min}}^3){\Omega }_1}{\int }_{{\Omega }_1}d\Omega {\int }_{r_{\text{min}}}^{r_{\text{max}}}{r}^{2}drn(r)\left[1+b(r)\sum _{\alpha }{\delta }^{\alpha }{\psi }_{\alpha }(r,\Omega )\right],$$

$$\frac{\partial \Delta n_{12}}{\partial {\delta }^{\alpha }}$$

$$=\frac{3}{r_{\text{max}}^3-r_{\text{min}}^3}{\int }_{r_{\text{min}}}^{r_{\text{max}}}{r}^{2}drn(r)b(r)\left[\frac{1}{{\Omega }_1}{\int }_{{\Omega }_1}d\Omega {\psi }_{\alpha }(r,\Omega )-\frac{1}{{\Omega }_2}{\int }_{{\Omega }_2}d\Omega {\psi }_{\alpha }(r,\Omega )\right].$$

$$⟨n⟩=\frac{3}{(r_{\text{max}}^3-r_{\text{min}}^3)}{\int }_{r_{\text{min}}}^{r_{\text{max}}}n(r)r^{2}dr.$$

$$N_{12}\equiv \mathrm{Var}\left(n_1-n_2\right)=⟨n⟩\left(\frac{1}{V_1}+\frac{1}{V_2}\right),$$

$$F_{\alpha \beta }^{\text{mit}}=\frac{\partial \Delta n_{12}}{\partial {\delta }^{\alpha }}(N_{12}{)}^{-1}\frac{\partial \Delta n_{12}}{\partial {\delta }^{\beta }}.$$

$$F_{\alpha \beta }^{\text{mit}}=({\mathbf{V}}^{\mathrm{T}}\mathbf{KV}{)}_{\alpha \beta }=V^a{}_{\alpha }{K}_{ab}{V}^b{}_{\beta },$$

$$C_{\text{SSC+mit}}^{\alpha \beta }=[({\mathbf{C}}_{\text{SSC}}^{-1}+{\mathbf{V}}^{\mathrm{T}}\mathbf{K}\mathbf{V}{)}^{-1}{]}^{\alpha \beta }={\mathbf{C}}_{\text{SSC}}^{\alpha \beta }-[{\mathbf{C}}_{\text{SSC}}{\mathbf{V}}^{\mathrm{T}}({\mathbf{K}}^{-1}+{\mathbf{VC}}_{\text{SSC}}{\mathbf{V}}^{\mathrm{T}}{)}^{-1}{\mathbf{VC}}_{\text{SSC}}{]}^{\alpha \beta }.$$

$$W^{aI}\equiv V^a{}_{\alpha }{C}_{\mathrm{SSC}}^{\alpha \beta }\frac{\partial O^I}{\partial {\delta }^{\beta }}$$

$$C_{\text{SSC+mit}}^{IJ}$$

$${\psi }_{\alpha }(r,\theta ,\varphi )=j_{l_{\alpha }}(k_{\alpha }r)Y_{l_{\alpha }{m}_{\alpha }}^{\text{R}}(\theta ,\varphi ),$$

$$Y_{lm}^{\text{R}}=\sqrt{\frac{2(l+1)}{4\pi }}\sqrt{\frac{(l-|m|)!}{(l+|m|)!}}(-1{)}^{|m|}{P}_l^{|m|}(\cos \theta )\{\begin{array}{ll}\sqrt{2}\sin |m|\varphi & \text{if }m<0\\ 1& \text{if }m=0\\ \sqrt{2}\cos |m|\varphi & \text{if }m>0\end{array}.$$

$${\delta }^{\alpha }(k_{\alpha })$$

$$=\frac{1}{N_{\alpha }}\int \frac{d^3\mathbf{k}}{(2\pi {)}^3}\delta (\mathbf{k}){\int }_0^{r_{\text{max}}}dr{r}^2{j}_{l_{\alpha }}(k_{\alpha }r)\int d^2\hat{\mathbf{r}}{e}^{ikr\hat{\mathbf{k}}\cdot \hat{\mathbf{r}}}{Y}_{l_{\alpha }{m}_{\alpha }}^{\text{R}}(\hat{\mathbf{r}})$$

$$=\frac{4\pi i^{l_{\alpha }}}{N_{\alpha }}\int \frac{d^3\mathbf{k}}{(2\pi {)}^3}\delta (\mathbf{k})Y_{l_{\alpha }{m}_{\alpha }}^{\text{R}}(\hat{k}){\int }_0^{r_{\text{max}}}dr{r}^2{j}_{l_{\alpha }}(k_{\alpha }r)j_{l_{\alpha }}(kr),$$

$${\int }_{S^2}{d}^2\hat{\mathbf{r}}{Y}_{l_{\alpha }{m}_{\alpha }}^{\text{R}}(\hat{r})e^{i\mathbf{k}\cdot \mathbf{r}}=4\pi i^{l_{\alpha }}{j}_{l_{\alpha }}(kr)Y_{l_{\alpha }{m}_{\alpha }}^{\text{R}}(\hat{k}).$$

$$N_{\alpha }$$

$$=I_{\alpha }(k_{\alpha },r_{\text{max}}){\int }_{\Omega }{Y}_{l_{\alpha }{m}_{\alpha }}^{\text{R}}(\hat{r})Y_{l_{\alpha }{m}_{\alpha }}^{\text{R}}(\hat{r})d^2\hat{\mathbf{r}}$$

$$=I_{\alpha }(k_{\alpha },r_{\text{max}}),$$

$$\begin{array}{rl}{I}_{\alpha }(k,r_{\text{max}})& ={\int }_0^{r_{\text{max}}}dr{r}^2{j}_{l_{\alpha }}(k_{\alpha }r)j_{l_{\alpha }}(kr)\\ & =\frac{\pi }{2}\sqrt{\frac{1}{k_{\alpha }k}}{\int }_0^{r_{\text{max}}}drr{J}_{l_{\alpha }+1/2}(k_{\alpha }r)J_{l_{\alpha }+1/2}(kr)\\ & =\frac{\pi }{2}\frac{r_{\text{max}}}{\sqrt{k_{\alpha }k}}\frac{\left[k_{\alpha }{J}_{l_{\alpha }+1/2}(k{r}_{\text{max}})J_{l_{\alpha }+1/2}^{\prime }(k_{\alpha }{r}_{\text{max}})-k{J}_{l_{\alpha }+1/2}(k_{\alpha }{r}_{\text{max}})J_{l_{\alpha }+1/2}^{\prime }(k{r}_{\text{max}})\right]}{k^2-k_{\alpha }^2}\\ & =\frac{\pi }{2}\frac{r_{\text{max}}}{\sqrt{k_{\alpha }k}}\frac{\left[k_{\alpha }{J}_{l_{\alpha }+1/2}(k{r}_{\text{max}})J_{l_{\alpha }-1/2}(k_{\alpha }{r}_{\text{max}})-k{J}_{l_{\alpha }+1/2}(k_{\alpha }{r}_{\text{max}})J_{l_{\alpha }-1/2}(k{r}_{\text{max}})\right]}{k^2-k_{\alpha }^2}\\ & =\frac{\pi }{2}\frac{r_{\text{max}}}{\sqrt{k_{\alpha }k}}\frac{k_{\alpha }{J}_{l_{\alpha }+1/2}(k{r}_{\text{max}})J_{l_{\alpha }-1/2}(k_{\alpha }{r}_{\text{max}})}{k^2-k_{\alpha }^2},\end{array}$$

$$N_{\alpha }=I_{\alpha }(k_{\alpha },r_{\text{max}})$$

$$=\underset{k\to k_{\alpha }}{\lim }\frac{\pi r_{\text{max}}}{2}\frac{J_{l_{\alpha }+1/2}(r_{\text{max}})J_{l_{\alpha }-1/2}(k_{\alpha }{r}_{\text{max}})}{k^2-k_{\alpha }^2}$$

$$=\underset{k\to k_{\alpha }}{\lim }\frac{\pi r_{\text{max}}}{2}\frac{r_{\text{max}}{J}_{l_{\alpha }+1/2}^{\prime }(k{r}_{\text{max}})J_{l_{\alpha }-1/2}(k_{\alpha }{r}_{\text{max}})}{2k}$$

$$=\underset{k\to k_{\alpha }}{\lim }\frac{\pi r_{\text{max}}^2}{2}\frac{\left[\frac{l_{\alpha }+1/2}{k{r}_{\text{max}}}{J}_{l_{\alpha }+1/2}(k{r}_{\text{max}})-J_{l_{\alpha }+3/2}(k_{\alpha }{r}_{\text{max}})\right]J_{l_{\alpha }-1/2}(k_{\alpha }{r}_{\text{max}})}{2k}$$

$$=-\frac{\pi r_{\text{max}}^2}{4k_{\alpha }}{J}_{l_{\alpha }+3/2}(k_{\alpha }{r}_{\text{max}})J_{l_{\alpha }-1/2}(k_{\alpha }{r}_{\text{max}}).$$

$$\begin{array}{rl}{C}_{\text{SSC}}^{\alpha \beta }=⟨{\delta }^{\alpha }{\delta }^{\beta }⟩& =\frac{(4\pi {)}^2}{N_{\alpha }{N}_{\beta }}\int \frac{d^3{\mathbf{k}}_1}{(2\pi {)}^3}\frac{d^3{\mathbf{k}}_2}{(2\pi {)}^3}⟨\delta ({\mathbf{k}}_1){\delta }^{\ast }({\mathbf{k}}_2)⟩Y_{l_{\alpha }{m}_{\alpha }}^{\text{R}}(\widehat{k_1})Y_{l_{\beta }{m}_{\beta }}^{\text{R}}(\widehat{k_2})\\ & \times {\int }_0^{r_{\text{max}}}dr{r}^2{j}_{l_{\alpha }}(k_{\alpha }r)j_{l_{\alpha }}(k_{1}r){\int }_0^{r_{\text{max}}}dr{r}^2{j}_{l_{\beta }}(k_{\beta }r)j_{l_{\beta }}(k_{2}r)\\ & =\frac{(4\pi {)}^2}{N_{\alpha }{N}_{\beta }}\int \frac{d^3\mathbf{k}}{(2\pi {)}^3}{P}^{\delta \delta }(k)I_{\alpha }(k,r_{\text{max}})I_{\beta }(k,r_{\text{max}})Y_{l_{\alpha }{m}_{\alpha }}^{\text{R}}(\widehat{k_1})Y_{l_{\beta }{m}_{\beta }}^R(-\widehat{k_1})\\ & =\frac{2}{\pi N_{\alpha }{N}_{\beta }}\int k^{2}dk{P}^{\delta \delta }(k){\delta }_{l_{\alpha },l_{\beta }}{\delta }_{m_{\alpha },m_{\beta }}{I}_{\alpha }(k,r_{\text{max}})I_{\beta }(k,r_{\text{max}}),\end{array}$$

Full text

1 Introduction

Observational cosmology has experienced a remarkable rebirth as a precision science in the last two decades. In 1998, two teams using Type Ia supernovae found that the expansion of the Universe is accelerating [1, 2]. Galaxy surveys such as the 2dF Galaxy Redshift Survey [3], SDSS [4], BOSS [5], and WiggleZ [6] have provided an independent means of exploring the recent expansion history of the Universe and have dramatically confirmed cosmic acceleration. Meanwhile, the high-redshift Universe has been probed to great precision using the full sky surveys of the cosmic microwave background by WMAP [7, 8] Planck [9, 10], ACT [11], and SPT [12].

The cause of the accelerated expansion of the universe is tied to the origin and ultimate fate of the universe. Depending on how the acceleration changes over time, the universe could expand forever, re-collapse, rip itself apart, or perhaps do something else entirely [13, 14]. The most common parametrization for the cause of cosmic acceleration introduces a fluid with positive energy density but negative pressure, called dark energy. Dark energy models are characterized by an equation of state parameter $w\equiv\frac{p}{\rho}$, which relates its pressure $p$ and the energy density $\rho$ [15]. In the simple case where $w=-1$, dark energy is a cosmological constant, which is the standard $\Lambda$CDM cosmological model. The most studied physical mechanisms which could produce a $w\neq-1$ would also allow $w$ to vary as a function of time, $w=w(z)$ [16]. Current cosmological observations do not give good constraints on $w(z)$ because errors on $w(z)$ at different times are highly correlated [17, 18].

Tomographic weak lensing surveys directly probe the statistical properties of the matter distribution in projection, using sources at a range of redshifts, and hence can be used to constrain $w(z)$ [17, 19]. Several completed and ongoing experiments have measured or will measure weak lensing, including CFHTLens [20], KiDS [21], DES [22], and HSC [23]. Ideally, the precision of weak lensing measurements would be limited primarily by fundamental statistical uncertainties due to their limited samples of galaxies and survey volumes. Such statistical uncertainties will be greatly reduced by the large volumes of future surveys, such as LSST [24], Euclid [25], and WFIRST [26]. However, future weak lensing surveys will also suffer from many sources of observational and astrophysical systematic errors, which will need to be properly understood to avoid systematic biases in the results [27, 28, 29, 30]. One source of statistical error for future experiments to consider is the coupling of measured results to matter density fluctuations outside the survey window, called super-sample covariance (SSC). The direct coupling between the amplitude of short-wavelength modes, which fit in the survey window, and long-wavelength modes, which do not, was described as beat-coupling in refs. [31, 32] and found to be the dominant contribution to the matter power spectrum covariance on small scales. Further studies have considered other aspects of the effect and described formalism for analyzing it [33, 34, 35, 36, 37, 38]. The SSC effect on lensing observations has been shown to depend strongly on the details of the survey geometry [39]. Survey geometry design is subject to many practical instrument constraints, such as observatory slew times, galactic plane avoidance, and calibration. Future experiments, such as WFIRST, should consider whether optimizing survey design to reduce super-sample covariance is a significant additional consideration. Additionally, missions should consider whether combining other probes of long-wavelength density fluctuations can be a useful means of mitigating their covariances [40, 41].

In this paper, we present a formalism for estimating the SSC contribution to the covariances of future experiments and investigating possible mitigation strategies in a realistic survey geometry. We implement our formalism in a public code, SuperSCRAM: Super-Sample Covariance Reduction and Mitigation.111SuperSCRAM is available at https://github.com/mcdigman/SuperSCRAM The code is intended to be modular and extensible to facilitate the addition of further observables, physical effects, and mitigation strategies. SuperSCRAM provides Fisher matrix forecasts of parameter constraints, which is useful for estimation of the relative magnitude of effects and effectiveness of mitigation strategies, but is not a substitute for a full Markov chain Monte Carlo likelihood code such as CosmoLike [40, 41]. We do not analyze other sources of non-Gaussian covariance, although SuperSCRAM could be extended to include them.

Refs. [42, 43, 44, 45] have applied a similar harmonic expansion formalism to calculate SSC effects on cluster counts including geometrical effects. This paper focuses on weak lensing observables, although the code’s structure is modular, so an extension to include cluster counts in the future would be straightforward.

The structure of the paper is as follows. In section 2, we describe the formalism and overall analysis workflow used throughout the paper. In section 3, we describe the specific implementation of our formalism in SuperSCRAM. In subsection 4.1, we present a practical demonstration of a useful application of SuperSCRAM, in which we assess the effectiveness of embedding the WFIRST weak lensing survey program in the LSST field of view. In subsection 4.2 and subsection 4.3, we discuss other applications and extensions of SuperSCRAM. Finally, in section 5 we discuss our conclusions and the outlook for future work.

2 Overview

2.1 Notation and Definitions

For our analysis, we define two types of observables, which must be handled in separate ways: short-wavelength (SW) observables, denoted $\{O^{I}\}$, and long-wavelength (LW) observables, denoted $\{O^{a}\}$. The cosmological parameters are denoted $\{\Theta^{i}\}$. Finally, the long-wavelength matter density modes are expanded in terms of a set of basis modes $\{\psi_{\alpha}({\bf r})\}$, $\delta({\bf r})=\sum_{\alpha}\delta^{\alpha}\psi_{\alpha}({\bf r})$. In this paper, we use indices $IJK...$ for SW observables, $abc...$ for LW observables, $ijk...$ for cosmological parameters, and $\alpha\beta\gamma...$ for the density mode basis.

Short-wavelength observables, discussed in subsection 3.2, are observables which can only directly probe fluctuations of shorter wavelength than the survey’s window, such as the shear-shear lensing power spectrum, or cluster counts. In this paper, they are the “primary” source of cosmological information, which a survey is designed to extract.

Long-wavelength matter density fluctuations, with wavelengths larger than the survey window, add covariance to estimators of the short-wavelength modes measured within the window, which ultimately degrades cosmological parameter estimation. Long-wavelength observables, discussed in subsection 3.3, are chosen to provide information about the amplitude of the long-wavelength matter density fluctuations, so the degradation of estimates of short-wavelength observables by super-sample covariance can be mitigated. Long-wavelength observables need not originate from the same survey that measures the short-wavelength observables.

2.2 Workflow

A sketch of the workflow of our mitigation procedure is provided in figure 1.

To provide cosmological parameter forecasts, we calculate the total Fisher information matrix for the cosmological parameters by combining prior constraints and the information obtained from the short-wavelength observables measured by the survey:

[TABLE]

where the indices $i$ and $j$ run over all possible combinations of cosmological parameters. For the purposes of this paper, we use the WMAP5 fiducial point with a set of forecast constraints from Planck for $F_{ij}^{\text{prior}}$, taken from ref. [46]. To calculate $F_{ij}^{\text{survey}}$, we use

[TABLE]

where here $O^{I}$ are the short-wavelength observable vectors. $F_{IJ}^{\text{tot}}=\left(C^{IJ}_{\text{tot}}\right)^{-1}$ is calculated as the inverse of the short-wavelength covariance matrix,

[TABLE]

In this paper, we do not consider the non-Gaussian covariance $C^{IJ}_{\text{ng}}$, as our analysis is primarily focused on estimating the contribution of the super-sample covariance $C^{IJ}_{\text{SSC}}$ term. We decompose the long-wavelength modes into a basis as described in Appendix A, and calculate the super-sample covariance in our basis using perturbation theory, such that we can write

[TABLE]

The response of the observable to a density fluctuation $\partial O^{I}/\partial\delta^{\alpha}$ is described in subsection A.1. $F_{\alpha\beta}^{\text{SSC}}$ is the Fisher information matrix for long-wavelength matter density fluctuations decomposed into our basis, which includes the SSC contribution from perturbation theory given in eq. (A.8), plus any information provided by mitigation strategies,

[TABLE]

The response to the long-wavelength fluctuations and the Fisher information matrices of long-wavelength observables contain explicit dependence on the geometry of the survey. Note that eq. (2.5) implicitly assumes that the long-wavelength observables are statistically independent on large scales. Statistical independence is generally a reasonable assumption if the two observables are coming from different types of measurements or non-overlapping surveys. However, care must be taken to avoid adding the same information twice if observables contain information from overlapping datasets.

It is worth noting that the SSC term is not purely an observational limitation from partial sky coverage, and $C^{IJ}_{\text{SSC}}$ would not vanish even for a full sky survey, because an important component of the SSC term is radial. In fact, the improvements to the Gaussian covariance from increasing the sky area may make the SSC term relatively more important. Therefore, the SSC term is an important component of the error budget even for large surveys [47]. The specific mitigation approach implemented in this paper — embedding a survey in a larger survey — would not work in this case, but it is possible that other approaches would.

2.3 Eigenvalue analysis

Inspecting the individual values of the elements of the covariance matrix of the cosmological parameters $C^{ij}_{\text{tot}}$ can give a rough idea of the magnitude of the SSC contamination or effectiveness of a mitigation strategy. However, a more systematic method is required to quantitatively assess the overall impact. The direction of parameter space most contaminated by SSC or improved by mitigation is very difficult to determine by visual inspection of two-dimensional slices of many-dimensional ellipsoids, such as figure 3.

To assess the overall impact, we use one covariance matrix, as a metric to identify the directions most changed in another perturbed covariance matrix. For our purposes, a useful metric is the purely Gaussian covariance in cosmological parameter space, $C^{ij}_{\text{g}}$. For this section, we will refer to the covariance in cosmological parameter space with the super-sample contribution but no mitigation as $C^{ij}_{\text{SSC}}$, and the covariance with all mitigation applied as $C^{ij}_{\text{mit}}$. We then consider the product matrix

[TABLE]

which has eigenvalues $\lambda^{(i)}$ and eigenvectors $v^{i}_{\;j}$ such that $\mathcal{M}^{i}_{\;j}v^{j(k)}=\lambda^{(k)}v^{i(k)}$. For this section, Latin indices $ijk...$ run over cosmological parameter space and must be raised or lowered by applying the $C_{\text{g}}^{ij}$ or its inverse $(C_{\text{g}}^{-1})_{ij}$. Indices in parentheses $(i)(j)(k)...$ are labels for the eigenvalues and eigenvectors, and are not summed over here. If $C^{ij}_{\text{SSC}}$ is a positive semidefinite perturbation to the metric, all of the eigenvalues satisfy $\lambda^{(i)}\geq 1$, and the largest eigenvalue corresponds to the direction in parameter space most contaminated by $C_{\text{SSC}}^{ij}$. The resulting eigenvalues are dimensionless, and can therefore be compared between cosmological parametrizations. Note that the matrix $\mathcal{M}^{i}_{\;j}$ is not symmetric in general, which is inconvenient for computing eigenvalues numerically, so we instead find the eigenvalues of the matrix

[TABLE]

where $L^{i\hat{j}}$ is the lower triangular Cholesky decomposition $C_{\text{g}}^{ij}\equiv L^{i\hat{k}}L^{j\hat{k}}$. $\widetilde{\mathcal{M}}_{\hat{i}\hat{j}}$ is a symmetric matrix with the same eigenvalues $\lambda^{(i)}$ as $\mathcal{M}^{i}_{j}$ and eigenvectors $u_{\hat{i}}^{\;(j)}$ related to the eigenvectors $v^{i(j)}$ of $\mathcal{M}^{i}_{\;j}$ by transforming $u_{\hat{i}}^{\;(j)}=(L^{-1})_{k\hat{i}}v^{k(j)}$. The transformed eigenvectors then exist in a new vector space where the metric is $\mathbb{1}$. Indices in this vector space, $\hat{i}\hat{j}\hat{k}...$, may therefore be raised and lowered freely.

Now, we would like to interpret the most contaminated directions in our physical parameter space. We define the one-form $\tilde{v}_{i}^{\;(j)}\equiv(C_{\text{g}}^{-1})_{ik}v^{k(j)}$. $\tilde{v}_{i}^{\;(j)}$ has the property $\tilde{v}_{j}^{\;(i)}C_{\text{SSC}}^{jk}\tilde{v}_{k}^{\;(l)}=u^{\hat{j}(i)}\widetilde{\mathcal{M}}_{\hat{j}\hat{k}}u^{\hat{k}(l)}$, such that $\tilde{v}_{i}^{\;(j)}$ gives the directions in the physical parameter space most contaminated by super-sample covariance relative to the Gaussian covariance. We require the normalization $u^{\hat{k}(i)}u_{\hat{k}}^{(j)}\equiv\mathbb{1}^{(i)(j)}$.

The covariance matrix in the basis of directions most contaminated by super-sample covariance is given:

[TABLE]

where $C_{\text{mit}}^{ij}$ is the covariance matrix with both super-sample covariance and mitigation applied. The elements of $\widetilde{C}^{(i)(j)}_{\text{SSC},\text{mit}}$ give the covariance relative to $C_{\text{g}}^{ij}$ in the directions most contaminated without mitigation. If mitigation is perfect, $C_{\text{mit}}^{ij}=C_{\text{g}}^{ij}$, we have $\widetilde{C}^{(i)(j)}_{\text{mit},C_{\text{g}}}=\mathbb{1}^{(i)(j)}$. If instead the mitigation had no effect, $C^{ij}_{\text{mit}}=C^{ij}_{\text{SSC}}$, we have $\widetilde{C}^{(i)(j)}_{\text{SSC},\text{SSC}}=\lambda^{(i)}\mathbb{1}^{(i)(j)}$.

The product of the eigenvalues is the ratio of the volumes of the uncertainty ellipsoid with and without contamination, $\prod_{(i)}\lambda^{(i)}={\rm det}({\bf C}_{\text{SSC}})/{\rm det}({\bf C}_{\text{g}})$. Note that because changing the cosmological parametrization affects both $C^{ij}_{\text{SSC}}$ and $C_{\text{g}}^{ij}$ in the same way, $\prod_{(i)}\lambda^{(i)}$ and the largest eigenvalue $\lambda^{\text{top}}$ are both relatively insensitive to the choice of parametrization, which makes them useful for examining the effect of super-sample covariance in isolation. The product of eigenvalues and top eigenvalue are two useful points of comparison between different survey and mitigation strategies.

In subsection 4.1, we examine the eigenvalues most contaminated by the addition of the SSC term before and after mitigation, $\lambda^{(i)}_{\text{SSC}}$, and $\lambda^{(i)}_{\text{mit}}$.

3 Methods and Code Organization

In this section, we present a framework and code for forecasting super-sample effects on future observations in in a survey with a given geometry. SuperSCRAM is written and tested in Python 2.7 and uses the SciPy [48], NumPy [49], and Astropy [50, 51] libraries. Convergence tests were run on the CCAPP condos of the Pitzer and Ruby clusters at the Ohio Supercomputer Center [52]. The overall structure of the important modules in SuperSCRAM is sketched out in figure 2.

3.1 Geometries

The base class for survey geometries is the Geo class, which contains basic functionality for setting up the tomographic lensing bins and the resolution bins to be integrated over. Classes that extend Geo for a specific angular window function must provide a method for calculating the angular area in steradians, angular_area(), and the spherical harmonic coefficients $a_{lm}$, a_lm(l,m).

We provide several possible methods of implementing survey geometries. Pixelated geometries extend the PixelGeo class, which allows for an arbitrary pixelation of the sky, though all current geometries use HEALPix [53]. Our PolygonGeo class allows a survey geometry to be input as an arbitrary bounding spherical polygon, with segments of great circles as edges. The primary advantage of this class is that the spherical harmonic decomposition $a_{lm}$ can be calculated analytically for such polygons, as described in Appendix C. For comparison with pixelated methods also provide the PolygonPixelGeo class, which implements spherical polygons using a HEALPix pixelation.

We have implemented several utility geometries which perform useful operations on other geometries. The AlmDifferenceGeo implements the difference between any two Geo objects, which is necessary for cutting a mask out of a survey window. AlmRotationGeo uses the Euler rotation method described in Appendix C to rotate any Geo to a specified position on the sky. PolygonUnionGeo and PolygonPixelUnionGeo calculate the union of PolygonGeo and PolygonPixelGeo objects respectively; we have not currently implemented a more general purpose class for calculating unions or intersections between two arbitrary geometries.

We also provide a number of demonstration and testing geometries in premade_geos.py. For our demonstration we use the LSSTGeo and WFIRSTGeo, which implement approximations of possible LSST and WFIRST window functions using PolygonGeo objects, combined using the various utility geometries. Various other geometries used for the demonstrations shown in figure 6 include StripeGeo, which describes a spherical rectangle, and CircleGeo, which approximates a circle as a many-sided spherical polygon. For testing purposes, we have implemented FullSkyGeo and HalfSkyGeo, which use the analytic $a_{lm}$ for the full sky or an entire hemisphere respectively.

Various utility functions, such as those needed to reconstruct the images of geometries in figure 6 from their $a_{lm}$ representations, are provided in alm_utils.py, ylm_utils.py, and geo_display_utils.py.

3.2 Short-Wavelength Observables

The SWObservable class provides an interface which implementations of specific short-wavelength observables must provide. A short-wavelength observable must implement the method get_dO_I_ddelta_bar, which should calculate ${\partial O^{I}}/{\partial\bar{\delta}}$ and return it as a NumPy array. All short-wavelength observables currently implemented in SuperSCRAM are subclasses of LensingObservable, which are projected power spectra calculated by ShearPower, as described in Appendix B. We treat lensing between each pair of tomographic bins as separate LensingObservable objects. To avoid unnecessarily recomputing power spectra in different tomographic bins, each SWSurvey object creates a single shared LensingPowerBase object to store all the ShearPower objects which its LensingObservable objects will need.

To create a specific lensing observable, a LensingObservable subclass, such as the shear-shear lensing power spectrum ShearShearLensingObservable, needs to be provided with two QWeight objects, in this case QShear, which are the functions $q_{A}(z)$ described in Appendix B.

3.3 Long-Wavelength Observables

In this section, we describe a sample mitigation strategy of embedding a weak lensing survey, such as WFIRST, in a wider but shallower survey of galaxy number counts, such as LSST. The difference between the number density measured by the larger survey inside the lensing survey window, $n_{1}$, and the number density outside, $n_{2}$, provides information about long-wavelength density fluctuations. Our observable is therefore $\Delta n_{12}\equiv n_{1}-n_{2}$. The number density can be written to linear order in perturbation theory as

[TABLE]

where $b(r)$ is the linear bias, $n(r)$ is the expected number density of galaxies at comoving distance $r$, $\Omega_{1}$ is the angular area of the window, and $\delta(r,\Omega)=\sum_{\alpha}{\delta^{\alpha}\psi_{\alpha}(r,\Omega)}$ is a density fluctuation expanded in our basis described in Appendix A.

For the number density of galaxies as a function of redshift in the mitigation survey, we use the expected number densities for the LSST 1-year survey to a depth of $i<24$ from ref. [54], which gives approximately $18$ galaxies/arcmin2. We have deliberately chosen a very conservative use of LSST. We obtain the bias by abundance matching to the Sheth-Tormen halo mass function, as described in Appendix D. Assuming the wider survey has uniform sensitivity as a function of sky position, we have the response to long-wavelength density perturbations:

[TABLE]

The “noise” in the long-wavelength observable is the shot noise, since for our application the clustering of the galaxies due to the underlying large-scale matter density field is the signal of interest. The mean number density of galaxies in this slice is

[TABLE]

We can calculate the variance under the assumption that the shot noise is Poissonian:

[TABLE]

where $V_{i}\equiv(r_{\text{max}}^{3}-r_{\text{min}}^{3})\Omega_{i}/3$. In real samples, the long-wavelength shot noise may deviate from the Poisson prediction due to non-linear clustering, exclusion, and satellite galaxies (see, e.g., ref. [55] for an extensive discussion). For our example application in section 4, the mitigation performance is only weakly dependent on $N_{12}$.222Increasing the galaxy number count survey’s limiting magnitude from $i<24.1$ to $i<25.1$ more than doubles the number of galaxies, but reduces the volume of the error ellipsoid in our 7D cosmological parameter space by $<1\%$.

Then, we can write the Fisher matrix for this mitigation strategy

[TABLE]

This long-wavelength observable is implemented in the DNumberDensityObservable class, which extends the general specifications for a long-wavelength observable, LWObservable. Subclasses of LWObservable should implement the method get_fisher, which DNumberDensityObservable obtains using eq. (3.5) by numerically integrating eq. (3.2) as described in subsection A.1. Alternatively, to use the more efficient numerical Fisher matrix manipulation method described in subsection 3.6.1, an LWObservable may instead implement a get_perturbing_vector method, which for DNumberDensityObservable returns $\frac{\partial\Delta n_{12}}{\partial\delta^{\alpha}}$ and $(N_{12})^{-1}$ for the perturbing vector ${\bf V}$ and inverse variance ${\bf K}$ respectively. For the current implementation $(N_{12})^{-1}$ is required to be diagonal. Which method should be used is toggled by setting the LWObservable object’s fisher_type variable to True to use get_fisher and False to use get_perturbing_vector.

3.4 Cosmological Parametrizations

In order to properly calculate the response of a set of short-wavelength observables $\{O^{I}\}$ to cosmological parameters $\{\Theta^{i}\}$, ${\partial O^{I}}/{\partial\Theta^{i}}$, we must select a cosmological parametrization. The parametrization is specified in the CosmoPie module, along with a set of rules for calculating derived parameters, so that the parametrization can easily be interchanged. Parameters relating to the dark energy equation of state $w(z)$ are handled separately from the rest. By default, SuperSCRAM uses $\{\Omega_{m}h^{2},\Omega_{b}h^{2},\Omega_{de}h^{2},n_{s},\ln{A_{s}}\}$ as its basic set of parameters, where $\Omega_{m},\Omega_{b},$, and $\Omega_{de}$ are the total matter, baryon, dark energy densities respectively. This parametrization is chosen to be compatible with the parametrization in [46]. For simplicity, we set $\Omega_{k}=0$. If we changed the fiducial model to a curved Universe, then we would have to build the basis modes $\psi_{\alpha}({\bf r})$ for a curved Universe, using the same ultraspherical Bessel functions that appear in CMB anisotropy calculations [56]. We also fix the neutrino masses to be $\sum m_{\nu}=0$ in the present implementation.

3.5 Matter Power Spectrum

SuperSCRAM provides several interchangeable sample implementations of the matter power spectrum accessible using a MatterPowerSpectrum object, including Python implementations of the Takahashi [57] and Casarini [58] revisions of the Halofit [59] model, linear matter power spectra from CAMB [60], and the FAST-PT [61, 62] implementation of the one-loop power spectrum. Additionally, we have extended both Halofit and FAST-PT to facilitate an arbitrary $w(z)$ using the same procedure used in refs. [58, 63] to extend the Halofit model as described in Appendix G, although we have not attempted to validate the resulting power spectra through simulations. In this paper, we report only results for the Halofit model.

3.6 Fisher Matrix Manipulation

Fisher and covariance matrix manipulations are primarily handled by the FisherMatrix class. The FisherMatrix class is optimized to efficiently perform the manipulations of Fisher matrices required by SuperSCRAM with sizes up to the memory limitations of the machine it runs on. For memory efficiency, it can internally store either a Fisher matrix, covariance matrix, or their Cholesky decompositions. The representations can be exchanged in place internally, and all external class methods intended for use by other modules can be used regardless of the internal representation. For some functions, the choice of internal representation can affect run time significantly, especially for large matrices. The overall logic for managing the Fisher matrices for the short and long-wavelength observables and cosmological parameters is handled by the MultiFisher class. A MultiFisher object takes an LWBasis object, an SWSurvey object, and a LWSurvey object. From these objects, it extracts the appropriate projection matrices for converting between the space of long-wavelength and short-wavelength observables, and the space of short-wavelength observables and cosmological parameters, ${\partial O^{I}}/{\partial\delta^{\alpha}}$ and ${\partial O^{i}}/{\partial\Theta^{i}}$ respectively. It then performs the necessary projections and applies priors internally. As an example, if we have already created a Fisher matrix object called multi_f, we can use

f_set = multi_f.get_fisher_set(include_priors=False)

to obtain a set of FisherMatrix objects in cosmological parameter space containing just the Gaussian covariance and total covariance with and without mitigation. To then obtain the eigensystems with and without mitigation using the Gaussian covariance as a metric, as described in subsection 2.3, we can use

eig_set = multi_f.get_eig_set(f_set)

3.6.1 Long-Wavelength Convergence

Obtaining well converged results for the super-sample covariance in our formalism may require a long-wavelength Fisher matrix larger than can be stored in memory. However, because the long-wavelength covariance matrix $C^{\alpha\beta}_{\text{SSC}}$ as written in Appendix A is block diagonal, it is possible to calculate the contribution to the short-wavelength covariance matrix $C^{IJ}_{\text{SSC}}$ from each block individually, which requires individual matrices only as large as $(2l_{\text{max}}+1)^{2}$. However, taking advantage of the block diagonal covariance matrix prevents directly perturbing the Fisher matrix with mitigating information as in eq. (2.5). Provided it is possible to decompose

[TABLE]

as in eq. (3.5),333The implementation in SuperSCRAM has a diagonal ${\bf K}$, but this is not essential to the formalism. we may instead use the Woodbury matrix identity

[TABLE]

Here $C^{\alpha\beta}_{\text{SSC+mit}}$ is not necessarily block diagonal, so storing it would require a full-size matrix. However, the matrix

[TABLE]

can be computed without storing the full $C^{\alpha\beta}_{\text{SSC}}$, as can $X^{a\alpha}\equiv V^{a}{{}_{\beta}}L^{\beta\alpha}_{\text{SSC}}$, using the fact that the Cholesky decomposition $L^{\alpha\beta}_{\text{SSC}}$ of a block diagonal matrix $C^{\alpha\beta}_{\text{SSC}}$ is also block diagonal. We can then write

[TABLE]

which can now be computed without ever storing a full-size matrix. In addition to allowing better-converged results due to the much larger number of practically usable basis elements, this procedure is also orders of magnitude faster and is more numerically stable than perturbing the Fisher matrix directly. This procedure is implemented in the SphBasisK module. To allow this procedure to be used, all LWObservable objects must set fisher_type=False and create a get_perturbing_vector method, which returns ${\bf V}$ and ${\bf K}$, for all long-wavelength observables in the survey. We implement this functionality for our DNumberDensityObservable class.

3.7 Long-Wavelength Basis

A basis for long-wavelength observations can be specified as an extension of the LWBasis class. A long-wavelength basis must at minimum implement get_fisher(), which returns a FisherMatrix object, and get_dO_I_ddelta_alpha(geo,integrand), which calculates ${\partial O^{I}}/{\partial\delta^{\alpha}}$ by integrating an observable input in a grid $[{\partial O^{I}}/{\partial\bar{\delta}}](z)$, where the grid should correspond to the fine $z$ grid of a given Geo object. The long-wavelength basis described in Appendix A is provided by the SphBasisK class.

4 Applications

4.1 Example Application: WFIRST+LSST

As a demonstration, we calculate the importance of the super-sample covariance in a simulated WFIRST weak lensing survey footprint, with mitigation from a simulated LSST footprint. For the LSST footprint, we use a simplified survey window from $+5^{\circ}$ to $-65^{\circ}$ equatorial latitude spanning $360^{\circ}$ in longitude, with a $\pm 20^{\circ}$ mask around the galactic plane, which gives about $13600$ deg2 of usable survey area. For the WFIRST footprint, we use an approximately $2100$ deg2 footprint, depicted along with the LSST footprint in the right panel of figure 6. Neither footprint is necessarily what will be implemented; for example, the LSST DESC has proposed an extended footprint going farther north [64].

For LSST, we use the forecast number density of galaxies for the LSST year 1 dataset with $i<24.1$ from ref. [54], which gives a total of about $18$ galaxies/arcmin2. We choose the $i<24.1$ cutoff to give a conservative estimate of the number of galaxies with good photometric redshifts; by the time of the WFIRST analysis, the photometric redshifts for these bright objects should be very good. We use 4 evenly spaced redshift bins of width $\Delta z=0.3$ from $z=0$ to $z=1.2$, which should be very large compared to the photometric redshift uncertainties. For the lensing source galaxies, we use the WFIRST weak lensing sources, using typical High Latitude Survey parameters of 5 exposures of 140 s each. The source densities were computed using v15 of the WFIRST Exposure Time Calculator [65] and the Phase A throughput tables.444See https://wfirst.gsfc.nasa.gov/science/WFIRST_Reference_Information.html The input galaxy catalog was from the CANDELS GOODS-S region [66], using the photometric redshifts from ref. [67]. Cuts were applied requiring detection S/N of the galaxy $>18$, ellipticity error per component $<0.2$, and resolution factor $>0.4$ in the convention of ref. [68]. These cuts leave about $34$ galaxies/arcmin2 for H band only. The actual source density will be larger in the combined images from all 4 filters, but, for simplicity, we focus only on the H band. The results are most sensitive to the number density chosen for the WFIRST lensing source galaxies because increasing the number of sources decreases the Gaussian covariance. These number densities should give an adequate approximation of the survey galaxy number densities for purposes of this demonstration.

The intent of this demonstration is to give a rough idea of the impact of super-sample covariance alone on a realistically shaped survey footprint and to give a rough idea of how much improvement could be expected from a simple mitigation strategy. It is not intended to be a full forecast of the actual constraints WFIRST will be able to achieve. In particular, we do not include any systematic or non-Gaussian effects beyond the super-sample covariance. This is because our philosophy is to compare each added term in the error budget (such as SSC) to the irreducible statistical error bars, and not to the combination of all other uncertainties where we have mitigation efforts in progress.

In Figures 3 and 4, we show traditional triangle plots of our constraints with different reasonable priors applied. It is apparent from these plots that super-sample covariance has an effect. However, the overall magnitude is difficult to read off, and different choices of priors change the apparent effect of super-sample covariance, which is undesirable for comparison purposes. Instead, we recommend comparing the magnitude of the effect in the most contaminated directions in parameter space, as described in subsection 2.3. The comparison is shown in figure 5, with the corresponding effect summarized numerically Table 1. From figure 5, it is much more apparent that there exist directions significantly affected by super-sample covariance. Additionally, our formalism removes the dependence on the choice of priors by isolating the effect of super-sample covariance on the constraining power of the survey. Figure 5 and Table 1 show more clearly that a mitigation strategy can significantly reduce the amount of constraining power lost due to super-sample covariance.

4.2 Sensitivity to Survey Geometry

One application of SuperSCRAM is to test the impact of varying the survey geometry. Previous analytic studies [40] have generally considered only a circular survey geometry. Results comparing several different survey geometries are shown in figure 6. The magnitude of the effect varies at approximately the $O(10\%)$ level between different contiguous, reasonably compact high latitude survey geometries. Because more compact geometries produce a smaller super-sample covariance, estimates using only a circular geometry may tend to understate the magnitude of the effect. The differences between geometries are small enough that it is possible that including additional super-sample effects beyond the density contrast effect considered here, such as tidal effects and redshift-space distortions, could change the relative ordering of the effect on different survey geometries. Because the mitigation reduces the super-sample contamination far more than the differences between any reasonable geometries, the possibility of applying mitigation strategies appears to be a more important survey-design consideration than the detailed survey geometry for controlling super-sample covariance.

4.3 Other Applications

A variety of possible applications are would be straightforward to implement but are not currently included. Additional short-wavelength observables, such as galaxy-galaxy lensing and galaxy power spectra, could be added to better understand the full constraining power of a weak lensing survey.

Many more mitigation strategies are possible. The current mean number density observable is only a simple demonstration of the information that can be extracted from a larger survey and does not include information about galaxy number density as function of sky position. Additionally, CMB lensing could be useful to constrain over-densities at higher redshifts. The mean tangential shear could also provide additional information.

It would be useful to compare the results of our Fisher matrix method to a Monte Carlo code. To achieve more realistic forecasts, photometric redshift uncertainties and other systematic and non-Gaussian effects would need to be implemented.

SuperSCRAM can be used to investigate the contribution of super-sample covariance to the tradeoff between a wider and a deeper survey. A more detailed model of the survey’s sensitivity as a function of redshift would be necessary to draw useful conclusions from this possible application.

5 Conclusions

In this paper, we have developed a Fisher matrix formalism and present the code SuperSCRAM for evaluating the impact of super-sample covariance on a realistically shaped weak lensing survey geometry. We have further applied our formalism to investigate the possibility of obtaining information from wider, shallower surveys to reduce the impact of super-sample covariance. Overall, we find that mitigation strategies can improve constraining power enough to merit serious consideration for next-generation weak lensing surveys.

Further work is needed to investigate other contributions to super-sample covariance beyond the density contrast effect discussed in this paper, such as tidal effects and redshift-space distortions. Including those effects would improve the ability to make detailed conclusions about the relative effectiveness of different survey geometries.

Additionally, it would be productive to investigate a wider array of mitigation strategies. In the simplest case, our mean number density observable does not exhaust the information available from a larger wider survey, as a more detailed analysis including observed number density as a function of sky position would likely improve that observable. Additionally, CMB lensing could potentially provide information about densities at much higher redshifts than can be obtained by a shallow galaxy survey, which could significantly improve constraining power. It may also be possible to use the tangential shear and other weak lensing related observables as further mitigating observables.

Acknowledgments

M.C.D. and C.M.H. were supported by the Simons Foundation, NASA, and the US Department of Energy during the preparation of this work.

We thank Elisabeth Krause for providing scripts to run CosmoLike for calibration of SuperSCRAM, as well as for general discussions. We thank Tim Eifler, Olivier Doré, Xiao Fang, and Benjamin Buckman for useful feedback.

The computations in this paper were run on the CCAPP condos of the Ruby and Pitzer Clusters at the Ohio Supercomputer Center.

This work is based in part on observations taken by the CANDELS Multi-Cycle Treasury Program with the NASA/ESA HST, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555.

Appendix A Basis Decomposition

In SuperSCRAM, we use a form of a spherical Fourier-Bessel basis, similar to the one described in ref. [69]. We use a spherical Bessel function for our radial basis and the real spherical harmonics for our angular basis, so that the elements of our basis can be written:

[TABLE]

where $Y^{\text{R}}_{lm}(\theta,\phi)$ is the real spherical harmonic, $k_{\alpha}$ is a solution of $j_{l_{\alpha}}(k_{\alpha}R_{max})=0$, and $\alpha$ is an index running over all possible sets $(k_{\alpha},l_{\alpha},m_{\alpha})$. To evaluate the covariance matrix with a finite number of modes, we choose a cutoff $k_{\text{max}}$ and take all combinations of $(k_{\alpha},l_{\alpha},m_{\alpha})$ which have $k_{\alpha}<k_{\text{max}}$. This cutoff procedure has the effect of gradually decreasing the number of modes included for a given $l_{\alpha}$ as $l_{\alpha}$ gets larger, and should converge faster than an arbitrary $l_{\text{max}}$ cutoff. The real spherical harmonics can be defined in terms of Legendre polynomials:

[TABLE]

We can then write a density fluctuation mode in this basis as

[TABLE]

where $N_{\alpha}$ is a normalization factor. In the second equality we use the Fourier transform relation $\delta(\mathbf{r})=(2\pi)^{-3}\int d^{3}\mathbf{k}\,e^{i\mathbf{k}\cdot\mathbf{r}}\delta(\mathbf{k})$, and in the third equality we use the identity

[TABLE]

The normalization factor $N_{\alpha}$ is defined

[TABLE]

where $\int_{\Omega}Y^{\text{R}}_{l_{\alpha}m_{\alpha}}(\hat{r})Y^{\text{R}}_{l_{\alpha}m_{\alpha}}(\hat{r})\,d^{2}\hat{\bf r}=1$. The integral $I_{\alpha}(k_{\alpha},r_{\text{max}})$ is defined as

[TABLE]

where in the last step we have use $J_{l_{\alpha}+1/2}(k_{\alpha}r_{\rm max})=0$, from the defining property of $k_{\alpha}$. In the special case $k=k_{\alpha}$, we can simplify $N_{\alpha}=I_{\alpha}(k,r_{\text{max}})$ using Bessel function identities:

[TABLE]

With these definitions, we can calculate the covariance matrix $C^{\alpha\beta}_{\text{SSC}}=\left(F^{\alpha\beta}_{\text{SSC}}\right)^{-1}$ in our basis:

[TABLE]

where $P^{\delta\delta}(k)$ is the matter power spectrum, and we used $(2\pi)^{3}\delta^{3}_{\text{D}}(\mathbf{k}+\mathbf{k}^{\prime})P^{\delta\delta}(k)=\langle\delta(\mathbf{k})\delta^{*}(\mathbf{k}^{\prime})\rangle$.

A.1 Response of Observables to Density Fluctuations

To calculate the response of an observable $O^{I}$ to a density fluctuation mode $\delta^{\alpha}$, we can use the chain rule:

[TABLE]

provided $[{\partial{O^{I}}}/{\partial\bar{\delta}}](z)$ can be calculated. In SuperSCRAM, the integral in eq. (A.9) is accomplished by calculating the integrand on a grid of $z$ values $\{z_{i}\}$ and using the trapezoidal rule. To calculate $[\partial\bar{\delta}/\partial\delta^{\alpha}](z_{i})$, we expand the mean density fluctuation $\bar{\delta}(z_{i})$, in our basis:

[TABLE]

where $r$ is the comoving distance in the range $r_{i}\leq r<r_{i+1}$, and $a_{l_{\alpha}m_{\alpha}}$ are the real spherical harmonic coefficients of a given survey window function $\Omega$. We can use eq. (A.10) to write the derivative:

[TABLE]

Appendix B Projected Power Spectra

The angular correlation function $w_{AB}(\hat{n}\cdot\hat{n}^{\prime})$ of the line of sight projections of two fields $A$ and $B$ can be expanded in terms of its angular power spectrum $C_{AB}\left(\ell\right)$ [70]:

[TABLE]

where $\hat{\bf n}$ and $\hat{\bf n}^{\prime}$ are unit vectors in the direction of observation and $P_{\ell}$ are Legendre polynomials. For a given field $A$, there is a weight function $q^{A}$ which relates the field to its line of sight projection $\tilde{A}(\hat{n})$, such that [70]:

[TABLE]

where $r$ is the comoving coordinate. In terms of $q_{A}$ and $q_{B}$, the angular cross-power spectrum can be written:

[TABLE]

where $j_{\ell}(kr)$ are spherical Bessel functions, and $k=(\ell+\frac{1}{2})/{r}$. In practice, this integral is inconvenient to compute numerically due to the rapid oscillations of the spherical Bessel functions for $kr\gtrsim\ell+\frac{1}{2}$. Therefore, most authors write $C^{AB}$ using the Limber approximation, which can be taken by expanding this expression to lowest order in $(\ell+\frac{1}{2})^{-1}\ll 1$,

[TABLE]

SuperSCRAM currently implements eq. (B.4) in the ShearPower class. Note that for $\ell\gg\frac{1}{2}$ many authors take $\ell+\frac{1}{2}\approx\ell$, although SuperSCRAM does not, both because the correction does not affect the code’s execution time, and for consistency with other implementations such as the one in CosmoSIS [71]. The next order correction is suppressed by a factor of $\mathcal{O}((l+\frac{1}{2})^{-2})$, which should be negligible for next generation weak lensing surveys.

For the weight functions $q_{A}(r)$ for shear-shear lensing, implemented as QShear, we use [41, 40]

[TABLE]

where

[TABLE]

where $p_{i}(r^{\prime})$ is the probability density function for source galaxies in tomographic bin $i$.

Appendix C Analytic Polygon Geometry

Using Stokes’s theorem, for a spherical polygon survey window, which has $N$ sides which are great circle arcs, we can write

[TABLE]

where $\partial\Omega_{n}$ denotes integration over the boundary of the $n$th arc and $\hat{z}_{n}$ is a unit vector orthogonal to the two unit vectors whose tips touch the ends of the $n$th arc, such that if ${\bf\hat{p}}_{n}$ is the unit vector at the start of the $n$th arc, ${\bf\hat{z}}_{n}\equiv{{\bf\hat{p}}_{n+1}\times{\bf\hat{p}}_{n}}/{|{\bf\hat{p}}_{n+1}\times{\bf\hat{p}}_{n}|}$, leaving ${\bf\hat{y}}_{n}\equiv{\bf\hat{z}}_{n}\times{\bf\hat{p}}_{n}$. The integral is most simple to evaluate if the great circle is along the equator at $\theta=\frac{\pi}{2}$, so for each side we rotate to a coordinate system where the side is along the equator, calculate $\Delta a^{\prime}_{lm,n}$ and rotate back to the global coordinate system using

[TABLE]

where $D^{n}_{lmm^{\prime}}$ is a spherical harmonic rotation matrix element.

The integrand restricted to the equator $\theta=\frac{\pi}{2}$ is given. Recalling eq. (A.2), we have

[TABLE]

where in the last step we have have used the recurrence relation

[TABLE]

Then, because the side has been rotated along the equator, we need only integrate $\phi$ from 0 to the side length $\beta_{n}$ and find

[TABLE]

where we have substituted $P_{(l-1)}^{m}(0)$ for $Y^{R}_{(l-1)m}(\frac{\pi}{2},0)$ using eq. (A.2). The side length is given $\beta_{n}\equiv\text{atan2}\left[|{\bf\hat{p}}_{n+1}\times{\bf\hat{p}}_{n}|,{\bf\hat{p}}_{n+1}\cdot{\bf\hat{p}}_{n}\right]$.

Now, to find the spherical harmonic rotation matrices, we obtain the Euler angles $\alpha_{n}$, $\psi_{n}$, and $\gamma_{n}$ for the $z-x-z$ Euler rotation necessary to rotate an arc from the global frame into a frame where the side lies along the equator in the $x-y$ plane. as depicted in figure 7. We define $\psi_{n}\equiv-\cos^{-1}\left({\bf\hat{z}}_{nz}\right)$, $\gamma_{n}\equiv\text{atan2}\left[{-{\bf\hat{z}}_{nx}},{{\bf\hat{z}}_{ny}}\right],$ $\alpha_{n}\equiv-\text{atan2}\left[-{\bf\hat{y}}_{nz},{\bf\hat{x}}_{nz}\right]$. Here are, by construction, the Euler angles for a $z-x-z$ Euler rotation which transforms from a local coordinate system where the $n$th arc lies in the $x-y$ plane to the global coordinate system, so that

[TABLE]

where the $z$-rotation matrices are given such that:

[TABLE]

We use an angle-doubling algorithm to compute the matrix ${\bf D}^{x}_{l}$. Writing ${\bf D}^{x}_{l}$ as a $(2l+1)\times(2l+1)$ matrix, ${\bf D}^{x}_{l}(\psi)={\bf E}_{l}(\psi)+\mathbb{1}$, and ${\bf E}_{l}(\psi)$ is computed by recursively applying the angle doubling formula ${\bf E}_{l}=2{\bf E}_{l}(\frac{\psi}{2})+[{\bf E}_{l}(\frac{\psi}{2})]^{2}$ to the infinitesimal rotation matrix ${\bf E}_{l}$. This matrix is constructed as ${\bf E}_{l}(\epsilon)=\mathcal{M}_{l}^{-1}\tilde{\bf E}_{l}\mathcal{M}_{l}$, where the infinitesimal complex-basis rotation around $x$ is

[TABLE]

and $\mathcal{M}_{l}$ transforms to the real spherical harmonic basis:

[TABLE]

In SuperSCRAM, the number of doublings can be specified by the user; for our WFIRST+LSST demo, 31 doublings are sufficient to recover $a_{lm}$ up to $l=100$ to within a maximum error of $\lesssim 0.1\%$, and approximately 82 doublings will converge all elements of the rotation matrices up to $l=200$ to within the limits of IEEE 745 64 bit precision. The repeated multiplications of $(2l+1)\times(2l+1)$ matrices are relatively time consuming for large $l$. An advantage of the analytic solution over a pixel based geometries is that the results are less vulnerable to pathological input survey geometries, such as a geometry with many narrow stripes (perhaps chip gaps) masked out, aligned so that none of the stripes contained any pixel centroids, which could produce poorer quality results in a pixel based geometry. Additionally the analytic solution can in principle be calculated for arbitrarily large $l$, while the pixelated solution is limited by the resolution of the pixelation scheme.

The various manipulations required to calculate eq. (C.2) are computed by the alm_utils module in SuperSCRAM. The PolygonGeo class provides a polygonal geometry with great circle arc sides which uses these analytic calculations, and the PolygonPixelGeo class provides the pixelated equivalent.

Appendix D Galaxy Number Density

To simulate a sample redshift distribution for a future weak lensing survey, we take all galaxies in the CANDELS [72, 73] GOODS-S catalogue with i band magnitude below a selected cutoff (i<24 is used as a default). We then calculate a simple smoothed number density as a function of redshift using a Gaussian smoothing kernel with user specified width $\sigma$ and reflecting boundary conditions at $z=0$:

[TABLE]

where $\Omega$ is the area of the CANDELS survey. Then we can calculate $n(z)$:

[TABLE]

Now, we want to use $n(z)$ to calculate $b(z)$, the bias. Because we are simply demonstrating a possible mitigation strategy and not attempting to model it in detail, we adopt a simple approximation in which we use this $n(z)$ to extract a $M_{\text{min}}$ for the halo mass function and use it to calculate the bias, neglecting any subtleties from the mass-luminosity relationship. We use the Sheth-Tormen halo mass function [74, 75, 76],

[TABLE]

where $\bar{\rho}$ is the background matter density in units of $h^{2}M_{\odot}\text{Mpc}^{-3}$, and $\sigma(M,z)$ is given by

[TABLE]

where $W(k,M)=W(k,R)=3j_{1}(kr)/(kr)$ is the Fourier transform of a spherical top hat window function, and $R=\frac{1}{h}[{3M}/(4\pi\bar{\rho})]^{1/3}$, $P_{\text{lin}}(k)$ is the linear power spectrum and $G(z)$ is the linear growth factor. The Sheth-Tormen mass function itself is

[TABLE]

where $A=0.3222$, $a=0.707$, and $p=0.3$ are empirically fit parameters, and we take $\delta_{c}=1.686$ to be the critical overdensity for spherical collapse to avoid considering any cosmology dependence. From here numerically solve $n(M_{\text{min}},z)=n(z)$ to obtain $M_{\text{min}}$. Then to calculate the linear bias, we use

[TABLE]

where $\sigma^{2}=\sigma^{2}(M_{\text{min}},z)$. SuperSCRAM could be extended to use improved fitting functions such as in refs. [77] or [78], although the Sheth-Tormen functions are sufficient for the demonstration in this paper.

In SuperSCRAM, eqs. (D.6) and (D.2) are calculated by the NZCandel class.

Appendix E Growth Factor

Simply modifying the calculation of the linear growth factor accounts for most of the correction to the matter power spectrum from dark energy with a variable equation of state. In general relativity, the linear growth factor $G(a)$ evolves according to

[TABLE]

where the primes denote derivatives with respect to a and the dots are derivatives with respect to cosmic time, and $H(a)$ is the Hubble rate [17]. In the presence of an arbitrary set of perfect fluids with densities $\{\Omega_{i}(a)\}$ and equations of state $\{w_{i}(a)\}$, the Hubble rate is

[TABLE]

and we can use $D(a)\equiv\frac{G(a)}{a}$ to arrive at a differential equation for $D(a)$:

[TABLE]

these equations take simple forms for constant w solutions, such as matter, with $w_{m}=0$, radiation with $w_{r}=\frac{1}{3}$, and curvature $w_{k}=-\frac{1}{3}$, and a cosmological constant $w_{\Lambda}=-1$. For dark energy, there are closed form solutions to the integral $e^{3\int_{a}^{1}(1+w(a^{\prime}))}$ for all dark energy parametrizations of $w(a)$ considered in this paper. SuperSCRAM can be extended to other parametrizations by simply providing functions returning $w(a)$ and $e^{3\int_{a}^{1}(1+w(a^{\prime}))}$. The CosmoPie class can then evaluate $D(a)$ using SciPy’s odeint to solve eq. (E.3).

Appendix F Code Tests

We have performed a number of tests to verify the various modules of SuperSCRAM perform as expected. Most of the code is covered by unit tests using pytest.

The FisherMatrix class and its associated algebra_utils module have unit tests for all their functions in fisher_tests.py and algebra_tests.py, which run every test with a suite of input matrices, including randomly generated matrices. For FisherMatrix, every test is also run for FisherMatrix objects with every possible internal state, to verify the external behavior of the functions is independent of internal state, as described in subsection 3.6.

The PolygonGeo, PolygonPixelGeo, and RectGeo classes have unit tests in polygon_geo_tests.py which verify that their results are both consistent between classes for the same geometries, and that their calculated $a_{lm}$ correctly describe the window function for the input survey geometry. Note that if results for $l>85$ are needed, SuperSCRAM requires the arbitrary precision mpmath package to avoid numerical overflows at double precision.

The power_response module’s predictions for $\frac{\partial P^{\delta\delta}}{\partial{\bar{\delta}}}$ are compared to the results of [79] in power_response_test.py, and qualitatively agree at the level expected given that they have convolved their power spectrum with a window function and we have not.

The response of a shear shear lensing observable $C^{\gamma\gamma}_{l}$ to a density perturbation as discussed in subsection A.1 and Appendix H and calculated by ShearPower should have functional dependences of the form $\frac{\partial C^{\gamma\gamma}_{l}}{\partial\bar{\delta}}(z_{s})\propto\frac{\Delta z_{s}}{z_{i}}C^{\gamma\gamma}_{l}$, where $\Delta z_{s}$ is the width of a resolution $z$ slice integrated over in (A.9) and $z_{i}$ is the average $z$ of the tomographic bin. The functional dependences of $\frac{\partial C^{\gamma\gamma}_{l}}{\partial\bar{\delta}}(z_{s})$ are checked by power_derivative_tests.py.

Our Python implementation of Halofit used by the MatterPowerSpectrum class described in subsection 3.5 is very similar to the Fortran implementation of the Takahashi Halofit prescription available in CAMB [60, 57], with some modifications to facilitate our prescriptions for arbitrary $w(z)$, and computing $\frac{\partial P^{\delta\delta}}{\partial\bar{\delta}}$. The agreement with the CAMB output is tested in power_comparison_tests.py. Our implementation agrees with CAMB’s implementation of the Takahashi Halofit prescription to within $0.2\%$ for $k\lesssim 10$ Mpc*-1*. Most of the residual difference is because we blend the transition between the linear and nonlinear power spectrum to avoid sharp spikes in $\frac{\partial P^{\delta\delta}}{\partial\bar{\delta}}$ in eq. (A.9), while CAMB’s implementation uses a sharp transition at $k=0.005h$ Mpc*-1*.

Our prescription for an arbitrary $w(a)$ prescription implemented in the WMatcher class, discussed in Appendix G is tested in w_matcher_tests.py. We test that our results for our w0wa model with $w(a)=w_{0}+(1-a)w_{a}$ match the results of [63], and that our 36 bin jdem model gives results consistent with the w0wa prescription for similar $w(a)$. We also test that both w0wa and jdem models correctly recover the known $w_{0}$ for constant $w(a)$ models.

The variance in a window with window function $W({\bf k})$ is given $\sigma^{2}=\int\frac{d^{3}{\bf k}}{(2\pi)^{3}}P({\bf k})W({\bf k})$ and can be written in our formalism $\sigma^{2}=\frac{\partial\bar{\delta}}{\partial\delta^{\alpha}}C^{\alpha\beta}\frac{\partial\bar{\delta}}{\partial\delta^{\beta}}$. The script variance_demo.py compares results from SuperSCRAM using PolygonGeo to a code which integrates the linear matter power spectrum directly for a rectangular window function. For an approximately cubic $500Mpc$ window, SuperSCRAM converges to within about $2\%$ of the $\sigma^{2}$ predicted by SuperSCRAM which directly integrates the matter power spectrum. The variance_demo.py module also enables convergence testing our basis decomposition. In general the decomposition converges faster when the volume of the geometry is increased, or when $z_{\text{max}}$ is decreased.

The script super_survey_tests.py does a full run of the sequence necessary to run SuperSCRAM, including various consistency checks on the calculated eigenvalues.

In addition to the testing modules described here, we have verified by conducting multiple runs that the results converged to ${\cal O}(0.01\%)$ with respect to all the various parameters, except for the number of basis elements.

Appendix G Modifying the Matter Power Spectrum

Motivated by the procedure in ref. [63], for a given $w(z)$, the w_matcher module calculates an effective equation of state for dark energy $\mathcal{W}(z)$ which represents, at a given $z$, the constant $w_{\text{eff}}$ which reproduces the same comoving distance to last scattering,

[TABLE]

where $E(z)=H(z)/H_{0}$ and $z_{\text{lss}}$ is the redshift of last scattering. Currently, the WMatcher class precomputes the left-hand side of eq. (G.1) on a grid of possible $w$ and $z$ values and interpolates to match the right-hand side. Then, the amplitude of the $z=0$ linear matter power spectrum must be rescaled using $P_{\text{lin}}(k,z=0)\rightarrow\mathcal{G}^{2}(z)P_{\text{lin}}(k,z=0)$ where $\mathcal{G}^{2}(z)$ is obtained from

[TABLE]

Note that this condition is equivalent to eq. (2.3) in ref. [63], but this form is clearer in parametrizations where $\sigma_{8}$ is not a fixed parameter. The w_matcher module matches this condition by precomputing a grid of possible $G(z,w)$ values and interpolating for a given $\mathcal{W}(z)$. Once $\mathcal{W}(z)$ is evaluated, the model is treated as having that constant $w$ value in all respects for calculations involving that $z$ value; for example, we must obtain a new linear power spectrum, because the small $k$ transfer function depends on $w$.

Appendix H Separate Universe Response

The response of a the matter power spectrum to a long-wavelength density fluctuation can be approximated by treating the overdense region as a ’separate universe’, which can be used to find the response of an observable to density fluctuations. The response of the linear matter power spectrum to a fluctuation $\bar{\delta}(t)$ is [79, 80, 81, 34]

[TABLE]

The one-loop correction simply picks up an additional factor of $\frac{26}{21}$,

[TABLE]

and for the Halofit power spectrum, we follow the prescription in ref. [79] by absorbing the factor of $\widetilde{G}(\tilde{a})=\left(1+\frac{13}{21}\bar{\delta}(t)\right)G(a)$ into $\tilde{\sigma}_{8}=\left(1+\frac{13}{21}\bar{\delta}(t_{0})\right)\sigma_{8}$, so that

[TABLE]

Because at present we only consider $w(z)$ models as perturbations around $w=$constant models, we can ignore any correction to these terms due to a variable $w(z)$. The response of the shear-shear lensing power spectrum can be approximated by plugging the response of the separate universe matter power spectrum into eq. (B.4):

[TABLE]

The response of the observable to fluctuations in our basis, $\frac{\partial C^{AB}(\ell)}{\partial\delta^{\alpha}}$ can then be calculated using the chain rule as described in subsection A.1, obtaining:

[TABLE]