© ESO, 2015

1. Introduction

Weak gravitational lensing (WL) probes matter structures in the Universe. It contains information from the linear growth of structures to the recent highly nonlinear evolution, going from scales of hundreds of Mpc down to sub-Mpc levels. Until now, most studies have focused on two-point-correlation functions, but the non-Gaussianity of WL cannot be ignored if one aims for a deep understanding of cosmology.

One simple way to extract higher order WL information is peak counting. Peaks are defined as local maxima of the projected mass measurement. They are particularly interesting for at least two reasons. First, peaks are tracers of high-density regions. While other tracers of halo mass such as optical richness, X-ray luminosity or temperature, or the SZ Compton-y parameter depend on scaling relations and often require assumptions about the dynamical state of galaxy clusters such as isothermal equilibrium and relaxedness, lensing does not. It therefore provides us with a direct way to study cosmology with the cluster mass function. Second, the lensing signal is highly non-Gaussian, and two-point-function-only studies deprive one of the information richness beyond second order. For example, Dietrich & Hartlap (2010) show that parameter constraints can be highly improved by joining peak counts and using second-order statistics, and Pires et al. (2012) find that peak counts capture more than the convergence skewness and kurtosis of the non-Gaussian information. Another advantage of WL peaks is information about the halo profile. Mainini & Romano (2014) showed that combining peak information with other cosmological probes provides an interesting way to study the mass-concentration relation.

For studies of the mass function via X-ray or the SZ effect, most works have adapted a reverse-fitting approach. This means that from diverse observables, one first establishes the observed mass function and then fit it with a theoretical model. To extract the mass function, this process needs to reverse the effect of selection functions, to use scaling relations, and to make further assumptions about sample properties. Alternatively, one can proceed with a forward-modeling approach: starting from an analytical mass function, we compute predicted values for observables and compare them to the data to carry out parameter fits (Fig. 1). The corresponding forward application of selection functions is typically much simpler than its reverse. Moreover, instrumental effects can be easily included, and model uncertainties can be marginalized over. Forward modeling requires well-motivated models of physical phenomena, which is challenging in the case of observables derived from baryonic physics, yet Clerc et al. (2012) still provide a forward analysis from X-ray observations. For WL peak counts, however, computing the observable prediction is more straightforward, as long as using some appropriate assumptions.

Fig. 1 Forward- and reverse-modeling diagram for the mass function studies. Two different approaches to establishing links between the theoretical mass function (the upper round rectangle) and the observables (the lower round rectangle) are to compare the observable mass function nobs with the theoretical one (reverse modeling), or to compare “predicted” observable values Xmodel with observations (forward modeling). In this paper, the forward modeling is adopted and we propose a new method to “predict” peak conuts.

One of the difficulties of predicting WL peak counts is that peaks can come from several mass overdensities at various redshifts due to projection effects (Jain & Van Waerbeke 2000; Hennawi & Spergel 2005; Kratochvil et al. 2010). This makes counting nonadditive even in the linear regime, and the prediction becomes less trivial. To overcome this ambiguous effect, some previous works have used N-body simulations, e.g., Dietrich & Hartlap (2010). They perform peak counts from N-body runs with different paremeter sets to obtain confidence contours for constraints. However, since N-body simulations are very costly in terms of computation time, input parameter sets should be carefully chosen, and an interpolation of results is needed. Thus the resolution in the parameter space is limited, and the Fisher matrix is only available for the fiducial parameters.

