CFHTLenS: Mapping the Large Scale Structure with Gravitational Lensing
Abstract
We present a quantitative analysis of the largest contiguous maps of projected mass density obtained from gravitational lensing shear. We use data from the deg covered by the CanadaFranceHawaii Telescope Lensing Survey (CFHTLenS). Our study is the first attempt to quantitatively characterize the scientific value of lensing maps, which could serve in the future as a complementary approach to the study of the dark universe with gravitational lensing. We show that mass maps contain unique cosmological information beyond that of traditional twopoints statistical analysis techniques.
Using a series of numerical simulations, we first show how, reproducing the CFHTLenS observing conditions, gravitational lensing inversion provides a reliable estimate of the projected matter distribution of large scale structure. We validate our analysis by quantifying the robustness of the maps with various statistical estimators. We then apply the same process to the CFHTLenS data. We find that the 2points correlation function of the projected mass is consistent with the cosmological analysis performed on the shear correlation function discussed in the CFHTLenS companion papers. The maps also lead to a significant measurement of the third order moment of the projected mass, which is in agreement with analytic predictions, and to a marginal detection of the fourth order moment. Tests for residual systematics are found to be consistent with zero for the statistical estimators we used. A new approach for the comparison of the reconstructed mass map to that predicted from the galaxy distribution reveals the existence of giant voids in the dark matter maps as large as degrees on the sky. Our analysis shows that lensing mass maps are not only consistent with the results obtained by the traditional shear approach, but they also appear promising for new techniques such as peak statistics and the morphological analysis of the projected dark matter distribution.
keywords:
dark matter  largescale structure of the Universe  gravitational lensing1 Introduction
Gravitational lensing is a powerful tool for the study of the dark matter distribution in the Universe. The statistical analysis of the distortion and magnification of lensed galaxies provides unique information on structure formation processes. For example, two point statistics of the shear or convergence can be used to constrain the dark matter power spectrum, the growth of structure, and cosmological parameters (see Munshi et al., 2008; Hoekstra & Jain, 2008, for recent reviews). Over the past ten years, multiple groups have reported improved constraints on the mass density parameter and the power spectrum normalisation using the shear correlation function. However, it is well known that gravitational lensing contains a lot more information than the amplitude and shape of the mass power spectrum. The distortion (shear), in particular, can be used to reconstruct the projected mass density, or mass maps, thus making the information on the distribution of dark matter available in a different form. The use of mass maps as a cosmological tool has received little attention so far.
The purpose of this paper is to explore the extent to which mass maps can access the cosmological information not captured by the two points statistics. Several groups have explored different theoretical routes beyond two points statistics. For instance, higher order shear measurements are a sensitive measure of the gravitational collapse process through mode coupling in the nonlinear regime (Bernardeau et al., 1997; Van Waerbeke et al., 2001; Takada & Jain, 2003; Kilbinger & Schneider, 2005), and can also be used as an indicator of nonGaussian nonlinearity in the primordial dark matter distribution (Takada & Jain, 2004; Valageas et al., 2005). More exotic statistical estimators involve global statistical tools, such as the Minkowski functional, as opposed to local measurements based on the shear correlation functions. Morphology of large scale structures (Mecke et al., 1994) is also a probe of the nonlinear processes in action during structure formation (Sato et al., 2001).
The peak statistics is also an estimator that can be used to probe the number of dark matter haloes as a function of redshift and mass (Van Waerbeke, 2000; Jain & Van Waerbeke, 2000; Kratochvil et al., 2010; Maturi et al., 2010), and it can potentially constrain the halo mass function (Yang et al., 2011). The effect of the large scale dark matter environment on the galaxygalaxy lensing signal could be investigated; for instance mass maps could be used to quantify to what extent specific galaxy properties (e.g. star formation rate, dust content, stellar population) depend on the global dark matter environment such as a supercluster or a giant void.
Most of these global statistics, like peaks and structure morphology, are nonlocal, and therefore cannot be expressed as a combination of a subset of moments of the shear, which are, by definition, local measurements. The study of global features of the mass distribution require the mapping of the full 2dimensional mass. Dark matter maps can also be used to understand global features in the baryon/dark matter relation via crosscorrelation with maps from different wavelengths surveys, whether SunyaevZel’dovich, Xray, or atomic hydrogen. The use of mass maps is therefore complementary to the analysis of the shear and magnification twopoint correlation functions. This motivates the need to have reliable mass maps from gravitational lensing data.
Mass reconstruction is a common tool for the study of galaxy clusters (Hoekstra et al., 2000; Clowe et al., 2006a, b; Heymans et al., 2008); it is mostly used to evaluate qualitatively the matching between the baryonic and dark matter distributions, and to explore the possible existence of dark clumps. A few rare cases suggest their possible existence (Erben et al., 2000; Mahdavi et al., 2007; Jee et al., 2012), although the significance of these detections is currently unclear. Mass reconstruction in cluster environments concern regions of relatively high lensing signal and small angular scale. Mass reconstruction of weaker lensing signal, beyond galaxy cluster scales, was performed in Massey et al. (2007) with the deg area of the COSMOS survey. However, our work shows that sampling variance is still important at even larger scales, implying that a relatively small survey like COSMOS does not provide a fair sample of the dark matter distribution.
In this paper we use the deg CFHTLenS lensing data (Erben et al., 2012; Heymans et al., 2012b; Miller et al., 2013; Hildebrandt et al., 2012) to perform lensing mass reconstruction over angular scales of several degrees. We do not attempt to perform a 3dimensional mass reconstruction (Massey et al., 2007; Simon et al., 2012; VanderPlas et al., 2011; Leonard et al., 2012) because such reconstruction using groundbased data would only yield meaningful results for galaxy clusters relatively close in redshift, , or for very massive and rare clusters at higher redshift, (Simon et al., 2009). We explore the measurement of convergence statistics from gravitational lensing maps and the connection between dark and baryonic matter up to a few degrees across. In particular, we focus on the importance of the effect of noise in the reconstruction, and how to account for it in the interpretation of the measurements. To this end, we use raytracing simulations to test our reconstruction and analysis procedure.
This paper is organized as follows. Notation and definitions are described in Section 2. The analysis of simulated reconstructed mass maps is discussed in Section 3. In Section 4, we present the results from the CFHTLenS data and our new approach for the comparison of dark matter and baryonic matter. We present our conclutions in Section 6.
2 Map making and cosmology
2.1 Mass reconstruction
Mass maps are proportional to the projected mass convolved with the lensing kernel. They can be constructed from the shear measurement. A galaxy at position on the sky is characterised by its redshift and its shear components , where . In this paper, we use the flat sky tangent plane approximation. The relation between and the convergence involves the gravitational deflection potential , defined as:
(1) 
where and . The function is given by a lineofsight integral of the 3dimensional matter gravitational potential via:
(2) 
The comoving radial distance at redshift is given by , and is the corresponding angular diameter distance. The comoving distance at the source redshift is . The 3dimensional gravitational potential depends on the mass density contrast via the Poisson equation:
(3) 
where is the Hubble constant, is the matter density, and is the scale factor. The convergence is defined as the lineofsight projection of the 2dimensional Laplacian of :
(4) 
where , and the derivative of the gravitational potential along the line of sight vanishes on average due to the Limber approximation (Kaiser, 1998). The observable quantity is the reduced shear :
(5) 
Due to the intrinsic ellipticity of galaxies, galaxy shapes provide a noisy estimate of . This estimator, , is given by:
(6) 
The intrinsic ellipticity is a generic term that contains a contribution from the intrinsic galaxy shape and from the measurement noise. The ensemble average is an unbiased estimator of the reduced shear (Kayser & Schramm, 1988). In the case of CFHTLenS, a complete description of the measurement noise is given in Heymans et al. (2012b) and Miller et al. (2013).
Weak lensing studies assume that both the shear and convergence amplitudes are much smaller than unity, () so that one can use the weak lensing approximation () for mass reconstruction. In that case, Eq.(6) reduces to , and it follows that by comparing Eq.(1) and Eq.(4), one can reconstruct the convergence from the observed shear in Fourier space (Kaiser & Squires, 1993). For the weak lensing approximation, we first write the relation between the shear and convergence Fourier components and , where we define the wavevector as the 2dimensional analogue of in Fourier space:
(7) 
Following Kaiser & Squires (1993), the shear data are first regularised with a smoothing window and then Eq.(7) evaluated from the smoothed Fourier components. Since real data also contain masked regions ( e.g. bright stars and sometimes area missing from the detector), the smoothing takes into account the number of pixels being masked within each smoothing window, so that the resulting averaged quantity in that window is not biased. The practical procedure for obtaining a mass map from reduced shear is as follows:

