Probing Cosmology with Weak Lensing Peak Counts
Abstract
We propose counting peaks in weak lensing (WL) maps, as a function of their height, to probe models of dark energy and to constrain cosmological parameters. Because peaks can be identified in two–dimensional WL maps directly, they can provide constraints that are free from potential selection effects and biases involved in identifying and determining the masses of galaxy clusters. We have run cosmological N–body simulations to produce WL convergence maps in three models with different constant values of the dark energy equation of state parameter, , and , with a fixed normalization of the primordial power spectrum (corresponding to present–day normalizations of , and , respectively). By comparing the number of WL peaks in 8 convergence bins in the range of , in multiple realizations of a single simulated degree field, we show that the first (last) pair of models can be distinguished at the 95% (85%) confidence level. A survey with depth and area (20,000 degree), comparable to those expected from LSST, should have a factor of better parameter sensitivity. We find that relatively low–amplitude peaks (), which typically do not correspond to a single collapsed halo along the line of sight, account for most of this sensitivity. We study a range of smoothing scales and source galaxy redshifts (). With a fixed source galaxy density of 15 arcmin, the best results are provided by the smallest scale we can reliably simulate, 1 arcminute, and provides substantially better sensitivity than .
pacs:
98.80.Cq, 11.25.w, 04.65.+eI Introduction
Weak gravitational lensing (WL) has emerged as one of the most
promising methods to constrain the parameters of both dark energy (DE)
and dark matter (DM) Albrecht et al. (2006). This method exploits the fact that
images of distant galaxies are distorted as their light is lensed by
inhomogeneities in the foreground gravitational potential. DE and DM
alter the expansion history of the universe, as well as the rate at
which cosmic dark matter structures grow from small initial
perturbations, and thereby modify the statistical features of the WL
signal (2). Several deep astronomical surveys plan to
cover large regions ( degree) of the sky in the future,
detecting millions of background galaxies, and measuring their shapes
precisely, in order to reveal the statistical features of the WL
distortion. Examples include the Panoramic Survey Telescope And Rapid
Response System (PanSTARRS) (3), the KiloDegree
Survey (KIDS)
Existing theoretical studies assessing the utility of such surveys to probe the properties of DE and of DM can be roughly divided into two types. The first type analyzes the WL distortion primarily on large scales, caused by fluctuations in the foreground mass distribution on correspondingly large scales. On these scales, the fluctuations, and hence the WL distortion patterns, are linear or quasilinear, and can be readily and reliably interpreted using existing semi–analytic tools (Kaiser, 1992; Jain and Seljak, 1997; Hu, 1999, 2002; Huterer, 2002; Refregier et al., 2004; Abazajian and Dodelson, 2003; Takada and Jain, 2004; Song and Knox, 2004). These studies (for example, those computing the two–dimensional power spectrum of the shear) have, in fact, been pushed into the non–linear regime, but the dependence of the non–linear power–spectrum on cosmological parameters has not been verified in simulations, to an accuracy necessary to take full advantage of the future surveys.
The second type analyzes the distortion on smaller scales, where it is dominated by the nonlinear mass fluctuations. In the limit of extremely non–linear, high– WL shear peaks, the individual peaks tend to be caused by a single, discrete, massive collapsed object – a galaxy cluster – in the foreground (White et al. (2002); (16); (17)). The abundance evolution and spatial distribution of galaxy clusters in different cosmologies can then be studied semi–analytically (calibrated by simulations), and has been shown to place very strong constraints on DE parameters, because it is exponentially sensitive to the growth rate of dark matter fluctuations. Several studies have explored constraints expected from samples of tens of thousands of galaxy clusters in Sunyaev–Zeldovich and X–ray surveys (e.g. (Wang and Steinhardt, 1998; Haiman et al., 2001; Weller et al., 2002)), and from the still larger samples of clusters expected to be detectable in large weak lensing surveys ((21); (22); (23)).
The amount of statistical information in these two regimes appears to be comparable ((24)); furthermore, the information from the two approaches is complementary, and using them in combination significantly improves DE constraints compared to using either method alone ((24); (25)).
The chief difficulty with using galaxy clusters identified in WL maps, however, is that the correspondence between peaks in WL maps and collapsed clusters is not one–to–one. A significant fraction of WL peaks are caused by chance alignments of multiple non–linear structures along the line–of–sight. Conversely, a non–negligible fraction of real clusters are missed because they do not show a significant WL peak, due to a chance cancellation of their shear signal by underdense regions in the foreground. In principle, the statistical correspondence between WL peaks and clusters can be precisely quantified ab–initio in numerical simulations, but in practice, this remains a challenge to the required level of accuracy. In particular, the completeness and purity of cluster selection are of the order of for peaks selected to lie above a convergence threshold, and both become progressively worse for lower S/N peaks (White et al. (2002); (21); Marian and Bernstein (2006); (27); (28)). Systematic selection errors are therefore likely to dominate over the statistical constraints for the majority of a WL–selected cluster sample.
A natural way to circumvent such selection effects is to study the statistics of non–linear peaks in WL maps directly. This method, however, does not lend itself to straightforward mathematical analysis – it requires numerical simulations and has been relatively much less well–explored (e.g. (Jain and Van Waerbeke (2000); (30))). In this paper, we study WL peaks in simulated maps, and their usefulness to constrain DE properties. More specifically, we use a set of Nbody simulations, combined with raytracing, to create convergence maps for three different values of the DE equation of state, . We then identify peaks in such maps, and tabulate them as a function of their height. We explicitly avoid a reference to clusters, cluster masses, or mass functions – we utilize only the directly observable height of the maximum of each convergence peak. This method should therefore provide one of the most robust cosmological probes, free from potential selection effects and mass measurements biases affecting cluster identifications.
For increasingly negative values of , DE starts to dominate at
increasingly later cosmic epochs, and one expects an excess of high
peaks in such models, compared to those with less negative . In
this first paper, we perform a feasibility study, to asses whether
changes in indeed have a measurable effect on the number of
peaks. To this end, all cosmological parameters, other than , are
held fixed at their fiducial values.
While the cosmological dependence of the cluster counts has been studied extensively, it is not a priori clear whether the full catalog of all WL peaks contains a comparable amount of information. A large number of peaks caused by largescale structure, with no –dependence, would constitute pure noise, and could wash out the theoretically known difference from peaks due to clusters. A strong dependence in the number of non–cluster peaks, in the sense opposite to that of the cluster peaks, could be even more detrimental, since this would explicitly cancel some of the signal contained in the cluster counts. Conversely, a strong dependence in the number of non–cluster peaks, in the same sense as the cluster counts, could significantly enhance the sensitivity. This remains a realistic possibility, since the dependence of the number of peaks caused by line of sight projections on cosmology has not yet been studied.
An intuitive expectation is that largescale structure filaments are expected to contribute more strongly to relatively low–height peaks, and will be more prominent for source galaxies at higher redshift. This is because most filaments are likely to line up only partially, thus creating a weaker signal than those caused by massive collapsed clusters, and because the probability to find large (Mpc) filaments that are lined up along the line of sight will increase significantly at larger distances. We perform our simulations at high mass–resolution (down to less than per particle), allowing us to resolve relatively low significance peaks in our convergence maps, where such projections of large scale structure are more common,. We also place the source galaxies out to redshifts as high as .
To keep the computational cost of the simulations reasonable, we sacrifice the angular size of our simulated field of view, limiting ourselves to creating WL maps covering degrees (for simplicity, we will hereafter refer to these as degree fields).
In these respects, our work differs from the recent related study by (22), which considered the projected mass function out to a lower redshift (to ), with simulations at significantly lower resolution, to demonstrate that the 2D projected mass function obeys a scaling with and that is similar to that of the three–dimensional halo mass function ((22) also do not vary ). Our work shows that the depth of the survey (at least to ) is vital for distinguishing between models with different , and that the main distinguishing power comes from medium–height convergence peaks, with convergence (on maps with 1 arcmin smoothing).
Another recent study, (30) explored probing cosmology with a related, simple statistic: the high tail of the one–point probability distribution of the convergence. The total area above a fixed convergence threshold includes a sum of the area of all peaks lying above this threshold, and we therefore expect it to capture less information than peak counts. Nevertheless, (30) have found that this statistic can already place constraints on DE properties, comparable to the traditional linear two point correlation function of shear or to the limited use of the nonlinear information from counting galaxy clusters ((30)). Our paper improves on this preliminary study, by employing direct simulations, rather than fitting functions, for the convergence and its dependence on cosmology, and also by examining a potentially more powerful statistic.
As our paper was being prepared for submission, we became aware of an independent preprint, addressing the dependence of WL peaks on the background cosmology (31). Their motivation is essentially identical to our own, and the conclusions reached are quite similar. However, the two studies are complementary, in that (31) study the peak dependence in the plane assuming a CDM cosmology, while we vary the dark energy equation of state parameter (with corresponding changes in ). There are additional differences, such as the specifications in the simulations that were used, and the statistical analysis methods, which we briefly discuss further in § IV.6 below.
The rest of this paper is organized as follows. In § II, we describe our calculational procedure, including the set–up and analysis of the simulations. In § III, we test our computations by comparing our numerical convergence power spectrum to semianalytic predictions. In § IV, we present and discuss our main results on the peak counts. Finally, in § V, we summarize our main conclusions and the implications of this work.
Ii Methodology
This section describes the procedure to create the simulated weak lensing maps, from generating the matter power spectrum with CAMB and scaling it to the starting redshift of the Nbody simulation, to generating the initial conditions for the simulation, raytracing the weak lensing maps from the simulation output, as well as post–processing and analyzing these maps.
ii.1 Nbody Simulations
For this study, we generated a series of cold dark matter (CDM) Nbody simulations with the code GADGET2 (32), which include DM only; baryons and neutrinos are both assumed to be absent. We have modified GADGET2 for use with arbitrary equation of state parameter , but in this work, we restrict ourselves to constant . As the reference model, we chose a CDM model with parameters close to the bestfit values from WMAP5, i.e. with cosmological constant , matter density parameter , Hubble constant . We fix the amplitude of curvature fluctuations at the pivot scale . This results in different values of the present–day normalization, , for the different dark energy models. In the three models we consider with , and , the values are , and , respectively.
For each of the three cosmologies, with all parameters other than and fixed, we ran two independent Nbody simulations, with two different random realizations of the initial particle positions and velocities. The simulations used particles, in a box with a comoving size of , starting at redshift . This corresponds to a mass resolution of per dark matter particle. The Plummerequivalent gravitational softening of the simulations was set to (comoving). Fixing the softening length in physical (rather than comoving) units is not warranted, since for a pure CDM simulation without gas, no further physical scale enters.
The initial conditions (ICs) for the simulations were generated with an IC generator kindly provided by Volker Springel. The matter power spectrum at , which is an input for the IC generator, was created with the EinsteinBoltzmann code CAMB (33), a modification/extension of CMBFAST (34). GADGET2 does not follow the evolution of the radiation density, which, starting with the correct power spectrum at , would accumulate to a error in its normalization by . To eliminate this small bias, we used CAMB to produce the power spectrum at and scaled it back to , with our own growth factor code, with the radiation component absent. This leads to a small deviation from the true power spectrum at the starting redshift, but the match will become better at the lower redshifts which are more important for our analysis.
The output of the simulations consists of particle positions at various redshifts between and . Similar to choices in most previous works (e.g. (37)), the output redshifts were chosen to span comoving intervals of along the line of sight in the CDM cosmology. For simplicity, the same redshifts were used in the other cosmologies. This results in a small difference in comoving separation between the planes in these cosmologies, but this has negligible impact on our results.
Most simulations and all post–processing analysis were run on the LSST Cluster at the Brookhaven National Laboratory, an 80CPU Linux cluster with 80 GB of memory. Some simulations were run on the Intel 64 Cluster Abe, with 9600 CPUs and 14.4 TB of memory at the National Center for Supercomputing Applications (NCSA), a part of the NSF TeraGrid facility.
ii.2 Weak Lensing Maps
We use our own codes to produce and analyze weak lensing maps from the Nbody simulation output. This subsection describes the procedures and algorithms, as well as the testing of the new codes.
Potential Planes
As mentioned above, the simulation outputs are produced at every Mpc, which results in an overlap between consecutive Mpc cubes. We prepare the Nbody simulation output for raytracing by truncating the snapshot cubes along the line of sight in order to remove this overlap, and then project the particle density of each cube on a 2D plane perpendicular to the line of sight at the center of the plane, located at the same redshift as the original cube. We use the triangular shaped cloud (TSC) scheme (35) to convert the individual particle positions in the simulations to a smoothed density field on the 2D plane. The projected density contrast field on plane , is calculated on a 2D lattice of points (here represents the coordinates in the plane perpendicular to the line of sight; is the radial coordinate, and is the position half–way between the and planes). We then use the 2D Poisson equation
(1) 
to calculate the 2D gravitational potential on plane . This equation is solved easily in Fourier space, in particular since the density contrast planes have periodic boundary conditions. The resulting potential plane is then stored in FITS format (36) for further use by our raytracing code. This multistep approach has proven useful to save memory, and for re–using intermediate steps of the lensing map construction.
RayTracing
To create weak lensing maps from the Nbody simulations, we implement the twodimensional raytracing algorithm described in (37). We follow light rays through a stack of potential planes, which, as mentioned above, are spaced apart in the case of the CDM cosmology. Starting at the observer, we follow the light rays backward to the source galaxies. At each lens plane , we calculate the distortion tensor and the gravitational lensing deflection from the equations:
(2)  
(3) 
where is the identity matrix, and denote the comoving angular diameter distance and the scale factor, evaluated at the comoving coordinate distance , and the so–called optical tidal matrix,
(4) 
contains the second spatial derivatives of the gravitational potential on plane . The derivatives are to be evaluated at the point on plane intersected by the light ray. Starting with , equation (2) can be solved to find for any .
The physical weak lensing quantities, the convergence and the two components of the shear, and , are extracted from the distortion tensor via
(5) 
The observable quantities, the magnification and the distortion , are given by
(6) 
where is the reduced shear. In the weak lensing limit these relations simplify to and .
For simplicity, for the purpose of this paper, we will work directly with the maps of the theoretically predicted convergence . In practice, in the absence of external size or magnification information, neither the shear, nor the convergence is directly observable, but only their combination, the reduced shear . One can, in principle, convert one quantity to the other; in practice, on small scales, where lensing surveys get much of their constraining power, the errors introduced by this conversion will not be negligible, and must be taken into account when extracting cosmological parameters from actual data (38).
ii.3 Source Galaxies, Smoothing and Noise
The point–particles in the N–body simulations represent lumps of matter with a finite spatial distribution. While gravitational softening is included in GADGET2 to prevent unrealistically close encounters, the simulations still follow the assembly of point particles. As a result, several additional types of smoothing are necessary to produce weak lensing maps. The TSC scheme, which is used to place particles on the density contrast planes, generates smooth two–dimensional density contrast distributions, thus avoiding unnaturally large deflections when a light ray would otherwise pass close to an individual particle. In the raytracing step, as light rays pass through arbitrary points on the potential planes, we employ the same TSC kernel to compute the potential at the ray intersection points.
Since the orientation of the disks of the galaxies with respect to the observer is unknown, this intrinsic ellipticity appears as noise in the maps. Following the interpretation of (39) in (24), we take the redshiftdependent r.m.s. of the noise in one component of the shear to be Song and Knox (2004)
(7) 
and perform the calculation analogously to (39) but with a square tophat function representing our pixels. The resulting convergence noise correlation function for our pixelized map where and are the discrete 2D pixel coordinates is:
(8) 
where is the Kronecker symbol, is the surface density of source galaxies, and is the solid angle of a pixel. For simplicity, in our analysis we assume that the source galaxies are located on one of three source planes, at fixed redshifts (), with on each source plane. In combination, this would yield a total surface density of . For comparison, we note that this is larger than the value (with a redshift–distribution that peaks at ), expected in LSST. However, we do not have sufficient simulation data to perform accurate tomography, and therefore, our results below utilize only a single source plane, with a conservative number of , rather than the full (see below). Thus, we add noise to the convergence in each pixel, drawn from a Gaussian distribution with the variance given in equation (8). Note that the noise between different pixels is uncorrelated.
In order to investigate the information content in peaks defined on different scales, once noise in each pixel is added, we apply an additional level of smoothing, with a finite 2D Gaussian kernel, on scales between arcminutes (see Sec. II.4 for exact values). The Gaussian kernel has the nice property of being decomposable, thus in practice, we can equivalently make two passes through the map with 1D Gaussian kernels, oriented once in each of the two orthogonal dimensions of the map. This saves computation time, especially for large smoothing scales.
The noise in the convergence introduced by intrinsic ellipticity, after the Gaussian smoothing, becomes
(9) 
which is the equivalent of (8), but with a Gaussian smoothing instead of a tophat function (39), and depends on the redshift distribution of the source galaxies, as well as on the smoothing scale applied to the maps. In general, the noise level should be computed as the weighted average of the redshift–dependent shear noise (7) over the redshift distribution of the source galaxies:
(10) 
Since in our case, for simplicity, we assume that the galaxies are located on three discrete source planes, is simply a delta function, at , or at (or a sum of delta functions, for the purpose of redshift tomography).
ii.4 Statistical Comparisons of Convergence Maps
The main goal of this paper is to discern between different cosmologies by counting the convergence peaks. Since the maps have been smoothed, as mentioned above, it is sufficient to identify peaks in the maps as pixels which correspond to local maxima, in the sense that all eight of their neighboring pixels have a lower value. For illustration, we show an example of the convergence maps generated in two of our cosmologies, before and after the addition of noise and smoothing, in Figure 1. The panels with exhibit visibly more structure at high convergence, consistent with the expectation that DE starts to dominate later in this model, allowing more time for growth.
Once the peaks are identified, their heights are recorded in a histogram. This histogram is divided into a reasonable number of bins (the choice in the number of bins will be discussed below), and a vector of values is formed from these bins. We will also attempt to utilize the source galaxy redshifts for tomography, and to combine different smoothing scales, so that each peak is further associated with a value of and . Finally, peak identification is repeated in each realization of a given cosmology. We therefore denote the number of peaks in bin with , where and indicate cosmology and realization, respectively, indicates the redshift of the sources and is the smoothing scale applied to the final map. To simplify notation we hereafter drop the last two lower indices and merge them into the single index .
In order to assess the statistical significance of any difference in peak counts, we construct the covariance matrix for cosmology of the counts in the different height, redshift–, and smoothing–scale bins and ,
(11) 
where is the average number of peaks in bin in the fiducial model , and the summation in equation (11) is over the total number of realizations constructed for every cosmology.
We use 200 such realizations. However, as noted above, we have run only two strictly independent realizations of each cosmology, i.e., two different simulations with different realizations of the initial conditions. Thus we produce additional pseudoindependent realizations by the standard procedure of randomly rotating, shifting, and slicing the data in the simulation output cubes at each redshift, prior to generating the 2D density planes. Our box–size, Mpc is much larger than the few comoving Mpc correlation length of large–scale structures. As mentioned above, the data cubes are truncated at each output, so that we utilize only a radial slice covering about 1/3rd of its length; the degree field furthermore utilizes only a fraction of this slice. As a result, at intermediate redshifts, where the lensing kernel peaks, each realization utilizes about a twelfth of the volume of the snapshot cube, allowing us to access, in each N–body run, effectively truly independent realizations of each lens plane. At the very highest redshift, the situation becomes worse, with up to a third of the volume of the data cubes used. However, due to the distance–dependence of the lensing kernel, the planes near , with these heavily sampled boxes, do not contribute much to the final lensing signal. The random rotations and shifts will create many more than 12 essentially independent random realizations of the incidental line–ups of multiple small clusters and filaments. On the other hand, the very largest clusters – which may create a substantial lensing peak by themselves – will tend to reappear in multiple realizations, and could lead to an underestimation of the variance. Fortunately, our main results come from medium–height peaks, which typically are not caused by a single cluster, so that any such underestimation at the high–convergence end contributes relatively little to our final distinguishing power.
We use the (inverse of) the covariance matrix to compute a value of for each realization of a “test” cosmology, labeled by , against the expectations in the fiducial cosmology, labeled by ,
(12)  
(13) 
where is the difference between the number of peaks in bin in realization of the test cosmology, and the average number of peaks in the same bin in the fiducial cosmology.
The resulting probability distributions for , using 200 realizations of each cosmology, fixed source redshift and smoothing scale , and 8 peak height bins, is shown in Figure 2. In the top left panel, the realizations were drawn from the simulation of the fiducial model itself, which was used to compute the mean number of peaks . Our simulated distribution is very close to a genuine distribution with degrees of freedom (shown by the black curve). In the bottom left panel, the realizations were again drawn from the cosmological model, but from the simulation with an independent realization of the initial conditions. The good agreement with the top left panel tests the accuracy of our covariance matrix (in particular, if the covariance had been underestimated, then this second realization would have produced larger values; see Appendix VI). The two right panels show distributions computed in the realizations of different cosmologies, with (top right panel), and (bottom right panel). The and CDM cosmologies are distinguished at the 68% (95%) confidence level in 84% (39%) of the realizations; the cosmology is distinguishable from CDM at the same confidence in 61% (24%) of the realizations.
For reference, the above confidence levels can be compared to a simpler estimate of the statistical distinction between the two pairs cosmologies, based on the value of , computed using the mean number of peaks in each bin (averaged over all realizations),
(14)  
(15) 
where is the difference between the average number of peaks in bin in cosmology and in cosmology . We find between and (when then same realization of the initial conditions are used in both cosmologies) and (when different realizations are used). The covariance matrix is calculated from the cosmology in both cases. Assuming that the likelihood distributions are Gaussians, these values would correspond to a and difference between the two cosmologies. This is in good agreement with the 95% exclusion of the mean of the distribution we find from the individual realizations in the model – indicating that interpreting the directly as a probability would only slightly overestimate the statistical difference between the two cosmologies. The corresponding value between the and models is (with the same realization of the initial conditions), again confirming that the true likelihood is only somewhat weaker than this simple estimate for the exclusion of the model.
In the above analysis, we have chosen convergence bins whose width was varied, so that they contain approximately equal number of counts ( peaks/bin at 1 arcmin smoothing). More specifically, the bin boundaries are located at . While equal–width bins would perhaps be a more natural choice, the number of peaks falls off exponentially with peak height, and this would have the undesirable result of having many almost empty bins, leading to large values excessively often (overestimate of improbability of unlikely events), and strong deviations from a true distribution (since the peak counts do not have a Gaussian distribution, especially for bins with low average counts).
Iii Weak Lensing Power Spectrum
Because our pipeline to produce the WL maps is newly assembled, GADGET2 was modified to allow , and the raytracing and WL map analysis codes were newly written, we performed a consistency check, by comparing the power spectrum of the convergence in our maps with the power spectrum derived from semi–analytical predictions. This comparison also determines the range of scales which is adequately represented by our simulations due to resolution and box size. In Figure 3, we show the power spectrum of convergence, averaged over 200 realizations of our simulations (solid curves), together with the theoretical expectations (dashed curves). The latter were obtained by direct line–of–sight integration using the Limber approximation (40), and using fitting formulas for the nonlinear 3D matter power spectrum from (41). From top to bottom, the three sets of curves correspond to (green curves), (red curves), and (blue curves).
We find significant deviations from the theoretical power spectra on both small and large scales, as expected. On small scales, with , we are missing power because of the finite resolution on the mass planes and on the raytracing grid. On large scales, with , the spectrum is suppressed because of the finite size of the simulation box. These two wave numbers thus set the range, from arcmin to degree, on which the absolute value of the converge power spectrum is accurately captured.
While we can not trust our results beyond these limits in the subsequent analysis, it should be noted that the ratio of the power spectra in different cosmologies is preserved accurately, even where the simulated power spectra start to deviate from the theoretical predictions. Since we are mainly interested in the comparison between the cosmologies, this makes our results even more robust.
Another consistency check we perform is a comparison between the number of peaks above a certain threshold to the number of clusters creating lensing peaks above that same threshold. This comparison will presented in Sec. IV.1 below.
Iv Results and Discussion
We have found several important quantitative and qualitative results in our study. We first discuss the difference in peak counts between the cosmologies – in the total number of peaks, as well as in the number of peaks in individual convergence bins. We then present our inferred sensitivity to distinguish these cosmologies, and extrapolate this sensitivity from our degree field to an LSST–like survey. We discuss the validity of such an extrapolation by studying the variance in peak counts as a function of angular scale. In the next subsection, we comment on correlations between peaks in different bins, with source galaxies at different redshifts, and with different angular smoothing scales, as well as crosscorrelations among these. Finally, we briefly contrast our results to other predictions for an LSST–like WL survey, namely constraints from the weak lensing shear power spectrum, the PDF of found in (30), as well as cluster number counts.
iv.1 Peak Counts
The fundamental quantity we are interested in is the number of peaks in different convergence bins. This quantity is shown in Figure 4 in the three different cosmologies, derived from convergence maps with a single source galaxy plane at , after the addition of ellipticity noise and Gaussian smoothing with arcmin. From top to bottom (at the peak of the distribution), the three curves correspond to (blue curves), (red curves), and (green curves).
As this figure demonstrates, the distribution of peak heights becomes narrower with increasing .
The total number of convergence peaks in the noisy, smoothed maps is , , and for , , and , respectively. Formally, this is a difference for and a somewhat weaker difference of for . The same numbers for the smoothed noiseless maps (still with 1 arcmin smoothing) are , , and , respectively. Interestingly, the addition of ellipticity noise boosts both the total number of peaks (by ; note that a similar increase was observed by (16)), and the significance in the difference in peak counts between pairs of cosmologies. Apparently, the number of new peaks introduced by the presence of noise is, itself, –dependent, in a way that helps in distinguishing the two cosmologies.
While this result may first sound surprising, the addition of noise is known boost the discriminating power of related observables. For example, scatter in the massobservable relation generally increases the number of detectable clusters and tightens constraints from cluster number counts (e.g Haiman et al. (2001)). In our case, one can image that the number of new peaks on 1 arcmin scales, at a given peak height, introduced by the presence of noise, depends on the amplitude of the –fluctuations on larger angular scales, which, itself, is cosmology–dependent. The precise origin of this result will be investigated in a followup paper.
Most of the peaks, of course, have a low amplitude. In order to compare our results with previous work, it is of interest to consider the number of peaks at higher significance. As an example, above the commonly used threshold of (for a source galaxy density of arcmin at redshift , the value of this threshold corresponds to [eq. 9]), the average number of peaks we identify in our 200 pseudoindependent realizations is , , and for , and respectively, which represent differences between the and the models. These numbers are with ellipticity noise included. Without ellipticity noise, the number of peaks above the same threshold is , , and , which, again, represent a differences. Although the number of these rare high– peaks is also boosted by the presence of noise, apparently this boost does not enhance the distinguishing power between cosmologies.
It is also useful to compare the number of peaks to the estimated number of clusters causing a convergence signal above the same threshold . To be closer to predictions in earlier work, for this comparison, we take arcmin, resulting in a threshold of for 1 arcmin smoothing. In our noiseless maps, we obtain and for and respectively. To find the corresponding number of clusters, we repeat the calculation in (24), which is based on the fitting formula for the DM halo mass function from (42), and assumes clusters follow the NavarroFrenkWhite profile (43), but we adjust the cosmological parameters (including the –dependent value of ) to the values adopted here. We find the number of clusters is , and for and , respectively (the errors bars here are approximately twice the Poisson errors, enhanced by cosmic variance, following (44)). Apparently, there are approximately twice as many convergence peaks in the maps without ellipticity noise as there are clusters, indicating a substantial number of projections. This is roughly consistent with the results of (White et al. (2002); (16); (17)), which self–consistently identify halos in their simulations to quantify this correspondence more reliably, although those works find a somewhat closer correspondence between clusters and peaks. Adding ellipticity noise increases the number of peaks by another factor of two, resulting in plenty of convergence peaks that are not due to a single cluster, which can carry additional information about the cosmology.
An important point to note here is that, as we emphasized earlier, we keep the primordial amplitude fixed as we vary , which leads to a different for each . The difference in the average number of clusters (48.8 versus 31.7) quoted above arises essentially entirely from this difference in . We expect that the difference in peaks counts will similarly be driven primarily by . Tomography should break some of the degeneracy between and ; this is an issue that will be clarified and explored in our followup work.
iv.2 Three Peak Types
It is instructive to next take a more detailed look at how the number of convergence peaks depends on the peak height, and, in particular, how the counts in different convergence bins vary with cosmology. In the top panel of Figure 5, we show the difference in peak numbers (averaged over the 200 realizations) in degree fields in convergence intervals , with between the (blue) and the CDM models, as well as CDM and (green). The middle panel in the same figure shows the standard deviation in each bin, and the bottom panel depicts the contribution to the and thus to the distinction power, derived from each bin (note that this bottom panel shows the unequal bins that were used in our calculation).
Based on the results shown in this figure, we identify three categories of peaks.