Alternatively, there have been several attempts at peak-count modeling. Maturi et al. (2010) propose to study contiguous areas of high-signal regions instead of peaks, and provide a model that predicts the amount of this alternative observable. Meanwhile, Fan et al. (2010, hereafter FSL10 propose a model for convergence peaks by supposing at most one halo exists on each line of sight. Both models are analytical and based on calculation from Gaussian random field theory. A comparison of the FSL10 model with observation has been shown by Shan et al. (2014), using the data from the Canada-France-Hawaii Telescope Stripe 82 Survey.

However, these models encounter difficulties for additional complications and subtleties. On one hand both models require Gaussian noise and linear filters, otherwise the Gaussian random field theory becomes invalid. As a result, nonlinear, optimized reconstruction methods of the projected overdensity are automatically excluded. On the other hand realistic scenarios, such as mask effects and intrinsic ellipticity alignment, introduce asymmetric changes into the peak counts. The impact of these additional effects are unpredictable in purely analytical models. This encourages us to propose a new model for WL peak counts.

In this paper, we adopt a probabilistic approach to forecasting peak counts. This can be handled by our Camelus algorithm (Counts of Amplified Mass Elevations from Lensing with Ultrafast Simulation). Unlike N-body simulations which are very time-consuming, we create “fast simulations” by sampling halos from the mass function. The only requirement is a cosmology with a known mass function and halo mass profiles. To validate this method and to justify various hypotheses that our model makes, we compare results from our fast simulations to those from N-body runs. This approach is similar to the sGL model of Kainulainen & Marra (2009, 2011a,b), where they show that the stochastic process provides a quick and accurate way to recover the lensing signal distribution.

The outline of this paper is as follows. In Sect. 2, we recall some of the WL formalism and theoretical support for our model. In Sect. 3, a full description of our model is given. In Sect. 4, we give the details concerning the N-body and the ray-tracing simulations. Finally, the results are presented in Sect. 5, before we summarize and conclude in Sect. 6.

2. Theoretical basics

In this section, we define the formalism necessary for our analysis. To model the convergence field lensed by halos, we need to specify their profile, projected mass, and distribution in mass and redshift, which is the mass function.

2.1. Weak-lensing convergence

Observationally, galaxy shape distortions can be displayed at linear order in the form of the lensing distortion matrix . For an angular position θ, is given by (1)which defines two WL observables: convergence κ and shear γ. The latter is a complex number given by γ = γ₁ + iγ₂. This linearization of the light distortion can be calculated explicitly in general relativity. Accordingly, the matrix elements are linked to second derivatives of the Newtonian gravitational potential φ by (2)where f_K is the comoving transverse distance and δ_ij the Kronecker delta. In particular, an explicit expression of κ is given as follows (see, e.g., Schneider et al. 1998), (3)where H₀ is the Hubble parameter, Ω_m the matter density, c the speed of light, a(w′) represents the scale factor at the epoch to which the comoving distance from now is w′, and δ is the matter density contrast.

2.2. Halo density profile and its projected mass

Consider now a dark matter (DM) halo with a Navarro-Frenk-White (NFW) density profile (Navarro et al. 1996, 1997), given by (4)where ρ_s and r_s are the characteristic mass density and the scale radius of the halo, respectively, and α is the inner slope parameter. The concentration parameter c_NFW is defined as the ratio of the virial radius to the scale radius, c_NFW = r_vir/r_s. We assume the following expression (proposed by Takada & Jain 2002): (5)where M is the halo mass and M_⋆ the pivot mass such that δ_c(z = 0) = σ(M_⋆), with δ_c the threshold overdensity for the spherical collapse model, and σ²(M) is the variance of the density contrast fluctuation smoothed with a top-hat sphere with radius R such that .

In this paper, we take c₀ = 8, α = 1, and β = 0.13. The value of α corresponds to the classical NFW profile. The value of β is provided by Bullock et al. (2001), and c₀ corresponds to the best-fit value, using r_vir, r_s, z, M derived from the N-body simulations that we use and fixing β. For δ_c, we use the fitting formula of Weinberg & Kamionkowski (2003) with (6)and (7)Lensing by an NFW halo is characterized by its projected mass. More precisely, defining the scale angle θ_s = r_s/D_ℓ as the ratio of the scale radius to the angular diameter distance D_ℓ between lens and observer, we get (following Bartelmann 1996; Takada & Jain 2003)^1,2(8)with (9)where the quantities D_s and D_ℓs are the angular diameter distances between source and observer, and lens and source, respectively, and (10)We have truncated the projected mass distribution at θ = c_NFWθ_s. Equation (8) is used and computed for the ray-tracing simulations with NFW halos.

2.3. Halo mass function

The halo mass function n(z,<M) indicates the halo number density with mass less than M at redshift z³, often characterized by a function f(σ,z) as (11)where is the current matter density, and σ(z,M) is defined as σ(M) multiplied by the growth factor D(z). In this study, we adopt the model proposed by Jenkins et al. (2001) in which a fit for f is given as (12)

3. A new model for WL peak counts

3.1. Probabilistic approach: fast simulations

Our model is based on the idea that we can replace N-body simulations with an alternative random process, such that the relevant observables are preserved, but the computation time is drastically reduced. We call this alternative process “fast simulations”, which are produced by the following steps:

• 1.

  generate halo masses by sampling from a mass function;

• 2.

  assign density profiles to the halos;

• 3.

  place the halos randomly on the field of view;

• 4.

  perform ray-tracing simulation.

One can notice that we have made two major hypotheses. First, we assume that diffuse, unbound matter, for example cosmological filaments, does not significantly contribute to peak counts. Second, we suppose that the spatial correlation of halos has a minor influence, since this correlation is broken down in fast simulations. Previous work has shown that correlated structures influence number and height of peaks by only a few percentage points (Marian et al. 2010). Furthermore, assuming a stochastical distribution of halos can lead to accurate predictions of the convergence probability distirbution function (Kainulainen & Marra 2009). One may also notice that halos can overlap in 3D space, and indeed we do not exclude this possibility. We test and validate these hypotheses in Sect. 5, and discuss possible improvements to our model in Sect. 6

Although we have chosen NFW profiles for the density of DM halos, using any halo profile model for which the projected mass is known is of course possible, such as triaxial halos or profiles offered by baryonic feedback (Yang et al. 2013). In addition, our prediction model is completely independent of the method by which peaks are extracted from the WL data. The same analysis can be applied to data (or N-body simulations + ray-tracing) and to fast simulations. Moreover, survey characteristics, such as masks, photometric redshift errors, PSF residuals, and other systematics, can be incorporated and forward-propagated as model uncertainties. Furthermore, the halo sampling technique is much faster than a full N-body run. For instance, it only takes a dozen seconds on a single-CPU desktop computer to generate a box that is large enough for our use (see specifications in Sect. 4.2).

This is a probabilistic approach to forecast peak counts, and we compare the convergence peaks obtained with those from full N-body runs in order to validate our forward model. This is described in Sect. 5.1.

3.2. Peak selection

In this paper, we focus on convergence peaks. We have followed a classical analysis used in former studies (e.g., Hamana et al. 2004; Wang et al. 2009; Fan et al. 2010; Yang et al. 2011) to extract peaks.

First, we should highlight that κ and κ_proj (respectively given by Eqs. (3)and (8)) do not follow the same definition. Actually, Eq. (8) can be recovered by replacing δ with in Eq. (3). This means that κ_proj does not take lensing by underdense regions into account and is shifted by a constant value, which corresponds to the mass-sheet degenerency. To obtain a model that is consistant with a zero-mean convergence field, we subtract the mean value of κ_proj from our convergence maps, so that (13)We use this approximation throughout this study when ray-tracing is done with projected mass.

Fig. 2 Comparison between an analytical mass function (blue line), halo mass histograms for N-body runs (green points), and sample histograms for realizing of the fast simulation (red crosses). The plots are drawn at 4 different redshift planes, and for each the thickness is dz = 0.1.

Consider now a reconstructed convergence field κ_n(θ) in the absence of intrinsic ellipticity alignment. The presence of galaxy shape noise leads to the true lensing field κ(θ) being contaminated by a linear additive noise field n(θ), such that (14)In general, κ is dominated by n, and one way to suppress the noise is to apply a smoothing: (15)where W(θ) is a window function, chosen to be Gaussian in this study as (16)which is specified by the smoothing scale θ_G. We denote K_N(θ), K(θ), and N(θ) as corresponding smoothed fields to Eq. (14), such that (17)and set θ_G = 1 arcmin in the following.

If intrinsic ellipticities are uncorrelated between source galaxies, N(θ) can be described as a Gaussian random field (Bardeen et al. 1986; Bond & Efstathiou 1987) for which the variance is related to the number of galaxies contained in the filter. This is given by Van Waerbeke (2000) as (18)Here, n_g is the source galaxy number density, and is the variance of the intrinsic ellipticity distribution. We then define the lensing signal-to-noise ratio (S/N) as (19)and the peaks are extracted from the ν field, defined as pixels that have a S/N value higher than their eight neighbors. This implies that peak analyses require S/N values on a well-defined grid (e.g., HEALPix grid). Furthermore, we suppose that source galaxies are uniformly distributed in this study, so σ_noise is a constant. However, this does not have to be true in general.

In summary, convergence peaks are selected by the following steps:

• 1.

  compute the projected mass κ_proj(θ) by ray-tracing;

• 2.

  subtract the mean to obtain κ(θ);

• 3.

  add the noise to obtain κ_n(θ);

• 4.

  smooth the field and acquire K_N(θ);

• 5.

  determine the S/N ν(θ); and

• 6.

  select local maxima and compute the density n_peak(ν).

Only positive peaks are selected, and the analysis is based on the abundance histograms from peak counts. From fast simulation, through ray-tracing, to peak selection, the calculation is carried out by our Camelus algorithm.

4. Simulations

4.1. N-body simulations

As provided by A. Evrard, the N-body simulations “Aardvark” have been used in this study. They were generated by LGadget-2, a DM-only version of Gadget-2 (Springel 2005). The Aardvark parameters had been chosen to represent a WMAP-like ΛCDM cosmology, with Ω_m = 0.23, Ω_Λ = 0.77, Ω_b = 0.047, σ₈ = 0.83, h = 0.73, n_s = 1.0, and w₀ = −1.0.

The DM halos in Aardvark were identified using the Rockstar friends-of-friends code (Behroozi et al. 2013). The field of view is 859 deg². This corresponds to a HEALPix patch with n_side = 2 (for HEALPix, see Górski et al. 2005).

4.2. Fast simulations

As described in Sect. 3.1, our model requires a mass function as input. We chose the model of Jenkins et al. (2001, see Sect. 2.3 to sample halos. This is done in ten redshift bins from z = 0 to 1. We set the sample mass range to the interval 10¹² and 10¹⁷M_⊙/h. For each halo, the NFW parameters were set to be (c₀,α,β) = (8.0,1.0,0.13). We suggest seeing Sect. 2.2 for their definitions.

Figure 2 shows an example of our halo samples, compared to the original mass function, and mass histograms established from the Aardvark simulations. Although halos with high mass can be 10³–10⁵ times less populated than low-mass halos, our sampling is still in a perfect agreement with the original mass function. One may notice a shift and a tilt in the Aardvark halo mass function for low and high redshifts, however, in these regimes, the lensing efficiency is low because of the distance weight term D_ℓD_ℓs/D_s, so this mismatch is not very large.

4.3. Ray-tracing simulations

For the Aardvark simulations, ray-tracing was performed with CALCLENS (Becker 2013). Galaxies were generated using ADDGALS (by M. Busha and R. Wechsler⁴). Ray-tracing information is available only on a subset of 53.7 deg² (a HEALPix patch with n_side = 8), which is 16 times smaller than the halo field. In this study, only galaxies at redshift between 0.9 and 1.1 were chosen for drawing the convergence map. It led to an irregular map, and in order to clearly define eight neighbors to identify peaks, we used a 2D-linear interpolation to obtain κ values on a grid. This was done after carrying out a projection to Cartesian coordinates.

For computational purposes, in order not to handle too many galaxies at a time, we split the field into four “ray-tracing patches”, the size of which is 13.4 deg² each (corresponding to n_side = 16). We then project the coordinates with regard to the center of each patch using the Gnomonic projection. The size lengths of the ray-tracing patches are between 3.5 and 6.2 deg, so small enough to retain a good approximation.

Fig. 3 Patch of a map projected with regard to its center, taken from a realization of fast simulations. The left, middle, and right panels give the fields κ1 / 3, , and , respectively. We have taken the cubic root to emphasize the contrast. It is clear that the signal is completely dominated by the noise. Even though the smoothed map is quite different from the original one, the structures, high-signal regions are still conserved and traced.

For the fast simulations and the two intermediate cases that we study in Sect. 5.1, source galaxies have a fixed redshift z_s = 1.0. They are regularly distributed on a HEALPix grid and placed at the center of pixels. Each ray-tracing pixel is a HEALPix patch with n_side = 16 384, for which the characteristic size is θ_pix ≈ 0.215 arcmin. Thus, the galaxy number density is arcmin^-2. Ray-tracing for fast simulations is carried out after splitting and projection to Cartesian coordinates. There are 64 ray-tracing patches in a halo field, and each patch contains 1024 × 1024 pixels. The convergence was computed using Eqs. (8)–(10). As a remark, no mask was applied in this study.

4.4. Adding noise

Shape noise n(θ) is added to each pixel after we obtain κ(θ) from N-body runs or fast simulations. It is modeled as a Gaussian random field with a top-hat filter with a size that corresponds to the pixel area A_pix. The variance of this is given by Van Waerbeke (2000) as (20)We choose σ_ϵ = 0.4 which corresponds to a CFHTLenS-like survey, and n_gA_pix is chosen to be 1 so that each pixel represents one galaxy. This leads to σ_pix ≈ 0.283. We can also estimate σ_noise with Eq. (18) and obtain σ_noise ≈ 0.024. This shows that a real map is in general dominated by the noise (Fig. 3). Even for a peak at ν = 5, the lensing signal is only on the order of κ = 0.12, less than half of the pixel noise amplitude.

5. Results

5.1. Validation of our model: comparison to N-body runs

To validate our model, we compare it to the N-body simulations. We compute peak abundance histograms from both simulations, together with two intermediate steps. This results in four cases in total:

• Case 1:

  full N-body runs;

• Case 2:

  replacing N-body halos with NFW profiles with the same masses;

• Case 3:

  randomizing angular positions of halos from Case 2;

• Case 4:

  fast simulations, corresponding to our model.

These cases form a progressive transition from full N-body runs toward our model. More precisely, Case 2 tests the hypothesis corresponding to the second step of our model (see Sect. 3.1); i.e., diffuse, unbound matter contributes little to peak counts. Case 3 additionally tests the assumption made in the third step. (Halo clustering plays a minor role.) Finally, Case 4 completes our model with the missing first step. As a result, the halo population and their redshifts are identical to N-body runs in Cases 2 and 3.

Fig. 4 Comparison of the peak abundance from different cases. Upper panel: blue solid line: full N-body runs (Case 1); green circles: replacement of halos by NFW profiles (Case 2); red squares: replacement of halos by NFW profiles and randomization of halo angular positions (Case 3); cyan diamonds: fast simulations, corresponding to our model (Case 4); magenta dashed line: peaks from noise-only maps. Lower panel: we draw the upper and lower limits of error bars shifted with regard to the N-body values. This refers to the standard deviation over 4 maps (green dash-dotted line for Case 2) or 16 maps (red dashed line for Case 3, cyan solid line for Case 4). The field of view is 53.7 deg2.

Figure 4 shows the peak abundance histograms for all four cases. In this section, the field of view is 53.7 deg², since we are limited by the available information of ray-tracing for the N-body runs. For Cases 1 and 2, we compute the average in each histogram bin for eight noise maps. For Cases 3 and 4, this is done with eight realizations (of randomization and of fast simulations, respectively) and eight noise maps, thus 64 maps in total. The error bars therefore refer to the combination of the statistical fluctuation due to the random process and the shape noise uncertainty.

For low peaks with ν ≤ 3.75, we observe that n_peak(ν) remains almost unchanged between the different cases. This is not suprising because in this regime, n_peak(ν) is mainly contributed by noise. This argument is supported by the noise-only peak histogram. The lower panel of Fig. 4 shows that there exist some systematic overcounts in this regime on the order of 10%. The cause of this bias is ambiguous. One possibility might be the use of NFW profiles for ray-tracing simulations. It might also come from the subtraction of the mean κ value from the maps. We leave this to future studies. Another observation in this regime is that by adding the signal to the noise field, the number of peaks with ν ≤ 2.75 decreases. This proves that the effect of noise is not additive for peak counts.

In the regime of ν ≥ 3.75, we observe that replacement by NFW profiles enhances the peak counts, while position randomization introduces an opposite effect of a similar order of magnitude. The enhancement from Case 2 may be explained by the halo triaxiality. A spherical profile such as the NFW model may lead to an overestimation of the projected mass at the center of halos if the major axis is not aligned with the line of sight, and this would probably be the case for most of the N-body halos. It could also be an effect of the M − c relation: we might overestimate c_NFW for large M.

Comparing Cases 2 and 3, we discover that position randomization decreases peak counts by 10% to 50%. Apparently, decorrelating angular positions breaks down the two-halo term, so that halos overlap less on the field of view and decreases high-peak counts. Yang et al. (2011) shows that high peaks with ν ≥ 4.8 are mainly contributed by one single halo, and about 12% of total high-peak counts are contributed by multiple halos. This number agrees with the undercount from our hypothesis of randomization. Combining this step with the former one, we confirm that considering lensing contribution from spatially decorrelated clusters is a good approximation for peak counts.

The impact of the mass function is shown by comparing Case 3 to Case 4. Peak counts are more numerous in our forward model based on the mass function of Jenkins et al. (2001). This excess compensates for the deficit from randomization. However, as shown by Fig. 2, the real mass function in N-body runs is coherent to the analytical model that we use, except for the low-mass deficit tails from N-body runs. To test the impact from this, we ran fast simulations with different lower limit for the halo sampling, and we discover that peak counts do not depend on the lower sampling limit M_min when M_min remains lower than 10¹³M_⊙/h. This proves that the deficit tails are not the cause of the peak count enhancement. Without this explanation, we may have to test with another N-body simulation set to understand the origin of this effect.

Fig. 5 Similar plot to Fig. 4, but in a larger field. Cases 2–4 are carried out for 859 deg2. Case 1 should only be taken as an indication, since its size of field is the same as in Fig. 4, and therefore 16 times smaller than cases 2–4. The fluctuation from high ν bins is much reduced compared to Fig. 4.

Figure 5 shows a similar study of Cases 2–4 for a larger field of 859 deg². One can recover the same effects: compensation of effects deriving from NFW profiles and randomization. Therefore, the difference between our model and N-body simulations is at the same order of magnitude of the one between the analytical and the N-body mass functions. We would like to point out that the Poisson fluctuation has been largely suppressed. A quick calculation shows that, for a given peak density n and a survey area A, the ratio of the Poisson noise to peak density is . The error bars for high peaks in both Fig. 4 and Fig. 5 stay within 50% of the values given by this formula. As a result, we argue that to reduce the Poisson fluctuation at the level of 10%, a survey of more than 150 deg² is preferable using WL peaks with ν ≲ 5.25 and 800 deg² using peaks with ν ≲ 6.25.

5.2. Comparison to an analytical model

In Fig. 6, we draw peak histograms obtained from the analytical model of FSL10 and from our model. The computation for the FSL10 model is done with the same halo profiles and parameters, and the same mass function. For our model, we use our large-field result as mentioned in the previous section. Both models are computed with the same parameter set as the Aardvark N-body simulation inputs. We observe that the FSL10 model is also in good agreement with N-body runs. The prediction from the FSL10 model is more consistent with N-body values for high-peak counts, whereas our model performs better in the low-peak regime. In general, the deviation of both models for ν ≤ 5.25 stays under 25%.

5.3. Sensitivity tests on cosmological parameters

Finally in this section, we show how our model depends on cosmological parameters. Weak lensing is particularly sensitive to Ω_m and σ₈, hence we carry out nine series of fast simulations for which (Ω_m,σ₈) is chosen from , where and are input from our N-body runs. The values of ΔΩ_m and Δσ₈ are chosen to be 0.03 and 0.05, respectively, and the remaining parameters are identical to the N-body simulations. Each scenario is the average over 16 combinations of four fast simulation realizations and four noise maps.

Figure 7 shows four plots that correspond to four variation directions on the Ω_m-σ₈ parameter plane, with regard to . Both upper panels show the variation of only one parameter. They reveal that our model performs a neat, progressive difference of peak abundance in every bin, ranging from ν = 4 to 6. We notice that the differences between cyan diamonds (higher value of Ω_m or σ₈) and red squares (N-body value) are always narrower than those between green circles (lower value of Ω_m or σ₈) and red squares. This is triggered by the banana-shape constraint on the Ω_m-σ₈ plane, from which a horizontal or a vertical cut will result in an asymmetric confidence level for a single parameter.

The two lower panels are variations in the diagonal and anti-diagonal directions. Like what we expect, the diagonal variation is the most efficient discriminant of Ω_m-σ₈. In contrast, peak counts for different parameter sets completely merge together in the lower right-hand panel, since the anti-diagonal direction corresponds roughly to the degenerency lines. Furthermore, all error bars (for 3.75 ≤ ν ≤ 6.25) remain smaller than 5%, which shows the robustness of our model. We recall that blue solid lines correspond to a small 53.7 deg² field, such that the Poisson noise might bias high peak counts, as explained in Sect. 5.1. At the end of the day, the performance of our model at distinguishing different cosmological models has been confirmed.

Fig. 6 Comparison of the FSL model (orange triangles) to our model (cyan diamonds). The full N-body peak histogram is shown as a blue line. Lower panel: we draw the difference between the FSL model and N-body data within an orange dashed line. The cyan-colored zone represents the error bars for our model. The field of view for fast simulations is 859 deg2. The N-body data is only indicative.

Figure 7 also shows that systematic biases of our model could lead to parameter biases. A simple interpolation for the bin ν = 5 shows that N-body peak counts correspond to a cosmology with Ω_m ≈ 0.212 if the knowledge of σ₈ is perfect. The bias is then ΔΩ_m ≈ 0.018. Similarly, the bias on σ₈ is Δσ₈ ≈ 0.030 if Ω_m is known. The origin of the biases of our model is complex. We discuss a list of possible improvements that reduce potential systematics in the following section.

6. Summary and perspectives

WL peaks probe cosmological structures in a straightforward way, since they are directly induced by total-mass gravitational effects, and they especially probe the high-mass part of the mass function. Unlike other tracers, WL peaks provide a forward-fitting approach to study the mass function and cosmology. This makes WL peaks a very competitive candidate for improving our knowledge about structure formation.

In this paper, we presented a new model that predicts weak-lensing peak counts. We generated fast simulations by sampling halos from analytical expressions. By assuming that halos in these simulations are randomly distributed on the sky, we count peaks from ray-tracing maps obtained from these simulations to predict number counts. In this model, we have supposed that unbound matter contributes little to the lensing and that halo clustering has little impact on peak counts.

Fig. 7 Sensitivity tests on (Ωm,σ8). The four plots indicate different variation directions in the Ωm-σ8 plane. The upper left, upper right, lower left, and lower right panels: variation in Ωm, σ8, diagonal, and anti-diagonal direction, respectively. Blue solid lines represent full N-body runs, while red squares are always our model with N-body input. The field of view is 859 deg2. The N-body data is only indicative.

We validated our approach by comparing number counts with N-body results. In particular, we focused on peaks with ν ≈ 4–6, since lower ν are dominated by shape noise, and higher ν are dominated by the Poisson fluctuation. We showed how the three steps corresponding to the main assumptions of our model influence convergence peak abundance. First, NFW profiles tend to shift some medium peaks to higher values, in spite of the lack of unbound objects. Second, the number of peaks decreases when halo positions are randomized. Last, the difference between the N-body mass function and the analytical one is observable in produced peak counts. In summary, our model is in good agreement with results from full N-body runs.

We also tested the dependence of our model on Ω_m and σ₈. For a 859 deg² sky area, the Poisson fluctuation is reduced to a reasonable level for peaks with ν ≲ 6. It turns out that different scenarios are discernable for ν ≳ 4, with a degerency direction corresponding roughly to the anti-diagonal in the plane Ω_m-σ₈. Tests on a large set of different parameters are feasible with our model thanks to the short computation time.

Our probabilistic model has other potential advantages. Repeated simulations for the same cosmological parameters generate the distribution of observables. This allows us to compare observations with our model without the need to define a likelihood function or to assume any Gaussian distribution. For example, model discrimination can be carried out using the false discovery rate method (FDR, Benjamini & Hochberg 1995, an application can be found in Pires et al. 2009), approximate Bayesian computation (ABC, see for example Cameron & Pettitt 2012; Weyant et al. 2013), or other statistical techniques. Another powerful advantage of our model is its flexibility. Additional effects such as intrinsic ellipticity alignment, alternative methods such as nonlinear filters, and realistic survey settings, such as mask effects, magnification bias (Liu et al. 2014), shape measurement errors (Bard et al. 2013), and photo-z errors, can all be modeled in this peak-counting framework. The forward-modeling approach allows for a straightforward inclusion and marginalization of model uncertainties and systematics.

Several improvements to our model are possible. Using perturbation theory, we may take halo clustering into account in fast simulations. This can be done by some fast algorithms, such as PTHalos (Scoccimarro & Sheth 2002), Pinocchio (Monaco et al. 2002; see also Heisenberg et al. 2011), and remapping LPT (Leclercq et al. 2013). In addition, we can go beyond the idealized setting considered in this work by including a realistic source distribution, intrinsic alignment, mask effects, etc. We also expect that nonlinear filters and tomography studies may bring some more refined results for cosmology from peak counting. Finally, peak counts can be supplemented with additional WL observables, such as magnification and flexion.

The Camelus algorithm is implemented in C language. It requires the Nicaea library for cosmological computations. The Camelus source code is released via the website⁵.

Acknowledgments

This study is supported by Région d’Île-de-France under grant DIM-ACAV and the French national program for cosmology and galaxies (PNCG). The authors acknowledge the anonymous referee for useful comments and suggestions. We would like to thank August Evrard for providing N-body simulations. We also thank Zuhui Fan, Xiangkun Liu, and Chuzhong Pan for giving constructive comments on the preprint. Chieh-An Lin is very grateful for inspiring discussions with François Lanusse, Yohan Dubois, and Michael Vespe.

References

All Figures

Fig. 1 Forward- and reverse-modeling diagram for the mass function studies. Two different approaches to establishing links between the theoretical mass function (the upper round rectangle) and the observables (the lower round rectangle) are to compare the observable mass function nobs with the theoretical one (reverse modeling), or to compare “predicted” observable values Xmodel with observations (forward modeling). In this paper, the forward modeling is adopted and we propose a new method to “predict” peak conuts.

In the text

	Fig. 2 Comparison between an analytical mass function (blue line), halo mass histograms for N-body runs (green points), and sample histograms for realizing of the fast simulation (red crosses). The plots are drawn at 4 different redshift planes, and for each the thickness is dz = 0.1.
In the text

Fig. 3 Patch of a map projected with regard to its center, taken from a realization of fast simulations. The left, middle, and right panels give the fields κ1 / 3, , and , respectively. We have taken the cubic root to emphasize the contrast. It is clear that the signal is completely dominated by the noise. Even though the smoothed map is quite different from the original one, the structures, high-signal regions are still conserved and traced.

In the text

Fig. 4 Comparison of the peak abundance from different cases. Upper panel: blue solid line: full N-body runs (Case 1); green circles: replacement of halos by NFW profiles (Case 2); red squares: replacement of halos by NFW profiles and randomization of halo angular positions (Case 3); cyan diamonds: fast simulations, corresponding to our model (Case 4); magenta dashed line: peaks from noise-only maps. Lower panel: we draw the upper and lower limits of error bars shifted with regard to the N-body values. This refers to the standard deviation over 4 maps (green dash-dotted line for Case 2) or 16 maps (red dashed line for Case 3, cyan solid line for Case 4). The field of view is 53.7 deg2.

In the text

	Fig. 5 Similar plot to Fig. 4, but in a larger field. Cases 2–4 are carried out for 859 deg2. Case 1 should only be taken as an indication, since its size of field is the same as in Fig. 4, and therefore 16 times smaller than cases 2–4. The fluctuation from high ν bins is much reduced compared to Fig. 4.
In the text

	Fig. 6 Comparison of the FSL model (orange triangles) to our model (cyan diamonds). The full N-body peak histogram is shown as a blue line. Lower panel: we draw the difference between the FSL model and N-body data within an orange dashed line. The cyan-colored zone represents the error bars for our model. The field of view for fast simulations is 859 deg2. The N-body data is only indicative.
In the text

Fig. 7 Sensitivity tests on (Ωm,σ8). The four plots indicate different variation directions in the Ωm-σ8 plane. The upper left, upper right, lower left, and lower right panels: variation in Ωm, σ8, diagonal, and anti-diagonal direction, respectively. Blue solid lines represent full N-body runs, while red squares are always our model with N-body input. The field of view is 859 deg2. The N-body data is only indicative.

In the text