text:
# Mining Weak Lensing Surveys

## 1 Introduction

Mapping the matter distribution of the universe is one of the principal aims of cosmology. The traditional approach to this problem has been the use of galaxy surveys, the 2dF Galaxy Redshift Survey \[Colless et al. (2001\] and the Sloan Digital Sky Survey \[York et al. (2000\] being the most relevant examples today. However, galaxy surveys only map the luminous matter in the universe; generalising to all forms of matter requires the additional assumption that the luminous matter faithfully traces the total matter distribution. On the other hand, weak lensing, or the coherent distortion of the shapes of background galaxies by intervening matter, requires no such assumption and is emerging as a powerful tool to map the matter in the universe.

The possibility of lensing by large scale structure (LSS) was first pointed out in pioneering work by ?) and has since been theoretically and numerically studied by a number of authors (\[Blandford et al. (1991\], \[Miralda-Escude (1991\], \[Kaiser (1992\], \[Bernardeau et al. (1997\], \[Kaiser (1998\], \[Jain & Seljak (1997\], \[Wittman et al. (2000\], \[Jain et al. (2000\], \[Wittman et al. (2000\]). However, detecting this “cosmic shear” had to wait for advances in imaging technology and has only recently become possible. There now are a number of detections (\[Bacon et al. (2000\], \[Van Waerbeke et al. (2000\], \[Rhodes et al. (2001\], \[Hoekstra et al. (2002\]) by various groups and with more observations in progress or planned, this number will continue to grow. The next generation of lensing observations such as the NOAO deep field, the CFHT legacy survey and the Deep Lens Survey will go beyond simply detecting cosmic shear but will map it over large areas of the sky. Such large scale surveys will herald in an era of precision cosmology for weak lensing. At the same time, these surveys bring with them the same computational challenges that current CMB experiments and galaxy surveys are facing. The number of background galaxies is $`N_{gal}`$ 105 to 106, making any brute force approach prohibitively expensive. Analysing these surveys therefore requires the development of algorithms that make the problem computationally tractable. Furthermore, the algorithms must be able to handle aspects of real data including noise, incomplete sampling, arbitrary cuts and so on.

This work has two principal goals; the first is the development of an algorithm that solves the computational challenges of large $`N`$ surveys, where $`N`$ can be number of galaxies or just pixelized intensity (as for the CMB). We then test this algorithm (described in detail in the Appendix) on various applications of weak lensing surveys using simulated data. The paper therefore serves the dual purpose of being a test of our algorithm as well as a survey of the potential of weak lensing surveys in cosmology.

We focus on three applications of weak lensing that we believe are the most interesting for cosmology. The first of these, reconstructing the matter distribution, is the basic goal of large lensing surveys. Lensing measures the integrated line of sight matter distribution, or the convergence ($`\kappa `$) map. Since this reconstruction includes all matter, luminous and dark, such a map gives us a glimpse of the true distribution of matter in the universe. Such maps could also be combined with the high quality imaging data of the lensing survey to understand the relationship between luminous and dark matter. Weak lensing reconstructions have been considered under two broad contexts, cluster mass reconstructions (\[Kaiser & Squires (1993\], \[Kaiser et al. (1995\], \[Seitz & Schneider (1995\]) and LSS reconstructions \[Seljak (1998\]. Since the surveys we are considering in this paper are large field surveys, we will focus on methods for the latter, including effects of non uniform sampling, irregular boundaries and noise.

The next scales of interest are the largest gravitationally bound systems in the universe, clusters of galaxies. As clusters are the rarest and most massive of all structures, they provide a sensitive probe of structure formation and the initial conditions of the universe. One such test is the number density of clusters \[Bahcall & Fan (1998\] as a function of mass and redshift, which depends sensitively both on the total matter in the universe as well as the presence of dark energy and its properties. In order to use such tests, one needs to have a large catalogue of clusters, with well defined selection criteria and completeness fractions. One approach is to use large photometric or spectroscopic surveys of galaxies \[White & Kochanek (2001\] to construct such catalogues; lensing provides us with an alternative approach (\[Hennawi et al. (2001\], \[White et al. (2001\]). An advantage to lensing is that it detects all massive structures, whether they are luminous or not, while other methods assume the presence of luminous matter. Since theory predicts the total number of such objects, lensing allows for the most direct comparison to predictions. Also, any possible existence of “dark clusters” (e.g. \[Miralles et al. (2002\]) would bias the latter methods, while lensing surveys would be unaffected. However, lensing has its own disadvantages that must be understood, if it is to be reliably used.

This paper develops a method to detect clusters in lensing surveys, as well as addressing the theoretical limitations of lensing for cluster searches. It then characterises the completeness and reliability (defined below) of such surveys. We also examine the potential of such surveys to constrain the profile of the dark matter haloes associated with such clusters. This work parallels and extends the work in ?), although we use different methods.

Finally, we address the measurement of the convergence power spectrum. The power spectrum is one of a host of statistics that can be measured, but it has emerged as the statistic of choice for cosmology, since many models predict a gaussian random field distribution on large scales, for which the power spectrum contains all the information. Also, the power spectrum has a natural interpretation from linear structure formation, and therefore is a powerful probe of cosmology, especially in conjuction with other astrophysical measurements (\[Jain & Seljak (1997\], \[Kaiser (1998\], \[Hu & Tegmark (1999\]). The techniques for estimating the power spectrum are similar to those already employed in CMB and galaxy redshift surveys. Here we adapt them to weak lensing.

The paper is organised as follows: We start by describing the features of our algorithm. In $`\mathrm{\S }3`$, we review the basic formalism of weak lensing, while $`\mathrm{\S }4`$ describes the simulations used in the paper. We then examine each of our chosen weak lensing applications, image reconstructions ($`\mathrm{\S }5`$), cluster detections ($`\mathrm{\S }6`$), and power spectrum estimation ($`\mathrm{\S }7`$ ). We conclude in $`\mathrm{\S }8`$. The numerical algorithms used in this paper are described in the Appendix.

## 2 The Algorithm: Features

As mentioned in the introduction, large weak lensing surveys have the same computational challenges that the next generation of CMB and galaxy survey analyses will face. Maximizing the weak lensing signal involves a high number density of background galaxies distributed over large areas of the sky, implying $`N_{gal}`$ 100,000 galaxies. All the algorithms described in this paper are written in terms of matrix operations111Most modern cosmology algorithms are cast in this form., and the computational challenge is in implementing these $`N_{gal}\times N_{gal}`$ matrix operations efficiently. The goal of this paper is to show that the challenge is tractable.

The actual numerical method is described in the Appendix; we discuss the various features of the algorithm here.

1. The algorithm works in real space: There are a number of reasons why a real space algorithm is preferable to a harmonic space approach:

1. The noise properties of data are typically trivially representable in real space.
    2. Real data is not uniform, but involves a number of cuts that arise from bad data, star removal, cosmic ray subtraction etc. In real space, these cuts are trivially included. In harmonic space, these cuts would introduce artifacts due to the non uniformity of coverage, and these are generically difficult to correct for.
    3. Since we do not work directly with Fourier modes there is no need to impose periodic boundary conditions. Thus, there are no issues associated with aliasing of power from the modes on the scale of the survey.
2. There is no pixelisation required: A particularly simple solution to the computational challenges of large data sets is to pixelise the data if possible. There are a number of methods to do this, and one can choose pixels that maximize the signal given a fixed number of pixels. Pixelisation however does involve some degree of data loss at all scales, and no information is preserved below the pixel scale. While this may be acceptable for image reconstruction and power spectrum measurements on large scales, it is undesirable for cluster searches, where small scale information is essential. Furthermore, pixelisation is particularly problematic if the survey geometry is complicated.
3. The time complexity is $`N_{gal}\mathrm{log}N_{gal}`$ : We have tested our algorithm in the range $`N_{gal}=`$ 30,000 to 180,000 using a single processor on a workstation. Within this range, the computational time scales (effectively) linearly with the number of galaxies; the processing time ranged from a few minutes for the image reconstructions and halo searches, to a few hours for the Fisher matrix estimation. We also note that the algorithms in this paper, especially the power spectrum estimation, are readily parallelized.
4. The algorithm is generalisable : The algorithm does not use any properties that are unique to weak lensing; most of the methods described here have either been borrowed from other applications, or else are trivially generalisable. Consequently, this algorithm can be used in other applications where large $`N`$ is a computational restriction.

## 3 Lensing Formalism

Let us begin by briefly reviewing the basic formalism of weak lensing to establish notation and our conventions. The reader is referred to ?) for more details; our treatment closely parallels the discussion in ?).

The gravitational deflection of light can be described as a mapping between a source (S) and image (I) plane. The mapping can be written as

$$\delta x_i^S=A_{ij}\delta x_j^I,$$
(1)
where $`\delta 𝐱`$ is the displacement vector between two points on a given plane. In the weak lensing regime, the mapping has the form

$$A_{ij}=(1\kappa )\delta _{ij}\gamma _1\sigma _3\gamma _2\sigma _1,$$
(2)
where $`\sigma _i`$ are the Pauli matrices, $`\kappa `$ is the convergence, and $`\gamma `$ is the two component shear. Naively, one might expect that this mapping depends on three independent parameters. However, the three components of $`A_{ij}`$ are not independent; the relation between them most easily expressed in Fourier space (we will assume a small enough patch of sky in this paper to ignore curvature effects),

$$\gamma _1(𝐥)=\kappa (𝐥)\mathrm{cos}(2\varphi _𝐥);\gamma _2(𝐥)=\kappa (𝐥)\mathrm{sin}(2\varphi _𝐥),$$
(3)
where $`\varphi _𝐥`$ is the direction of the $`𝐥`$ mode. Within the weak lensing approximation, the expectation value of the ellipticity is proportional to the shear. The proportionality constant depends on the definition of the ellipticity; we adopt

$$e=\gamma .$$
(4)
We now compute the various lensing two point statistics. To do so, we Fourier decompose the shear field into the so-called $`E`$ and $`B`$ modes,

$`\gamma _1(\widehat{𝐧})={\displaystyle \frac{d^2l}{(2\pi )^2}[E(𝐥)\mathrm{cos}(2\varphi _𝐥)B(𝐥)\mathrm{sin}(2\varphi _𝐥)]e^{i𝐥.\widehat{𝐧}}}`$
$`\gamma _2(\widehat{𝐧})={\displaystyle \frac{d^2l}{(2\pi )^2}[E(𝐥)\mathrm{cos}(2\varphi _𝐥)+B(𝐥)\mathrm{sin}(2\varphi _𝐥)]e^{i𝐥.\widehat{𝐧}}}`$ (5)
The two point statistics of these quantities are specified by their power spectra,

$`E(𝐥)E(𝐥^{})=(2\pi )^2\delta (𝐥𝐥^{})P_l^{EE}`$
$`B(𝐥)B(𝐥^{})=(2\pi )^2\delta (𝐥𝐥^{})P_l^{BB}`$
$`E(𝐥)B(𝐥^{})=(2\pi )^2\delta (𝐥𝐥^{})P_l^{EB}.`$ (6)
We note that weak lensing by density perturbations only produces $`E`$ modes, while shot noise and systematic effects can produce both. In terms of power spectra, $`P_l^{EE}=P_l^{\kappa \kappa },P_l^{EB}=P_l^{BB}=0`$ for weak lensing, while $`P_l^{EE}=P_l^{BB},P_l^{EB}=0`$ for shot noise. Systematic effects can produce all three power spectra. For the rest of this paper, we will ignore $`P_l^{EB}`$ since it is parity violating; we refer to $`P_l^{EE}=P_l^{\kappa \kappa }`$ as the $`E`$ mode power, while $`P_l^{BB}`$ is the $`B`$ mode power.

Using the above expectation values and standard trigonometric and Fourier identities, we calculate the various covariance matrices,

1. $$S_{(ij)(ab)}^{\gamma \gamma }=\gamma _{ia}\gamma _{jb}=\frac{ldl}{4\pi }\underset{X=ϵϵ,\beta \beta }{}P_l^XI_{ab}^X,$$
(7)
with

$`I_{ab}^{ϵϵ}=\left[\begin{array}{cc}J_0+c_4J_4& s_4J_4\\ s_4J_4& J_0c_4J_4\end{array}\right]`$ (10)
$`I_{ab}^{\beta \beta }=\left[\begin{array}{cc}J_0c_4J_4& s_4J_4\\ s_4J_4& J_0+c_4J_4\end{array}\right].`$ (13)
2. $$S_{(ij)(a)}^{\kappa \gamma }=\frac{ldl}{4\pi }P_l^{\kappa \kappa }I_a,$$
(14)
with

$$I_a=(2c_2J_2,2s_2J_2).$$
(15)
3. $$S_{ij}^{\kappa \kappa }=\frac{ldl}{4\pi }P_l^{\kappa \kappa }J_0.$$
(16)

For convenience, we have introduced the shorthand notation $`c_n=\mathrm{cos}(n\varphi ),s_n=\mathrm{sin}(n\varphi ),J_n=J_n(l\theta )`$. The displacement vector between two points $`(i,j)`$ is described in polar coordinates by $`(\theta ,\varphi )`$; the components of the shear are $`(a,b)`$. We will often omit the superscripts for notational simplicity when referring to these matrices if they are implicit in the context.

## 4 Simulations

The simulated convergence maps used in this paper are derived from the Virgo N-body simulations \[Jenkins et al. (1998\], run at the Edinburgh Parallel Computing Center and the Computing Centre of the Max Planck Society at Garching. These simulations use a parallel AP3M code (\[Couchman et al. (1995\],\[Pearce & Couchman (1997\]) to evolve 2563 particles from z=50 to z=0. The cosmological parameters used in the simulation are in Table 1. The dark matter distribution is then projected onto a series of planes in the redshift range z=2 to z=0. Dark matter haloes with a mass greater than $`10^{14}`$ M are identified with a FOF group finder with b = 0.2, forming the halo catalogue we use in the rest of the paper.

We compute the convergence map from these planes by using discrete multiple plane lensing \[Schneider et al. (1992\],

$$\kappa =\frac{3H_0^2}{2}\mathrm{\Omega }_M\underset{i=1}{\overset{N}{}}g_i\frac{\delta _i}{a_i}\mathrm{\Delta }\chi _i,$$
(17)
where $`\delta _i`$ is the density perturbation at plane i and $`g_i`$ is the geometrical distance ratio, given by

$$g(\chi ,\chi ^{})=\frac{r(\chi ^{})r(\chi \chi ^{})}{r(\chi )},$$
(18)
with $`\chi `$ the comoving distance and $`r(\chi )=\mathrm{sin}_K(\chi )`$ the usual curvature distance. Two sets of convergence maps are created, 3.5$`{}_{}{}^{}\times `$ 3.5 maps assuming all the sources are at z=1 and 2$`{}_{}{}^{}\times `$ 2 maps using the redshift distribution,

$$n(z)z^2\mathrm{exp}(z/z_0),z_0=0.4,$$
(19)
to weight the planes in the sum. A representative $`1^{}\times 1^{}`$ subsection is shown in Fig. 1. In order to create “independent” maps, the origin of the planes is randomly chosen and we restrict ourselves to using $`1^{}\times 1^{}`$ subsections to create the galaxy catalogues. While the maps are independent on small scales, this breaks down at larger scales, which manifests itself in the high precision numerical experiments of Sec. 6.

In order to generate galaxy catalogues, we create shear maps by Fourier transforming the $`\kappa `$ maps and using Eq. 3. In order to eliminate artifacts from the periodic boundary conditions of the FFT, only the central $`1^{}\times 1^{}`$ of the shear map is used. We randomly assign galaxies positions, and bilinearly interpolate the shear map and use Eq. 4 to compute their ellipticities. In addition, each galaxy is given a random ellipticity drawn from a Gaussian distribution with $`\sigma =0.4`$. For the rest of the paper, we refer to the intrinsic ellipticities of galaxies as “noise”. In general, the noise would also contain measurement errors, but we do not simulate these. Note that our method simulates the effect of non-uniform coverage due to the clustering of background sources, although we limit ourselves to Poisson clustering. Table 2 summarizes the different catalogues used in the various sections of the paper.

## 5 Image Reconstruction

One of the basic aims of a weak lensing survey is the reconstruction of the 2D convergence field. A particularly simple nonparametric method for this is Wiener Filtering (WF) \[Seljak (1998\]. In the case that the data are gaussian distributed, the WF estimator coincides with the maximum posterior probability estimator \[Zaroubi et al. (1995\] and is therefore, optimal. Even when the data are not gaussian, this estimator still minimizes the variance (as defined below) among all linear estimators; however, it is no longer guaranteed to be optimal. For LSS reconstructions, the deviations from gaussianity on large scales are expected to be small, and so WF is, in an appropriate sense, the best that one can do.

### 5.1 Theory and Implementation

To derive the estimator, let us organise the ellipticities of N galaxies into a 2N - vector, $`𝐞`$. The estimator for the convergence $`\kappa `$, at M points can be written as $`\widehat{\kappa }=𝚽𝐞`$. We wish to minimize the variance,

$$(\kappa 𝚽𝐞)^t(\kappa 𝚽𝐞)$$
(20)
with respect to $`𝚽`$. This gives

$$𝚽=\kappa 𝐞^𝐭\mathrm{𝐞𝐞}^𝐭^1=S^{\kappa \gamma }(S^{\gamma \gamma }+N)^1$$
(21)
where $`S^{ab}`$ are the covariance matrices of Sec. 2, and $`N`$ is the noise covariance matrix. We note that $`S^{\gamma \gamma }+N`$ is a 2N $`\times `$ 2N matrix while $`S^{\kappa \gamma }`$ is an M $`\times `$ 2N matrix.

In addition, it is essential to generate an error estimate for the reconstruction. We do this by generating mock catalogues created by randomizing the galaxies’ ellipticities and reconstructing the $`\kappa `$ map for each of these catalogues. The variance over $`N_{MC}`$ (=200 in this paper) maps is a measure of the error of the reconstructed convergence field. It is important to note that in generating the mock catalogues, the galaxy positions are not altered, thereby explicitly taking into account the non-uniform (here Poisson) sampling of the data.

Using WF requires knowing the exact power spectrum (in the covariance matrices) for it to be optimal. For this one can just use the estimate of the power spectrum from the data using the methods of Sec. 6 \[Seljak (1998\]. However, empirical tests demonstrate that the reconstruction is relatively insensitive to the power spectrum used if it is roughly approximates the true one, even though the estimator is no longer strictly optimal.

Finally, on the implementation of this algorithm and others in the paper; first, it is not necessary to estimate $`\kappa `$ at the positions of the galaxies. Indeed, it is more useful to estimate it on a uniform grid; all the figures in this paper are reconstructions on a 2562 or 5122 grid. Second, while the estimator is theoretically simple, implementating the required matrix operations numerically for $`10^5`$ galaxies and greater, is more of a challenge. There are a number of possible solutions to this; we refer the reader to the appendix for our implementation.

### 5.2 Results

The WF reconstructions of the field (Fig.1) for a variety of background densities of galaxies and realistic noise are in Fig. 2, while Fig. 1 shows the reconstruction in the case where the noise has been reduced by a factor of 100. The reconstructions resemble smoothed versions of the original map. This can be understood by considering the WF operator in Fourier space assuming diagonal noise. In Fourier space, the covariance matrices are diagonal with $`S^{\kappa \gamma }S^{\gamma \gamma }P_l^{\kappa \kappa }`$. The WF operator then simply weights each Fourier mode by $`P_l^{\kappa \kappa }/(P_l^{\kappa \kappa }+\widehat{\sigma }^2)`$, where $`\widehat{\sigma }^2`$ is the amplitude of the noise power spectrum; i.e. every mode is simply weighted by signal/(signal+noise). Considering Fig.3, we see that the principal effect of the WF is to act as a low pass filter, removing modes greater than $`l`$ 1000. Reducing the level of the noise increases the cutoff frequency of this low pass filter, as can be seen both from the images as well as the power spectra.

The advantages of implementing the WF in real space are evident from the lack of artifacts at the edges of the reconstructions. The presence of an edge is irrelevant to a real space algorithm, and it correctly reconstructs structures at the edge of a field. This is not true for harmonic space approaches which are sensitive to the entire field, and generically produce artifacts at edges.

For completeness, we also considered the reconstructions of fields assuming a redshift distribution of sources and verified that within our implementation of the redshift distribution, there is no effect on the reconstruction. This will continue to be true for more realistic implementations of the redshift distribution except in the limit of sparse sampling (low number of background galaxies). However, WF has only limited value in the sparse sampling limit as is evident from the panels of Fig.2.

## 6 Clusters

We turn to the construction and characterization of halo catalogues from weak lensing data. Although weak lensing, being sensitive only to the lensing mass, would appear to be the ideal method for constructing such catalogues, it has its own theoretical limitations. We examine these first, then turn to the construction of an optimal filter to detect clusters, and test it. Finally, we address the issue of constraining cluster profiles from weak lensing data.

### 6.1 Theoretical Estimates

In order to understand the limitations of weak lensing halo searches, we measure the expected signal from our catalogue of dark matter haloes, with no noise added. We define the signal as the ratio $`M_{est}/M_{true}`$, where $`M_{est}`$ is the lensing mass estimated from the convergence map, while $`M_{true}`$ is the halo finder mass. The lensing mass is estimated by integrating $`\kappa `$ outward from the halo centre upto a fixed fraction of the virial radius222To get the mass, one must multiply this by $`\mathrm{\Sigma }_{\mathrm{crit}}`$. Equivalently, one could scale the halo finder mass by this same ratio., while the halo mass is determined by scaling the halo finder mass assuming an isothermal density profile. The results are shown in Fig. 4. The slight trend in the mass ratio decreasing with radius is a result of clusters being steeper than isothermal in the outer parts of the cluster, such as in NFW profile.

A feature of this plot is the presence of negative masses, i.e. approximately 5% of the haloes are undetectable. These negative masses result from the fact that lensing measures the line of sight integrated mass; low mass haloes can therefore be masked by underdensities. The converse, low mass haloes being masked by heavier haloes, also occurs and is responsible for the long tail of Fig. 4. Removing line of sight contamination (Fig. 5) by considering each lensing plane individually (Eq. 17) removes both the negative masses and the tail, verifying our interpretation. This figure shows that the lensing masses are systematically greater than the halo finder masses. This is the combination of two effects; the true profiles are not isothermal, and so using an isothermal profile to scale the halo finder masses will underestimate the mass. Also, the lensing masses are 2D projected masses, while the halo finder measures the 3D distribution.

Another related concern are filaments oriented along the line of sight that could masquerade as haloes. In order to determine the degree of contamination, we consider all pixels in our $`\kappa `$ maps that exceed a particular threshold (we tried a few different thresholds, but our results were insensitive to the particular choice). We then compute the distance, $`\mathrm{\Delta }\theta `$, between the pixel and the closest halo. This is a measure of the correlation between the $`\kappa `$ overdensities and the haloes. We find that all but a negligible fraction of overdensities did not correspond ($`\mathrm{\Delta }\theta >1`$ arcmin) to real haloes, implying that contamination due to filaments is insignificant.

### 6.2 Cluster Detection: matched filters

The problem of halo detection can be formally stated as follows - given an input signal with a spatial distribution $`f(𝐱)`$ and amplitude $`A`$, the measured signal can be written as

$$f^{}(𝐱)=Af(𝐱)+N(𝐱)$$
(22)
where $`N(𝐱)`$ is the noise. Assuming that the noise has zero mean and and its statistical properties are spatially independent, the minimum variance estimator, $`\mathrm{\Phi }`$, has the form (in Fourier space) \[Haehnelt & Tegmark (1996\],

$$\mathrm{\Phi }(𝐤)\frac{f(𝐤)}{N(𝐤)}.$$
(23)
Specializing to the case of lensing, the estimator is given by

$$\mathrm{\Phi }\frac{S^{\kappa \gamma }}{N}.$$
(24)
The normalization is usually determined by requiring that it be unbiased. However, since we are interested in the ratio of the signal to the noise, the normalization cancels and we leave it arbitrary. The signal matrix, $`S^{\kappa \gamma }`$, has the same form as the unprimed matrix of equation 14, except that the input power spectrum is replaced by the Fourier transform of the halo profile.

We now describe our halo detection algorithm. We start by convolving the shear map with the matched filter; in practice, this involves multiplying the data (organised into a vector) by the matrix in Eq. 24. Let us denote this convolved map by $`\stackrel{~}{\gamma }(𝐱)`$ \- we must now determine whether the value of $`\stackrel{~}{\gamma }`$ at a given point determines a cluster or not. In order to do this, we start with the null hypothesis, i.e. there is no halo at $`𝐱`$ and ask whether $`\stackrel{~}{\gamma }(𝐱)`$ is consistent with this. Under the null hypothesis, the expectation value of $`\stackrel{~}{\gamma }(𝐱)`$ is zero with variance $`\sigma ^2`$. We define a cluster detection if $`\stackrel{~}{\gamma }(𝐱)>n\sigma `$, where $`n`$ is some chosen threshold. We compute $`\sigma `$ as in the previous section, by randomising the galaxy ellipticities while keeping their positions fixed to create 200 noise maps, and measuring the variance of the resulting convolved maps.

We must emphasize a number of points at this stage. The first is that for any given threshold, we expect to have points not associated with a halo to have $`\stackrel{~}{\gamma }(𝐱)>n\sigma `$; this fraction of false detections will decrease with increasing $`n`$. If the noise is Gaussian, then the fraction is known analytically; however, as we will see in the following sections, the noise is not consistent with being Gaussian. Indeed, the principal sources of noise come from extraneous structures not associated with the halo; these structures do not constitute a Gaussian random field. It is however possible to calibrate the expected false fraction from simulations and include it in theoretical analyses. We also re-emphasize that when computing the variance, it is important to leave the galaxy positions unchanged and only randomise the ellipticities. Only by leaving the galaxies’ positions unchanged will the signal, which is dependent on the distribution of background galaxies, be properly estimated.

The cluster profile is still a free parameter. A particularly simple choice is the singular isothermal sphere, which in projection is $`\theta ^1`$. This choice has the disadvantage of excessively weighting the outer parts of the halo (the integrated profile is logarithmically divergent), which are both noise dominated and contaminated by external structures. Numerical experiments with it also verify that it is suboptimal.

A simulation motivated choice is the NFW profile in projection (\[Navarro et al. (1997\], \[Bartelmann (1996\]). However, the NFW profile has a $`r^1`$ inner cusp, which is not resolved by our simulations. Considering the NFW profile, convolved with the simulation pixels, suggests a profile of the form \[White & Kochanek (2001\],

$$f(\theta )=\frac{1}{(1+\theta /\theta _s)^2},$$
(25)
where $`\theta `$ measures angular seperation on the sky. This profile has the same aysmptotic $`\theta ^2`$ behaviour of the NFW profile; within the scale radius $`\theta _s`$, it possesses a core. Note that we have introduced an angular scale into the filter, and that there is no way to theoretically choose this scale for all clusters, since they will be at different distances. When analyzing real data one may of course choose a projected NFW profile instead, which however still has a physical scale that cannot correspond to a single angular scale. We resolve this problem by simply using a series of scale radii to perform the reconstructions.

The results for $`\theta _s=`$ 5.0, 1.0, 0.5 and 0.1 arcminutes are shown in Figs. 6 and 7. As is evident from the figures, the scale radius corresponds to a smoothing scale. A larger scaling length results in fewer false detections, but tends to coalesce seperate haloes into single structures and misses a greater fraction of clusters. Decreasing the scale radius resolves seperate structures better, but is also more susceptible to noise and extraneous structures, as can be seen from the number of detected structures that do not correspond to haloes.

### 6.3 Completeness/Reliability

It is useful to recall at this stage that we have two lists of data, a list of haloes from the halo finder and a list of pixels whose halo signal was above a certain threshold. There are two obvious questions one can ask - what fraction of haloes have associated pixels, and what fraction of pixels have associated haloes. We define the completeness as the fraction of haloes that have their closest detected pixel a specified distance away from them. Conversely, the reliability is defined as the fraction of pixels that have their closest halo a specified distance away. We note that the natural measure of distance for completeness is in units of the virial radius, whereas the reliability distance is measured as a physical angular distance.

The completeness fractions, for three different thresholds, are shown in Figs. 8 and 9. We make the following observations about these results:

1. As expected, increasing the threshold reduces the completeness fraction. Note however that this happens at 2-3$`\sigma `$ level, so one cannot choose a high level of significance (thereby rejecting spurious detections with high confidence) and have a high level of completeness at the same time.
2. The width of the completeness distribution is also seen to widen as the scale length decreases. This can be understood heuristically by noting that as the maps become more noise dominated, the point closest to the halo may be unrelated to it, thereby broadening the distribution. A crude modelling of the distribution in the limit of purely random, independently distributed points obtains,

$$f(x)x(1x^2)^n;x=r/r_{vir},$$
(26)
which, for $`n5`$ resembles the distribution for the smallest scale radius333The distribution assumed $`n`$ points distributed within an area of $`\pi r_{vir}^2`$.. Note however that this is meant to be illustrative, and not a close approximation to the actual distribution.
3. A redshift distribution of the sources decreases the completeness fraction. This is the result of two competing effects - using a redshift distribution allows one to probe a larger volume of space, and one would hope to detect more clusters in that region. On the other hand, the redshift distribution results in a more non uniform signal which reduces the signal.

The corresponding reliability fractions are shown in Figs. 10 and 11, and the fraction of detected pixels that had no corresponding halo within 5 arcminutes is in Table 3. Note that we have only shown the fractions for the simulations without redshift distributions; including a redshift distribution produces similar results. We note the following:

1. The reliability fraction becomes more peaked as the threshold increases.
2. The above trend becomes more pronounced as the scale radius decreases, implying that a higher threshold is necessary to eliminate noise at these scales. However, the number of detected pixels at both high thresholds and small scales also steeply decreases, again pointing to intermediate scales for optimal cluster detection.
3. Unlike the completeness, the reliability appears to be insensitive to a redshift distribution.
4. Finally, the dropoff near the origin is a purely geometrical effect, due to the fact that the differential area grows linearly with radius.

It is also useful to consider completeness as a function of halo mass, we do this in Figs. 12 and 13. The completeness plotted is an integrated quantity - the fraction of haloes that have a detected pixel within 0.2 $`r_{vir}`$ of them. The completeness is a steep function of halo mass, with the surveys being nearly 100% complete at the high mass end, but less than 5% complete for the lowest masses. The noisy nature of these plots is simply an artifact of the fact that our simulations have few massive haloes.

An important observation to make is that the achieved completeness is significantly lower than what one might have expected theoretically. This discrepancy can be traced back to the assumption of a uniform signal in our theoretical estimates. The Poisson clustering of the background galaxies does not satisfy this assumption; indeed, configurations of background galaxies make certain haloes undetectable. This would be even worse if one included clustering of the background galaxies, and must be taken into account when computing predicted cluster finding efficiency as a function of halo mass.

### 6.4 Cluster Profiles

There has been considerable interest recently \[Sheldon et al. (2001\] in using weak lensing to statistically constrain cluster halo profiles. We examine this by considering both the integrated mass profile (Fig. 14) and the tangential shear profile (Fig. 15). The integrated mass profile is simply given by summing our convergence map with no noise added. For the shear profile, we use our simulated catalogues. The ellipticity of each galaxy is decomposed, relative to the halo center, into a tangential and radial component, which we then average in radial bins to get the profile. Unlike the integrated mass profiles, the shear profiles are estimated in the presence of noise.

We divide clusters into three mass bins, and average the profiles within those bins assuming we know their central positions. This simulates the process of “stacking” clusters with centers known from X-ray or optical data. The profiles are seen to be well constrained for the highest mass bin and poorly constrained at low masses.

We can also compare the measured shear profiles to what one would expect if one assumed a convergence profile of the form in Eq. 25. The expected profiles for the four scale radii we use are shown in Fig.16. As we might have anticipated from the results of the two previous subsections, the best agreement is when $`\theta _s`$ is between 0.5 and 1.0 arcminutes. It is also important to emphasise that the dropoff seen in Fig.15 is an effect of the pixelisation of the map; a similar effect is seen in Fig.16 where the core radius mimics the effects of pixelisation. If clusters had significant cores, then just such an effect would be also physically expected.

## 7 Power Spectrum Estimation

In this section we consider the measurement of the convergence power spectrum from weak lensing data. The subject of the optimal measurement of the power spectrum from noisy data has received a lot of attention, with regards to the CMB and galaxy redshift surveys, and we will simply import the techniques that have been developed and validated there to weak lensing. It should be emphasized that weak lensing, unlike galaxy redshift surveys, measures the matter power spectrum directly, eliminating the complications of bias.

Since the formalism for optimal weak lensing power spectrum measurement has been discussed in detail elsewhere (see for eg. \[Seljak (1998\], \[Tegmark et al. (1998\]), we will limit ourselves to a brief discussion (following the notation in \[Padmanabhan et al. (2001\]) to establish the formalism we will be using. A similar discussion that uses a different, although related, approach is in ?).

### 7.1 Theory

Let us parametrize the power spectrum by Np step functions such that $`P(l)=p_i`$ for $`l_{i1}ll_i`$, where $`i`$ ranges from 1 to Np. We can now arrange the $`p_i`$ into an Np-vector, $`𝐩`$. The problem now reduces to estimating $`𝐩`$ from the data $`𝐱`$, where $`𝐱`$ is a 2N vector consisting of the galaxy ellipticities. Define the covariance matrix of the data, $`C=S^{\gamma \gamma }+N`$ where $`S^{\gamma \gamma }`$ is the signal covariance matrix (see Sec. 2) while $`N`$ is the noise matrix. Recall that $`N=N_{meas}+N_{int}`$ where $`N_{meas}`$ is the measurement noise, while $`N_{int}`$ is the intrinsic noise due to galaxy ellipticities. This intrinsic noise is not known a priori and must be estimated from the data, an issue we address at the end of this section.

Since the power spectrum is the sum of step functions, we rewrite the covariance matrix,

$$C=\underset{i=1}{\overset{N_p}{}}p_iC_i+N=\underset{i=0}{\overset{N_p}{}}C_i,$$
(27)
where $`C_i=S^{\gamma \gamma }/p_i`$ for $`1iN_p`$, and we have defined $`C_0=N`$ and introduced a dummy parameter, $`p_0=1`$ for notational convenience.

We now form the minimum variance quadratic estimators (\[Hamilton (1997a\], \[Hamilton (1997b\], \[Tegmark (1997\]),

$$q_i=\frac{1}{2}𝐱^tC^1C_iC^1𝐱$$
(28)
and group them into an Np+1 vector, $`𝐪`$. The quadratic estimators have the following properties,

$`𝐪=F𝐩`$
$`𝐪q^t𝐪𝐪^t=F`$ (29)
where $`F`$ is the Fisher information matrix (\[Tegmark et al. (1997\] and refs. therein),

$$F_{ij}=\frac{1}{2}tr\left[C^1C_iC^1C_j\right].$$
(30)
The best estimate of the power spectrum can be written as

$$\widehat{𝐩}=M(𝐪𝐛)$$
(31)
where $`𝐛=F_{i0}`$, and we restrict all indices to run from 1 to Np. The matrix $`M`$ is an arbitrary N$`{}_{p}{}^{}\times `$ Np matrix with the property that the rows of $`MF`$ sum to unity and is chosen such that our power spectrum estimates have certain desirable properties. The choices for $`M`$ that we consider and their properties are listed in Table 4. Note that there is no gain in information in this final step; the information content is the same as what is in the quadratic estimators and the entire Fisher matrix. However, if one only presents the diagonal of the Fisher matrix (eg. errorbars in a plot), then it is important to ensure that the errorbars are uncorrelated (the third choice). Our estimates $`\widehat{𝐩}`$ now have the following properties,

$`\widehat{𝐩}=W𝐩=MF𝐩`$
$`\widehat{𝐩}\widehat{𝐩}^t𝐩𝐩^t=MFM^t,`$ (32)
where we note that the rows of $`W`$ have the obvious interpretation as window functions.

#### 7.1.1 Choosing a prior

In the preceding discussion, we have implicitly assumed a prior power spectrum, and all the properties have rested on this prior being the true power spectrum. This begs the question - how does one choose the prior? Any reasonable guess to the power spectrum works as an initial guess. The power spectrum estimated then can be used as the prior, and the process can be iterated. Several authors (\[Bunn (1995\], \[Padmanabhan et al. (2001\] ) have explicitly verified that the choice of prior does not bias the result, and therefore, in practice, only one or two iterations are required.

#### 7.1.2 Measuring the intrinsic ellipticity scatter

The intrinsic scatter in the ellipticities can be estimated by simply computing the r.m.s. value of the ellipticity of the data. This will produce the correct answer so long as the shear correlation length is smaller than the size of the field considered. However, if the scatter is underestimated, then the extra power will manifest itself as excess power in the $`E`$ and $`B`$ modes. In order to correct for this, we follow ?) and introduce a shot noise parameter, $`p_{noise}`$, to the power spectrum, whose contribution to the power spectrum is $`C_{noise}=p_{noise}\delta _{ij}`$. As with the power spectrum, the excess power measured this way can be corrected in the next iteration.

### 7.2 Results

The Fisher matrix formalism is only exact for Gaussian random fields, and so, in order to test it, we simulate 200 gaussian random convergence fields with a known $`E`$ mode power spectrum, and no $`B`$ mode power. Figs. 17 and 18 summarize the results. The triangles are the input power spectrum, while the heavily shaded boxes are the input power spectrum convolved with the window functions. The width of the boxes are given by the Fisher matrix error estimates. The lightly shaded boxes are the measured power spectrum, averaged over 200 realizations, with the width of the box given by the standard deviation over the realizations. If the method is correct, then both the power spectrum and the errors should agree. The agreement is best for the $`E`$ mode power, where the input power was non-zero. We note that although the input $`B`$ mode power was zero, we expect to see power measured due to leakage from the $`E`$ modes through the window functions. This is significant only for the leftmost band, that probes scales larger than the survey444The modes were included to prevent aliasing of power.. We observe that the remaining $`B`$ mode power is consistent with zero, as expected.

Figs. 19 and 20 show the expected and measured power for 200 realizations drawn from N-body simulations, which allow us to test whether the breakdown of the gaussian approximation biases our results. We find that the Fisher matrix underestimates the true errors in the power spectrum. A similar effect was found by ?). Although this is not a large effect, N-body simulations must be used to calibrate the Fisher matrix errors when using the power spectrum estimates to extract cosmological information. The second observation is the discrepancy between the expected and measured $`E`$ mode power in the leftmost bands. This is due to the fact that all our realizations were derived from the same N-body simulations, and although we expect the realizations to be independent on small scales, this will break down at large scales. Fig. 17 verifies this expectation. Note that these bands probe scales larger than those simulated; it is interesting that we correctly estimate the power to within a factor of two in these bins.

All of the above results use the third choice of Table 4 for $`M`$. Figs. 21 and 22 contrast the first and third choices for $`M`$. As expected, the errorbars for the first are smaller, but are correlated, as compared to the third. The effects of decorrelating are pronounced for the $`B`$ mode power; the points are randomly distributed above and below zero, while in Fig. 21, they are clearly correlated. We therefore reemphasize the importance of decorrelating points when visually presenting data. We do not present the results for the second choice since the errorbars are too large to be useful in this case. This is a direct result of the fact that the window functions in this case are delta functions.

## 8 Discussion

Weak lensing is emerging as an important tool in cosmology. One of its principal advantages is that it probes the matter distribution directly, making no assumptions about the dynamical state of the matter. This is desirable both because it eliminates complications of interpretation, but also because it gives us an opportunity to study the physical processes underlying those assumptions.

In this paper, we considered three applications of large weak lensing surveys. We summarise the results here, comparing it to previous work.

### 8.1 Image Reconstruction

The simplest goal of a weak lensing survey is to produce a map of the distribution of matter in the universe. Wiener filtering provides a simple, and almost optimal reconstruction of the matter distribution, and is our method of choice. An important characteristic of Wiener filtering is that it suppresses power on scales with low signal to noise. For the weak lensing maps we consider here, this scale is at $`l1000`$, while there is no power at scales larger than $`l360`$, the scale of the survey. This suggests the obvious design strategy of large area surveys to probe large scales, with high background number densities $`n25`$ to resolve features on small scales. However, Wiener filtering has limited cosmological applicability because of its inability to resolve smaller structures.

### 8.2 Cluster Detection

Clusters, as the most extreme structures in the universe, are a sensitive probe of cosmology. Weak lensing has the advantage that it searches for clusters directly as mass enhancements, independent of the presence of luminous matter. However, in order to compare with theoretical predictions, one must, in addition to compiling a catalogue of clusters, understand the selection criteria and the completeness of the catalogue.

A disadvantage of weak lensing is that it measures the projected mass distribution, and is therefore susceptible to contamination from uncorrelated haloes in the line of sight. This is worst for the low mass haloes, since they can be obscured by either more massive haloes or underdensities. More massive haloes are less sensitive to this effect and the theoretical maximum completeness approaches 100%.

The philosophy that we adopt in this paper is that weak lensing cluster searches will be done in conjunction with other measurements such as galaxy, X-ray or Sunyaev-Zel’dovich surveys. One could either imagine using weak lensing to identify candidate clusters and verifying them with follow up measurements, or correlating different measurements to remove false detections. With this in mind, we develop an optimal filter approach to detecting clusters. Detections are defined as those regions where the measured signal to noise ratio exceeds a certain threshold. We then compute the completeness (the number of haloes with a corresponding detection) and the reliability (the number of detections with a corresponding halo). We observe that completeness is strongly dependent on the mass of the cluster; for cluster masses $`>`$ 6 $`\times `$ 1014 M, our cluster samples were virtually 100% complete independent of the threshold used, implying that one could construct a uncontaminated, yet complete sample, by choosing a high enough threshold. At lower masses, the completeness drops drastically as the threshold in increased. At low thresholds however, the false detection rate rises to $``$ 25%. This is a direct result of the fact that our estimator is sensitive to spurious structure masquerading as a halo at low thresholds. We should however note that even for low thresholds, the completeness for low mass haloes is less than 40%.

The optimal filter defines a scale length that physically corresponds to a smoothing scale. Excessive smoothing reduces the number of false detections, since the noise is reduced but it also smooths away or coalesces small sclae structures. As we reduce the smoothing scale, we start to resolve smaller structures but are more susceptible to spurious detections from extraneous structures that don’t belong to any halo. The optimal scale length is $`1^{}`$. However, some caution should be exercised when interpreting this value since it may be affected by pixelisation in the N-body maps.

This work parallels a similar study by ?), although we use different methods to detect clusters. ?) use both a simple Gaussian smoothing with a smoothing scale of 1’ - 2’, and aperture mass measures \[Schneider (1996\] with a scale of 1’-5’ and a signal to noise threshold varying from $`S>1`$ to $`S>5`$. As with our methods, the completeness drops with increasing scale radius, with the maximum at $`1^{}`$, although ?) do not consider smaller scale lengths. In addition, they conclude that their catalogues are complete for masses $`>5\times 10^{14}`$ M, consistent with our results.

An extension of this is the inclusion of a redshift distribution of the background sources. Including this distribution does not qualitatively change any of the trends one observes for, but reduces the overall completeness. The halo catalogues, constructed assuming a redshift distribution, are complete for masses $`>8\times 10^{14}`$ M . This is what one might naively expect - the most massive haloes are unaffected by the redshift distribution, while the signal from the less massive haloes is reduced and a larger fraction of them are missed. The reliability is however unaffected by a redshift distribution, because false detections are due to noise and extraneous structure, and one would not expect those to be affected by a redshift distribution.

There are two cautionary lessons that we draw from this work. The first is that the optimal nature of matched filters is easily affected by deviations from ideality. These are caused here by deviations from circular symmetry due to extended haloes, as well as internal and external substructure. The second is that the completeness of cluster surveys are strongly affected by clustering of the background galaxies. Certain haloes are rendered undetectable by the Poisson clustering of our simulated background galaxies. Real background galaxies are not Poisson distributed, but are correlated, and we expect this to have a worse effect. Therefore, all analyses must take into account the clustering of source galaxies to get a representative measure of the expected completeness of the survey.

We note that we have restricted ourselves here to consider the information from lensing only. However, lensing surveys will produce high resolution multi-colour images of the survey region. One might imagine using both the lensing and imaging data simultaneously to detect clusters; such an approach would naturally reduce the number of false detections. Another advantage to using the imaging data jointly with the lensing is that each complements the other - lensing would identify haloes with no optical counterpart, while optical searches could detect objects with low lensing signals.

### 8.3 Power Spectrum

The convergence power spectrum is a measure of the clustering of the recent universe, a regime that was accessible until now only to redshift surveys. In contrast with the galaxy power spectrum, the convergence power spectrum is not affected by the complications of biasing, and provides a direct measure of the matter power spectrum.

The formalism for estimating the power spectrum, borrowed from the CMB and galaxy surveys, assumes that the convergence field is Gaussian distributed. While this is not true for the true convergence field on small scales, it is approximately true on the scales that lensing probes. We have explicitly shown, using N-body simulations, that the Gaussian approximation does not bias the estimators. However, the Fisher matrix underestimates the error bars and although this is not a large effect, it means that the error bars must be calibrated against N-body simulations. Our results are similar to those obtained by ?) although our methods differ in details. The main difference is that our method does not pixelise and so extracts all the information from the data.

In concluding, we observe that the methods presented here are not unique in their application to weak lensing. We have used methods developed for analysis of other data sets and adapted them to weak lensing. New numerical solutions presented here may be adapted to other similar problems in cosmology, particularly those where brute force evaluations are prohibitely expensive. In our application we have managed to reduce an instrinsic $`O(N^3)`$ numerical problem to $`O(N\mathrm{log}N)`$. The same methods can be used in other applications, such as the analysis of CMB temperature and polarization anisotropies and galaxy surveys. We expect that both lensing and other areas of cosmology will benefit from the growing synergy in the field.

We would like to thank SISSA, Trieste, where this work was begun, for its hospitality. We thank VIRGO collaboration for making N-body simulations available to the community. We would also like to thank Joe Hennawi, Catherine Heymans, Christopher Hirata, Kevin Huffenberger and David Spergel for useful discussions. N.P. was supported by a Centennial Graduate Fellowship from Princeton University. U.S. is supported by NASA, NSF CAREER, David and Lucille Packard Foundation and Alfred P. Sloan Foundation. U.P. is supported by NSERC grant 72013704. Computing infrastructure was provided in part by the Canada Foundation for Innnovation PScinet alphacluster.

## Appendix A Numerical Implementations

An important feature of all the algorithms presented in this paper is that they are explicitly written as linear algebra operations. The basic building block of any implementation is therefore a routine to perform matrix-vector multiplication. Unfortunately, the dimensionality of the vector space, $`N`$, is given by the number of data points. A naive application of the algorithms yields $`NN_{gal}100,000`$ for which straight matrix-vector multiplication, an $`N^2`$ process, becomes computationally impractical.

One can approach this problem in two ways, either by reducing $`N`$, or by using properties of the matrices that appear to speed up the vector multiplications. The former approach is equivalent to pixelising the data and a number of pixelisation schemes have been suggested, ranging from direct binning to “optimal” Karhunen-Loéve pixelisations of the data.

We propose and implement an approach that is based on the latter approach. We start by observing that all the matrix operations that we require are of the form $`C𝐱`$ or $`C^1𝐱`$ where $`C`$ is a correlation matrix of $`\mathrm{\S }3`$. The latter operations can be recast as direct matrix-vector multiplications by performing the matrix inversion iteratively (we specify the exact algorithm below). We then observe that the multiplication $`C𝐱`$ is simply a convolution of the data by the appropriate correlation function, $`\zeta `$,

$$(C𝐱)_i=\underset{j}{}C_{ij}x_j=\underset{j}{}\zeta (𝐫_i𝐫_j)x_j,$$
(33)
where the last identity follows from the fact that the correlation function only depends on the seperation between the two points. The problem is now explicitly translationally invariant and one can readily apply the Fourier convolution theorem to perform the vector multiplications efficiently. The asymptotic scaling is now $`O(N\mathrm{log}N)`$ instead of $`O(N^2)`$, making $`N10^5`$ tractable on a workstation.

The careful reader will no doubt point out that the use of Fourier methods requires the data to be uniformly sampled, which our data is not. We solve this problem by resampling the data onto a grid of $`N_{grid}4N_{gal}`$ points where the additional factor of 4 comes from the need to zero pad in 2 dimensions. The exact scaling of the algorithm is therefore $`N_{grid}\mathrm{log}N_{grid}`$. We emphasise that although an auxiliary grid is used, this grid is only an intermediate step which does not impose periodic boundary conditions, and whose discreteness effects can be compensated using multigrid and direct summation methods as described below. Our approach is thus conceptually very different from pixelisation approaches. It is also useful at this stage to point out the obvious analogy between our approach and gravity solvers of PM (particle - mesh) N-body simulations.

While the above approach works well for large scale modes, the pixelisation introduces inaccuracies at scales comparable to the pixel resolution. Our analogy with N-body simulations come to the rescue here; PM simulations have similar inaccuracies on small scales that can be corrected by introducing a direct summation between pairs of particles at small seperations (P3M - particle-particle particle-mesh simulations). We start by splitting the convolution kernel into short and long range pieces,

$$\zeta =\zeta _{long}+\zeta _{short}=f(r)\zeta +g(r)\zeta ,$$
(34)
where $`f`$ and $`g`$ are filter functions with the properties that

$`g=1(r<r_{min}),`$
$`f=1(r>r_{max}),`$
$`f=0(r<r_{min}),`$
$`g=0(r>r_{max}),`$
$`f+g=1(r_{min}<r<r_{max}),`$ (35)
with the definitions of long and short range are determined by $`r_{min}`$ and $`r_{max}`$. The multiplications are now done in two stages; the first is to do the long range piece by the Fourier method described above, while the short range correlations are done by direct summation. Possible filter functions have the form $`1\mathrm{cos}^2(\theta )`$ in the regime $`r_{min}<r<r_{max}`$, where $`\theta `$ is an appropriately scaled length. This form is chosen to minimise the inaccuracies that result from the truncation of the correlation function.

A second scheme to improve the resolution is to use a multi-grid approach. The single-grid scheme described above is only the simplest such implementation. One can trivially generalise it to multiple scales by introducing a series of filters, $`f_i`$. In such implementations, the direct summation is performed only for the shortest range; all other convolutions are done by the Fourier method with the grid becoming coarser for larger scales. Our codes used 3 scales - the innermost scale for direct summation, an intermediate scale with twice the resolution, and a coarse grid for the largest scales.

Note that while the $`O(N\mathrm{log}N)`$ scaling is dependent on a flat space assumption, it can be generalized to the sphere with no significant loss of efficiency. This is because on a sphere only the coarse grid needs to be done with the slower $`O(N_{\mathrm{coarse}}^{3/2})`$ scaling, while the subsequent grids can still be $`O(N_{\mathrm{grid}}\mathrm{log}N_{\mathrm{grid}})`$ as long as the flat sky approximation holds on these grids.

### A.1 Matrix Inversions

The problem we focus on is of the form

$$𝐲=C𝐱=(S+N)𝐱,$$
(36)
where $`𝐲`$ is known, $`S`$ is the signal matrix and $`N`$ is the noise. We will assume that $`N`$ is diagonal in real space here, although the algorithm can be modified for the case of sparse $`N`$. Also, we assume that multiplication by $`S`$ can be efficiently performed by the methods above. The matrix inversion is then performed by an underrelaxed Jacobi iteration,

$$𝐱^{n+1}=𝐱^n+(\overline{S}I+N)^1\left[𝐲(S+N)𝐱^n\right],$$
(37)
where $`I`$ is the identity matrix, and $`\overline{S}`$, the underrelaxation parameter, is equal to $`e_{max}`$, the maximum eigenvalue of $`S`$. This can be estimated using the iterative scheme,

$`𝐱^{n+1}={\displaystyle \frac{S𝐱^n}{𝐱^n}}`$
$`e_{max}=𝐱^n,`$ (38)
which can intuitively be verified by considering the case of $`S`$ diagonal.

The underrelaxed Jacobi iteration is only one of a number of possible approaches to computing the matrix inverses. For our matrices, we obtained a fractional precision of $`10^4`$ in $`100`$ iterations and were therefore not limited by it. However, for ill-conditioned matrices, it might be necessary to go to conjugate gradient or multigrid methods.

### A.2 Trace estimation

The only operation that cannot be trivially written in terms of matrix-vector operations is the computation of the Fisher matrix (Eq. 30) which involves the trace of the product of four matrices. This is an intrinsically $`O(N^3)`$ process. If we assume a random ensemble of vectors, $`𝐱`$, with the property that,

$$\mathrm{𝐱𝐱}^t=I,$$
(39)
we can use the following statistical identity to estimate the trace,

$$tr(A)=tr(AI)=tr(A\mathrm{𝐱𝐱}^t)=𝐱A𝐱^t.$$
(40)
The estimator is completely determined by specifying the ensemble $`𝐱`$. Note that if we choose an ensemble of $`dim(𝐱)`$ orthogonal vectors, then we can exactly recover the trace; however, taking the trace would then be an $`O(N^2\mathrm{log}N)`$ operation. This is already a gain from $`O(N^3)`$ and is achieved by the fast $`O(N\mathrm{log}N)`$ convolution methods of $`C𝐱`$ and $`C^1𝐱`$ operations. However, $`O(N^2\mathrm{log}N)`$ is still slow compared to other operations. One might attempt to modify this by choosing a smaller subset of random, but orthogonal, vectors; we however find that this is slow to converge to the correct value. We obtain best convergence by using the real stochastic $`Z_2`$ esimator, where $`𝐱`$ is a random vector consisting of 1’s and -1’s. For $`N=90,000`$, we measure the trace with a fractional precision of $`10^5`$ with an ensemble of 400 random vectors. This means that we have reduced the scaling from $`O(N^3)`$ to $`O(N\mathrm{log}N)`$ only.