High Peaks (): As is known from previous work (White et al. (2002); (16); (17)), a significant fraction (and perhaps the majority) of these peaks are caused by individual collapsed clusters, which can cause strong peaks by themselves. These peaks are more numerous in cosmologies with more negative . This is as expected: as mentioned above, the more negative is, the later dark energy starts to dominate the energy density, and the more time massive clusters had to form.

Medium Peaks (): The number of these peaks significantly exceeds the expected number of clusters that would cause, by themselves, peaks with a similar height. We examined the contributions of individual lens planes along the line of sight (LOS) to the total convergence in these peaks. We found that for the majority of these peaks, multiple planes contribute significantly, and we therefore suspect these peaks are typically caused by projection along the LOS. In contrast with the high peaks, these medium peaks are more numerous in cosmologies with less negative . This was unexpected. Furthermore, we find that the difference in the counts of these medium peaks dominates the cosmological sensitivity. Clarifying the physical nature of these peaks is important, and will be presented in a followup paper (45).

Low Peaks (): Similarly to the high peaks, these peaks are more numerous in cosmologies with more negative . However, their number difference compared to their variance is too low to make an appreciable contribution to the distinction between the cosmologies (the number difference of the low peaks is only a fifth of that of the medium peaks, yet their variance still about half). We speculate that these peaks typically reside in “convergence voids” (large areas of the map with very low values of convergence). We have indeed verified that the maps in cosmologies with more negative show a larger fraction in voids, which gives very low–amplitude peaks a larger area to reside in.
iv.3 Distinguishing Models from Cdm
We use the above difference in peak numbers in bins with different values of to quantify the statistical distinction between the three models with , and . We find that the number of peak height bins used has a minor effect, as long as their number is kept between 4–9. Below this value, the binning is too coarse to capture all the information of the rich structure in the top panel of Figure 5, while higher numbers of bins show signs of overbinning, especially when several source redshifts and smoothing scales are combined (see Appendix VI).
As our fiducial set–up, we use a single source plane at and arcmin smoothing, as in the previous section. Redshift tomography and using different smoothing scales are addressed in Sec. IV.5 below. The for each of our degree convergence fields (with ellipticity noise and smoothing), compared to the reference cosmology CDM, is calculated as described in § II.4. This leads to the probability distribution of the shown in Figure 2.
We can derive our statistical conclusions from these distributions. For example, the mean of the distribution for lies at the 95% tail of the distribution for . This suggests that the distinction from a single degree field corresponds to . The corresponding result for the case of is an exclusion at the 85% confidence level.
iv.4 Extrapolation to an All–Sky Survey
While we are limited by computational costs to studying single degree fields, the forthcoming surveys listed in § I will cover much larger areas of the sky  up to 20,000 degree. Therefore, first we would like to extrapolate our results to a larger field of view – or equivalently, translate it to the tighter sensitivity to expected from larger fields.
We will use two methods of extrapolation, an optimistic and a pessimistic one. In the optimistic case, we neglect largescale correlations and assume that, with increasing survey size, the variance per area–squared in numbers of peaks just falls off as the number of peaks, i.e. . In the pessimistic case, we extrapolate the scaling of the variance with survey size from the scaling behavior in subfields of our degree weak lensing maps. We note that correlations on large scales generally fall off faster than the Poisson–like extrapolation would suggest ((30)), and in reality, both cases may therefore be conservative estimates.
We illustrate these two methods of extrapolations in Figure 6, where we plot the variance in peak numbers for progressively larger fields, starting with a arcminute field and going all the way to degrees (always obtained from our 200 pseudoindependent realizations).
The optimistic extrapolation, shown by the black line in Figure 6, scales to a stronger sensitivity for a 20,000 degree field of view. The difference in the dark energy equation of state parameters in our study is , which we discern at . Therefore, roughly, a 20,000 degree survey will be able to discern a difference of at the same level. The somewhat more pessimistic case of extrapolating from existing subfields, as depicted by the green line in Figure 6, yields a variance that is larger by a factor of 1.3. Therefore, for this case, our estimated sensitivity from an LSST size survey is .
While this is an interesting level of sensitivity, one has to keep in mind that this is a single–parameter constraint, and is not directly comparable to marginalized constraints usually quoted in the literature which fold in degeneracies with other parameters (see § IV.6 below).
iv.5 Different Redshifts and Smoothing Scales
In our fiducial set–up so far, we have used a single source plane () and smoothing scale ( arcmin). Here we examine using other source planes and smoothing scales, including the possibility of using multiple scales in combination.
Single Redshifts and Smoothing Scale Combinations
We have re–computed distributions, following the procedure described above, for a variety of different smoothing scales ( arcmin) and source redshifts (). We quote the results by indicating, in each case, the fraction of the realizations of the model which can be excluded from the cosmology at the 68% (95%) confidence level. These fractions are to be compared to the values of 32% and 5%, expected in the absence of any signal (i.e., for the same cosmology, compared to itself). These fractions are presented in Table 1. We list in the table as well, but we caution that, as discussed above, our simulations do not have sufficient resolution to reproduce the convergence power spectrum on these scales, and therefore these results may not be reliable. Unless noted otherwise, we exclude these values from the discussion in the rest of this paper.
% exclusion chance  