The data are first placed on a regular grid and then smoothed. At pixel location , the smoothed ellipticity is given by:
(8) where is a normalised Gaussian smoothing window:
(9) and is the weight associated with the measurement noise of . Details on the weight specific to CFHTLenS data are given in Section 4.2.

The reconstructed mass is obtained from Eq.(7). An absence of galaxies should result in a pixel value of zero. However, note that since the Gaussian filter has infinite spatial extension, there are no pixels in the final grid with a value of zero. There are, however, pixels with higher noise (fewer objects) than others, especially when a mask overlaps with the central region of the filter. The noisiest pixels in the reconstructed map are removed if the effective filling factor within the Gaussian window is below percent. This point is discussed later in Section 3.2, where the noise in the reconstructed maps plays an important role in the cosmic statistics measurement from the convergence.
The alternative to the Kaiser & Squires (1993) mass reconstruction approach is the full nonlinear mass reconstruction as described in Bartelmann et al. (1996). A detailed comparison between the two methods is left for a future study. For the purposes of this paper, the KS93 reconstruction is sufficient because we are not trying to recover the lensing signal in the nonlinear regime at subarcminute scales. We note, however, that the ability of the nonlinear mass reconstruction to recover large scale structure statistics was demonstrated in Van Waerbeke et al. (1999). Therefore, moving to a full nonlinear reconstruction does not represent a conceptual challenge.
As shown in Crittenden et al. (2002), the shear can be split into and modes, which correspond to the curlfree and curl shear components respectively. In the absence of residual systematics, the scalar nature of the gravitational potential leads to a vanishing mode. As shown in Schneider et al. (1998), the split between and modes is performed by applying the transformation to , which is the same as rotating each galaxy by degrees. Our mass maps are reconstructed using the and modes and the corresponding convergence is called and . As we will see later on, it will be necessary to distinguish between the reconstructed (noisy) map, and the underlying true convergence value. We therefore introduce three additional convergence terms. We call the reconstructed convergence such that:
(10) 
where is the true underlying signal and is the reconstruction noise. This equation is the convergence equivalent of Eq.(6) in the weak lensing regime. We also introduce the convergence , reconstructed from the galaxies rotated by degrees such that:
(11) 
where should be consistent with zero if the residual systematics are negligible. The validity of this statement will be verified statistically.
2.2 Cosmology
The noisefree convergence map , smoothed with a Gaussian filter of size , can be expressed in terms of the 3dimensional mass density contrast :
(13)  
where is the normalised Gaussian window given by Eq.(9) and represents the radial distance at infinite redshift. The function accounts for the distribution of source redshifts:
(14) 
where is the source redshift distribution. Eq.(13) represents the convergence map obtained from the reconstruction process after the shear data have been pixelated and smoothed as described in Section 2. Therefore, the positions on the sky are pixels on the regular grid where the map has been computed. The remainder of this paper will focus on subdegree scales where the flat sky approximation applies. However, we will use the exact full sky formalism for all the lensing quantities, which has the advantage of providing a robust numerical integration via a direct, one pass, summation of the wavevectors. Furthermore, a conservative cut at is used for both the second and third order statistics predictions. The second order moment of the convergence can be expressed from the convergence power spectrum :
(15) 
where are the multipole moments of the smoothing window. For the Gaussian window given by Eq.(9), the multipole moments are given by . The power spectrum can be derived from Eq.(13), assuming that the small angle and Limber approximations apply (MiraldaEscude, 1991; Kaiser, 1992):
(16)  
(17) 
The 3dimensional mass power spectrum is defined as:
(18) 
The third order moment can also be expressed analytically using the convergence angular bispectrum , the smoothing filter multipoles , and the Wigner 3j symbols:
(19)  
(20)  
with the angular bispectrum given by:
3 Simulations
3.1 The CFHTLenS data
The purpose of the following Section is to validate our mass reconstruction approach using mock catalogues that replicate the true CFHTLenS observing conditions. CFHTLenS spans a total survey area of deg, covered with a mosaic of individual pointings observed by the one square degree imager at the CanadaFranceHawaii Telescope. The survey consists of four compact regions called W1, W2, W3 and W4, which cover approximately , , and deg respectively. Details on the data reduction are described in Erben et al. (2012). The effective area is reduced to deg by the masking of bright stars, artificial and natural moving objects, and faulty CCD rows. The observations in the five bands of the survey allow for the precise measurement of photometric redshifts (Hildebrandt et al., 2012). The shape measurement with lensfit is described in detail in Miller et al. (2013).
3.2 Mock catalogues
The mock catalogues are constructed from a mixture of real data (the noise, galaxy position, and masking structure of CFHTLenS) and simulated data (the shear signal from Nbody simulations). The Nbody simulations used in this work are described in HarnoisDéraps et al. (2012). The procedure for generating mock catalogues is as follows:

