Unfolding the matter distribution using 3-D weak gravitational lensing

Combining redshift and galaxy shape information offers new exciting ways of exploiting the gravitational lensing effect for studying the large scales of the cosmos. One application is the three-dimensional reconstruction of the matter density distribution which is explored in this paper. We give a generalisation of an already known minimum-variance estimator of the 3-D matter density distribution that facilitates the combination of thin redshift slices of sources with samples of broad redshift distributions for an optimal reconstruction. We show how, in principle, intrinsic alignments of source ellipticities or shear/intrinsic alignment correlations can be accommodated, albeit these effects are not the focus of this paper. We describe an efficient and fast way to implement the estimator on a contemporary desktop computer. Analytic estimates for the noise and biases in the reconstruction are given. The bias -- a spread and shift of structures in radial direction -- can be expressed in terms of a radial PSF, comprising the limitations of the reconstruction due to source shot-noise and the unavoidably broad lensing kernel. We conclude that a 3-D mass-density reconstruction on galaxy cluster scales is feasible but, for foreseeable surveys, a map with a S/N>~3 threshold is limited to structures with M200>~10^14 Msol/h, or 7x10^14 Msol/h, at low to moderate redshifts (z=0.1 or 0.6). However, we find that a heavily smoothed full-sky map of the very large-scale density field may also be possible as the S/N of reconstructed modes increases towards larger scales. Future improvements of the method can be obtained by including higher-order lensing information (flexion) which could also be implemented into our algorithm. [ABRIDGED]


INTRODUCTION
Since its detection about a decade ago (Bacon et al. 2000;Kaiser et al. 2000;Van Waerbeke et al. 2000;Wittman et al. 2000) weak gravitational lensing, or cosmic shear -the small, coherent distortion of distant galaxy images due to the intervening large-scale matter distribution in the cosmos -has become a well-established tool for studying the dark Universe (for recent reviews see: Hoekstra & Jain 2008;Munshi et al. 2006;Van Waerbeke & Mellier 2003). The phenomena of gravitationally scattering of light by massive structures is well described by General Relativity. Hence lensing provides us with a direct probe of the mass-distribution in the Universe.
The spectrum of applications for exploiting this pool of information includes the interpretation of correlations of the cosmic shear signal to constrain cosmological parameters, the cross-correlation of cosmic shear with galaxy positions for studying the typical matter environment of galaxies or galaxy biasing, and the fascinating possibility of mapping the large-scale mass distribution in the Universe.
Provided that General Relativity itself is accurate (Huterer & Linder 2007;Uzan & Bernardeau 2001) the statistical analysis of the weak shear signal for cosmological parameter estimation provides us with ample indirect evidence that most of the matter of the Universe is not made up of familiar baryonic matter, but of so-called Dark Mat-ter. Moreover the largest piece in the energy budget in the present Universe is occupied by the so-called Dark Energy which seems to be accelerating the expansion of the Universe (see e.g. Peacock et al. 2006;Albrecht et al. 2006, for a recent review). As the area, A, of weak lensing surveys grows, the statistical uncertainty on the parameters of this cosmological model will decrease as 1/ √ A, allowing even General Relativity itself to be tested. The main issue for lensing surveys for extracting cosmological parameters will be the control of systematics at the 1%-level (Huterer et al. 2006).
As well as statistical analysis, weak lensing can also be used to directly visualise the matter distribution in the Universe. Kaiser & Squires (1993) first showed how to map the projected 2-D mass distribution from weak lensing. Taylor (2001), Bacon & Taylor (2003) and Taylor et al. (2004) showed how to extend this to 3-D mapping by combining shear data with redshift distances (3-D weak lensing). Visualisation of the Dark Matter distribution has the clear advantage that one can directly compare the matter distribution with that of galaxies and gas (e.g. Clowe et al. 2006).
For a given resolution scale, mass mapping does not benefit from a larger area. The noise of a mass map scales as σǫ/ √ N (Kaiser & Squires 1993;Hu & Keeton 2002) where N = nθs is the surface number of lensing relevant galaxies (sources) within a resolved area, θs, the mean surface density of galaxies is n, and σǫ is the intrinsic random ellipticity distribution of sources. It is unlikely that n and σǫ will change compared to current surveys as dramatically as the survey area, while the uncertainty in redshifts must also be added for 3-D mass mapping.
Ambitious survey campaigns are either under way (e.g. CFHTLS, RCS2, COSMOS), about to commence (e.g. Pan-STARRS, VST-KIDS, DES) to cover larger areas of sky, while surveys being planned for the near future (e.g. LSST, Euclid) will map the full sky to great depths and high resolution. This will provide the lensing community with large galaxy catalogues with billions of galaxies including multicolour information for photometric redshifts and accurate galaxy image measurements. These surveys mark the beginning of a new era for lensing as the availability of a large set of sources with redshifts and multi-colour data will allow us to accurately probe the evolution of mass-clustering as a function of cosmic time.
This paper addresses in detail the question of how to optimally reconstruct the 3-D matter density solely based on a lensing catalogue with redshift information. Applying Wiener filtering to cosmic shear follows in the tradition of other successful cosmological applications, initially outlined for astrophysical applications by Rybicki & Press (1992) and applied to galaxy redshift surveys Zaroubi et al. 1995;Webster et al. 1997), CMB maps (Bunn et al. 1994;Bouchet et al. 1999) and peculiar velocity fields (Zaroubi et al. 1999). Here we describe a linear minimum-variance reconstruction operator which is a generalisation of the Wiener filter method proposed by Hu & Keeton (2002). Bacon & Taylor (2003) also proposed a Wiener filtering method, which was first applied in Taylor et al. (2004) for a 3-D reconstruction of the COMBO-17 survey ) and in  to the COSMOS survey. However this was to reconstruct the gravitational potential, rather than the mass-density field, which improved the noise properties of the 3-D reconstruction on cluster scales (see Bacon & Taylor 2003). Hu & Keeton (2002) pointed out that on larger scales the signal-to-noise improves allowing density-reconstruction. Both Hu & Keeton (2002) and Bacon & Taylor (2003) proposed reconstructing the 3-D convergence field in a first step, and then applying the Wiener filter along each individual line-of-sight, ignoring the correlation of errors between different line-of-sights (radial mapmaking). In the application to the COMBO-17 and COSMOS surveys, all of the redshift bins where assumed to have the same noise level, which is sub-optimal.
In this paper, the transition from 3-D shear to lensplane density contrast is done in one step, along all line-ofsights simultaneously, including all correlations between different lines-of-sight and lens planes. Another extension consists in introducing an additional parameter that regulates the impact of the Wiener regularisation. With the approach outlined here we show how all available redshift information, i.e. redshift slices and sub-samples with broad redshift distributions and individual statistical weights of sources, can be combined in one framework. The primary focus of this study will be the signal-to-noise and biases that are introduced in the 3-D reconstruction if we regularise the process with a prior. Regularisation will be unavoidable due to the poor signal-to-noise of a 3-D density map based on lensing alone. Besides the theoretical considerations, we also describe how the reconstruction algorithm can be efficiently implemented in practise.
This paper is laid out as follows. We start off in Sect. 2 with a very general discussion of inverse estimators for linear problems. Sect. 3 continues by specifying the problem of 3-D reconstructions in two steps: the convergence tomography and the lens-plane matter density contrast from the cosmic shear tomography. Both problems fall into the category of linear inversion problems. Sect. 4 discusses the signal-to-noise of mass density reconstructions, the main focus of this paper, and biases therein that are introduced by a prior. As prior we consider expected density correlations in radial or transverse direction. The practical implementation of the estimators, Appendix A, requires some optimisation, because a naive implementation would force us to handle extremely large matrices even on relatively small grids. Finally in Sect. 5, we apply the algorithm to a fictive, simulated survey with specification we should expect for a weak lensing surveys within the next twenty years.
We use as fiducial cosmology a ΛCDM model (adiabatic fluctuations) with Ωm = 0.24, Ω b = 0.0416, ΩΛ = 1 − Ωm and H0 = h 100 km s −1 Mpc −1 with h = 0.732. The normalisation of the matter fluctuations within a sphere of radius 8 h −1 Mpc at redshift zero is assumed to be σ8 = 0.76. For the spectral index of the primordial matter power spectrum we use ns = 0.96. These values are consistent with the thirdyear WMAP results (Spergel et al. 2007).