@ 68 (95)% CL  
48 (13) %  76 (32) %  90 (46) %  
56 (16) %  68 (18) %  84 (34) %  
53 (12) %  66 (13) %  50 (18) %  
41 (10) %  48 (14) %  40 (6) %  
39 (6) %  32 (7) %  37 (12) %  
38 (6) %  28 (8) %  39 (10) %  
26 (2) %  28 (5) %  25 (4) % 
Table 1: The probability of distinguishing the model from CDM, based on peak counts in a single degree field, at the 68% (95%) confidence level for various source redshifts and smoothing scales , using 8 peak height bins.
There are several trends worth noting in Table 1. First, using source galaxies at high redshift () provides substantial advantage over using sources at redshift or even (recall that the surface density is fixed in each case to be 15 arcmin). Dark energy constraints from Type Ia supernovae, and from most other techniques also benefit from deep surveys and high redshifts. For weak lensing, going to higher redshifts is particularly important, as the best lensing efficiency comes from approximately halfway in distance between sources and observer.
Second, as is well known, weak lensing derives most of its power from small scales, and indeed, the smaller the smoothing scale in our maps, the stronger the distinction. This is true down to the smallest angular scales that we can reliable simulate, . Smoothing scales larger than do not provide very strong constraints by themselves, but may still be useful when combined with smaller ones (especially for real surveys with much larger fields of view – the number of peaks at these large scales becomes very small for our small degree fields).
Correlations
There are several types of correlations that would be useful to investigate: i.e., correlations in the number of peaks in different height bins, identified for different source redshifts, and with different smoothing scales, as well as the crosscorrelations between these quantities. In order to study the effect of these correlations, we computed the distributions from covariance matrices in which only diagonal elements, or diagonal blocks were retained, and the crossterms outside these were set to zero.
We find that the peak height bins are correlated. Neglecting all correlations, by retaining only the diagonal elements of the covariance matrix, lead to histograms that resemble distributions with fewer degrees of freedom, and an excessive amount of outliers, both typical signs that the dataset contains more correlations than is assumed in the analysis. In particular, as expected, the distributions (analogous to those shown in Figure 2) develop wider tails toward high values, with their maximum moving to somewhat smaller values. The net result is that the overlap of the distributions between a pair of cosmologies is enhanced, and the distinction between the cosmologies is degraded. Table 2 illustrates this in numbers.
% exclusion chance  