Projected shear and convergence maps are constructed from the combination of redshift slices for all source redshifts available in the simulation
^{3} . The maps are sampled on a grid and cover deg each, which corresponds to a native pixel size of arcmin. These maps represent individual ‘tiles’. The CFHTLenS fields W1, W2, W3, and W4 are covered with the maximum number of nonoverlapping tiles that fit within their respective areas. Given the area of the four CFHTLenS mosaic fields, one can only cover W1 and W3 with tiles each, and W2 and W4 with one tile each. The final simulated area is therefore deg; which is similar to the deg covered by CFHTLenS. Note that since each of the tiles used in our simulations is an independent line of sight, we expect sampling variance to affect the mock catalogues less than the four CFHTLenS fields, where only four fields are statistically independent. 
The 2D positions on the sky of all galaxies in CFHTLenS are preserved. The orientation of each galaxy’s ellipticity is randomized, but its amplitude is held fixed as this is used as the only source of shape noise in the mock catalogue. The redshift of each galaxy is resampled using redshift probability distribution obtained from the photometric redshift distribution function. We refer the reader to Benjamin et al. (2012) for the details. The galaxies are then placed on the mosaic tiles and each galaxy’s redshift determines uniquely the simulated shear and convergence. This is calculated using Eqs.(1) and (4) with the reduced shear computed in Eq.(5). The final ‘observed’ ellipticity in the mock catalogue is obtained by Eq.(6). The masks are applied to the mock catalogue (i.e. no ‘missing galaxy’ has been added to fill in the gaps), as are all the other characteristics of the survey (e.g. the weight associated with galaxy shape measurements is also preserved).

Following the procedure described in Section 2, the mass reconstruction is performed on a regular grid, which is a rebinning of the simulations native pixel grid. The mass reconstruction is performed for different smoothing scales, with a radius ranging from approximately to arcminutes.