Estimators for linear problems
All problems discussed in this paper are linear inversion problems. For such problems, a vector of observables, d, is known to be a linear transformation of a certain vector of underlying parameters, s: where R is a linear projection from s-space to d-space. The vector n denotes a set of random variables that introduces noise to the observables. The statistical uncertainties of the observables are thought to be vanishing on average, n = 0, and may be a function of s. All involved matrices and vectors may be complex numbers, the dagger † denotes the complex conjugate transpose. The inversion of Eq. 1 with respect to s for a given d can, in general, only be done in a probabilistic way because of a) the noise in the observation and b) because R is usually not regular or not even a square matrix. One approach, the maximum-likelihood (ML) approach, uses the posterior likelihood function for s given the observation d and known probability distributions of the involved random variables to find the most likely solution for s (e.g. Kitaura & Enßlin 2008): The Bayesian evidence P (d) is a constant with no interest for the scope of this paper, P (s) is a prior on s and P (d|s) the likelihood function of d given some parameters s. In a weaker approach to the inversion problem, one can demand to find a linear operator, Hd say, acting on the observation d that minimises the variance of residuals between the guessed signal and true signal, (s − Hd)(s − Hd) † , without the strong assumption of P (s) and P (s|d) (Seljak 1998). All that is now required, are the statistical means of s and n, both assumed to be zero in our case (or can be made zero by definition), and their covariances N ≡ nn † and S ≡ ss † . Minimising the covariance of the residuals for all conceivable noise and signal configurations results in a minimum variance (MV) estimator: where 1 is a unity matrix of the dimensions of S. We call this filter H a Wiener filter. Here, we have explicitly assumed no correlation between signal and noise, i.e. sn † = 0 (see Appendix B). If we would not like or cannot make strong assumptions on the signal covariance, we can assume it to have infinite variance, S −1 = 0, avoiding any preference in parameter space and only using the noise covariance as weight (inverse variance filter): In the Bayesian interpretation of the MV-estimator we would not have any prior on s. For multivariate Gaussian distributions P (s) and P (d|s) the maximum-likelihood and minimum-variance solution are identical, so that in the Gaussian case the most likely solution is the minimum-variance solution. For this paper, the reconstructed fields will be cosmological density fields smoothed to large scales which makes them roughly Gaussian fields.
Commonly, the sizes of the vectors and matrices can become large so that multiplications between vectors and matrices, a O(N 2 ) process, or matrix inversion, a O(N 3 ) process (Press et al. 1992), can take a very long time or require much computer memory if done without optimisations. The last step in Eq. 4 ensures that we have one matrix inversion less.
It may be that in the estimators Eq. 4 or, in particular, Eq. 5 the matrices in the square brackets are not regular and hence cannot be inverted. This happens in reconstructions of the lensing convergence due to the mass-sheet degeneracy. This problem can be solved by enforcing another regularisation that asserts the mean over all elements of sMV to be exactly a particular constant s0. This makes the minimisation problem of finding H a minimisation with additional constraint s † s = s0 breaking the degeneracy (Hu & Keeton 2002).
Alternatively, we can assert s † s to be close to zero (or another constant value) confined by a Gaussian probability distribution as prior ∝ e − 1 2 λs † s ; λ is a constant constraining the allowed regime of s † s. This is achieved by adding another equation to the minimisation of the aforementioned residual matrix for H and results in the so-called Tikhonov regularisation (Wright 1996) also breaking the degeneracy: Note that using the Wiener filter, and therefore a constraint on the allowed variance of s † s via S, usually already does the job and breaks, in particular, the mass-sheet degeneracy. A fusion of a Wiener filter and a filter without prior, i.e. S −1 = 0, is introduced by a tuning parameter ("α-tuning"), α ∈ [0, ∞[, for the Wiener filter that changes the impact of the Wiener filter by scaling up the expected signal-to-noise in the data by S −1 → αS −1 (Saskatoon filter : Tegmark 1997;Tegmark et al. 1997): For α = 1, the Saskatoon filter is identical to the normal Wiener filter, and for α = 0 we would apply the filter without any prior. The Saskatoon filter will be of central importance for the reconstruction technique outlined in this paper.

Covariance
The (noise) covariance of the estimator Eq. 7 is found by working out with . . . being the ensemble average over all realisations of the noise n but leaving the embedded signal unchanged. This results in the expression Here, X ≡ˆR † N −1 R˜− 1 is the covariance of sMV for the no-prior case and W ≡ S[S + αX] −1 the Wiener filter in sspace. The covariance matrix has two limiting regimes. In the regime where the data is noise dominated, |SX −1 | ≪ 1, it is approximately SX −1 S. On the other hand, for data where the noise is small compared to the signal, |SX −1 | ≫ 1, one finds X, which is exactly the covariance of the MV-estimator without any Wiener regularisation. The power of the signal that gets through the filter (on average), obtained by taking the ensemble average . . . ′ over all possible signal realisations while setting the noise n = 0, is 1 In the signal dominated regime, this is roughly S, but in the noise dominated regime it is well described by SX −1 SX −1 S. This shows that a Wiener filter suppresses the signal if the signal-to-noise becomes low.

Bias
We briefly address here the question if the MV-estimators are unbiased estimators, i.e. whether sMV ′ = s. Using the above equations and d ′ = Rs we find for the Wiener filter. For sufficiently small noise, N −1 → ∞ the Wiener filer asymptotically approaches XX −1 , which is the identity matrix if R † N −1 R = X −1 is invertible and undefined otherwise. We may call problems in which X −1 is invertible well-posed. Note that in the forgoing discussion of the asymptotic estimator covariances we assumed a wellposed problem. Therefore, for well-posed problems the Wiener filter is asymptotically unbiased for an increasingly higher signalto-noise, and the non-Wiener filter is always unbiased (no prior). As we will see later, the Wiener filter applies a smoothing to the reconstruction by rescaling of low S/Nmodes to some degree which explains why it is not an unbiased estimator in general.
Non-well posed problems loose information when applying R to the original signal. Therefore, even with no noise one cannot recover the exact signal. Those problems can be made treatable, by introducing the Tikhonov regularisation, if we make X −1 + λ1 invertible. This, however, will also to some extent bias the estimator. One finds W = (X −1 + λ1) −1 X −1 is the bias matrix in this case.
For a Wiener filter an original signal in one of the components of s is, on average, in the reconstruction rescaled and spread over the other vector components sMV. The ith column of the (bias) matrix W tells us how the original signal si is spread over the other vector elements of sMV in the reconstruction.
The apparently obvious thing to remove the bias in the estimate would be to multiply the MV-estimator by the inverse, if it exists, of the bias matrix: This, however, yields exactly the estimator Eq. 5, removing the benefits of the Wiener filter. The noise covariance is now X again as can be seen from the second term in the r.h.s. of the last equation.
The Saskatoon filter, Eq. 7, is a compromise between the advantages of the Wiener filter, Eq. 4 (improved signalto-noise but biased), and the no-prior filter, Eq. 5 (noisy but unbiased). However, less Wiener filtering means less smoothing in the, then noisier, reconstruction which here will require an additional (non-optimal) smoothing after filtering. Zaroubi (2002) derive another linear filter, a specific quadratic estimator, for unbiased reconstruction of the large-scale matter distribution from galaxy redshift surveys. Principally, this filter cold be employed in the context of this paper as well. However, as it turns out, the UMV-estimator of Zaroubi (2002) is identical to Eq. 7 with setting N = 1 (equal weighing of all data points) and α = 0, which is the no-prior filter, Eq. 5 with N = 1 and is therefore automatically included in the following discussion. Moreover, as already pointed out, the no-prior filter is singular for lensing applications so that some regularisation, such as α > 0 or λ > 0, would be required for the UMV-filter, then being a biased as well.

Lensing convergence tomography
Let us see how this fits into the context of three-dimensional mass reconstructions based on cosmic shear tomography. In weak gravitational lensing (e.g. , the observables are complex ellipticities, E (i) k with k being a galaxy index, of galaxy sources that, for our purposes, belong to one out of i = 1 . . . Nz different redshift bins. The source ellipticity is, in the weak lensing regime, an unbiased estimator of the cosmic shear, γ, at the position, θ k , of the source. An angular position on the sky is denoted by θ which, for convenience, is a two-dimensional Cartesian vector expressed as complex number.
For the following, we assume that the vector contains the complex ellipticities binned on a grid, Ng is the number of grid points. Only sources belonging to the same redshift bin are binned together onto an angular grid. All source ellipticities, E (i) j , contained inside the kth grid pixel and the same redshift bin, i, are combined to one average ellipticity where E (i) j and w (i) j are the ellipticity and statistical weight of the sources associated with the kth grid pixel, respectively. If there are no sources inside the pixel or the sum of all w (i) j is vanishing, set ǫ (i) k = 0. We will further assume that sources of all redshift bins are binned onto the same angular grid so that same elements ǫ (i) k and ǫ (j) k of different bins i and j lie along the same line-of-sight, θ k . By A we denote the solid angle of the grid pixels. Naturally, the number of grid points is identical for all redshift bins.
The affiliation of a source to the ith redshift bin means that we have prior knowledge on its redshift. We denote this prior probability on the redshift, z, by p (i) z (z). Although we are using the term "redshift bin" for p (i) z (z) in this context, we would like to point out that the "bins" may be arbitrary, including overlapping, distributions of sources in redshift. It may be conceivable, for example, that only a small fraction of sources have redshift information and can thereby be ordered into non-overlapping z-bins. The cosmic shear information of all the other sources could then be added as well by putting them into one "bin" with a wide redshift distribution, which is overlapping with all the other bins. This way, we can combine all available shear catalogues to acquire the best possible reconstruction. 2 In the weak lensing regime, the lensing convergence κ (i) -also arranged as Ng-dimensional vector of convergences at position θ k analogues to γ (i) -is related to the source ellipticities by a linear transformationPγκ, a complex Ng × Ng matrix (Hu & Keeton 2002): where n (i) γ is a vector of random intrinsic source ellipticities. The linear transformationPγκ is for k = l (Kaiser & Squires 1993) [ where θ * kl is the complex conjugate of the separation vector θ l − θ k . For k = l, we set [Pγκ] kl = 0.
Note that, so far, we have carefully separated the different redshift bins from each other by using the superscript " (i) ". For each bin i there is one Eq. 16. We can combine everything to one compact system of equations, by putting together all information into the vectors: Pγκ ≡ diag{Pγκ,Pγκ, . . . ,Pγκ} .
The round brackets, grouping together the vector arguments, should be understood as big vectors obtained by piling up all embraced vectors on top of each other. Employing the linear relation between shear and convergence and plugging them into the Saskatoon filter, Eq. 7, with R = Pγκ, S = κκ t ′ and N = nγ n † γ gives an estimator, κMV, for the lensing convergence seen by the sources within the various redshift bins.
This estimator is ideally, only for exactly vanishing Bmodes in the cosmic shear signal, a purely real vector. The imaginary part of κMV is at best pure noise, but contains in general a significant signal owing to a possible B-mode signal in the cosmic shear field. Thus, the imaginary part of κMV can be used as parallel B-mode reconstruction, whereas the real part is the sought E-mode of the convergence tomography. E-mode and B-mode are interchangeable by rotating the source ellipticities by 45 • .
For the statistical uncertainty, needed for N, of the source ellipticity attached to the kth grid pixel and source bin i, we take Here the sum for the overall ellipticity variance, [σ (i) ǫ ] 2 , runs over all source ellipticities belonging to the ith source bin, and the sum for the effective number of sources inside the kth pixel, N (i) k , only over the statistical weights of sources belonging to that pixel. For pixels containing no sources, we set σ (i) ǫ = σ∞ with σ∞ = 10 3 . The noise covariance is a diagonal matrix, δ K kl and δ K ij denote Kronecker symbols, since cross-correlations of the intrinsic ellipticities of sources are neglected.
Although it is assumed for the scope of this paper that intrinsic ellipticities of sources are uncorrelated to any other relevant quantity, it may be possible to accommodate also cases where this is not true. For example, if we expect intrinsic alignments between sources of the same or different zbins, i.e. correlations of intrinsic ellipticities (Heavens et al. 2000;Croft & Metzler 2000;Crittenden et al. 2001). This case can be treated by adding the two-point correlation of intrinsic alignments, ξ (ij) ia (θ) as additional noise contaminant to N: Another possible known contaminant of lensing applications is the possible correlation between the intrinsic alignment of sources and the lensing signal (Hirata & Seljak 2004). This would give rise to a correlation between noise and signal, sn † = 0, which was excluded here from the very beginning of the derivation of the linear filters and is being assumed to be negligible for this work. We defer the treatment of this contaminant to a future paper but already outline in Appendix B how this can in principle be included in the estimator.
The signal covariance, S, is the two-point correlation function of the lensing convergence at the separation of two grid points, which is where ξ (ij) + is a two-point correlation function of the cosmic shear (e.g.  between the ith and jth source bin, γ (i) (θ k ) and γ (j) (θ l ).

Density contrast on lens planes
We can go one step further than the lensing convergence by recognising that the convergence κ (i) itself, as seen in the sources of the ith redshift bin, is within the framework 0.3 0.5 0.7 0.9 1.1 1.3 1.5 Figure 1. Lensing efficiency, Q ij , for j = 1 . . . 15 lensplanes, equally spaced between redshifts z = 0, 1.5(∆z = 0.1), and i = 1 . . . 30 source redshift bins, equally spaced between z = 0, 1.5(∆z = 0.05). The numbers labelling some of the lines denote the mean source redshift; the source redshift increases from left to right. The galaxy distribution within the redshift bins is assumed to be flat.
of general relativity a linear transformation of the matter density contrast, δ, along the line-of-sight: .
The function fK(w) is the angular diameter distance as function of the comoving radial distance w and the curvature, K, of the fiducial cosmological model: Furthermore, we denote the vacuum speed of light and the cosmological scale factor by c and a(w), respectively. For the density contrast, we have chosen the comoving coordinate frame such that the first argument x ⊥ in δ(x ⊥ , x ) denotes a position in a plane perpendicular to the line-of-sight, and the second argument, x , the distance of the plane along the line-of-sight. Note that we are assuming a flat sky.
For a practical approximation of Eq. 28, we postulate δ to be constant over l = 1 . . . N lp comoving distance ranges x ∈ [w l , w l+1 ]: This means, we assume that the 3-D density contrast stays constant along the line-of-sight as long as we stay inside the same matter slice; δ (i) k is here the average 3-D density contrast, and not the total, over the width of the lens-slice integrated 3-D contrast as can also be found in the literature. The coefficients Q il (lensing efficiency) define the response of the convergence tomography to the matter density on the individual lens-planes. A mean density contrastδ, over the radial width of the ith slice, induces on average shear polarisation of sources belonging to the jth redshift bin of amplitude Qjiδ. The coefficients Qij are independent of the line-of-sight direction, Fig. 1 displays an example for Q il .
Due to the approximation in Eq. 31, the mean matter density contrast within a slice is: hence a mean weighed with p (i) (w). For relatively narrow lens-planes, the weights p (i) (w) will be roughly constant over the slice width so that they can overall be approximated by top-hat functions with width ∆wi = wi+1 − wi. This is what we will assume for the remainder of the paper. As before, we can cast the relation between mean matter density and convergence tomography into a compact notation by writing: and where 1 is the Ng × Ng unity matrix. The relation between the observable source ellipticities and the lens-plane density contrast is therefore: As before for the lensing convergence tomography, an estimator for the average density contrast, δMV, on the lensplanes is given by Eq. 7 with R = PγκQ, N = nγ n † γ and S = S δ ≡ δδ t ′ . For an efficient, practical implementation of the filter see Appendix A1.
Analogous to the estimator κMV, δMV will be a complex vector although the density contrast is by definition a real number. The real part of δMV will be the mass reconstruction due to the cosmic shear E-mode, whereas the imaginary part has to be attributed to the cosmic shear B-mode.

Transverse Wiener filter
How can a prior on the signal covariance S δ be computed for a given fiducial cosmological model? By exploiting the Limber equation in Fourier space , the relation between 3-D matter density and projected density can be used to work out the second-order correlations of δ Here P 3d (k, w) is the fluctuation power spectrum of the 3-D matter distribution for a spatial scale k at a comoving distance w. For the scope of this paper, the approximation of Smith et al. (2003) will be used. The function P (i) δ (ℓ) quantifies the fluctuation power spectrum of δ (i) (θ), the projected density contrast on the lens-planes. Note that crosscorrelations between δ (i) 's belonging to different lens-planes vanish in the Limber approximation. Fig. 2 gives an example for ω (i) (θ). The window kernel F (ℓ) accounts for the binning on the grid, it is the Fourier transform of a grid pixel, which we assume as an approximation to be a circular top-hat function, A = πΘ 2 s , with radius Θs, hence: Jn(x) is the Bessel function of first kind. We call this prior on the correlations of density contrasts in different directions a transverse Wiener filter.

Radial Wiener filter
Hu & Keeton (2002) discuss another Wiener filter, also applied in Bacon & Taylor (2003) for the 3-D potential, that puts priors only on the correlation between structures (average density contrast inside pixel needle light cones) in the radial direction, i.e.
fi δ wherê The solid angle Θs denotes the radius of a circular pixel in radians andwi = wi + ∆wi/2 is the comoving (radial) distance of the centre of the ith lens-plane.
In the total absence of correlations between pixels of different line-of-sights, the angular power spectra are flat power spectra (white spectra), which for a circular top-hat smoothing window equals 3 We use this prior for the radial Wiener filter. Obviously, the amplitude of the flat signal power depends on the radius, Θs, of the pixels since Eq. 43 is also smoothing density fluctuations in transverse direction. The following discussion will adopt Θs = 1 ′ . This may appear too specific, but a change of the radius will only rescale the overall amplitude of Pij which, in turn, will solely act as a rescaling of the tuning parameter α. The essential feature of the radial filter, however, namely its scale-independence, will remain unaffected by the choice of Θs. We find that Pij is close to diagonal with negative offdiagonal elements roughly two orders of magnitude or more smaller than the diagonal elements. Thus, the main difference between our radial and transverse Wiener filter is the scale-independence of the radial filter, the off-diagonal elements are not of much importance which we checked by setting them to zero.
The reader may be reminded that the two different Wiener priors discussed here -the radial and the transverse filter -differ only in the way assumptions about correlations between pixels on the lens planes are made. The terms do not mean that radial and transverse modes are reconstructed separately. All modes are determined simultaneously during the process of reconstruction, no decomposition into radial and transverse modes is done.

Fourier space reconstructions
For a discussion of the noise properties of the 3-D mass reconstruction in the following section a Fourier representation of the problem is extremely helpful. Moreover, it allows with some idealistic assumptions to perform a quick-and-dirty reconstruction in Fourier space which is of practical interest, see Appendix A2.
An idealised survey is considered. Such a survey has either periodic boundary conditions or is infinite in extent (flat sky), and the statistical noise n (i) γ is homogeneous for all source redshift bins. The latter means that the statistical uncertainty in the shear estimate is equal for each grid cell of the same redshift bin.
In this case the matter density contrast, and all other fields, can equivalently be expressed in terms of its angular modes 3 As a related example, take the intrinsic shear ellipticity noise which has the flat power σ 2 ǫn −1 . The intrinsic noise inside a circular pixel is then σ 2 ǫ /nπΘ 2 s = σ 2 ǫN −1 , whereN is the mean number of sources inside a pixel. Hence, we would expect the intrinsic noise variance to go down as ∝ 1/ √N (Poisson shot noise).
which defines the density contrast on a typical scale 2π/|ℓ|. The advantage is that now only angular modes of the same ℓ couple to produce the angular modes of the shear tomography because the conversion from lensing convergence to shear is a convolution, where D(ℓ) ≡ ℓ/ℓ * (Kaiser & Squires 1993) and ℓ = ℓ1 + iℓ2 with ℓ1 and ℓ2 being the Cartesian components of the angular mode.
Devising the filter Eq. 7 of Sect. 2 one gets for the MVdensity modẽ Note thatδ andγ (the Fourier transform of the gridded shear catalogues) are written here in a compact vectorial form with mere N lp /Nz elements, grouping all lensplanes/source bins together, as outlined in the foregoing sections; α is the usual tuning parameter. Here we furthermore use which -for the sake of a simplistic noise model used in the next section -assumes a shot-noise term only 4 ;ni is the mean number density of sources belonging to the ith redshift bin. The transverse Wiener filter constraints the expected signal power by Noise and signal are subject to the same smoothing if we bin our data on a grid. Therefore, the effect of a smoothing window cancels out in the filter, sinceS δ andÑ appear always as ratios in the filter. The (expected) signalS δ is a diagonal matrix for transverse filtering in the Fourier space representation. This may change, if we choose to apply radial Wiener filtering, where only correlations between matter densities along the same line-of-sight are regularised. Since correlations between different directions are neglected, the power spectra used for regularisation are flat power spectra (Sect. 3.2.2), independent of ℓ:

Gravitational potential on lens planes
We may ask the question how the MV-estimator of δ is related to the MV-estimator of that is thought to be a linear, invertible transformation of δ. Such cases could be (lens-plane-) smoothed density reconstructions, for which F (i) would be a (invertible) smoothing operator acting on δ (i) , or the 2-D gravitational potential on the lens-planes which can be pictured as a special type of smoothing of δ (i) (Bacon & Taylor 2003): or in a form for a discrete grid assuming that the grid sizes for φ and δ are equal: The density contrast is constant over the size of one (round) pixel. For θ l = θ k one finds therefore [F (i) ] kk = −A/2π. Note that Eq. 56 gives only one possible solution for the 2-D gravitational potential because the potential for a given Following the arguments in Sect. 2 it is easy to show that the MV-estimator for φ is This means, once we have got the MV-solution of δ we can simply apply F to it in order to acquire the MV-solution of φ. We do not have to go through the process of another full MV-reconstruction, if we want to further smooth a MV-solution. That is, provided the smoothing is invertible.

BIAS AND NOISE PROPERTIES OF RECONSTRUCTIONS
As shown in Hu & Keeton (2002) a mass reconstruction, even with as many sources as 100 arcmin −2 and with heavy smoothing, is very noisy. Therefore, an additional regularisation of the reconstruction, by means of Wiener filtering for instance, is an absolute necessity. With regularisation, however, a careful analysis of the remaining statistical uncertainties and biases in the reconstruction will be required. This will be the focus of this section.

Average signal-to-noise of reconstructions
Devising the MV-filter in Fourier space, Eq. 51, offers a convenient way to estimate the covariance of the noise in the reconstruction. The following conclusions, however, will be based on the simplifying assumption that the noise pattern is homogeneous and that there are no gaps in the data, in particular the shear pattern of the full (infinite) flat sky is known.
We start by considering the signal-to-noise for individual ℓ-modes and assuming that the matter fluctuations of the large-scale structure that provide the signal in the shear catalogue have cosmic average amplitudes.
The covariance of the statistical uncertainties in Eq. 51 is according to Eq. 9: The latter is the noise covariance without a Wiener prior or if the prior is unimportant. How much signal-to-noise can we expect on average for a 3-D mass reconstruction as function of angular scale? To answer this question, we compare, for a ΛCDM fiducial cosmological model, the expected signal power as function of scale ℓ,S, getting through the filter (Eq. 11), to the noise power: As expected density fluctuation power we take the prediction based upon the Limber equation.
Since the noise covariance,Ñ, scales as ∝ σ 2 ǫ /n, we can derive useful scaling relations for the signal and noise power with a fixed redshift distribution of sources. One finds that the tuning parameter inside the filter scales as if we change the intrinsic shape variance from σǫ to σ ′ ǫ or the source number density fromn ton ′ . Thus, increasing the number density ton ′ , for instance, yields a filter corresponding ton but with a different, smaller α. Similarly, we get: so that the signal-to-noise of modes scales as This property is useful for preparing plots of the estimated signal-to-noise spanning a wide range of fiducial surveys. Fig. 3 is a particular example of the signal-to-noise of reconstructed density modes using realistic survey parameters. For the source distribution with redshift a function with z0 = 0.57 was used; the mean redshift isz = 0.85. The distribution was truncated beyond a redshift of z = 1.5. The redshift distribution is somewhat shallow but not too unrealistic for contemporary surveys if we take into account that source redshifts are required. Obviously, even for a very optimistic source density of 100 arcmin −2 , Wiener filtering is utterly necessary to improve the signal-to-noise in the reconstructions. Simple transverse smoothing of the reconstruction obtained without Wiener prior, i.e. suppressing high-ℓ modes, does not suffice as even the lowest ℓ-modes have a signal-to-noise less than unity.
The transverse or radial Wiener filter bring the signal up to a signal-to-noise of roughly unity, somewhat higher on very large angular scales. If we tune down the Wiener filter by lowering the α-parameter, we get a signal-to-noise that lies somewhere in between of the full Wiener filter and the very noisy reconstruction with no prior (not shown). These figures also demonstrate the attenuation of low-S/N modes by the Wiener filter, i.e. the ratio . This tells us that the Wiener filter will, on average, not recover the original amplitude of a fluctuation mode. Instead, the amplitude will, depending on the redshift of the lens plane and the scale ℓ, be attenuated. This is an unavoidable sideeffect of Wiener filtering, which, on the other hand, thereby improves the signal-to-noise in the reconstruction. Another side-effect concerning the radial distribution of densities will be discussed in the next section.
For the transverse filter especially, small scale modes are damped, which thereby acts as a low-pass filter performing an automatic smoothing on the lens-planes. In case of a radial Wiener filter, the damping is roughly the same on all scales due to the scale-independence of the regularisation.
Notice that in Fig. 3 the tuning parameters for the transverse, α, and radial filter, β, are not completely equivalent, although they have the same effect: For our transverse filter, the theoretical signal power is identical to the expected signal power plugged into the filter asS δ =S (both based on Limber's equation), whereas for the radial Wiener filter, the expected powerS -still the same as before to have a comparison assuming the same input signal -is different from S δ used inside the Wiener filter (flat power spectrum).
The improvement of the signal-to-noise due to Wiener filtering comes with another price that has to be paid. The statistical errors between the modes of different lens-planes, which are weakly correlated if no Wiener filtering is applied, become correlated. Fig. 4 shows the matrix for the previous survey parameters with 30 arcmin −2 source density. With no prior at all one finds that the errors between different lens-planes quickly decorrelate in a oscillatory manner (bottom right). As pointed out in Hu & Keeton (2002) this particular form of error correlations is owed to the finite difference approximation of a fourth-derivative in radial direction. By means of the Wiener filter, however, the errors become correlated over a wider range (the filter also smoothes in radial direction). Errors between different ℓ's remain uncorrelated as they were. For a signal-to-noise close to one, the error correlations limit our ability to pin down the redshift of structures along the line-of-sight, that is the redshift resolution of the reconstruction. The correlations have also a positive side: the radial oscillations become stretched out by the Wiener filter.
The foregoing way of estimating the signal-to-noise in a reconstruction takes a quite pessimistic view. Namely, it presumes that the matter fluctuations of the structure we    would like to recover are actually of an amplitude we expect for a cosmic average, i.e. we expect |δ(ℓ)| 2 ∼ P δ (ℓ). In fact, in individual regions of space we will find matter concentration exceeding the average (large galaxy clusters), or we will find less matter (voids).
Let's quickly make an rough estimate of how much power above the average we may can expect. For a Gaussian random field -probably being a fair approximation of the mass density field on large scales and, thus, the regime we can realistically hope to reconstruct -the expected variance of the fluctuation power is (exploiting Wick's theorem): Therefore, we may expect, owing to the natural variance in a Gaussian random field, to find regions where the signal power is actually up to four (2σ) times the average fluctuation power, enhancing the signal-to-noise in Fig. 3 by a factor of two -or even more in the non-linear regime. In Sect. 4.4 we will work out the signal-to-noise of specific matter haloes.

Response in Fourier space
There will be a fundamental limit in recovering the radial matter distribution, even if we have peaks with a high signalto-noise. As discussed in Sect. 2, a Wiener based reconstruction is unavoidably biased, unless the signal-to-noise of the data is infinite. This bias is quantified by the bias matrix f W that relates the original vector of ℓ-mode coefficients to the (statistical average) vector of modes in the reconstruction. It is a convolution applied to the true matter distribution.
Specifically, the ith column of f W expresses how the Wiener reconstruction, on average, responds in the reconstructed density contrast field to a valueδ (i) (ℓ) in the true matter distribution. Ideally, we would like to have a unity bias matrix (unbiased estimator), but we find that, in the reconstruction, the Fourier coefficient becomes rescaled (attenuation), shifted and spread out in the radial direction, contributing to otherδ  Figure 5. Examples for the response of the Wiener filtered (transverse) matter reconstruction algorithm to peaks in angular Fourier space of coherence scale ℓ = 800 in the true matter distribution at low, z ∼ 0.05, and high, z ∼ 1.05, redshift (downward arrows). The fiducial survey has the parameters as outlined in the caption of Fig. 1 assuming a source density of 30 arcmin −2 with α = 1 (thick lines). In the 3-D reconstruction, the localised peaks are shifted in redshift and smeared out, independent of the underlying signal, which limits our ability to reconstruct the matter distribution in radial direction. The thin lines correspond to a Wiener filter withn = 100 arcmin −2 and α = 1. Increasingn or decreasing α reduces the shift and spread of the response.
certain ℓ-mode is demonstrated by Fig. 5. In this particular but realistic example, an angular mode of a lens-plane at redshift z = 1.05 (z = 0.05) is shifted to z ∼ 0.5 (z = 0.15) and smeared out by σz ≈ 0.3 (1σ-width; σz ≈ 0.1) to adjacent lens planes. The smearing of the peak is the reason for the strong correlation of modes of neighbouring lens-planes, which has been discussed in the forgoing section.
Especially structures at the upper redshift end of the survey a prone to the shift and smearing. Both effects are solely a function of angular scale, ℓ, larger coherent structures (smaller ℓ's) are easier to resolve than smaller structures. For small enough ℓ's all redshifts are eventually shifted to one redshift and similar response functions.

Point spread function
The ℓ-dependence of the response in the reconstruction as described above is hard to relate to the question how well we can recover the redshift matter clump as function of mass and redshift. Therefore, the following translates the Fourier picture of the reconstructions back to real (angular) space, where we further study the Wiener filter impact.
Assume we would like to calculate the effect of the Wiener filter on a halo sitting at some redshift, or equivalently on some lens plane i, with projected density contrast profile Furthermore, we consider for simplicity only profiles that are rotationally invariant, i.e.
A single "point" would be just a Dirac delta function in direction of θ.
The prefactor ∆w denotes the range of comoving distances spanned by the matter slice the halo is located in. In our framework, δ 2d (θ) is the average 3-D density contrast, constant inside the slice represented by the lens-plane, and not the integrated 3-D density contrast over the width of a slice, see Sect. 3.2.
If we place the halo at the centre of our coordinate system, θ = 0, and smooth the overdensity over the solid angle, A = πΘ 2 s , of the circular pixel with radius Θs (our grid binning) centred on θ = 0, we would find on average after Wiener filtering, performing the integration over all modes ℓ, a smoothed over-density of whereδ(ℓ) is the 2-D Fourier transform of the radially symmetric halo profile: The window F (x) defines the smoothing window (also radially symmetric) of our lens-plane pixels, which for circular pixels is defined in Eq. 42. For the scope of this analysis, the original halo signal, δ(ℓ), is a vector of mode amplitudes on the lens-planes, being zero except for the plane on which the halo is located. After filtering, l.h.s. of Eq. 76, the signal will be spread along the line-of-sight, i.e. spread over the output vectorδ(Θs, α), as already seen for the Fourier modes, Fig. 5. Moreover, as everything is linear, the response of a sum of haloes is just the sum of the individual responses.
There is also a transverse (spread inside same lensplane) point spread function (PSF) associated with the filter, which can be evaluated at separation ∆θ from the halo centre by offsetting the position of δ(θ) in Fourier space, and reassessing Eq. 76 (for a radially symmetric profile): Obviously, the transverse PSF is symmetric (only a function of the modulus of ∆θ) and will, therefore, not bias the angular position of a halo. Moreover, usually a strongly downtuned Wiener filter will be used, α ≪ 1, which will make the transverse spread relatively small. As a matter of fact, we will have to apply an additional transverse smoothing to get a reasonable signal-to-noise data cube. Conversely, the radial PSF will be of more concern. Fig. 6 gives a toy model example of the response of the Wiener filter to some halo with profile δ(θ, z h ), sitting at a singular redshift z h = 0.5. If the response is a directionindependent function C(z h , z) -in a toy model fashion a Gaussian with mean z h and width σz = 0.2 -then we will find in the reconstruction the 3-D map δ(θ, z h )C(z h , z). The Figure 6. Toy model example of the response to a halo with some density profile sitting at a redshift of z = 0.5, originally with no radial extension. Plotted are four different iso-density surfaces with decreasing iso-densities, growing bigger in size for smaller iso-density values. Clearly, the surfaces become stretched -the more, the smaller the iso-density -in radial direction in a cigarlike manner (we loose redshift information) due to the response of the Wiener filter. As shown in Fig. 8, the effect becomes smaller for smaller α (less Wiener prior) which, however, will increase the noise in the reconstruction. The iso-surfaces responding to structures situated at the redshift boundaries of the survey will be radially stretched in a asymmetric manner (not shown).
general effect is that the iso-density surfaces of the profile are stretched radially. Although all δ(θ, z h ) are being stretched with the same C(z h , z), surfaces with smaller iso-densities are stretched more than higher iso-density values. Qualitatively, this toy model describes quite well what we find in the real reconstructions. Note that the surfaces are centred about the true halo redshift since no z-shift due to the response C(z h , z) was assumed for this illustration.
Under the idealised circumstances the radial point spread function will be the same in all directions. Since, in reality, we will have a varying source number density as function of direction and a finite survey area we can expect a more complex, direction dependent response of the Wiener filter which accounts for inhomogeneous source ellipticity noise patterns.
Assuming as a lowest-order approximation a homogeneous noise pattern of the sources, we can work out the expected radial PSF of the Wiener filter, i.e. the radial convolution C, by means of Eq. 77: Here, an original pixel on the lens plane has the same shape and size as a reconstructed pixel, i.e.δ(ℓ) ∝ F (ℓΘs). A pixel on the ith lens plane is spread over the other lens planes according to Cei, ei is a unity vector non-zero only at the original position of the pixel, or put in another way: Cij is the signal-leakage of pixels on the ith lens-plane into pixels on the jth lens-plane. In an ideal, unbiased reconstruction, f W = 1, the matrix kernel of the integral would be the unity matrix so that for a properly normalised pixel window function F . Thus, the pixelised reconstruction would yield exactly the original matter distribution on the lens-planes smoothed with the pixel window. In reality, however, we will find a situation more like the one shown in Fig. 7. Now, Fig. 8 uses the results of the previous discussion to quantify the quality of a 3-D reconstruction for a particular fiducial cosmology and a realistic fiducial survey employing a Wiener filtering scheme with α-tuning. The plots can be scaled to different σǫ orn as explained in the caption. We can see in the plots that, the shift and spread of the signal in the reconstruction due to Wiener filtering decreases if the filter is tuned down, which, on the other hand, will decrease the signal-to-noise. We also find that the redshifts of objects near the "edges" of the survey (here: z ∼ 0 and z ∼ 1.0) are mostly biased, high redshift objects appear to be more affected. The bias tends to shift low-z objects to apparently higher redshifts, while high-z objects are moved in direction of small redshifts.
In analogy to these plots, Fig. 9 shows the same parameters but this time obtained with a radial Wiener filter. We find a similar behaviour compared to the transverse filter, keeping in mind that the tuning parameters of the radial and transverse filter are not comparable: an β of the radial filter ten times smaller than the α of the transverse filter yields roughly the same signal-to-noise (see next section). The expected signal based on the radial filter gives, in contrast to the transverse our filter, a different prediction as the "true" signal which is calculated from Limber's equation. Looking more closely, one finds that a radial filter might be used to slightly decrease the z-shift due to the PSF (β ∼ 10 −3 compared to α ∼ 10 −2 ), although the spread in radial direction may be larger. As this spread is partly into the negative regime, see right Fig. 7, it may not become immediately obvious but can have the effect of diminishing structures along the line-of-sight -or generating structures as a response to an underdensity.
We noticed during the course of this work that a very simple radial prior that gives all lens-planes the same prior, Pij = πΘ 2 s , gives often somewhat less biased signal-to-noise maps so that it is worthwhile to compare in a concrete case the benchmarks of this simple radial filter to the previous radial filter.
As already pointed out in the introduction, we could in principle deconvolve out reconstruction with the known PSF. This, unfortunately, would completely remove the effect of the Wiener filter leaving us with the extremely noisy and radially oscillating reconstruction of the no-prior estimator.

Signal-to-noise of matter haloes
Here we would like to compare the expected pixelised signal in a reconstruction to the expected noise level within the same pixel.
If the density modes,δ(ℓ), of Sect. 4.3 belong to a random field -such as the noise in the reconstruction or the random density fluctuations in a particular fiducial cosmological model -the variance of the density contrast inside the a pixel will turn out to be for a power spectrum, P (ℓ), quantifying the fluctuations iñ δ(ℓ). This follows from Eq. 75. If we apply this integral transformation to the power spectra (matrices) Eqs. 62 and 64, we can transform the Fourier space representation to signal-to-noise estimates inside pixels as function of lens-plane redshift. From that one can infer that the average signal-to-noise inside a circular pixel will be approximately the signal-to-noise of one effective ℓ-mode, Fig. 3. As the integration kernel, x|F (x)| 2 , in Eq. 83 peaks at about ℓΘs ≈ 1.36, the effective mode frequency will be near ℓ eff ≈ 4666 (Θs/1 ′ ) −1 .
The expected covariance of noise inside lens-plane pixels arranged along the line-of-sight will be: Cov(δ(Θs, α)) = 1 2π This can be directly compared to the blurred signal from a single halo (Eq. 77), sitting at some redshift, being filtered by our Wiener filter: Note that we only consider the smoothed signal about the halo centre here, where we expect to see the highest signalto-noise of the halo reconstruction. The Wiener filtering does not bias the angular position of a halo (the phases of the Fourier modes), in particular the centre position. Therefore, we will, on average, find the halo centre in the reconstruction in the original direction as before filtering. For a singular isothermal sphere (SIS; e.g.  we have where σv is the velocity dispersion inside the SIS and a(w h ) the scale factor at the radial comoving distance of the halo centre, w h . The Fourier transform of the SIS profile is As the signal inside a pixel, Eq. 77, is linear in the original underlying amplitude of δ, whereas the noise inside the pixel is independent of δ (Eq. 84) we can already infer that the signal-to-noise of the SIS in the reconstruction has to scale with S/N ∝ σ 2 v , roughly the mass of the halo. We can expect this to be approximately true for more realistic halo profiles as well.
More specifically,  give for the virial mass 5 of the SIS where E(z) ≡ p Ωm(1 + z) 3 + (1 − Ωm − ΩΛ)(1 + z) 2 + ΩΛ (90) is The signal-to-noise of the central pixel is highest for small halo redshifts and decreases rapidly towards higher redshifts. This is because the number of sources behind the halo decreases as we move towards higher halo redshifts. Therefore, the technique identifies haloes most efficiently at low redshifts and becomes increasingly ineffective towards higher redshifts. We also found by varying the redshift distribution of the sources, while keeping their number density and intrinsic ellipticity distribution constant, that the signal-to-noise fall-off is shallower for a deeper survey, although the z-spread increases owing to the same number of sources being distributed over a larger z-range.
We also gather from the results, that the signal-to-noise decreases as we tune to smaller values of α, thereby reducing the Wiener regularisation. For a transverse filter and the given fiducial survey, values of α 10 −2 are necessary to reduce the bias of the Wiener filter (1σ-z-spread and zshift less than ∼ 0.2) to an acceptable level, Fig. 8. For our radial filter, a value of β 10 −3 should be used. We excluded values smaller than α 10 −3 (radial: β 10 −4 ) as this brings us too close to the no-prior filter which exhibits heavy, undesired oscillations in radial direction and very poor signal-to-noise. Note that α also depends on the redshift distribution of the sources.
An example reconstruction (transverse filter) of a set of SISs planted into a one-square-degree field-of-view between redshift z = 0...1 can be found in Fig. 11. This example underlines the theoretical expectation for our fiducial survey that the signal-to-noise of the haloes in the reconstruction drops quickly if we move towards higher redshifts; the most distant haloes (F-I) are barely or not distinguishable from the background noise. Note that the signal-to-noise levels cited in the figure caption were computed based on 500 FFT-reconstructed noise realisations, see Appendix A3. One can also observe the effect of the PSF that spreads out Figure 10. Signal-to-noise of a reconstructed SIS, positioned at different redshifts, x-axis, with σv = 10 3 kms −1 ; this corresponds to M 200 = 6.6 × 10 14 M ⊙ h −1 (1.6 × 10 15 M ⊙ h −1 ) at z = 0 (1). Plotted is as function of tuning parameter (line styles) the signal-to-noise at maximum response. The assumed fiducial survey has the parameters outlined in Fig. 8. To obtain the corresponding plots for different " 1/2 , as described in Sect. 4.1. Left figure: transverse Wiener filter. Right figure: radial Wiener filter. Figure 11. A set of nine SISs (σv = 10 3 km s −1 ) with redshifts z = 0.1, 0.2, . . . , 0.9 (A-I in alphabetical order) has been planted into a particular noise realisation of the fiducial survey defined in Fig. 8. The left panel is a 2-D projection of the reconstruction onto the sky while the right panel is a projection showing the extension of the volume in redshift. The differently coloured contours are constant signal-to-noise levels corresponding to S/N = 8, 6, 5, 3, 2.5 (darkest to brightest colour). The reconstruction is smoothed with a kernel of radius 1 ′ . The employed transverse Wiener filter has α = 0.01. Statistical errors can dissolve peaks as can be seen for D and E or shift peaks towards wrong redshifts, see D for instance. Note that all redshifts are biased to some extend, left panel of Fig. 8. Figure 12. Reconstruction of the lensing convergence on a 60 ′ × 60 ′ patch falsely assuming that all shear signal originates from matter inside the patch (left: B-mode, right: E-mode). Actually, the signal is solely produced by a SIS, σv = 10 3 km s −1 and z = 0.2, sitting in a separation of 6 ′ off the right edge of the patch. The noise in the images is due to the finite number density of (intrinsically circular) sources,n = 30 arcmin −2 , used for the fiducial survey here.
the peaked matter distributions and moves, for instance, the maximum of the response of low-redshift haloes towards higher redshifts. Due to the noise in the data, some peaks (C and D) fragment into different parts and turn up at clearly wrong redshifts. An analogous reconstruction (β ∼ 10 −3 ) employing a radial filter looks essentially the same.

Edge effects
The reconstruction algorithm assumes that the cosmic shear signal within the patch of the survey originates solely from matter within the field-of-view. This is the same for all Kaiser-Squires type algorithms (KS; Kaiser & Squires 1993). This introduces, especially at the boundaries of the field, a bias if the actual source of the shear within the field is outside the patch. This problem would not exist for a full-sky implementation of the algorithm.
To have an qualitative assessment of the impact of this bias, we consider one particular, extreme case which is shown in Fig. 12. For this mock survey, a square patch with 60 ′ on each side and the previous redshift distribution of sources, we made up the case that there are no sources of shear at all inside the patch area. The total shear signal observed stems from one massive singular isothermal sphere located close to the right edge of the patch. The figure shows a KS reconstruction (left: B-mode, right: E-mode) of the convergence inside the patch ignoring an outside source as possible origin of the shear. The 3-D reconstruction of the density contrast for a lens plane at about the redshift of the SIS exhibits exactly the same pattern. Lens planes well separated from the SIS redshift show no or little reconstructed mass density while neighboured lens planes may be affected due to the radial spread of the reconstruction PSF.
It can be seen from this test that a source outside the observation area can indeed produce mass density inside the reconstruction patch. However, this "ghost signal" seems to be mainly focused on the edges. The largest inferences can be seen on the edge closest to the SIS, namely at the middle and at the corners, and on the edge opposite to that with roughly the same but inverted signal pattern.
Furthermore, it can be seen that this ghost signal is associated with a B-mode signal of roughly the same amplitude than the E-mode signal, also most prominent at the edges (closest and opposite) but leaking somewhat more into the patch. Therefore, if an off-patch source produces a significant signal inside the patch area it may reveal itself by an equally significant B-mode signal.
Note that the response pattern looks different if we use FFTs for the reconstruction, as they assume periodic boundary conditions. With FFT the SIS can actually be reconstructed as mirror image near the edge opposite to the edge closest to the true position of the SIS. Yet, we find that even with periodic boundary conditions a B-mode signal is generated inside the patch which has the same order of magnitude than the E-mode signal.
This bias can be controlled by making the reconstruction area larger than the patch area, attributing to this additional patch frame infinite noise. By doing this, we explicitly allow the reconstruction algorithm to place sources of shear outside the observed patch area.

Reduced shear
Throughout this paper we assumed the weak lensing approximation to be perfectly valid, i.e. ellipticities of galaxy images are unbiased estimators of the cosmic shear and not the reduced shear, γ/(1 − κ), as they actually are.
If the distinction between reduced shear and shear is ignored, the convergence will in the vicinity of density peaks biased towards smaller values and hence the density will be biased towards smaller values as well.
Presumably, the algorithm can be modified in this respect by running the linear algorithm on ǫ , where is κ = Qδ is the minimum variance convergence reconstruction of the previous run; for the first run one starts off with κ (i) j = 0. For a 2-D convergence reconstruction this iterative approach is known to converge quickly to the right solution (Seljak 1998, Sect. 4.3). The noise of the nonlinear reconstruction would have to be estimated by a series of noise realisations. We will postpone a test of this approach for a 3-D reconstruction to a future paper.

OUTLOOK
For a demonstration and an outlook of what may be achieved with this technique in the near future we apply the reconstruction method to a simulated lensing survey patch of four square degree size. The simulated survey has a source redshift distribution as in Eq. 70 but now with z0 = 1.0, a source number density ofn = 100 arcmin −2 , corresponding to 1.6 × 10 6 sources, and an average intrinsic ellipticity variance of σǫ = 0.3 for all sources. It is further assumed that we have (photometric) redshifts for all individual galaxies and that we are unable to obtain (photometric) redshifts beyond a redshift of z = 2, whence we truncate p(z) at a redshift of two. The mean redshift of the sources is thereforez = 1.1. These survey parameters roughly reflect a space-based lensing survey which we consider as the best possible choice for Figure 13. 3-D mass density reconstruction of a 2 deg × 2 deg patch of a mock lensing survey with 100 arcmin −2 source number density and mean source redshiftz = 1.1 with z 2.0. The patch was chosen to have a large number of high-density peaks. The intrinsic shape variance is chosen to be σǫ = 0.3. The reconstruction was done on a 128 × 128 pixel 2 grid with 25 source z-bins (∆z = 0.08), and 20 lens planes (∆z = 0.1). The reconstruction was done on grid larger than the observed area, allowing for a frame with 7 ′ .5 width. A transverse Wiener filter with α = 0.05 was employed. For the signal-to-noise maps, which are shown here, the data was randomised 100 times. The upper left (projection on the sky; x vs y), lower left (x vs. redshift) and lower right panel (y vs. redshift) show the signal-to-noise data cube of the reconstruction with 3-D iso-levels of S/N = 6, 5, 4.5, 4, 3.5, 3, 2.5. The upper right panel shows for comparison true high over-densities (not over-densities due to projection) on the lens-planes used for the ray-tracing. The most significant peaks, S/N 4, can be identified and matched with true matter peaks (solid circles). Several less significant peaks can be matched as well (not all encircled). But there are also peaks in the reconstruction with no apparent real counterpart (dotted circles) or high density peaks in the true distribution with no counterpart in the reconstruction (dashed circles); not all cases have been highlighted. The numbers denote the estimated (left) or true redshifts (right) of the structures. The most significant peaks in the upper left panel are labelled with alphabetic letters, A-G, and can be identified in the radial projections in the lower row. Label X denotes a case discussed in the text. No attempt has been made to correct for the z-shift bias which is most evident for z 0.15. mass reconstructions due to the small ratio σǫ/ √n , which is the main factor in the signal-to-noise of the maps.
The mock survey was generated by ray-tracing through an N -body simulation run with the publicly available version of GADGET-2 (Springel 2005). We have used 256 3 Dark Matter particles in a box with a side length of 150 h −1 Mpc, which together with the adopted ΛCDM cosmology (Ωm = 0.25, Ω b = 0.04, ΩΛ = 0.75, σ8 = 0.78) results in a particle mass of 1.2 × 10 10 h −1 M⊙. We have traced 2048 2 light rays trough 25 matter slices with a thickness of 150 h −1 Mpc up to z = 2 using the standard multiple-lens-plane algorithm (e.g. Hilbert et al. 2009, and references therein). This yields the Jacobian matrix of the lens mapping to the back side of each matter slice. We have then created the mock survey by randomly sampling the slice boundaries with galaxies according to the desired redshift distribution, onto which we interpolate the shear from the light ray positions.
For the reconstruction, we subdivide, equally in redshift, the source catalogue into 25 sub-samples with width ∆z = 0.08. This number of source bins is fixed here by the number of lens planes originally used for ray-tracing to obtain the simulated data. In reality, ∆z reflects the average uncertainty in the redshift estimates which is probably somewhat too pessimistic with 0.08. However, we would like to stress that a much finer binning does not make much difference in the reconstruction -at least if the redshift estimates are not biased and the redshift distribution inside the bins is accurately known -since the lensing efficiency is but a slowly changing function with source redshift. This is related to the known fact that increasing the number of source bins beyond a few, in constraining cosmological parameters with lensing tomography, does not significantly increase the constraints. Hence, we do not expect a notable improvement by using a larger number of source bins.
The number of lens planes, on which the matter density contrast is to be estimated, is set to twenty ranging between z = 0 . . . 2 with ∆z = 0.1. We do not use more lens planes, as a radial resolution exceeding ∆z = 0.1 should not be expected according to the foregoing discussion. In the following plots, we interpolated between the lens planes. Fig. 13 shows the reconstruction of the simulated patch as signal-to-noise map either projected onto the sky or in two different radial projections. Note that the sky-projections are actually 3-D signal-to-noise maps. A signal-to-noise map of the projected matter distribution, as commonly extracted from lensing surveys (smoothed convergence maps), can be anticipated to have higher significance. For comparison also a sky-projection of the most massive density peaks on the original lens planes is displayed. For this demonstration a patch particularly rich in over-density peaks was chosen, on average the number of peaks in a four square degree patch in a WMAP3-like universe should be expected to be smaller.
We observe that although we have a quite optimistic lensing survey, the significant features in the reconstruction mainly correspond to the most massive haloes in the Dark Matter density field. A detection, S/N 3, is restricted to the regime z 1. Below S/N ≈ 3.5 we find cases in the reconstruction which do not seem to have a real counterpart. Those can be explained as pure statistical flukes of the noise pattern which are still to be expected on this signal-to-noise level. It may be that 100 noise realisations, as performed here, do not suffice to estimate the noise pattern accurately enough everywhere. We also seem to have cases where we have no significant detection of arguably high peaks in the true Dark Matter density. The problem may be here that we are locally lacking the required sampling of background sources. Peaks along the same or close to the same lineof-sight are not always distinguished from each other and merged into a mediate redshift instead (see E), although there are cases where the algorithm succeeds in doing so (see C).
Quite interesting is the case X, upper right quadrant, where we have a large peak at very low redshift z = 0.02 which however is completely missed by the reconstruction or might be hinted at by a very weak detection at a completely wrong redshift. The explanation for this case probably is that the lensing efficiency for a structure at a redshift that low is either too small for sources at higher z or the number density of sources which in principle would be sensitive enough to light-deflections at that low redshifts is too small (low z). In this sense, the lensing approach appears to be blind towards structure at very low z.
Seemingly, the fiducial survey does not perform a lot better than the inferior survey in Fig. 11, where we assumed n = 30 arcmin −2 and a somewhat shallower mean redshift. An increase inn from the latter to the former survey, although observationally rather challenging, only improves the average signal-to-noise by a mere factor of p 10/3 ∼ 1.8. This again highlights the problem of all lensing mapping schemes.
The signal-to-noise of the reconstruction can be strongly increased by using a more conservative, larger, tuning parameter (here: α = 0.05). For instance, structure D can be boosted to a signal-to-noise of approximately ten by using α = 1, however on the expense of loosing information on the redshifts, revealed by a significantly larger spread in radial direction, and more z-shift bias.
Note that the z-shift bias is somewhat different in a signal-to-noise data cube as opposed to a density contrast data cube which was discussed in the foregoing section. Moreover, we have not tried to correct this bias in our redshift estimates. In the plot, therefore, particularly structures at low redshift, z 0.15, are artificially shifted to larger redshifts, z ∼ 0.2.
The redshift estimates of structures more significant than S/N ≈ 4 are all in all reliable within ±0.1, provided they are not unresolved structures at different redshifts, but become increasingly inaccurate below that signal-to-noise level.
The scenario discussed here contemplates a case where we have a small but deep patch of which we would like to obtain a 3-D map with a preferably good spatial resolution. Another scenario would be a (almost) full-sky survey to which heavy smoothing is applied to study the very largescale density modes. Fig. 3 implies that for those modes the reconstruction may become easier than for the modes on a small patch as the signal-to-noise improves for lower ℓ and the z-shift bias becomes less severe.
The situation may be further improved if the shear data is combined with higher-order lensing information, such as gravitational flexion, that may be available in the future (Bacon et al. 2006;Okura et al. 2007). In the weak lensing regime, higher-order moments of lensing distortions can be modelled as linear functions of the lensing convergence, as derivatives of κ. Therefore, higher-order moments are easily implemented into the algorithm as new extra set of redshift bins with a new linear projection operator "PFκ" as opposed to Pγκ, Q remains unchanged. Very crudely -same number of sources as for cosmic shear and same (uncorrelated) signal-to-noise provided -this will be equivalent to effectively increasing the number density of sources in a conventional cosmic shear survey by a factor of N + 1, where N is the number of new lensing distortion moments that is included; the signal-to-noise in the reconstruction consequently goes up by √ N + 1. Actually, as already known from first 2-D applications (Bacon et al. 2006;Leonard et al. 2007), the 3-D maps will become better on small angular scales with less improvement on larger scales as the signal-to-noise of flexion-sampled mass density fields is scale-dependent with most information on small scales. The radial resolution of the maps, however, will only moderately benefit from the new information, i.e. no more than from an extra cosmic shear catalogue with same size, as the dependence of higher-order lensing on redshift is identical to that of cosmic shear.

SUMMARY AND CONCLUSIONS
In this paper we presented a linear algorithm -and two practical ways of implementing it (Appendix A1,A2) -that allows to reconstruct the three dimensional distribution of matter based on a lensing survey with redshift information. All available redshift information, accurate source redshifts binned into thin slices and broad redshift information such as wide redshift distributions for faint, high-z source subsamples, can be combined to find a matter distribution that should be on average closest to the true distribution. The statistical uncertainties of the source redshifts, attached as PDF's of the redshifts to the corresponding source subsample, and the statistical uncertainty of the shape measurement can be properly factored into the reconstruction. Principally, as shown in the paper, it is also possible to account for intrinsic alignments of the sources, that have an impact on the expected shape noise covariance entering the filter, and shear-intrinsic alignments.
The presented algorithm, being a generalisation of the algorithm proposed by Hu & Keeton (2002), yields the most probable matter distribution, if the matter density field obeys Gaussian statistics, or a minimum-variance solution otherwise. Due to the heavy smoothing required in a 3-D matter distribution reconstruction based on lensing tomography, the reconstruction is probably close to the maximumlikelihood solution.
The algorithm is a Wiener filtering of the lensing tomography data. The strength of the Wiener filter prior, i.e. the assumed second-order correlations between densities in the reconstruction, can be tuned with an additional parameter, α. We looked into two different cases for a regularisation, namely one case for which only correlations between densities at the same redshift are given (transverse filtering) and a case where only correlations along the line-of-sight are specified (radial filter). Both priors show similar properties, it is not possible to single out one filter as superior filter.
We showed that a reconstruction with virtually no prior, α ∼ 0, gives an unbiased albeit extremely noisy reconstruction. Moreover, a reconstruction with small α exhibits strong radial oscillations. The oscillations vanish and the signal-tonoise improves considerably when one increases α, thereby enhancing the effect of the Wiener prior. However, a Wiener filtered reconstruction represents, as all Wiener filter estimators, a biased image of the true underlying signal. This bias can in our context be expressed as a radial PSF or response of the reconstruction for which we give an analytic expression in Eq. 81. The effect of the PSF is that a peak in the original density field will, on average, be stretched out and shifted in radial direction, see Fig. 7. The stretch and shift depend on the fiducial parameters of the survey and the tuning of the filter. The bias essentially reflects our inability to pin down exactly the redshift of a density peak due to the noise in the data.
Generally, the larger the α-tuning, the better the signalto-noise of the reconstruction but the more biased becomes the reconstruction. In the extreme case of overregularisation, one obtains essentially a 2-D reconstruction on the sky stretched out over the entire radial range lacking any radial information, providing the best signal-to-noise of a map of the projected matter distribution though.
As strategy, we suggest to choose a parameter α with moderate z-shift bias, see e.g. Figs. 8 and 9, and to relabel the redshifts of the reconstruction lens-planes according to the expected statistical average of the z-shift bias. At low redshift the radial PSF flattens out, as can be seen in these figures (there: z 0.3). A similar behaviour can be observed at the high redshift end of the survey. In these regions, we are hence unable to uniquely relabel the lens plane redshifts. At the low redshift range of that fiducial survey, all we can say is that the structure generating a peak in the reconstruction has to be located somewhere below z 0.3.
The reconstruction will be relatively noisy so that a final smoothing in 2-D (lens planes only) or 3-D, for example with a Gaussian kernel, has to be applied. The thereby obtained smoothed reconstruction is still a minimum-variance reconstruction of the true density field subject to the same smoothing, an invertible smoothing provided.
As can be seen from the implementation strategy, the algorithm is in essence a series of linear operations -projections, pixel-rescalings and convolutions -that have to be applied step by step to the pixelised input data. We would like to point out here that for the derivation of this algorithm we can equally imagine the data to be gridded on spherical shells in a spherical coordinate system; Fourier modes would be spherical harmonics on the unity sphere instead of waves on a plane. In fact, by doing so we will get exactly the same algorithm that, however, now tells us to linearly combine different lens-shells or to perform convolutions on the unity sphere. Therefore, the recipe described in Sect. A1 is a recipe for a full-sky reconstruction as well, for which the discussed optimisations and biases apply too.
We also gave an analytic estimate for the noise covariance of density pixels in the reconstruction, Eq. 84, for an idealised survey with homogeneous noise, no gaps and infinite extent. Several numerical algorithms for the covariance of a realistic survey are outlined in Appendix A3. A comparison of the noise to the signal of a SIS, Eq. 77, shows that the signal-to-noise of a matter halo quickly decreases towards increasing redshifts in a reconstruction (Fig. 10).
We noted that the algorithm (if implemented on a flat sky), as other similar reconstruction methods known in lensing, implicitly assumes that all cosmic shear signal originates from matter fluctuations inside the reconstruction area. If a shear of galaxy images inside the patch actually stems from a source outside the patch, this will introduce ghost images near the edges of the patch. In principle, the generation of significant ghost images can be detected by significant features in the mass density reconstruction based on the Bmodes of the shear field, which come along with the E-mode ghost structures, and have the same amplitudes. These edge effects can be reduced by putting an additional frame around the reconstruction patch with infinite noise.
This redshift limit increases for deeper surveys which we demonstrated by doing a reconstruction of a mock survey that broadly mimics the characteristics of a future spacebased lensing survey. Still, even boosting the number density of sources to 100 arcmin −2 and going deeper toz = 1.1 only moderately improves the performance of the reconstruction. For such a survey, structures in the map of S/N ≈ 4 or greater have accurate redshifts within ±0.1, below this level we have to expect false positives. Generally, we have to expect to miss out on higher density peaks, even for a powerful lensing survey as the one assumed here.
The situation can be expected to become better if a fullsky reconstruction with heavy smoothing is done to study -in contrast to structures on non-linear scales -density modes on very large scales: the signal-to-noise increases towards smaller ℓ, the bias effects in the reconstruction ease. We expect the inclusion of higher-order lensing distortions into the 3-D reconstruction to further improve the reconstruction, very roughly increasing the effective number of "cosmic shear sources" by a factor of N + 1 for N extra distortion moments. This, however, needs to be explored in more detail elsewhere.
The authors see cosmography as main application of the reconstruction technique. With 3-D lensing conventional 2-D lensing convergence maps can be supplemented to some degree by redshift information, mainly restricted to larger mass overdensities or to large density modes. This may not give us the full 3-D appearance of the mass distribution of, say, individual galaxy clusters -rather the 2-D projection stretched out in radial direction about the maximum-likelihood radial distance -but may enable us to disentangle physically unconnected structures along the line-of-sight providing a better understanding of the matter distribution inside a 3-D volume. This can be vital for a more sophisticated modelling later on.
As an obvious application, the 3-D data cubes may be used to make catalogues of matter overdensities endowed with redshifts without any assumptions about density profiles. Estimating the density profiles of individual structures, essentially projected to 2-D, in the 3-D lensing map is conceivable, too, but only if the full reconstruction PSF, including the transverse direction, is accounted for. Owing to the map biases an estimate for the mass of structures should be obtained by directly fitting a density model to the shear data once they have been identified in the non-parametric 3-D reconstruction.
In principle the 3-D maps can be compared to the 3-D distribution of galaxies as derived from galaxy redshift surveys. In this case, again, care has to be taken for the biases and limited radial resolution in the lensing map. For a fair comparison, this may be done by convolving the map of galaxies with the radial PSF of the lensing maps. For a statistical comparison of galaxy and (lensing) matter maps required for studying the galaxy biasing (e.g. Simon et al. 2007), however, we would suggest to forward and fit a 3-D model of matter density fluctuations and galaxy biasing to the 3-D shear tomography and its cross-correlation with galaxy positions, rather than producing a biased Wiener filtered 3-D matter map that is cross-correlated with the 3-D distribution of galaxies. In the best case, with all biases in the 3-D mass map being fully understood, both ways should give same results anyway.
The same holds for attempts to estimate matter density power spectra from the 3-D lensing maps. As shown in this paper, due to low signal-to-noise the power spectra on the lens planes will be strongly influenced by the Wiener prior which has to be taken into consideration for an estimate of the density fluctuation power. Even if this is done properly, it should only give a result equivalent to fitting a 3-D model for the shear tomography correlation functions which is straightforward in comparison .
Casting the estimator into this form is actually more convenient from a numerical point of view if we have to consider a non-diagonal noise covariance, N, e.g. in the case of intrinsic alignments. Here, we split the noise covariance, N = N d +No, into the sum of diagonal, N d , and off-diagonal, No, elements; the diagonal matrix N d is easily inverted and multiplied with a vector.
For an implementation of this estimator we proceed in three steps. For the first step, see the underbraces, the binned shear data γ is multiplied by the inverse of the diagonal elements of N. This yields the intermediate result "x" that is actually a set of shear grids in which each grid point is weighed by its noise variance. For the second step, x has to be multiplied by the matrix M −1 . We employ the method of conjugate gradients to invert the equation My = x with respect to y (Press et al. 1992). This essentially boils the matrix inversion down to a series of multiplications with M until the solution converges. Again, multiplying by the matrix M is treated as a one-by-one multiplication by the series of individual matrices contained within the square brackets. As initial guess for y for the iterative algorithm one commonly uses y = 0 which appears to work fine in the context of this paper. Finally, for the last step y is multiplied by S δ Q t P t γκ . Multiplying a matrix with a vector means the application of a linear operation on the input vector. Here, each matrix has a different effect on the vector that is hard-wired as vector-in/vector-out sub-routine individually. The multiplication with Q (or its transpose Q t ) is easily implemented as it is just summing k along every line-of-sight θ k . Multiplications with inverses of diagonal matrices are also easy as mentioned above.
The time consuming matrix operations are connected to Pγκ, S δ and No, if considered, which denote convolutions of the input vector. Only in the situation where we are using radial Wiener filtering, S δ x will become an easy operation requiring only sums along line-of-sights analogous to Qx. The most effective way of doing convolutions uses Fast-Fourier-Transforms (FFTs) as an efficient application of the convolution theorem (Press et al. 1992). For that we need to FFT the different layers of x (either lens planes or source bins) separately and multiply the Fourier coefficients with the power spectrum of S δ , the Kaiser-Squires kernel D(ℓ), in case of Pγκ, or the power of intrinsic alignments; the power spectra are functions of lens planes or source bins. The window function of the grid cells has to be taken into account as well. The transpose of Pγκ corresponds to D * (ℓ). After that we re-FFT to real space. Note that the whole chain of operations PγκQS δ Q t P t γκ can be done completely in Fourier space with no back and forth FFTing in between.
The FFT method has two drawbacks: i) It requires grids sizes which are powers of two so that not any arbitrary fieldgeometry can be processed. ii) Also, as well known, we will have edge effects as FFTs are assuming periodic boundary conditions (aliasing). Regarding i) we have to find the smallest bounding box to enclose the whole field of view at the desired resolution. Resolutionwise a number of grid cells of 128 or 64 for each coordinate axis is usually enough. With respect to ii) we can use the zero-padding technique or downweighting towards the edges (see next section) to mitigate the aliasing.
There is an alternative technique known in the literature that speeds up convolutions but works solely in real space (Padmanabhan et al. 2003), does not suffer from FFTaliasing and the geometry constraint and, in fact, does not require griding. This approach is a real alternative to FFT but considerably slower as we found. As both techniques yield essentially the same reconstructions, apart from the edges, we propose to stick to the FFT technique in this context.
Processing 25 lens planes and 40 source redshift bins, all with a grid size of 128 × 128 pixel 2 (the data vector, complex, and matter density vector, complex (E-, B-mode combined), have sizes of 1.3×10 6 and 8.2×10 5 elements, respectively) requires for the algorithm ∼ 60 MByte of computer memory (double precession). The computation time for the second step on a AMD Athlon 64 processor with 2.4GHz is approximately five seconds for every iteration of step 2, which is the main bottleneck of the algorithm. For noisy data and (full, i.e. α = 1) Wiener filtering usually about five iteration steps are needed, whereas a reconstruction with no prior or heavily tuned-down Wiener filter (α ∼ 0.01) requires up to ∼ 150 iteration steps, depending on the desired accuracy, because less regularisation is applied to the data so that more details are fitted. The benchmark parameters used here are usually, with contemporary real data, unnecessary. For that 64 2 -grids, N lp ∼ 10 lens planes and Nz ∼ 20 source bins are absolutely sufficient bringing down the computation time of one full reconstruction to roughly ten seconds.
After the algorithm has converged, we run a final (transverse) smoothing of the lens-plane grids with a Gaussian kernel. The smoothing is usually required as a down-tuned Wiener filter is applied that produces noisy images. The blurring can be done efficiently by again employing the FFT convolution or a kdtree implementation.

A2 Fourier space reconstruction
If an idealised survey is a good assumption for particular data, the Fourier space filter Eq. 51 is, due to the small matrices involved, a very fast and efficient way to perform a 3-D reconstruction completely in Fourier space. In particular, the matrix inversions can be done explicitly. To apply this filter, we FFT the gridded source z-bins, arrange the shear Fourier coefficients for each ℓ inside one vector,γ, and perform Eq. 51 for every angular mode. The matrices involved have only sizes of the order of Nz or/and N lp , which are easily and quickly inverted. The Fourier reconstruction is eventually back-FFTed to acquire the real space mass density field. Finally, the lens-planes need to be smoothed, since we are usually using a tuned down Wiener filter (α = 0.01, for instance). This smoothing can be done in real space or, more efficiently, directly in Fourier space immediately after the Fourier space reconstruction.
Usually for ℓ = 0 a regularisation (λ = 0), especially for the no-prior MV-estimator, is not necessary in a Fourier space reconstruction because the reason for the singularity of the filter is the mass-sheet degeneracy that affects solely the ℓ = 0 modes. Here, we lift this degeneracy by setting all ℓ = 0 Fourier coefficients to zero which amounts to saying that the mean density contrast on each lens plane is exactly zero. For a large enough field-of-view this should be a reasonable assumption.
In order to reduce the FFT edge effects, the real space grids should, initially, be downscaled by sin (π∆/2∆0) within a strip of width ∆0, where ∆ is the shortest distance of a pixel to the grid edge.
The noise covariance is estimated by means of the wellknown shot-noise formula [σ (i) ǫ ] 2 /ni (Sect. 3.3), ideally assuming that the grid cells are homogeneously filled with sources of roughly equal intrinsic ellipticity distribution.
A complete reconstruction on a Fourier grid for the same parameters (but 128 × 128 grids) as in the foregoing section takes less than a second.

A3 Covariance of statistical uncertainties
An estimation of the statistical errors, or more general the covariance, of pixel densities in the reconstruction is needed to convert the density reconstruction into a signal-to-noise 3-D map. As a lowest-order estimate we suggest to use the analytical solution, Eq. 84, that is based on homogeneous noise of the sources. It offers the r.m.s.-error for each lensplane and the correlation of errors along each line-of-sight.
As a final smoothing of each lens-plane will usually be necessary, the estimate has to be rectified slightly. Lets us assume that the errors, σ δ , of the unsmoothed map are uncorrelated, which will be roughly correct for a strongly downtuned (small α) filter. By smoothing a lens-plane with a normalised kernel, K(x), one will combine an effective number of N eff ≡ˆR R 2 d 2 θK(|θ|)˜2 R R 2 d 2 θ[K(|θ|)] 2 (A4) pixels with uncorrelated errors σ δ . The combined error of the weighted average (every smoothed pixel) will therefore be σ δ / √ N eff . For a Gaussian kernel with width σ one finds N eff = π 2 σ 2 . The truncation of the kernel due to the edges of the map can be accounted for using the formula for N eff and setting K(x) = 0 for areas outside the map. As the pixels one is smoothed over are slightly correlated, σ δ / √ N eff will give only a lower limit for the true error (upper limit for the signal-to-noise), whereas σ δ would be the noise upper limit being reached if all pixels were 100% correlated.
In reality, one encounters surveys with varying source number densities and intrinsic source shot noise. One may desire to get a better estimate for the signal-to-noise of a pixel than provided by the foregoing formula, which is given for all pixels by where X −1 = Q t P † γκ N −1 PγκQ, W = F † S δ [X −1 S δ + α1] −1 and F is the smoothing matrix (unity matrix in absence of final smoothing). However, the full covariance imposes a huge, if not almost impossible, computational task because of the enormous size of the involved matrices.
What can easily be done, though, is to compute single matrix elements of the covariance: where Fi is defined as the ith column of the matrix F, which is the smoothing weight of any grid pixel relative to a fixed grid pixel i. This vector is filled with many zero because we are only smoothing within the same lens-plane. As indicated, the recipe is to compute the vectors v i/j in a first step and then v † i N −1 vj in a second step. In order to attach error bars to single δMV one requires only one vector vi = vj .
The definition of v i/j is almost identical to δMV in Eq. A1 so that the algorithm discussed in Sect. A1 applies without much change. On the other hand, this means the effort to obtain a matrix element of NMV is roughly twice the effort (Only off-diagonals; the same effort for diagonals) to remake a whole 3-D reconstruction, so we may only wish to work out the exact signal-to-noise or correlation of errors for a selected small number of pixels such as the peaks in the mass map that are suspected to belong to galaxy clusters.
A third and most effective alternative for obtaining a signal-to-noise map is to randomise the ellipticities of the sources and to make a full reconstruction. With the algorithm outlined in Sect. A1 this can be done hundreds of times in a reasonable time, maybe running parallelly on different computers. The variance in the density contrast between the different noise realisations would serve as estimator for the noise level in the reconstruction.

APPENDIX B: SHEAR-INTRINSIC ALIGNMENT CORRELATIONS
The estimator Eq. A2 is the minimum variance solution, or the most likely solution for s in the Bayesian sense for pure Gaussian statistics, for a given data vector d provided that there is no correlation between noise and signal, i.e. sn † = 0. In the context of gravitational lensing, however, these correlations can occur if, for example, the shear signal itself is correlated to the intrinsic alignment of sources (Hirata & Seljak 2004). We ignore this effect for 3-D mass reconstructions in this paper but would like to sketch how it can, in principle, be included by extra terms in the minimum-variance estimator for s.
For that we assume that the signal/noise-correlations (here: correlations between the lens plane matter density and the intrinsic ellipticities of sources) between the various source z-bins are known either from measurement or from a physical model, U ≡ sn † . The minimum-variance Ansatz seeks to find a linear transformation, H, that minimises the residual matrix fi (Hd − s)(Hd − s) † fl (B1) for all possible noise and signal configurations. For the simpler case that U vanishes, this residual matrix is, leaving out terms independent of H, on average which is minimised for If we now drop the condition of a vanishing U, the variance B1 becomes more generally If we substitute in the last expression we recover exactly the variance B2 with the already known minimising solution Eq. B3. Undoing the former substitution in B3 (the substitution is independent of H) gives the more general minimum-variance that now also accounts for signalnoise correlations: Here, we already included the rescaling of S → α −1 S needed for the α-tuning (Saskatoon filter) supplemented by U → √ α −1 U to be consistent with the α-tuning (the signalnoise correlations are unchanged by changing α). We find that including the signal-noise correlation adds three new terms to Eq. B3. In practice, this estimator can be tackled by the same numerical techniques as outlined in Appendix A.
This paper has been typeset from a T E X/ L a T E X file prepared by the author.