@ 68 (95)% CL  
32 (1) %  46 (2) %  52 (2) %  
34 (2) %  54 (2) %  65 (8) %  
40 (3) %  51 (5) %  44 (8) %  
39 (2) %  44 (6) %  38 (4) %  
36 (2) %  30 (6) %  32 (5) %  
32 (4) %  32 (3) %  36 (6) %  
24 (2) %  32 (4) %  26 (3) % 
Table 2: Same as Table 1, but with correlations between peak height bins neglected. The likelihood of distinguishing from CDM are much lower, and approach the nosignificance values of 32 (5) %, in particular for the 95% confidence level numbers, indicating a large number of outliers.
We also attempted to study crosscorrelations between different redshifts and smoothing scales. When studying several redshifts and smoothing scales simultaneously, the number of elements in the covariance matrix grows quickly. Because of the low number of realizations and the corresponding underestimation of the variance, results obtained by using more than bins become unreliable (Appendix VI). Unfortunately, this prohibits us from studying the effect of correlations between redshifts and smoothing scales reliably: the number of different peak height bins still allowed becomes too small. In particular, when multiple source redshifts or smoothing scales are simultaneously included, only 2–3 peak heights are allowed by the low number of realizations. We found that such a low number of convergence bins does not sufficiently capture the peak number difference distributions in Figure 5, and no meaningful constraints are possible. We defer a study of these cross–correlations to a follow–up paper, with a larger number of independent realizations of each cosmology.
iv.6 Comparison to Other Forecasts for LSST
It is useful to place our results in the context of constraints expected to be available from other uses of large WL maps. Three methods we will briefly mention here are (i) the shear power spectrum, (ii) cluster number counts, and (iii) the one–point function of the convergence PDF.
(i) From an 11parameter fit to the tomographic shear power spectrum from LSST and constraints from Planck, Song and Knox (2004) obtain .
(ii) From a 7parameter fit to the number counts of 200,000 shearselected clusters, combined with Planck constraints, (21) obtain .
(iii) Using the fractional area statistic of areas of high convergence, (30) found for LSST and priors from Planck.
The above are all constraints, marginalized over a large number of parameters. We do not attempt a fair comparison of the three methods, which would require that identical assumptions are made about the survey specifications and about the set of free parameters (see Albrecht et al. (2006) for a similar exercise covering i. and ii.). Nevertheless, these numbers demonstrate, first, that the statistical information contained in the non–linear features in WL maps – captured by (ii) and (iii) – is at least roughly comparable to the information encoded in the power spectrum in (i). Second, our single–parameter sensitivity to can be compared meaningfully to those from clusters counts. In particular, as mentioned in § IV.1, the counts of clusters above the convergence threshold differs in the and cosmologies by . We find that the number of convergence peaks above the same threshold differs in these two cosmologies at a similar, significance level. We find further that considering lower height peaks, with convergence well below this threshold, significantly improves the distinction. While followup work, with simulations that address degeneracies between , , and other parameters is needed, these results suggest that the parameter sensitivity from WL peak counts will be competitive with those from the above three methods.
As mentioned in § I, in a recent preprint, (31) independently addressed the dependence of WL peaks on the background cosmology. (31) used the same simulation box size as in our work, but a smaller number of particles (, rather than ), which allowed them to sample a 2D parameter space and initial conditions, at the cost of mass and angular resolution. They compute the aperture mass for the peaks as a weighted integral of the tangential shear, rather than the convergence. Despite these differences, the conclusions of the two papers appear to be remarkably similar. In particular, if we attribute the distinction we find between our cosmologies with different values entirely due to the differences in in these cosmologies, we find a very similar sensitivity to the latter parameter as their Figure 3: (for a fixed ). We find, additionally, a subdominant type of peaks – low–amplitude peaks, residing in the voids – the distinction between the cosmologies is however dominated by medium peaks, with a strong contribution from high peaks.
V Summary
We propose the new method of using the counts of peaks, identified in weak gravitational lensing maps, to constrain dark energy and other cosmological parameters. The method makes no reference to cluster counts or their mass function, and so, by definition, is free from selection effects inherent in methods that target clusters.
We used N–body simulations to create theoretical convergence maps in three cosmologies with , , and . We identify three different peak types: high, medium, and low peaks, whose numbers depend differently on cosmology. The high () and low () peaks both become more numerous with increasingly negative , while the opposite trend is true for the medium peaks. For the high peaks, this dependence confirms the expectation that when is more negative, dark energy begins to dominate later, and more time is allowed for the growth of non–linear structures. The opposite behavior of the medium peaks is a surprise  which is all the more important since we find that these peaks dominate the distinction between the cosmologies. The majority of these peaks are not typically caused by individual collapsed clusters, and are more likely due to projection of large scale structure along the line of sight. The low peaks form a subdominant contribution. These peaks are typically due to random ellipticity noise, and we speculate that the larger surface area of voids in cosmologies with more negative , which we have detected in the maps, offers more opportunities for these small fluctuations to arise.
We have shown that based on these peak counts, cosmologies with and can be distinguished from CDM with a single degree convergence field at a 95% and 85% confidence level, respectively. Extrapolating these results to the 20,000 degree field of view of a large survey, such as LSST, we obtain a projected sensitivity to of –, depending whether the extrapolation is based on the variance measured in subfields of our maps, or whether a simple Poisson–like scaling is used.
In this work, we have explored three constant values of the dark energy equation of state, , , and , representing variations around the best fit CDM model to the 5year WMAP data. The most significant concern is that there will be strong degeneracies between and other parameters; our results may already be driven primarily by the sensitivity in the number of peaks. In a followup paper, we will study such degeneracies; it will also be necessary to perform a larger number of realizations of each cosmology, in order to study the utility of combining several source galaxy redshifts and to use information possibly contained in the angular size distribution of the peaks. In a separate publication (45), we also plan to investigate in more detail the physical origin of the three kinds of peaks detected in this work.
Acknowledgements.
We would like express our deep thanks to Lam Hui for numerous helpful discussions, and to Greg Bryan, Francesco Pace, and Matthias Bartelmann and his group for invaluable help during code development. We are also grateful to Christof Wetterich for kindly supporting an extended visit by JK at the University of Heidelberg. We also thank Puneet Batra, Wenjuan Fang, Eugene Lim and Sarah Shandera for useful discussions about statistics, CAMB, and lensing, and Volker Springel for his help with GADGET2 and for providing us with his parallelized initial conditions generator. JK is supported by ISCAP and the Columbia Academic Quality Fund. This work was supported in part by the Polányi Program of the Hungarian National Office for Research and Technology (NKTH), by NSF grant AST0507161, by the U.S. Department of Energy under contract No. DEAC0298CH10886, and by the Initiatives in Science and Engineering (ISE) program at Columbia University. The computational work for this paper was performed at the LSST Cluster at Brookhaven National Laboratory and with the NSF TeraGrid advanced computing resources provided by NCSA.Vi Appendix: Initial Conditions and the Number of Realizations
In this Appendix, we show our results as a function of the number of peak height bins used, and justify the choice for the number of bins we made in the body of the paper.
Ideally – in the limit that the covariance matrix was perfectly computed – the distinguishing power would asymptotically approach a constant value as one increases the number of bins, at which point using more bins will not provide any more benefit. Unfortunately, it is computationally expensive to generate a large ensemble of initial conditions from which to construct several strictly independent realizations of a given cosmology (a single realization takes approximately 4800 CPU hours, i.e. ca. 2.5 days on 80 CPUs). This, together with the storage required for the large amount of data produced in each run, ultimately limits the number of bins we can utilize.
In particular, using too few realizations, together with too many bins, will generally result in an underestimate of the variance, and an artificial boost in the values. The possible concern, then, is that the differences we find in the convergence maps between two cosmologies may be caused by such an underestimation of the variance.
In order to guard against this possibility, and to verify that the differences in the convergence maps are caused by genuine differences in cosmology, we generated a second set of independent initial conditions for the cosmologies with and with . The simple Ansatz we adopt for the maximum number of bins that can be safely used, without underestimating variances, is the following: when two different realizations of the same CDM cosmology are compared, we should recover the distribution, i.e., we must not be able to rule out the fiducial model itself.
In Figure 7, we show the percentage of the realizations of the test cosmology that lies beyond the 68% (left panel) and 95% (right panel) tail of the distribution (solid curves), together with the same quantities for the fiducial CDM cosmology itself (dashed curves). Our Ansatz requires that the dashed curves correspond to the confidence level at which the exclusion is to be achieved (e.g. no more than 32 / 5% of maps excluded at the 68 / 95% confidence level). We see that this condition is satisfied when 49 bins are used, but for 10 or more bins, there is a steadily rising difference between the two realizations of the CDM model. We therefore take 9 bins as the maximum acceptable total number of bins.
Figure 7 also shows that the distinction between and the CDM model is relatively steady for 49 bins, and is insensitive to the choice of the realization of the initial condition. However, for a larger number of bins, the distinction starts to be either less significant (when different realizations of the initial conditions are used in the two cosmologies; solid curves), or begin a modest rise in significance (when the same realizations are used; dotted curves).
While the above analysis is based only on the number of convergence bins, we find similar conclusions when the bins include multiple source galaxy redshifts and/or multiple smoothing scales. We conclude that the accuracy of our covariance matrix limits us to a total of 9 bins.
Footnotes
 preprint: astroph/yymmnnn, March 6, 2018
 http://www.strw.leidenuniv.nl/ kuijken/KIDS
 https://www.darkenergysurvey.org
 It is important to note here, however, that we fix the primordial amplitude of the density perturbations; this results in a –dependent value of the present–day normalization . In the rest of this paper, whenever we compare cosmologies with two different values of , it should be remembered, even when not explicitly stated, that these two cosmologies also have different values of .