For some of the following tests performed on simulations, we will be using the noisefree convergence map, obtained from the stacking of convergence maps at difference source redshifts using the redshift distribution from (ii) above. These noisefree maps are obtained directly from the simulated light cone, prior to the construction of the mock catalogues. For this reason, and to make a clear distinction with the quantity , we call the stacked noisefree convergence map .
Figure 1 shows an example of a mass reconstruction of one the tiles. The background image with the white contours shows a noisy mass reconstruction; the white contours represent the levels, where is the convergence rms determined from the map. The black contours represent the levels on the noisefree convergence map . This figure illustrates qualitatively how realistic shot noise (ellipticity noise) noticeably affects the position, amplitude, and sometimes even the presence or absence of peaks in the reconstructed map. Many reconstructed peaks do not match a real mass peak and the converse is also true. Although a quantitative analysis of peaks is left for another paper, this illustration is consistent with earlier work (Van Waerbeke, 2000; Athreya et al., 2002) showing that individual peaks are relatively noisy objects, with a higher chance of being a coincidence ( per cent for a peak) than what Gaussian statistics would predict in the field. Masks are also shown in Figure 1 as the black areas with sharp boundaries, which shows that our mass reconstruction procedure does not generate catastrophic edge effects near the masks. The rest of this Section explores the reliability of the mass reconstruction quantitatively.
3.3 Analysis of the mock catalogues
In this Section, we use the simulations and their reconstructed lensing maps to verify that the level of shape noise and masks in the CFHTLenS data will not introduce systematic errors in the statistical properties of the reconstructed convergence map. This will be performed by running two distinct tests.
For the first test, we focus on the crosscorrelation between two maps and , where and can be any of the signal , noisefree , galaxy rotated , or randomized maps. The crosscorrelation map is given by:
(24) 
where denotes the zero lag value of the autocorrelation map. This definition guarantees that is normalized, i.e. the central pixel of is equal to one if . is then azimuthally averaged within annuli as a function of the distance from the central pixel; this quantity is called .
Figure 2 shows for different combinations of reconstructed and noisefree maps for the Gaussian smoothing scale arcmin. The filled triangles show , the crosscorrelation profile between the reconstructed mass map and the galaxy rotated map, averaged over the lines of sight of the mock catalogue (see Section 3.2). The error bars, which represent the error on the average profile, show that is consistent with zero. The solid area shows the scatter of averaged over random realizations of . This quantity is also consistent with zero, which confirms that the field boundary, edges, and masks do not generate systematic effects, even for pure noise reconstructions. The filled circles show the crosscorrelation profile between the mass reconstruction and the noisefree map . The error bars illustrate the dispersion of the average over the lines of sight; the zerolag crosscorrelation coefficient is only , which is another indication that mass maps are noisy, especially for small smoothing scales. Larger smoothing scales, not shown here, lead to a stronger crosscorrelation amplitude, but also a larger correlation length. This first test demonstrates the absence of systematic effects (e.g. spurious peaks) around masks and lowdensity areas; if these were present, we would expect a significant nonzero residual crosscorrelation between the reconstructed mass map and or the pure noise reconstruction.
The second test consists of looking at various higherorder moments of the convergence. Like the traditional measurement of the second order moment of the shear (here ), we want to measure the true moments of the convergence , with . A measurement of the convergence moments from noisy reconstructed mass maps has never previously been reported. However, there is a complication with the measurement of convergence moments: the convergence can only be measured within a smoothing window, unlike the shear, which is given for each individual galaxy. Therefore, the convergence noise is correlated for different points on the grid, with the result that the moments of and are not equal. A shot noise removal, or denoising, is therefore required in order to remove this bias from the observed convergence moments .
We choose a direct approach for the denoising procedure, which consists of removing the correlated noise contribution from the moments measured on . Fortunately, two independent residual systematics tests can be performed; one on the ‘mass map’ of the rotated galaxies , and the other by comparing the denoised moments of to the moments of the noisefree maps . Figure 2 shows the absence of crosscorrelation between and the noise which means that, to first approximation, the signal and noise can be treated as statistically independent (Van Waerbeke, 2000). We use this lack of correlation to derive relatively straightforward relations between the moments of and . For instance, the observed second order moment is the quadratic sum of and . The extension of this relation to higher order moments defines our denoising procedure. Moments of the noise maps , needed for this step, are measured from a large number of pure noise reconstructions. The relations below show how, for each smoothing scale , the observed moments of are related to the true moments of and the moments of the noise map :
(25)  
(26)  
(27)  
(28) 
For clarity, denotes the moment average over one map, is averaged over several noise maps , and is the denoised moment. We limit our analysis to higher order moments up to the fifth order. Beyond the fifth order, we find that the measurements become too noisy with the CFHTLenS data. The Appendix includes the expressions of the residual systematics moments, i.e. all the combinations of and modes that we expect to vanish if residual systematics are negligible for . As the measurements will show, the residual systematics are indeed consistent with zero. Our denoising approach automatically takes care of varying noise across the field because the noise rms for different positions on a map is preserved for the different noise realizations. This is due to the fact that we keep the galaxys position and shape noise unchanged. A further improvement of the denoising technique was implemented by taking into account the fraction of pixels not masked within each smoothing window. We achieve this by constructing a filling factor map computed from the smoothed grid with pixels set to one, and then masking and smoothing with the same window that was applied to the lensing data for the mass reconstruction. Using this filling factor map, regions around masks in the lensing reconstructed maps are downweighted without biasing the signal. In order to eliminate the noisiest regions, we apply a cutoff of per cent such that pixels below this threshold are excluded from the moments analysis. We verified that the application of this cutoff has a negligible effect on the measurement with a cutoff varying from to per cent. Of course, high cutoff values boost the noise and it becomes difficult to evaluate changes in the signal. This procedure guarantees that pixels in the mass map contribute proportionally to the unmasked area located in a smoothing window. We will compare the denoised moments of with the noisefree moments measured on .
Figure 3 shows the convergence moments as measured from the mock catalogues using the denoising technique described earlier. We see that the denoised moments (filled circles) are consistent with the noisefree moments . Note that the noisefree moments represent the true answer in the sense that they are directly computed from the Nbody simulations, without mass reconstruction, and they therefore contain no noise, no discrete sampling and no mask. On the other hand, the moments are denoised moments obtained from the noisy mass reconstruction maps that include masks. The two moments agree, which demonstrates that our mass reconstruction procedure leads to reliable mass maps that preserve the statistical and cosmological information. In Figure 3, the denoised moments errors were computed from the variance between the lines of sight and divided by . The noisefree moments errors are shown by the filled regions around the solid lines, and are also divided by . The open symbols show, for the different moments, the residual systematics obtained from all the possible combinations of the reconstructed maps and , where the denoised moments were computed using the expressions derived in the Appendix. All moments involving one or more of the rotated galaxies maps are consistent with zero, showing that the mode is consistent with pure noise once it has been denoised.
These results from mock catalogues validate our approach for the CFHTLenS data and demonstrate that for the deg of the survey, our reconstruction process is stable. Furthermore, we have shown that realistic masking geometry does not alter the reconstruction, and therefore we are able to reliably quantify some of the most basic statistics of the projected mass density.
3.4 Beyond the moments: Convergence PDF and peak statistics
The first few convergence moments give only a partial description of the histogram of the convergence, which is also called the point Probability Distribution Function (PDF). However, the point PDF, along with its extensions to higher order, contains additional information about the moment hierarchy, and therefore about the structure formation process. For instance, Juszkiewicz et al. (1995) demonstrated that a combination of different moments of the PDF, through the Edgeworth expansion, probes different aspects of the gravitational collapse. General characteristics of the PDF, such as the height and the minimum of the convergence histogram, are important features that are not easily captured by convergence moments. Therefore, it would be interesting to use the convergence PDF itself as a cosmological probe. This possibility has been theoretically explored (Valageas & Munshi, 2004; Munshi et al., 2004) using the aperture filtered shear, the same filter which transforms the shear into a local, scalar quantity (Schneider, 1996). Doing this analysis on the convergence PDF, instead of the aperture filtered shear, would exploit the fact that the long wavelength modes are preserved; this analysis with the CFHTLenS mass maps is left for a forthcoming study.
Nevertheless, we can already illustrate with our simulations the expected level of point PDF signaltonoise. Figure 4 shows the measured point PDF for compared to the point PDF. The average point PDF obtained from pure noise reconstructions (with error bars) is also shown, and the mode PDF is consistent with pure noise. This result is in agreement with Figure 3, which shows a negligible mode for the convergence moments. One can also clearly see that the cosmological signal broadens the PDF compared to pure noise (or mode) PDFs. This point PDF study opens up a new statistical analysis opportunity for gravitational lensing surveys, which remains to be exploited. Alternative probes for cosmology using the convergence PDF could be constructed, such as using the minimum value of the convergence found in voids (modulo the convolution by the noise) as a direct measurement of the mass density parameter , which is relatively insensitive to the nonlinear clustering.
Another interesting use of convergence maps is the study of peak statistics. The idea of using peaks in dark matter maps to constrain the halo mass function was first proposed in Van Waerbeke (2000); Jain & Van Waerbeke (2000). This analysis relies heavily on the properties of peak statistics in Gaussian noise developed for the cosmic microwave background (Bond & Efstathiou, 1987). The peak statistic theory has recently been pursued further (Maturi et al., 2010; Marian et al., 2012). Peak statistics provides a higher level of statistical analysis of maps, where we expect the sensitivity to cosmology and residual systematics to be different from the statistical analysis using moments of the shear or convergence. It is important to distinguish between the maps constructed in this paper and the maps being discussed in Marian et al. (2012), which use the aperture mass filter, originally defined in Schneider (1996). The aperture mass filter is a passband filter, as opposed to the Gaussian (or tophat) smoothing window used in our study, which is a lowpass filter. The lowpass filter preserves the large scale modes, hence keeping visible the voids and large overdensities, unlike what happens with a passband filter (see Section 4.4 below). The peak statistics of a lowpass map allows the study of the dark matter halo mass function as a function of the larger scale dark matter environment, which is essentially not possible with a bandpass filter. The two approaches are complementary and it would be beneficial in the future to unify them into a single approach to making mass maps. Peak statistics, convergence PDF, and the morphological analysis of large scale structure from the CFHTLenS data are left for future studies.
4 Results from CFHTLenS
Residual systematics in the shear signal, based on the selected fields that passed the systematics tests (the ”good” fields) are given in Heymans et al. (2012b). In total, MegaCam pointings (out of ) passed the residual systematics tests. The convergence moments presented in this paper are based on these good fields, while the mass reconstruction is performed on all the fields. Each galaxy in the CFHTLenS catalogue has a shear estimate and a weight (Miller et al., 2013), a calibration factor (Heymans et al., 2012b), and a photometricredshift probability distribution function (Hildebrandt et al., 2012). In the following sections we describe how these elements are integrated into the mapmaking process.
For the lensing mass reconstruction, source galaxies are selected in the redshift range , which is sufficiently broad to avoid the complication of intrinsic alignment (Heymans et al., 2013), but narrow enough to guarantee that the 3point signal is detectable at a few sigma (Vafaei et al., 2010).
4.1 Mass reconstruction
Mass reconstruction on the CFHTLenS data follows the same procedure as described in Section 2 using the KS93 algorithm (Kaiser & Squires, 1993). The regular grid over which the mass reconstruction and smoothing are performed consists of square pixels with an area of approximately arcmin. The average of the galaxy ellipticities is first calculated within each pixel, using the calibration correction implemented in Miller et al. (2013) in order to account for the constant correction (Heymans et al., 2012b), the multiplicative shear measurement bias , and the weighting for each galaxy located at within that pixel. The shear estimate per galaxy is not divided by , but instead the average in each pixel is given by the weighted sum over the galaxies in that pixel:
(29) 
The smoothing of this gridaveraged ellipticity map is then performed using the Gaussian window function given by Eq.(9). Mass maps are reconstructed with the following smoothing scales: , , , , and arcmin. The final maps have a size of , , , and for W1, W2, W3, and W4 respectively, which is chosen to be close to the grid size that was used on the simulations. The fact that the angular resolution is different for the different fields has no effect on our results. Mass reconstruction is also performed with the galaxies rotated by degrees in order to probe the modes. As a sanity check, we computed the crosscorrelation map between and , following Eq.(24), which shows that the two maps are uncorrelated, hence supporting the conclusion that the mode in the data is consistent with zero.
4.2 Cosmic statistics on CFHTLenS mass maps
Following the procedure described in Section 3.2, independent pure noise maps are reconstructed for each CFHTLenS field and each smoothing scale. The pure noise maps are necessary for the denoising of the convergence moments. This paper is the first to measure cosmic shear statistics directly from reconstructed mass maps and it is therefore important to quantify the reliability of the measurements against predictions. Theoretical predictions for the second and third order moments are calculated using Eqs.(15) and (20). The key ingredient for these predictions is an accurate redshift distribution of the galaxies. In order to estimate , we use the probability distribution function of the photometric redshifts (Hildebrandt et al., 2012), the robustness of which has been thoroughly tested in Benjamin et al. (2012). The implementation of this approach in other CFHTLenS papers, where the second order moment of the shear has been used to constrain cosmological parameters, shows that the redshift distribution derived from CFHTLenS photometric redshifts provides a robust and consistent interpretation of the weak lensing data as a function of source redshift (Benjamin et al., 2012; Heymans et al., 2012a; Kilbinger et al., 2013; Simpson et al., 2013).
The data points in Figure 5 show the stacked redshift distribution PDF for our galaxy selection with photometric redshifts , the same galaxy selection used for the lensing mass reconstruction. In order to make predictions for the lensing statistics, we fit the observed redshift distribution with the following fourparameter doubleGaussian model:
(30) 
Figure 5 shows that this is a reasonable fit to the data, in particular in capturing the long, slowlydecreasing tail at high redshifts. The best fit values are , , , and .
Figure 6 shows the second to fifth order moments of the convergence measured on CFHTLenS lensing data. The solid lines show the predictions for the second and third order moments using the redshift distribution derived above. The predictions are not a fit to the data. They are computed using the CAMB nonlinear power spectrum
4.3 The connection between large scale dark matter and baryons
Construction of the mass map predicted from baryons
Sections 3 and 4 describe the mass reconstruction process in the presence of realistic noise and masks. A cosmological signal is clearly measured on the CFHTLenS data and the level of residual systematics is shown to be consistent with zero. In this section, we are interested in the comparison between the reconstructed dark matter and the matter distribution of the stellar content. The connection between dark matter and baryons is quantified in the CFHTLenS companion papers (Gillis et al., 2013; Velander et al., 2013) at the galaxy and galaxygroup scales. Here we are interested in the connection between dark matter and baryons at much larger scales. For example, we would like to explore whether there are large scale features common to both the dark matter and the baryons, such as voids. To this end we develop a new approach, which consists of first predicting the dark matter map from the stellar mass distribution (Wilson et al., 2001), and then comparing the peak distribution from the predicted map to the lensing map.
In order to construct a predicted map, a NavarroFrenkWhite dark matter density profile (Navarro et al., 1997) is assigned to each galaxy in the CFHTLenS catalogue; the procedure is detailed below. There is no need to separate the galaxy population into sources and lenses because most galaxies act as both due to our broad redshift selection. Galaxies are assigned a dark matter mass following the relation between stellar mass and halo mass from Leauthaud et al. (2012), which provides a relation for redshifts between and , covering almost the entire redshift range of galaxies in our work. The only complication is that this relation is only provided for central galaxies, while in our sample a large fraction are satellites. Unfortunately, it is impossible to separate the central galaxies from the satellites on an individual basis; it can only be done statistically (Velander et al., 2013). We therefore proceed by assuming that all galaxies are central galaxies and all follow the stellar mass to halo mass relation in Leauthaud et al. (2012). This will overestimate our predicted total mass, but to a first approximation it should not dramatically affect the relative distribution of mass. Based on Figure 4 from Leauthaud et al. (2011), we anticipate that, on average, there are roughly one to two satellite galaxies for every central galaxy. This should lead us to overestimate the total mass by roughly a factor of to . The exact calculation is not needed as we are only interested in an order of magnitude estimate of how wrong our predicted convergence can be. In a future work, the same strategy will be applied to clusters instead of individual galaxies, which should mitigate this effect. In order to complete our convergence prediction from the galaxies, we need to assign a concentration to each halo. To this end, we use the massconcentration relation calibrated from numerical simulations in MuñozCuartas et al. (2011).
At this stage, each galaxy in the CFHTLenS catalogue is associated with a dark matter halo of known concentration and mass. The last step is to apply the lensing kernel in order to predict the convergence based on the redshifts of the lenses and sources. Every galaxy is simultaneously both a lens and a source, depending on whether it is in the background or the foreground relative to other galaxies. One can then compute for each galaxy a predicted convergence based on the foreground mass distribution coming from all galaxies located at lower redshift. For a source galaxy at location on the sky with redshift and foreground lenses at redshifts the total convergence predicted from the baryonic distribution is given by:
(31) 
where is the projected halo mass of lens centred at , and is the critical density given by:
(32) 
Note that the average predicted convergence is calculated only after all haloes have been assigned to the galaxies. The critical density depends on the observerlens, lenssource, and observersource angular diameter distances , , and respectively. The skyaverage predicted convergence is set to zero by subtracting the mean in Eq.(31). We assume that the weak lensing approximation applies, which means that the convergence from the different lenses can be added linearly. It is important to emphasize that the lens redshift goes down to the lowest reliable value , and that the sources only cover redshifts and higher (in order to be consistent with our source galaxy selection outlined in Section 4.2).
The convergence predicted from the baryonic content is assigned to each galaxy in the CFHTLenS catalogue. Following Section 4.2, the are placed on the same regular grid that is used for the lensing mass reconstruction. The lensfit weighting is used to determine the average within each pixel, and the same Gaussian smoothing is applied. For most galaxies, the size of a single galactic halo is comparable to the size of a pixel, which is of the order of half an arcminute, in the final resolution map. The galaxies in the masked regions therefore do not impact on our analysis as the extension of their halo is small
Comparing the lensing with the predicted dark matter maps
In this Section, we perform a comparison of the and maps. Following an approach similar to Section 3.3, when comparing noisefree and noisy simulated data, we first compute the crosscorrelation profile between the lensing reconstructed mass map and the predicted map and between the rotated galaxies map and . These profiles are called and respectively, as defined by Eq.(24). Figure 7 shows that the correlation increases for larger smoothing windows, and that the overall correlation level always remains below . The relatively low crosscorrelation coefficient for small smoothing scales is due to the fact that the noise level in the reconstructed lensing mass maps is high, and for small smoothing windows, many of the peaks and structures we see in mass maps are the result of noise, which is consistent with the description given in Section 3.3. However, the predicted maps , based on real galaxies that have been detected by CFHTLenS, have a lower noise level than the lensing maps. The sources of noise in the predicted maps are likely dominated by intrinsic stochastic biasing between dark matter and baryons and the undetected galaxies in the survey. For this reason, we adopt the predicted maps as the reference from which the peaks will be detected and then compared to the lensing maps. Using the lensing maps as the reference is formally equivalent, but we found that the level of noise in the comparison is reduced when the predicted maps are used instead.
Next, we want to compare the 2D spatial distribution of peaks between the maps. The peak distribution is a powerful tool that helps visually identify the large scales structures. We will see that this comparison reveals the existence of large underdensities (voids) that cannot be identified with a statistical analysis using moments. Given that for a fixed smoothing scale, the noise in lensing maps is higher than the noise the predicted maps, we decided to detect peaks in the predicted map using the arcmin smoothing scale and compare it to the lensing map using the smoothing scale of arcmin. A peak location is defined as a pixel where all surrounding pixels have a lower amplitude. Figure 8 shows the location of peaks for W1 (shown as white circles) superimposed on the reconstructed lensing map shown as the continuous coloured background. Contours are shown for the lensing reconstruction map at 1,2,3, and 4 sigma levels, which is a common way of indicating the significance of structures in lensing maps. On average, the distribution of peaks matches the lensing mass overdensities. A quantitative comparison between the predicted convergence and the lensing convergence is shown in the left panel of Figure 9; the small dots in this figure show, for each peak detected on the map, the corresponding value of the lensing map at the same location. Note that for Figure 9 we have used a smoothing scale of arcmin for the lensing map as well, hence the high noise rms for the lensing peak amplitude. Although only peaks have been used in this plot, it does not mean that the lensing convergence is positive. In fact, there are a substantial number of peaks inside low density regions with negative convergence, despite the fact that, on average, peaks live in positive convergence regions. A comparison of the binned peak values to the averaged lensing convergence in the bin shows that the reconstructed and predicted convergence are strongly correlated, consistent with the fact that baryons trace dark matter to first approximation. The lensing convergence is roughly times lower than the predicted convergence; this is expected, as discussed in Section 4.4, because we assumed that all galaxies are central galaxies of the host dark matter halo. Following the numerical calculations from Leauthaud et al. (2011), this erroneous assumption would lead to an overestimate of the predicted mass by a factor of . Note that for an individual peak, the noise has rms , consistent with the fact that there are approximately lensed galaxies within each smoothing window. Figure 9 also shows the peak analysis performed with the mode lensing map , and one can see that the correlation between and peaks and troughs vanishes.
The analysis of the 2D distribution of peaks provides a consistent picture where the galaxy distribution is correlated with the lensing mass, as one would expect. We can perform a similar analysis on troughs (local minima). A trough location is defined as a pixel where all surrounding pixels have a higher amplitude. A high number density of troughs would be an indication of cosmic voids. Figure 10 shows such an analysis for W1, where white triangles represent troughs in . One can see that the triangles preferentially populate underdense regions of the lensing map extending over a few degrees; this is a convincing illustration of the detection of giant projected voids in large scale mass maps based on lensing data. A comparison of Figures 8 and 10 reveals how the distribution of peaks and troughs provides a mapping of large scale structures and voids, respectively. The right panel on Figure 9 shows how troughs compare to the lensing mass reconstruction , with troughs populating large void regions preferentially. This shows that the detection of voids is significant, with residual systematics also being consistent with zero in the underdense regions. Troughs are located primarily in negative convergence regions, as one would expect. The predicted convergence is closer to the reconstructed convergence for the troughs than for the peaks (Figure 9). A possible explanation is that there are less satellite galaxies in voids than in overdense regions.
Figures 11 and 12 show the same analysis for the three other CFHTLenS fields W2, W3, and W4. The peak and void statistics shown in Figure 9 include all four fields. The key message from the comparison of the baryon and dark matter maps is that dark matter maps trace large scale structures reliably, but evaluating the significance of most of the individual mass peaks is challenging for two reasons; first, the shot noise is high and very few peaks are above the secure threshold. Second, the maps clearly show that the large scale structures are a source of noise comparable to the shot noise (Hoekstra, 2001; Hoekstra et al., 2011).
5 Conclusion
This paper is the first quantitative cosmological analysis of mass maps reconstructed from lensed galaxy shapes. We validated our approach by using Nbody simulations and then applied it to the CFHTLenS data. We find that convergence maps contain reliable cosmological information that has the potential to go beyond the traditional analysis using 2point statistics. Peak statistics and morphological analysis are the next studies to be performed on mass maps.
Using Nbody simulations, we have shown that the reconstruction process is stable and that the input cosmological signal is recovered accurately despite the presence of masks, the relatively high level of shape noise, and the nonPoissonian spatial distribution of the background sources. Mass reconstruction was performed and tested with the traditional KS93 algorithm (Kaiser & Squires, 1993). The testing of the full nonlinear mass reconstruction at the same level of precision is left for a future study.
The application to the CFHTLenS data shows for the first time that windowed statistics of the convergence can be measured on mass maps. We find excellent agreement between the second and third order moments measured on the reconstructed mass maps and the predictions for a cosmological model determined from a shear correlation function analysis of the same data (Benjamin et al., 2012; Heymans et al., 2012a; Kilbinger et al., 2013; Simpson et al., 2013). Our attempt at measuring higher order statistics shows a marginal detection of the fourth order moment. For all the moments of the convergence, the residual systematics are found to be consistent with zero. In Section 4.4, we compared the reconstructed convergence with the predicted convergence using galaxies as tracers of dark matter haloes, where a halo is assigned to each galaxy. We have shown that the predicted and reconstructed mass maps are strongly correlated with each other. The maps reveal the existence of large voids in the projected dark matter distribution, which span regions as large as degrees on the sky.
We consider this paper to be a feasibility study that strongly suggests that future precision cosmology on mass maps is possible. Cosmology with mass maps enables studies that are currently not possible with only the shear (shape) information: global structure morphology, peak statistics with long wavelength modes included, cross correlation with other cosmology maps, and statistics of convergence probability distribution functions. However, before mass maps can become a completely reliable cosmological probe, several technical issues must be addressed:

A correct nonlinear reconstruction will have to be implemented to account for the most massive structures, such as clusters of galaxies.

Our approach consisted of removing the noise bias from the moments measured on the reconstructed mass. A completely different strategy is to denoise the map itself before any measurement (Starck et al., 2006). It remains to be tested whether the latter is a robust and better approach for precision measurements.

The photometric redshift uncertainty will have to be included in the comparison of dark matter and baryons and in the map making process itself. How to construct a predicted mass map from a distribution of galaxies with photometric redshift errors is an open problem.

A potential fundamental limitation with maps is that the current mass reconstruction process does not distinguish between intrinsic galaxy alignment and gravitational lensing, which both lead to correlated galaxy shapes. From our particular choice of redshift range, and from the residual systematics studies in CFHTLenS companion papers, we know that our measured signal is mainly caused by gravitational lensing (Heymans et al., 2013; Kitching et al., 2013), but future precision cosmology using mass maps will require a clear identification and separation of the different causes leading to galaxy shapes correlations. It is not clear whether this is possible with mass maps, and further theoretical studies will be needed.
Acknowledgments
We would like to thank P. Schneider and B. Gillis for insightful comments on the manuscript, and Jennifer Flood for a careful editorial work on the manuscript. This work is based on observations obtained with MegaPrime/MegaCam, a joint project of the CanadaFranceHawaii Telescope (CFHT) and CEA/Irfu, at CFHT, which is operated by the National Research Council (NRC) of Canada, the Institut National des Sciences de l’Univers (INSU) at the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. This research used the facilities of the Canadian Astronomy Data Centre operated by the NRC of Canada with the support of the Canadian Space Agency. We thank the CFHT staff, in particular J.C. Cuillandre and E. Magnier, for the observations, data processing and continuous improvement of the instrument calibration. We also thank TERAPIX for quality assessment, and E. Bertin for developing some of the software used in this study. CFHTLenS data processing was made possible thanks to support from the Natural Sciences and Engineering Research Council of Canada (NSERC) and HPC specialist O. Toader. The Nbody simulations were performed on the TCS supercomputer at the SciNet HPC Consortium. The early stages of the CFHTLenS project were made possible thanks to the European Commissions Marie Curie Research Training Network DUEL (MRTNCT2006036133) and its support of CFHTLenS team members LF, HHi, and BR.
The Nbody simulations used in this analysis were performed on the TCS supercomputer at the SciNet HPC Consortium. SciNet is funded by: the Canada Foundation for Innovation under the auspices of Compute Canada; the Government of Ontario; Ontario Research Fund  Research Excellence; and the University of Toronto.
LF acknowledges support from NSFC grants 11103012 & 10878003, Innovation Program 12ZZ134 and Chen Guang project 10CG46 of SMEC, and STCSM grant 11290706600 & Pujiang Program 12PJ1406700. CH and FS acknowledge support from the European Research Council (ERC) through grant 240185. TE is supported by the DFG through project ER 327/31 and the Transregional Collaborative Research Centre TR 33. HHo acknowledges support from Marie Curie IRG grant 230924, the Netherlands Organisation for Scientific Research (NWO) through grant 639.042.814 and from the ERC through grant 279396. HHi is supported by the Marie Curie IOF 252760 and by a CITA National Fellowship. TDK is supported by a Royal Society University Research Fellowship. YM acknowledges support from CNRS/INSU and the Programme National Galaxies et Cosmologie (PNCG). LVW and MJH acknowledge support from NSERC. LVW also acknowledges support from the Canadian Institute for Advanced Research (CIfAR, Cosmology and Gravity program). BR acknowledges support from the ERC through grant 24067, and the Jet Propulsion Laboratory, California Institute of Technology (NASA). TS acknowledges support from NSF through grant AST0444059001, SAO through grant GO011147A, and NWO. ES acknowledges support from the NWO grant 639.042.814 and support from ERC under grant 279396. MV acknowledges support from NWO and from the Beecroft Institute for Particle Astrophysics and Cosmology.
Author Contributions: All authors contributed to the development and writing of this paper. The authorship list reflects the lead author of this paper (LVW) followed by two alphabetical groups. The first alphabetical group includes key contributors to the science analysis and interpretation in this paper, the founding core team and those whose longterm significant effort produced the final CFHTLenS data product. The second group covers members of the CFHTLenS team who made a significant contribution to either the project, this paper, or both. The CFHTLenS collaboration was coled by CH and LVW.
6 Appendix
In this Appendix, we derive the relation between the denoised moments , and the observed moments , , where is the Gaussian smoothing scale. The observed convergence map is the sum of the true signal and an uncorrelated noise component . We are repeating Eqs.(10) and (11):
(33) 
The observed mode convergence map , which is the reconstructed mass from the degree rotated galaxies, is similarly related to the true mode and the noise contribution :
(34) 
Note that the true mode should be zero if there is no residual type systematics. Since we are testing this hypothesis precisely, we keep as a quantity to be measured from the data instead of making it equal to zero. The noise contribution to each convergence moment measured in this work are then given by:
(35)  
(36)  
(37)  
(38)  
(39)  
(40)  
(41)  
(42)  
(43)  
(44)  
(45)  
(46) 
The lefthand side is the denoised (unbiased) moment and the righthand side shows the combination of moments leading to this unbiased measurement. The measured signal, mode, and random noise maps are given by , , and respectively. As a reminder from the definition given in Section 4.2, a few ensemble averages have been introduced: denotes the moment average over one map, is averaged over several noise maps , and is the denoised moment. One can see from these equations that there are as many systematics moments as there are permutations of and for a given moment order. A total of one hundred random maps have been generated for each smoothing scale. In order to derive Eqs.(46), we assumed that all the terms with odd order moments of the noise maps can be neglected, as they only introduce noise in the equation above.
Footnotes
 pagerange: CFHTLenS: Mapping the Large Scale Structure with Gravitational Lensing–References
 pubyear: 2012
 The source redshifts are .
 http://camb.info/
 Note that the haloes from low redshift galaxies do extend over several pixels, but their lensing efficiency is small due to their proximity to the observer