References
 Albrecht, A., G. Bernstein, R. Cahn, W. L. Freedman, J. Hewitt, W. Hu, J. Huth, M. Kamionkowski, E. W. Kolb, L. Knox, J. C. Mather, S. Staggs, and N. B. Suntzeff, 2006: Report of the Dark Energy Task Force. eprints astroph/0609591 (2006).
 A. Refregier, Ann. Rev. Astron. Astrophys. 41, 645 (2003). [arXiv:astroph/0307212].
 N. Kaiser, J. L. Tonry and G. A. Luppino, arXiv:astroph/9912181 (1999).
 J. A. Tyson, D. M. Wittman, J. F. Hennawi and D. N. Spergel, Nucl. Phys. Proc. Suppl. 124, 21 (2003). [arXiv:astroph/0209632].
 J. A. Tyson and the LSST Collaboration, Proc. SPIE Int. Soc. Opt. Eng. 4836, 10–20 (2002). See also http://www.lsst.org
 N. Kaiser, “Weak Gravitational Lensing Of Distant Galaxies,”
 B. Jain and U. Seljak, Astrophys. J. 484, 560 (1997).
 W. Hu, Astrophys. J. 522, L21 (1999).
 W. Hu, Phys. Rev. D 66, 083515 (2002).
 D. Huterer, Phys. Rev. D 65, 063001 (2002).
 A. Refregier et al., Astron. J. 127, 3102 (2004).
 K. N. Abazajian and S. Dodelson, Phys. Rev. Lett. 91, 041301 (2003).
 M. Takada and B. Jain, Mon. Not. Roy. Astron. Soc. 348, 897 (2004).
 Y. S. Song and L. Knox, Phys. Rev. D 70, 063510 (2004).
 M. J. . White, L. van Waerbeke and J. Mackey, Astrophys. J. 575, 640 (2002).
 T. Hamana, M. Takada and N. Yoshida, Mon. Not. Roy. Astron. Soc. 350, 893 (2004).
 J. F. Hennawi and D. N. Spergel, Astrophys. J. 624, 59–79 (2005).
 L. M. Wang and P. J. Steinhardt, Astrophys. J. 508, 483 (1998).
 Z. Haiman, J. J. Mohr and G. P. Holder, Astrophys. J. 553, 545 (2001).
 J. Weller, R. Battye and R. Kneissl, Phys. Rev. Lett. 88, 231301 (2002).
 S. Wang, J. Khoury, Z. Haiman and M. May, Phys. Rev. D 70, 123008 (2004).
 L. Marian, R. E. Smith and G. M. Bernstein, Astrophys. J. 698, L33 (2009).
 T. Chantavat, C. Gordon, and J. Silk, Phys. Rev. D, 79, 083508, (2009).
 W. J. Fang and Z. Haiman, Phys. Rev. D 75, 043010 (2007).
 M. Takada, and S. Bridle, New Journal of Physics, 9, 446 (2007).
 L. Marian and G. M. Bernstein, Phys. Rev. D 73, 123525 (2006).
 M. Schirmer, T. Erben, M. Hetterscheidt and P. Schneider, Astron. Astrophys. 462, 875 (2007).
 J. P. Dietrich, T. Erben, G. Lamer, P. Schneider, A. Schwope, J. Hartlap and M. Maturi, Astron. Astrophys. 470, 821(2007).
 B. Jain and L. V. Van Waerbeke, Astrophys. J. 530, L1 (2000).
 S. Wang, Z. Haiman and M. May, Astrophys. J. 691, 547 (2009).
 J. P. Dietrich and J. Hartlap, eprint arXiv:0906.3512.
 V. Springel, Mon. Not. Roy. Astron. Soc. 364, 1105 (2005).
 A. Lewis, A. Challinor and A. Lasenby, Astrophys. J. 538, 473 (2000).
 U. Seljak and M. Zaldarriaga, Astrophys. J. 469, 437 (1996).
 R. W. Hockney, J. W. Eastwood, Adam Hilger, Bristol (1988).
 D. C. Wells, E. W. Greisen and R. H. Harten, Astronomy and Astrophysics Supplement Series 44, 363370 (1981).
 T. Hamana and Y. Mellier, Mon. Not. Roy. Astron. Soc. 327, 169 (2001).
 S. Dodelson, C. Shapiro and M. J. . White, Phys. Rev. D 73, 023009 (2006).
 L. Van Waerbeke, Mon. Not. Roy. Astron. Soc. 313, 524 (2000).
 D. N. Limber, Astrophys. J. 117, 134 (1953).
 Smith, R. E., J. A. Peacock, A. Jenkins, S. D. M. White, C. S. Frenk, F. R. Pearce, P. A. Thomas, G. Efstathiou, and H. M. P. Couchman, MNRAS, 341, 1311 (2003).
 A. Jenkins et al., Mon. Not. R. Astron. Soc. 321, 372 (2001).
 J. F. Navarro, C. S. Frenk and S. D. M. White, Astrophys. J. 490, 493 (1997).
 W. Hu, & A. V. Kravtsov, Astrophys. J. 584, 702 (2003).
 X. Yang, J. .M. Kratochvil, Z. Haiman, M. May, in preparation.
 S. Pires, J. L. Starck, A. Amara, A. Refregier and R. Teyssier, eprint arXiv:0904.2995