References
 Athreya R. M., Mellier Y., Van Waerbeke L., Pelló R., Fort B., DantelFort M., 2002, A&A, 384, 743
 Bartelmann M., Narayan R., Seitz S., Schneider P., 1996, ApJ, 464, L115
 Benjamin J. et al., 2012, ArXiv1212.3327, accepted
 Bernardeau F., Van Waerbeke L., Mellier Y., 1997, A&A, 322, 1
 Bond J. R., Efstathiou G., 1987, MNRAS, 226, 655
 Clowe D., Bradač M., Gonzalez A. H., Markevitch M., Randall S. W., Jones C., Zaritsky D., 2006a, ApJ, 648, L109
 Clowe D. et al., 2006b, A&A, 451, 395
 Crittenden R. G., Natarajan P., Pen U. L., Theuns T., 2002, ApJ, 568, 20
 Erben T. et al., 2012, ArXiv eprints: 1210.8156
 Erben T., Van Waerbeke L., Mellier Y., Schneider P., Cuillandre J. C., Castander F. J., DantelFort M., 2000, A&A, 355, 23
 Gillis B. R. et al., 2013, ArXiv1301.7421, accepted
 HarnoisDéraps J., Vafaei S., Van Waerbeke L., 2012, MNRAS, 426, 1262
 Heymans C. et al., 2008, MNRAS, 385, 1431
 Heymans C. et al., 2012a, MNRAS submitted
 Heymans C. et al., 2012b, MNRAS, 427, 146
 Heymans et al., 2013, MNRAS, submitted
 Hildebrandt H. et al., 2012, MNRAS, 421, 2355
 Hoekstra H., 2001, A&A, 370, 743
 Hoekstra H., Franx M., Kuijken K., 2000, ApJ, 532, 88
 Hoekstra H., Hartlap J., Hilbert S., van Uitert E., 2011, MNRAS, 412, 2095
 Hoekstra H., Jain B., 2008, Annual Review of Nuclear and Particle Science, 58, 99
 Jain B., Van Waerbeke L., 2000, ApJ, 530, L1
 Jee M. J., Mahdavi A., Hoekstra H., Babul A., Dalcanton J. J., Carroll P., Capak P., 2012, ApJ, 747, 96
 Juszkiewicz R., Weinberg D. H., Amsterdamski P., Chodorowski M., Bouchet F., 1995, ApJ, 442, 39
 Kaiser N., 1992, ApJ, 388, 272
 Kaiser N., 1998, ApJ, 498, 26
 Kaiser N., Squires G., 1993, ApJ, 404, 441
 Kayser R., Schramm T., 1988, A&A, 191, 39
 Kilbinger M. et al., 2013, MNRAS, 735
 Kilbinger M., Schneider P., 2005, A&A, 442, 69
 Kitching et al., 2013, MNRAS, submitted
 Komatsu E., et al., 2011, ApJS, 192, 18
 Kratochvil J. M., Haiman Z., May M., 2010, Phys. Rev. D, 81, 043519
 Leauthaud A., Tinker J., Behroozi P. S., Busha M. T., Wechsler R. H., 2011, ApJ, 738, 45
 Leauthaud A. et al., 2012, ApJ, 744, 159
 Leonard A., Dupé F.X., Starck J.L., 2012, A&A, 539, A85
 Mahdavi A., Hoekstra H., Babul A., Balam D. D., Capak P. L., 2007, ApJ, 668, 806
 Marian L., Smith R. E., Hilbert S., Schneider P., 2012, MNRAS, 423, 1711
 Massey R. et al., 2007, Nature, 445, 286
 Maturi M., Angrick C., Pace F., Bartelmann M., 2010, A&A, 519, A23
 Mecke K. R., Buchert T., Wagner H., 1994, A&A, 288, 697
 Miller L. et al., 2013, MNRAS, 429, 2858
 MiraldaEscude J., 1991, ApJ, 380, 1
 MuñozCuartas J. C., Macciò A. V., Gottlöber S., Dutton A. A., 2011, MNRAS, 411, 584
 Munshi D., Valageas P., Barber A. J., 2004, MNRAS, 350, 77
 Munshi D., Valageas P., Van Waerbeke L., Heavens A., 2008, Phys. Rep., 462, 67
 Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493
 Sato J., Takada M., Jing Y. P., Futamase T., 2001, ApJ, 551, L5
 Schneider P., 1996, MNRAS, 283, 837
 Schneider P., Van Waerbeke L., Jain B., Kruse G., 1998, MNRAS, 296, 873
 Scoccimarro R., Couchman H. M. P., 2001, MNRAS, 325, 1312
 Simon P. et al., 2012, MNRAS, 419, 998
 Simon P., Taylor A. N., Hartlap J., 2009, MNRAS, 399, 48
 Simpson F. et al., 2013, MNRAS, 429, 2249
 Starck J. L., Pires S., Réfrégier A., 2006, A&A, 451, 1139
 Takada M., Jain B., 2003, MNRAS, 344, 857
 Takada M., Jain B., 2004, MNRAS, 348, 897
 Vafaei S., Lu T., Van Waerbeke L., Semboloni E., Heymans C., Pen U. L., 2010, Astroparticle Physics, 32, 340
 Valageas P., Munshi D., 2004, MNRAS, 354, 1146
 Valageas P., Munshi D., Barber A. J., 2005, MNRAS, 356, 386
 Van Waerbeke L., 2000, MNRAS, 313, 524
 Van Waerbeke L., Bernardeau F., Mellier Y., 1999, A&A, 342, 15
 Van Waerbeke L., Hamana T., Scoccimarro R., Colombi S., Bernardeau F., 2001, MNRAS, 322, 918
 VanderPlas J. T., Connolly A. J., Jain B., Jarvis M., 2011, ApJ, 727, 118
 Velander et al., 2013, MNRAS, submitted
 Wilson G., Kaiser N., Luppino G. A., 2001, ApJ, 556, 601
 Yang X., Kratochvil J. M., Wang S., Lim E. A., Haiman Z., May M., 2011, Phys. Rev. D, 84, 043529