<figure><img src="image_1.png" /></figure>

# Cosmological constraints from weak lensing peak statistics with Canada-France-Hawaii Telescope Stripe 82 Survey

Xiangkun Liu, Chuzhong Pan, Ran Li, Huanyuan Shan, Qiao Wang, Liping Fu, Zuhui Fan, Jean-Paul Kneib, Alexie Leauthaud, Ludovic van Waerbeke, et al.

► To cite this version:

Xiangkun Liu, Chuzhong Pan, Ran Li, Huanyuan Shan, Qiao Wang, et al. Cosmological constraints from weak lensing peak statistics with Canada-France-Hawaii Telescope Stripe 82 Survey. Monthly Notices of the Royal Astronomical Society, 2015, 450 (3), pp.2888-2902. 10.1093/mnras/stv784. hal-01439964

HAL Id: hal-01439964

https://hal.science/hal-01439964

Submitted on 15 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

# Cosmological constraints from weak lensing peak statistics with Canada-France-Hawaii Telescope Stripe 82 Survey

Xiangkun Liu,<sup>1★</sup> Chuzhong Pan,<sup>1</sup> Ran Li,<sup>2</sup> Huanyuan Shan,<sup>3</sup> Qiao Wang,<sup>2</sup>

Liping Fu,<sup>4</sup> Zuhui Fan,<sup>1,5★</sup> Jean-Paul Kneib,<sup>3,6</sup> Alexie Leauthaud,<sup>7</sup>

Ludovic Van Waerbeke,<sup>8</sup> Martin Makler,<sup>9</sup> Bruno Moraes,<sup>10,11</sup> Thomas Erben<sup>12</sup>
and Aldée Charbonnier<sup>13,14</sup>

Accepted 2015 April 3. Received 2015 March 6; in original form 2014 December 18

## ABSTRACT

We derived constraints on cosmological parameters using weak lensing peak statistics measured on the ~130 deg<fnref n="2" /> of the Canada-France-Hawaii Telescope Stripe 82 Survey. This analysis demonstrates the feasibility of using peak statistics in cosmological studies. For our measurements, we considered peaks with signal-to-noise ratio in the range of v = [3, 6]. For a flat Λ cold dark matter model with only (Ω_m, σ_8) as free parameters, we constrained the parameters of the following relation Σ_8 = σ_8(Ω_m/0.27)α to be Σ_8 = 0.82 ± 0.03 and α = 0.43 ± 0.02. The α value found is considerably smaller than the one measured in two-point and three-point cosmic shear correlation analyses, showing a significant complement of peak statistics to standard weak lensing cosmological studies. The derived constraints on (Ω_m, σ_8) are fully consistent with the ones from either WMAP9 or Planck. From the weak lensing peak abundances alone, we obtained marginalized mean values of Ω_m = 0.38<sup>+0.27</sup><sub>-0.24</sub> and σ_8 = 0.81 ± 0.26. Finally, we also explored the potential of using weak lensing peak statistics to constrain the mass-concentration relation of dark matter haloes simultaneously with cosmological parameters.

**Key words:** gravitational lensing: weak–dark matter–large-scale structure of Universe.

## 1 INTRODUCTION

Large-scale structures in the Universe perturb the propagation of light rays from background sources causing small shape distortions and luminosity changes for their observed images (e.g. Bartelmann

& Schneider 2001). Such effects, namely weak lensing effects, are closely related to the formation and evolution of foreground structures and the global expansion history of the Universe, and therefore are known to be one of the most promising probes in cosmological studies (e.g. Albrecht et al. 2006; LSST Dark Energy Science Collaboration 2012; Amendola et al. 2013; Weinberg et al. 2013). The cosmic shear two-point (2-pt) correlation analysis has been demonstrated to be a powerful statistics in extracting weak

<footnote n="1"><sup>1</sup>Department of Astronomy, Peking University, Beijing 100871, China</footnote>

<footnote n="2"><sup>2</sup>Key Laboratory for Computational Astrophysics, The Partner Group of Max Planck Institute for Astrophysics, National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100012, China</footnote>

<footnote n="3"><sup>3</sup>Laboratoire d'astrophysique (LASTRO), Ecole Polytechnique Fédérale de Lausanne (EPFL), Observatoire de Sauverny, CH-1290 Versoix, Switzerland</footnote>

<footnote n="4"><sup>4</sup>Shanghai Key Lab for Astrophysics, Shanghai Normal University, Shanghai 200234, China</footnote>

<footnote n="5"><sup>5</sup>Collaborative Innovation Centre of Modern Astronomy and Space Exploration, Nanjing 210093, China</footnote>

<footnote n="6"><sup>6</sup>Aix Marseille Université, CNRS, LAM (Laboratoire d'Astrophysique de Marseille) UMR 7326, F-13388 Marseille, France</footnote>

<footnote n="7"><sup>7</sup>Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU, WPI), The University of Tokyo, Kashiwa, Chiba 277-8582, Japan</footnote>

<footnote n="8"><sup>8</sup>Department of Physics and Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, BC V6T 1Z1, Canada</footnote>

<footnote n="9"><sup>9</sup>Centro Brasileiro de Pesquisas Físicas, Rua Dr. Xavier Sigaud 150, Rio de Janeiro, RJ 22290-180, Brazil</footnote>

<footnote n="10"><sup>10</sup>Department of Physics and Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>

<footnote n="11"><sup>11</sup>CAPES Foundation, Ministry of Education of Brazil, Brasilia/DF 70040-020, Brazil</footnote>

<footnote n="12"><sup>12</sup>Argelander Institute for Astronomy, University of Bonn, Auf dem Hügel 71, D-53121 Bonn, Germany</footnote>

<footnote n="13"><sup>13</sup>Observatório do Valongo, Universidade Federal do Rio de Janeiro, Ladeira do Pedro Antônio 43, Saúde, Rio de Janeiro, RJ 20080-090, Brazil</footnote>

<footnote n="14"><sup>14</sup>Centro Brasileiro de Pesquisas Físicas, Rua Dr. Xavier Sigaud 150, Rio de Janeiro, RJ 22290-180, Brazil</footnote>

<footnote>* E-mail: lxlk98479@pku.edu.cn (LX); fanzuhui@pku.edu.cn (ZF)</footnote>
lensing signals from shape measurements of background galaxies (e.g. Fu et al. 2008; Heymans et al. 2012; Jee et al. 2013; Kilbinger et al. 2013; Kitching et al. 2014). On the other hand, 2-pt correlations can only reveal part of the cosmological information embedded in weak lensing signals given the non-linearity of the structure formation. To overcome this limitation, higher order cosmic shear correlation analyses are a natural extension (e.g. Semboloni et al. 2011; van Waerbeke et al. 2013; Fu et al. 2014). Weak lensing peak statistics, i.e. concentrating on high signal regions, is another way to probe efficiently the non-linear regime of the structure formation, and thus can provide important complements to the cosmic shear 2-pt correlation analysis (e.g. White, van Waerbeke & Mackey 2002; Hamana, Takada & Yoshida 2004; Hennawi & Spergel 2005; Tang & Fan 2005; Dietrich & Hartlap 2010; Kratochvil, Haiman & May 2010; Yang et al. 2011; Marian et al. 2012; Hilbert et al. 2012; Bard et al. 2013; Lin & Kilbinger 2014)

Observationally, different analyses have proved the feasibility of performing weak lensing peak searches from data (e.g. Wittman et al. 2006; Gavazzi & Soucail 2007; Miyazaki et al. 2007; Geller et al. 2010). However, up to now, few cosmological constraints are derived from weak lensing peak statistics in real observations. There are two main reasons for lack of such analyses. First, weak lensing observations are just starting to reach significantly large survey areas to provide reasonable statistics for peak abundances (e.g. Shan et al. 2012, 2014; van Waerbeke et al. 2013). The second reason is the theoretical difficulty to calculate the cosmology dependence of peak abundances.

Theoretically, high signal-to-noise (S/N) weak lensing signal regions are expected to be associated closely with massive structures along lines of sight (e.g. White et al. 2002). Therefore, in principle, weak lensing peak abundances should reflect the underlying mass function of dark matter haloes weighted by the lensing efficiency kernel (e.g. Hamana et al. 2004). In practice, however, the correspondence between weak lensing peaks and the massive dark matter haloes is influenced by various effects, such as the noise from the intrinsic ellipticities of source galaxies, the projection effect of large-scale structures, and the hierarchical mass distribution of dark matter haloes (e.g. van Waerbeke 2000; Tang & Fan 2005; Hamana et al. 2012; Yang et al. 2013). Thus, it is not straightforward to predict the cosmology dependence of weak lensing peak abundances. One possible solution is to create a large number of simulation templates for weak lensing peak statistics densely sampled in cosmological-parameter space. By comparing the observational measurements with the templates, we can derive cosmological constraints (e.g. Dietrich & Hartlap 2010; Liu et al. 2015, hereafter LPH2015). Considering the large number of cosmological parameters and different physical and observational effects, such an approach can be numerically expensive. Another efficient way is to build theoretical models, which take into account the impact of the different effects.

Based on simulation studies, Marian, Smith & Bernstein (2009, 2010) developed a phenomenological model for hierarchically detected weak lensing peak abundances in which the 2D peak mass function is scaled to the 3D mass function of dark matter haloes. Hamana et al. (2004, 2012) derived a fitting formula for weak lensing peak abundances by incorporating a probability function in relating peak heights and underlying dark matter haloes at a given mass and redshift. Calibrated with numerical simulations, such a probability function tends to include the effects of noise from intrinsic ellipticities of source galaxies, the projection effects of large-scale structures and the non-spherical matter distributions of dark matter haloes. Assuming Gaussian random fields for both

the projected field of large-scale structures and the shape noise, Maturi et al. (2010) proposed a theoretical model to calculate the number of contiguous areas above a given threshold in the filtered convergence field. This is equivalent to the genus in Minkowski functionals. When the threshold is high, this statistics corresponds well to the number of peaks. By comparing with simulations, it was shown that the model can predict well the number distribution for relatively low thresholds, but underestimates the high threshold regions that are mostly related to individual massive haloes (Maturi et al. 2010; Petri et al. 2013).

In Fan, Shan & Liu (2010, hereafter F10), we have presented a theoretical model taking into account the shape noise effects. In this model, we divide a given area into halo regions occupied by dark matter haloes with the size limited by their virial radii, and the regions outside dark matter haloes. We first calculate the weak lensing peak abundances in a halo region by assuming a density profile for the halo and the Gaussianity of the shape noise. By employing the mass function of dark matter haloes with a lower mass cut representing the halo mass above which single haloes dominantly contribute to weak lensing peaks along their lines of sight, we can then calculate statistically the peak distribution in regions occupied by massive haloes. For the rest of the regions, we assume that the peaks are purely noise peaks. Our model has been tested extensively by comparing with numerical simulations (F10; Liu et al. 2014, hereafter LWPF2014). It was shown that the model results are in very good agreement with simulation results. It is noted that the model in its present form does not contain the projection effects of large-scale structures. For current generation of weak lensing surveys with the surface number density of lensing-usable galaxies around $n_g \sim 10 \text{ arcmin}^{-2}$, the shape noise is dominant over the projection effects for the smoothing scale $\sim 1$ arcmin. Thus neglecting the projection effects should not affect the model prediction significantly. For future surveys with much improved statistics, the projection effects need to be considered carefully, and we have started to look into this problem.

For the Canada-France-Hawaii Telescope (CFHT) Stripe 82 Survey (CS82; e.g. Li et al. 2014; Shan et al. 2014; Hand et al. 2015), the number density of galaxies used in weak lensing studies is $\sim 10 \text{ arcmin}^{-2}$, and the survey area excluding the masked regions is $\sim 130 \text{ deg}^2$. For this survey, the shape noise is the dominant source of contaminations on weak lensing peak analyses. We thus expect that our model can work well in predicting theoretically the peak abundances. This in turn allows us to perform cosmological constraints from observational weak lensing peak abundances.

As this work was being completed, we became aware of the study by LPH2015. They also analysed the cosmological application of weak lensing peak statistics using CFHTLenS data. Their studies are based on interpolations from a suite of simulation templates on a grid of 91 cosmological models in the parameter space of ($\Omega_m, \sigma_8, w$), where $\Omega_m, \sigma_8$ and $w$ are, respectively, the dimensionless matter density of the Universe, the amplitude of the extrapolated linear matter density fluctuations smoothed over a top-hat scale of $8h^{-1}\text{Mpc}$, and the equation of state of dark energy. Using different and independent approaches, LPH2015 and our work both showed the promising potential of weak lensing peak statistics in cosmological studies.

The paper is organized as follows. In Section 2, we describe briefly the CS82 survey. In Section 3, we present the procedures of weak lensing peak analyses. In Section 4, we show the cosmological constraints derived from peak abundances. Summary and discussion are given in Section 5.
## 2 CFHT STRIPE 82 SURVEY AND WEAK LENSING CATALOGUES

The CS82 survey was conducted thanks to the collaboration between the Canadian, French and Brazilian CFHT communities. CS82 covers a large fraction of the Sloan Digital Sky Survey (SDSS) Stripe 82 with high-quality i-band imaging under excellent seeing conditions in the range of 0.4–0.8 arcsec with an average of 0.59 arcsec. The survey contains a total of 173 tiles, 165 of which from CS82 observations and 8 from CFHT-LS Wide (Erben et al. 2013). Each CS82 tile was obtained from four consecutive dithered observations each with an exposure time of 410 s. The derived $5\sigma$ limiting magnitude in a 2 arcsec diameter aperture is $i_{AB} \sim 24$. After removing overlapping regions and applying all the masks across the entire survey, the effective survey area is reduced from ~173 to ~130 deg$^2$.

The same forward modelling LENSFIT pipeline (Miller et al. 2007, 2013) as that for CFHTLenS was used for the CS82 shape measurements. As described in Miller et al. (2013), the LENSFIT algorithm applied to CFHTLenS was calibrated using different sets of simulated images with different observing conditions and point spread functions (PSFs). For the shape measurement errors written in the form of $\epsilon = (1+m)\epsilon_{true} + c$, it is found that the multiplicative bias factor $m$ can be well modelled as a function of galaxy S/N ratio and size. For the additive bias $c$, the simulation calibration shows that it is consistent with zero. However for real data, the additive bias can occur. Heymans et al. (2012) show that for the LENSFIT measurements of CFHTLenS, the $c_1$ component is consistent with zero, but there are small residues for $c_2$. Similarly to $m$, the $c_2$ term also depends on galaxy S/N ratio and size. It is noted that CFHTLenS has a wide range of seeing, number of exposures, noise and depth, and their influences on shape measurements can all be encoded into the two parameters of galaxy S/N ratio and size. Therefore for CS82, although the observing conditions are different from CFHTLenS, the LENSFIT pipeline is well applicable. We should emphasize that although we use the same LENSFIT pipeline as CFHTLenS, we measure the PSFs and calculate the S/N ratio and size and the corresponding bias terms $m$ and $c$ for CS82 source galaxies ourselves.

In our weak lensing analyses, the selection criteria of source galaxies are weight $w > 0$, FITCLASS = 0, MASK $\le 1$. Here the weight factor is the inverse variance weight accorded to each source galaxy given by LENSFIT. The FITCLASS is an index for star/galaxy classification provided by LENSFIT with FITCLASS = 0 for galaxies. The index MASK describes the mask information at an object's position. Objects with MASK $\le 1$ can safely be used for most weak lensing analyses (Erben et al. 2013). No magnitude cut is applied for the catalogue as fainter galaxies have lower weights. These criteria result in a total number of source galaxies of 9281 681. The total effective number of galaxies taking into account their weights is 5475 318, and the corresponding average effective number density is ~11.8 galaxies arcmin$^{-2}$.

Not all the source galaxies have redshift information. In our theoretical calculations, we therefore adopt a redshift distribution derived for the whole population of source galaxies, which is obtained by magnitude matching of Cosmic Evolution Survey (COSMOS) galaxies with the CS82 source galaxies (Shan et al. 2014; Hand et al. 2015). It is given by

$$p_z(z) \propto \frac{z^a + z^{ab}}{z^b + c}, \qquad (1)$$

where $a = 0.531$, $b = 7.810$ and $c = 0.517$. The median redshift is $z_m = 0.76$ and the mean redshift is $z = 0.83$. The normalized redshift distribution is shown in Fig. 1.

<figure><img src="image_4.png" /><figcaption>Figure 1. Normalized redshift distribution of CS82 source galaxies.</figcaption></figure>

Because the COSMOS field is small, the sample variance can be significant. There are also errors in the photometric redshift estimations for COSMOS galaxies. Thus the CS82 redshift distribution derived from COSMOS can have uncertainties. We will discuss the impact of such uncertainties on peak analyses in Section 4.3.

## 3 WEAK LENSING PEAK ANALYSIS

### 3.1 Theoretical aspects

In the weak lensing theory under the framework of general relativity, the deflection of light rays from a source can be written as the gradient of a lensing potential $\psi(\theta)$. The induced observational effects can be described by the Jacobian matrix $A$, which is given by (e.g. Bartelmann & Schneider 2001)

$$A = \left( \delta_{ij} - \frac{\partial^2 \psi(\theta)}{\partial \theta_i \partial \theta_j} \right) = \begin{pmatrix} 1 - \kappa - \gamma_1 & -\gamma_2 \\ -\gamma_2 & 1 - \kappa + \gamma_1 \end{pmatrix}, \quad (2)$$

where the convergence $\kappa$ and the shear $\gamma_i$ lead to the isotropic change and the elliptical shape distortion of the observed image, respectively, with respect to the un lensed image. They are related to the potential by

$$\kappa = \frac{1}{2} \nabla^2 \psi, \quad \gamma_1 = \frac{1}{2} \left( \frac{\partial^2 \psi}{\partial \theta_1^2} - \frac{\partial^2 \psi}{\partial \theta_2^2} \right), \quad \gamma_2 = \frac{\partial^2 \psi}{\partial \theta_1 \partial \theta_2}. \qquad (3)$$

In the weak lensing regime under the Born approximation, we have

$$\kappa = \frac{3 H_0^2 \Omega_m}{2c^2} \int_0^\chi d\chi' \frac{f_K(\chi') f_K(\chi - \chi') \delta[f_K(\chi')\theta, \chi']}{f_K(\chi) a(\chi')}, \qquad (4)$$

where $H_0$ is the Hubble constant, $\chi$ is the comoving radial distance, $f_K$ is the comoving angular-diameter distance, $a$ is the scale factor of the universe and $\delta$ is the density perturbation along the line of sight.

The convergence $\kappa$ is directly related to the projection of line-of-sight density fluctuations weighted by the lensing efficiency factor. Physically, massive structures generate large weak lensing signals along their lines of sight. These peak signals are best seen visually in the weak lensing convergence field. On the other hand, weak lensing signals directly extracted from galaxy shape measurements are the shears, or more precisely the reduced shears defined as $g_i = \gamma_i/(1-\kappa)$, rather than the convergence. Therefore, weak lensing peak analysis usually involves procedures to construct quantities representing the projected mass distribution from lensing shears
based on the relation between the two quantities shown in equation (3).

Specifically, assuming precise shape measurements, the observed ellipticity of a galaxy at redshift z located at the sky position θ can be written in a complex form given by (e.g. Seitz & Schneider 1997)

$$ \epsilon(\theta, z) = \begin{cases} \frac{\epsilon_s(\theta, z) + g(\theta, z)}{1 + g^*(\theta, z)\epsilon_s(\theta, z)} & \text{for } |g(\theta, z)| \le 1 \\ \frac{1 + g(\theta, z)\epsilon_s^*(\theta, z)}{\epsilon_s^*(\theta, z) + g^*(\theta, z)} & \text{for } |g(\theta, z)| > 1 \end{cases} \quad (5) $$

where '*' represents the complex conjugate operation, g = g₁ + ig₂ is the complex reduced shear, and εs is the intrinsic ellipticity of the galaxy. Here the complex ellipticity is defined as ε = (a - b)/(a + b) exp(2iφ) with a, b and φ being the length of the major and minor axes and the orientation of the approximate ellipse of the observed image, respectively. For source galaxies at a fixed redshift z, it has been shown that the average of ε over a large number of galaxies near θ gives rise to an unbiased estimate of g(θ, z) or 1/g(θ, z) assuming εs = 0. For galaxies with a redshift distribution, the average of ε over galaxies near a given sky position may have a complicated relation with the lensing signal we are interested in if both |g| < 1 and |g| > 1 occur in the region for galaxies at different redshifts. On the other hand, for subcritical regions with |g| < 1 for all the redshifts, the average of the observed ellipticity 〈ε〉 gives rise to an estimate of 〈g〉 weighted by the redshift distribution of source galaxies. In the case κ ≪ 1, we have 〈ε〉 ≈ 〈γ〉. For the weak lensing peak analysis, we therefore need to construct a field closely related to the matter distribution from the (reduced) shear estimate 〈ε〉.

The aperture-mass peak analysis, also referred to as the shear peak analysis, is to study peaks in the aperture mass M<sub>ap</sub> field constructed from the tangential shear component with respect to the point of interest with a filtering function Q (e.g. Schneider et al. 1998; Marian et al. 2012; Bard et al. 2013). Theoretically, M<sub>ap</sub> corresponds to the convergence field filtered with a compensated window function U where U and Q are related. One of the advantages of M<sub>ap</sub> studies is that because of the compensated nature of U, M<sub>ap</sub> is independent of the lensing mass-sheet degeneracy. Furthermore, in the case with κ ≪ 1 and g ≈ γ, M<sub>ap</sub> can be obtained directly from the tangential component of the observed 〈εt〉. In the peak regions where the lensing signals are high, the difference between g and γ is not negligible. Therefore, noting that the M<sub>ap</sub> constructed from 〈εt〉 itself carries cosmological information, it is not the same as the filtered convergence field.

Another approach for weak lensing peak studies is to reconstruct the convergence field from 〈ε〉 taking into account the non-linear relation between g and γ (e.g. Kaiser & Squires 1993; Bartelmann 1995; Kaiser, Squires & Broadhurst 1995; Seitz & Schneider 1995, 1997; Squires & Kaiser 1996; Jauzac, Jullo & Kneib 2012; Jullo et al. 2014). In this approach, to avoid unphysical results, it is important to filter ε first and then to proceed with convergence reconstruction using the filtered 〈ε〉. Different reconstruction schemes have been studied. For the classical Kaiser–Squires (KS) reconstruction (Kaiser & Squires 1993; Squire & Kaiser 1996), boundary effects and the mass-sheet degeneracy problem can exist. However, for a field of view of about 1 deg × 1 deg and larger, such effects are expected to be insignificant. In this paper, we reconstruct the convergence field from the filtered 〈ε〉 with the non-linear KS method (e.g. Bartelmann 1995). From numerical simulations, we find that the regions with |g| > 1 are negligible, and therefore assuming subcriticality for all the regions is an excellent approximation.

## 3.2 The convergence reconstruction and the peak identification

Our convergence reconstruction procedures are described below.

For a source galaxy used in weak lensing analyses in the CS82 catalogue, we first correct the additive errors

$$ \epsilon_1^c = \epsilon_1, \quad \epsilon_2^c = \epsilon_2 - c_2, \quad (6) $$

where εi and εc are the uncorrected and corrected ellipticity components, respectively, and c2 is the additive bias given by CS82. With εc, we apply smoothing and obtain a smoothed field of ε on regular 1024 × 1024 grids over the field of view of one pointing. With the multiplicative errors m taken into account statistically (van Waerbeke et al. 2013), we have

$$ \langle \epsilon \rangle (\theta) = \frac{\sum_j W_{\theta_G}(\theta_j - \theta) w(\theta_j) \epsilon^c(\theta_j)}{\sum_j W_{\theta_G}(\theta_j - \theta) w((\theta_j)(1+m_j))}, \quad (7) $$

where θ and θj are for the grid position and the galaxy position, respectively, W<sub>θG</sub> is the normalized smoothing function, w is the weight for source galaxy shape measurements given by CS82. The summation is over all the source galaxies. In the subcritical approximation, 〈ε〉(θ) is an unbiased estimate of 〈g〉 smoothed over the window function W<sub>θG</sub> and weighted by the source redshift distribution. We use the Gaussian smoothing function W<sub>θG</sub> with

$$ W_{\theta_G}(\theta) = \frac{1}{\pi \theta_G^2} \exp \left( - \frac{|\theta|^2}{\theta_G^2} \right). \quad (8) $$

The smoothing scale θG is chosen to be θG = 1.5 arcmin, suitable for cluster-scale structures that are closely related to high weak lensing peaks. Within the smoothing kernel, the number of galaxies is about n<sub>G</sub>θ<sub>G</sub><sup>2</sup> ~ 20, for n<sub>G</sub> ~ 10 arcmin<sup>-2</sup>. We expect that the statistics of the residual shape noise after smoothing is approximately Gaussian from the central limit theorem (e.g. van Waerbeke 2000).

With 〈ε〉(θ), we perform the convergence reconstruction iteratively by using the relation between κ and γ in equation (3). Particularly, we use their relation in Fourier space with

$$ \hat{\gamma}(k) = \pi^{-1} \hat{D}(k) \hat{\kappa}(k), \quad (9) $$

and

$$ \hat{D}(k) = \pi \frac{k_1^2 - k_2^2 + 2ik_1k_2}{k_1^2 + k_2^2}. \quad (10) $$

We start by assuming κ⁽⁰⁾ = 0 everywhere, and thus γ⁽⁰⁾ = 〈ε〉 (Bartelmann 1995). At n-th step, we obtain κ⁽ⁿ⁾ from γ⁽ⁿ⁻¹⁾ via equation (9) and the subsequent inverse Fourier transformation. We then update γ to γ⁽ⁿ⁾ = (1 − κ⁽ⁿ⁾)〈ε〉 for next iteration. The reconstruction process is stopped when the converging accuracy of 10⁻⁶ (the maximum difference of the reconstructed convergence between the two sequential iterations) is reached. For CS82, the reconstruction is done pointing by pointing each with the field of view of about 1 deg × 1 deg.

To evaluate the shape noise level in each pointing for subsequent peak analyses, we randomly rotate the corrected ellipticity of each galaxy. Then the same procedure is applied to obtain the reconstructed random noise convergence field.

It is noted that there are regions with no reliable shape measurements for galaxies indicated with the index MASK > 1 in the CS82 catalogue. These galaxies are excluded in the weak lensing analyses. The existence of these masked regions can affect the weak lensing peak abundances significantly if they are not treated properly. In LWPF2014, we study in detail the mask effects. With the S/N ratio defined by the average noise level, the number of high
peaks increases around the masked regions, which can lead to considerable bias in cosmological-parameter constraints derived from weak lensing peak abundances. To reduce the mask effects, regions around masks should be excluded in peak counting. If they are kept, the noise effects in these regions should be considered separately from the regions away from masks (LWPF2014).

To quantify the mask effects on the number of usable galaxies in the convergence reconstruction, we calculate the galaxy filling factor at each grid point similar to that done in van Waerbeke et al. (2013). Summing over galaxies outside masked regions, we define the galaxy filling factor as

$$f(\theta) = \frac{\sum_j W_{\theta_G}(\theta_j - \theta)w(\theta_j)}{f_0}, \qquad (11)$$

where $f_0$ is calculated by randomly populating galaxies over the full area of a tile with

$$f_0 = \langle \sum_n W_{\theta_G}(\theta_n - \theta) \tilde{w}(\theta_n) \rangle. \qquad (12)$$

Here $\langle \rangle$ is for the average over $\theta$. Specifically, for each tile, we calculate the average number density of galaxies in the area excluding the masked regions. With this number density, we then randomly populate galaxies over the full area of the tile including the masked regions. For each galaxy, we randomly assign it a weight $\tilde{w}$ according to the weight distribution of the real observed galaxies. We then calculate the quantity $\sum_n W_{\theta_G}(\theta_n - \theta) \tilde{w}(\theta_n)$ at each grid point $\theta$ where the summation is over all the populated galaxies. The average value over all the grid points gives rise to $f_0$.

Hence, for each pointing, we obtain the reconstructed lensing convergence, noise and filling factor maps, respectively. We have a total of 173 sets of maps corresponding to the 173 pointings. In Fig. 2, we show an example of the reconstructed convergence map for one pointing and the corresponding map of the filling factor. In each map, the dark blue regions are regions with the filling factor $f < 0.5$. To avoid the mask effects on weak lensing peak analyses, we exclude these regions in peak counting (van Waerbeke et al. 2013; LWPF2014). The black circles in the plots show the clusters in the field detected using the red sequence Matched-filter Probabilistic Percolation (redMaPPer) algorithm (Rykoff et al. 2014). It is seen that the clusters have a good association with weak lensing convergence peaks, but the correspondence is not one to one due

to the existence of noise and the projection effects of large-scale structures (e.g. Shan et al. 2014).

For weak lensing peak analyses, we detect peaks from the reconstructed lensing convergence maps as follows. Considering a pixel on a reconstructed convergence map (1024 × 1024 pixels), if its convergence value is the highest among its nearest eight neighbouring pixels, it is identified as a peak. We only count peaks in regions with the filling factor $f > 0.5$. To reduce the boundary effects, we also exclude the outer most 50 pixels (corresponding to ~3 arcmin ~ 2θ_G) in each of the four sides of a map in our peak counting. The effective area of CS82 for weak lensing peak studies is then reduced to ~114 deg² after mask-region and boundary exclusions. The S/N ratio of a peak is defined by

$$\nu = \frac{K}{\sigma_0}, \qquad (13)$$

where K is the reconstructed convergence value of the peak and $\sigma_0$ is the mean rms of the noise from the 173 noise maps. It is known that $\sigma_0$ depends on the number density of source galaxies and the smoothing scale of the window function used in obtaining the smoothed ellipticity field $\langle \epsilon \rangle$. It can vary somewhat from one pointing to another. In our study, the mean $\sigma_0$ is evaluated from all the noise maps considering only regions with the filling factor $f > 0.5$. For the smoothing scale $\theta_G = 1.5$ arcmin, we have $\sigma_0 \approx 0.022$. It is noted that we calculate $\sigma_0$ directly from the rotated galaxies in noise maps, and therefore we do not need to know explicitly the galaxy intrinsic ellipticity dispersion $\sigma_\epsilon$. On the other hand, we find that for CS82 galaxies, $\sigma_\epsilon \sim 0.4$ for the total of the two components.

We note that in our peak identification scheme, we do not group peaks together as some other methods do (e.g. Hamana et al. 2012). We will see in Section 3.3 that our theoretical model for peak abundances takes into account the noise peaks in halo regions and counts them as independent ones. Accordingly, we therefore do not need the peak grouping that may give rise to some artificial effects.

### 3.3 Theoretical model for weak lensing peak abundances

To derive cosmological constraints from observed weak lensing peak abundances, their dependence on cosmological models needs to be understood and quantified. Known to be closely related to

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Left-hand panel: the reconstructed convergence map for one specific tile with the Gaussian smoothing scale $\theta_G = 1.5$ arcmin. Regions with filling factor <0.5 are masked out in dark blue. Right-hand panel: the corresponding filling factor map. The redMaPPer clusters in the field are indicated by the black circles with the size indicating the richness of the clusters.</figcaption></figure>
line-of-sight matter concentrations, the existence of galaxy shape noise, the projection effects of large-scale structures, the complex mass distribution of dark matter haloes, etc., complicates the relation between weak lensing peak abundances and the underlying mass function of dark matter haloes. While building a large set of templates from numerical simulations densely sampled over the multidimensional cosmological-parameter space can be very useful and important, it is computationally expensive. On the other hand, theoretical modelling based on our physical understandings can be very insightful and valuable in disentangling different effects on weak lensing peak abundances. With explicit dependences on cosmological parameters and other physical parameters, it can be used to perform cosmological constraints efficiently. Depending on the assumptions and approximations employed in the modelling, their results can be less accurate than those from full simulations. Similarly to the calculation of the non-linear matter power spectrum based on the halo model but calibrated with simulations, which is widely used in weak lensing 2-pt correlation analyses, the combination of the two, that is, testing and calibrating a model with simulations, can be a very effective way to study the cosmological dependence of weak lensing peak statistics efficiently with high precision.

In this paper, we focus mainly on high peaks and adopt the model of F10 for the weak lensing peak abundances, which takes into account the dominant shape noise in the modelling. The model has been tested extensively by comparing with full ray-tracing simulations (F10; LWPF2014), and has also be confronted with observational studies (Shan et al. 2012, 2014).

Here we describe the important ingredients of the model. More details can be found in F10 and LWPF2014.

In F10, we assume that the smoothed convergence field can be written as $K_N = K + N$, where K represents the true lensing convergence and N is for the residual shape noise. The field N results from the contribution of the intrinsic ellipticities of different galaxies in the smoothing kernel. Without considering the intrinsic alignments of source galaxies, if the number of galaxies within the smoothing kernel is large enough, it has been shown that N is approximately a Gaussian random field from the central limit theorem (e.g. van Waerbeke 2000). As discussed in Section 3.2, for a smoothing scale $\theta_G = 1.5$ arcmin, the number of galaxies within the smoothing window in CS82 is ~20. Therefore, N can be well approximated as a Gaussian random field.

Concentrating on high peaks, it is expected that signals of true peaks mainly come from individual massive haloes (e.g. Hamana et al. 2004; Tang & Fan 2005; Yang et al. 2011). We therefore divide a given area into halo regions and field regions. Inside the region of an individual halo, we have $K_N = K + N$, where K is regarded as a known field from the halo convergence, and N is a Gaussian random field. Therefore $K_N$ itself is also a Gaussian random field modulated by the halo surface mass distribution K. Then the peak number distribution for $K_N$ is readily calculable using Gaussian statistics. The modulation effects from K involve K itself, its first derivatives $K^i = \partial K / \partial x_i$, ($i = 1, 2$) and its second derivatives $K^{ij} = \partial^2 K / \partial x_i \partial x_j$. The total number of peaks in halo regions can thus be obtained by the summation of the peaks over all the halo regions weighted by the halo mass function. In the field region, the numbers of peaks are directly computed from the noise field N.

Specifically, the total surface number density of peaks can then be written as

$$n_{\text{peak}}(\nu) d\nu = n_{\text{peak}}^{c}(\nu) d\nu + n_{\text{peak}}^{n}(\nu) d\nu, \quad (14)$$

where $\nu = K_N/\sigma_0$ is the S/N ratio of a peak. The term $n_{\text{peak}}^c(\nu)$ is for peaks in halo regions including not only the true peaks corresponding to real haloes but also the noise peaks within the halo regions. The second term $n_{\text{peak}}^n(\nu)$ is for pure noise peaks in field regions.

For $n_{\text{peak}}^c(\nu)$ in halo regions, it can be written as

$$n_{\text{peak}}^{c}(\nu) = \int dz \frac{dV(z)}{dz d\Omega} \int_{M_{\text{lim}}} dM n(M, z) f_{\text{p}}(\nu, M, z), \quad (15)$$

where $dV(z)$ is the cosmological volume element at redshift $z$, $d\Omega$ is the solid angle element, $n(M, z)$ is the mass function of dark matter haloes. Here we adopt the Sheth-Tormen mass function in the calculation (Sheth & Tormen 1999). The mass limit $M_{\text{lim}}$ is for the mass above which individual haloes contribute dominantly to the weak lensing peak signals along their lines of sight. From simulation analyses, we find that $M_{\text{lim}} = 10^{13.7} h^{-1} M_{\odot}$ is a suitable choice. The factor $f_{\text{p}}$ is for the number of peaks in the area within the virial radius of a halo of mass M at redshift $z$, and is given by

$$f_{\text{p}}(\nu, M, z) = \int_{0}^{\theta_{\text{vir}}} d\theta (2\pi\theta) \hat{n}_{\text{peak}}^{c}(\nu, \theta, M, z), \quad (16)$$

where $\theta_{\text{vir}} = R_{\text{vir}}(M, z)/D_{\text{A}}(z)$ and $D_{\text{A}}(z)$ is the angular-diameter distance. The physical virial radius is calculated by

$$R_{\text{vir}}(M, z) = \left[ \frac{3M}{4\pi\rho(z)\Delta_{\text{vir}}(z)} \right]^{1/3}, \quad (17)$$

where $\rho(z)$ is the background matter density of the Universe at redshift $z$ and the overdensity $\Delta_{\text{vir}}$ is taken from Henry (2000).

The function $\hat{n}_{\text{peak}}^c(\nu, \theta, M, z)$ in equation (16) describes the surface number density of peaks at the location of $\theta$ from the centre of the halo, which depends on the convergence profile of the halo. On the basis of the theory of Gaussian random fields, it can be derived explicitly as (F10)

$$\begin{aligned} \hat{n}_{\text{peak}}^{c}(\nu, \theta, M, z) = & \exp \left[ - \frac{(K^1)^2 + (K^2)^2}{\sigma_1^2} \right] \\ & \times \left[ \frac{1}{2\pi\theta_*^2} \frac{1}{(2\pi)^{1/2}} \right] \exp \left[ - \frac{1}{2} \left( \nu - \frac{K}{\sigma_0} \right)^2 \right] \\ & \times \int_0^\infty dx_N \left\{ \frac{1}{[2\pi(1 - \gamma_N^2)]^{1/2}} \right. \\ & \qquad \left. \times \exp \left[ - \frac{[x_N + (K^{11} + K^{22})/\sigma_2 - \gamma_N(\nu_0 - K/\sigma_0)]^2}{2(1 - \gamma_N^2)} \right] \times F(x_N) \right\}, \end{aligned} \quad (18)$$

where $\theta_*^2 = 2\sigma_1^2/\sigma_2^2$, $\gamma_N = \sigma_1^2/(\sigma_0\sigma_2)$, $K^i = \partial_i K$, and $K^{ij} = \partial_{ij} K$. Here the quantities $\sigma_i$ are the moments of the noise field N given by (e.g. van Waerbeke 2000)

$$\sigma_i^2 = \int dk k^{2i} \langle |N(k)|^2 \rangle, \quad (19)$$

where $N(k)$ is the Fourier transform of the noise field N. For $K(\theta)$, $K^i(\theta)$, and $K^{ij}(\theta)$ of a halo with mass M at redshift $z$, they depend on the mass profile of the halo and source redshift distribution. Here we assume the spherical Navarro-Frenk-White (NFW) mass distribution for dark matter haloes (Navarro, Frenk & White 1996, 1997). The concentration parameter $c_{\text{vir}}(M) = R_{\text{vir}}/r_s$ is calculated from the mass-concentration relation given in Bhattacharya et al. (2013), where $r_s$ is the characteristic scale of an NFW halo. For
the redshift distribution of source galaxies, we take equation (1) for CS82 data.

The function $F(x_N)$ in equation (18) is given by (F10)

$$F(x_N) = \exp \left[ - \frac{(K^{11} - K^{22})^2}{\sigma_2^2} \right] \\ \times \int_0^{1/2} d e_N 8 \left( x_N^2 e_N \right) x_N^2 (1 - 4 e_N^2) \exp (-4 x_N^2 e_N^2) \\ \times \int_0^{\pi} \frac{d \theta_N}{\pi} \exp \left[ -4 x_N e_N \cos(2 \theta_N) \frac{(K^{11} - K^{22})}{\sigma_2} \right], \quad (20)$$

where $x_N = (\lambda_{N1} + \lambda_{N2})/\sigma_2$ and $e_N = (\lambda_{N1} - \lambda_{N2})/(2\sigma_2 x_N)$. $\lambda_{N1}$ and $\lambda_{N2}$ are the two eigenvalues ($\lambda_{N1} \ge \lambda_{N2}$) and $\theta_N$ is the rotation angle in the range $[0, \pi]$ with the diagonalization of $(-K_N^{ij})$ (e.g. Bardeen et al. 1986; Bond & Efstathiou 1987).

As for the field term $n_{\text{peak}}^n(\nu)$ in equation (14), it is given by

$$n_{\text{peak}}^n(\nu) = \frac{1}{d\Omega} \left\{ n_{\text{ran}}(\nu) \left[ d\Omega - \int dz \frac{dV(z)}{dz} \times \int_{M_{\text{lim}}} dM n(M, z) (\pi \theta_{\text{vir}}^2) \right] \right\}, \quad (21)$$

where $n_{\text{ran}}(\nu)$ is the surface number density of pure noise peaks without foreground haloes. It can be calculated by equation (18) with $K=0$, $K^i=0$ and $K^{ij}=0$.

It is seen that in this model, the cosmological information is contained in the halo mass function, lensing kernel, cosmic volume element, and the density profile of dark matter haloes. We note that although we use the NFW density profile and the mean mass-concentration relation derived by Bhattacharya et al. (2013) for the full sample of dark matter haloes in our fiducial model calculations, in principle, the density profile parameters can be treated as free parameters. Therefore from weak lensing peak abundances, it is possible to constrain these structural parameters simultaneously with cosmological parameters.

For our model calculation of $n_{\text{peak}}$, multidimensional integrations are needed. With great efforts numerically and applying multiple parallel techniques, such as OpenMP and GPU programming, we have developed a fast and high-precision model calculation algorithm, which makes it possible for us to perform cosmological constraints from weak lensing peak abundances. An outline of our programming structures is given in the appendix.

# 4 COSMOLOGICAL CONSTRAINTS FROM CS82 WEAK LENSING PEAK COUNTS

## 4.1 Fitting method

As described in Section 3.2, from the reconstructed convergence maps, we identify and count peaks in regions where the galaxy filling factor is $f > 0.5$. This effectively excludes the masked regions in the peak counting to avoid the mask effects. The useful survey area is approximately 114 deg$^2$. We only consider high peaks with the S/N ratio $\nu \ge 3$. With the noise level $\sigma_0 \sim 0.022$ under the smoothing scale $\theta_G = 1.5$ arcmin for CS82, the high peaks have smoothed signals $K \ge 0.066$. For the considered survey area, there are few peaks with $\nu > 6$. We therefore concentrate on the peaks in the range of $3 \le \nu \le 6$.

We divide the peaks into five bins. We consider both equal bins with $\Delta\nu = 0.5$ and unequal bins with the number of peaks comparable in different bins. For equal bins, we do not include the number of peaks in the bin of $\nu = (5.5, 6]$ in cosmological studies because it

is only ~1 with large expected statistical fluctuations. In this case, the number of peaks is ~500 for the first bin with $\nu = [3, 3.5]$ and is ~10 for the last bin with $\nu = (5, 5.5]$. For unequal bins, we have the number of peaks in the range of (~160, ~80) for different bins with $\nu = [3, 3.1], (3.1, 3.25], (3.25, 3.5], (3.5, 4], (4, 6]$, respectively.

To derive cosmological-parameter constraints from weak lensing peak counts, we calculate the $\chi^2$ defined as follows (LWPF2014):

$$\chi_{p'}^2 = dN^{(p')} (\widehat{\mathbf{C}^{-1}}) dN^{(p')} = \sum_{i,j=1,...,5} dN_i^{(p')} (\widehat{C_{ij}^{-1}}) dN_j^{(p')}, \quad (22)$$

where $dN_i^{(p')} = N_{\text{peak}}^{(p')}(v_i) - N_{\text{peak}}^{(d)}(v_i)$ with $N_{\text{peak}}^{(p')}(v_i)$ being the prediction for the cosmological model $p'$ from F10 and $N_{\text{peak}}^{(d)}(v_i)$ being the observed data for the peak counts. This effectively assumes that the number fluctuation in each bin can be approximated by a Gaussian distribution. With the number of peaks in the equal bin case being larger than ~10 per bin and being larger than about 80 per bin in the unequal bin case, the Gaussian error distribution is expected to be a good approximation. The corresponding likelihood function is given by

$$L \propto \exp \left( -\frac{1}{2} \chi_{p'}^2 \right). \quad (23)$$

The matrix $C_{ij}$ is the covariance matrix of the peak counts including the error correlations between different $\nu$ bins. Here we apply bootstrap analyses using the CS82 observational data themselves to obtain an estimate of $C_{ij}$. We realize that such an estimate cannot reveal the cosmic variance over the full survey area. Ideally, $C_{ij}$ should be constructed by generating a large number of CS82 mocks from ray-tracing simulations following exactly the same galaxy distribution, mask distribution, survey geometry, etc., as CS82 data. It is noted that the CS82 survey covers a long stripe of ~90 deg $\times$ 2 deg. Therefore to fully mimic the survey geometry, we need very large simulations to cover the ~90 deg extension, which are difficult to achieve at the moment.

For $C_{ij}$, we thus generate 10 000 bootstrap samples by resampling the 173 tiles from real observation data sets. The covariance matrix $C_{ij}$ is then calculated from the bootstrap samples by

$$C_{ij} = \frac{1}{R-1} \sum_{r=1}^{R} \left[ N_{\text{peak}}^{r}(v_i) - N_{\text{peak}}^{(d)}(v_i) \right] \left[ N_{\text{peak}}^{r}(v_j) - N_{\text{peak}}^{(d)}(v_j) \right], \quad (24)$$

where $r$ denotes for different samples with the total number of samples $R = 10000$, and $N_{\text{peak}}^r(v_i)$ is for the peak count in the bin centred on $v_i$ from the sample $r$. The inverse of the covariance matrix is then calculated by (Hartlap, Simon & Schneider 2007)

$$\widehat{\mathbf{C}^{-1}} = \frac{R - N_{\text{bin}} - 2}{R - 1} (\widehat{\mathbf{C}^{-1}}), \quad N_{\text{bin}} < R - 2, \quad (25)$$

where $N_{\text{bin}}$ is the number of bins used for peak counting.

We note that such bootstrap analyses implicitly assume the independence of the peak distribution between different tiles each with an area of about 1 deg$^2$. We have made a test by dividing the survey area into units each containing four adjacent tiles (i.e. ~2 deg $\times$ 2 deg) and then performing bootstrap analyses by resampling these units. The resulted inverse covariance, considering the diagonal elements which are much larger than the off-diagonal terms, shows less than ~8 per cent differences from the one using one tile as a unit. We also perform tests using 2 $\times$ 3 and 2 $\times$ 10 adjacent tiles as units, separately. The diagonal terms of the inverse covariance differ from those using one tile as a unit by $\le 8$ per cent and $\le 11$ per cent, respectively. We have carried out
another test to use our 15 sets of mocked data to calculate the covariance. However, our mock simulations do not mimic the long stripe geometry of the CS82 survey, and the independent unit has an area of $4 \times 3.5 \times 3.5 = 49 \text{ deg}^2$, as will be described in Section 4.2. We find that for the diagonal elements of the inverse covariance, the differences between the results from mock sets and that of using one tile as a unit are generally less than 10 per cent.

We also carry out Jackknife resampling for error estimations. For using one tile, $2 \times 2$ and $2 \times 3$ tiles as independent units, respectively, the diagonal elements of the resulted inverse covariance differ from the corresponding bootstrap resampling by less than 5 per cent.

We therefore adopt the bootstrap covariance estimated from CS82 data using one tile as an independent unit in our following fiducial analyses.

We use COSMOMC (Lewis & Bridle 2002) modified to include our likelihood function for weak lensing peak counts to perform cosmological constraints. In this paper, we mainly consider constraints on the two cosmological parameters ($\Omega_m, \sigma_8$) under the flat $\Lambda$ cold dark matter (ACDM) assumption. We adopt flat priors in the range of [0.05, 0.95] and [0.2, 1.6] for $\Omega_m$ and $\sigma_8$, respectively. We take the Hubble constant $h = 0.7$ in units of $100 \text{ km s}^{-1} \text{ Mpc}^{-1}$, the power index of the initial density perturbation spectrum $n_s = 0.96$, and the present baryonic matter density $\Omega_b = 0.046$. As discussed above, our fiducial constraints use the covariance estimated from bootstrap resampling of individual tiles. As a quantitative comparison, we also perform ($\Omega_m, \sigma_8$) constraints by using the covariance estimated from resampling $2 \times 2$ adjacent tiles. The constraint contours are nearly overlapped with those of our fiducial analyses with the area of $1\sigma$ region larger only by $\sim 1$ per cent.

To further show the potential of weak lensing peak statistics, we also perform constraints on the mass-concentration relation of dark matter haloes, assuming it follows a power-law relation, simultaneously with the cosmological parameters ($\Omega_m, \sigma_8$).

## 4.2 Mock CS82 analyses

Before presenting the results from CS82 observational data, in this part, we first show our mock CS82 analyses using ray-tracing simulations. As discussed above, the full mock of CS82 taking into account the long stripe geometry demands very large simulations, which are yet to be realized. Therefore in the current mock analyses presented here, we do not attempt to mimic the survey geometry of CS82, but are limited to the tests of our peak analyses procedures, including the convergence reconstruction, mask effects exclusion, cosmological-parameter fitting, etc. These mocks also serve as a further test of our theoretical model on the peak abundances.

To construct the CS82 mocks, we carry out dark-matter-only N-body simulations in the flat $\Lambda$CDM framework. The cosmological parameters are taken to be $\Omega_m = 0.28$, $\Omega_\Lambda = 0.72$, $\Omega_b = 0.046$, $\sigma_8 = 0.82$, $n_s = 0.96$ and $h = 0.7$. Our ray-tracing procedures largely follow our previous studies of LWPF2014, but with different box paddings.

In accord with the redshift distribution of CS82 galaxies, we perform ray-tracing calculations up to $z = 3$. For our considered cosmological model, the comoving distance to $z_s = 3$ is approximately $4.54 h^{-1} \text{ Gpc}$. We design to pad 12 independent simulation boxes to $z = 3$, with eight simulation boxes each with a size of $320 h^{-1} \text{ Mpc}$ to redshift $z = 1$ and 4 larger simulations each with a size of $600 h^{-1} \text{ Mpc}$ from $z = 1$ to $z = 3$. With these independent simulations, we can perform ray-tracing calculations straightforwardly with no repeated use of same structures.

<table><caption>**Table 1.** Redshifts of the lens planes. The planes at $z_1 > 1.0$ are produced from four independent $L = 600 h^{-1} \text{ Mpc}$ simulations, while those at lower $z_1$ are obtained from eight independent simulations with $L = 320 h^{-1} \text{ Mpc}$.</caption><tr><td>0.0107</td><td>0.0322</td><td>0.0540</td><td>0.0759</td><td>0.0981</td><td>0.1205</td></tr><tr><td>0.1432</td><td>0.1661</td><td>0.1893</td><td>0.2127</td><td>0.2364</td><td>0.2604</td></tr><tr><td>0.2847</td><td>0.3094</td><td>0.3343</td><td>0.3596</td><td>0.3853</td><td>0.4113</td></tr><tr><td>0.4377</td><td>0.4645</td><td>0.4917</td><td>0.5193</td><td>0.5474</td><td>0.5759</td></tr><tr><td>0.6049</td><td>0.6344</td><td>0.6645</td><td>0.6950</td><td>0.7261</td><td>0.7578</td></tr><tr><td>0.7900</td><td>0.8229</td><td>0.8564</td><td>0.8906</td><td>0.9254</td><td>0.9610</td></tr><tr><td>0.9895</td><td>1.0289</td><td>1.0882</td><td>1.1496</td><td>1.2131</td><td>1.2789</td></tr><tr><td>1.3472</td><td>1.4180</td><td>1.4915</td><td>1.5680</td><td>1.6475</td><td>1.7303</td></tr><tr><td>1.8166</td><td>1.9066</td><td>2.0005</td><td>2.0987</td><td>2.2013</td><td>2.3087</td></tr><tr><td>2.4213</td><td>2.5393</td><td>2.6632</td><td>2.7934</td><td>2.9296</td><td></td></tr></table>

The N-body simulations are done with GADGET-2 (Springel 2005). For both small and large sized simulations, we use 640³ particles. The mass resolution is $\sim 9.7 \times 10^9 h^{-1} \text{ Mpc}$ and $\sim 6.4 \times 10^{10} h^{-1} \text{ Mpc}$ for small boxes and large boxes, respectively. The simulations start from $z = 50$ and the initial conditions are set by using 2LPTIC (Crocce, Pueblas & Scoccimarro 2006). The initial density perturbation spectrum is generated by CAMB (Lewis, Challinor & Lasenby 2000). The force softening length is about $\sim 20 h^{-1} \text{ kpc}$, which is good enough for our studies concerning mainly high weak lensing peaks corresponding to massive dark matter haloes.

For ray-tracing simulations, we use 59 lens planes up to $z = 3$. The corresponding redshifts of the planes are listed in Table 1. We follow closely the method of Hilbert et al. (2009). The detailed descriptions for ray-tracing calculations can be found in LWPF2014. In order to generate mock data for CS82 galaxies, we calculate shear and convergence maps at the far edge of each of the 59 lensing planes using the lower redshift planes. For a set of simulations with 12 independent boxes, we then can generate four sets of lensing maps each with an area of $3.5 \times 3.5 \text{ deg}^2$ sampled on 1024 $\times$ 1024 pixels. In each set, we have 59 shear and 59 convergence maps at 59 different redshifts corresponding to the far edges of the 59 lens planes. We run 12 sets of simulations, and generate lensing maps with the total area of $12 \times 4 \times (3.5 \times 3.5) = 588 \text{ deg}^2$. This allows us to generate three nearly independent mocks for CS82, and each mock is constructed from three sets of simulations with the area of $\sim 3 \times 49 = 147 \text{ deg}^2$.

For each mock, the generating procedure is as follows.

(i) With the $3 \times 49 \text{ deg}^2$ lensing maps, we place the tiles of CS82 observed galaxies behind. We note again that we do not attempt to mimic the true CS82 long stripe survey geometry here. Therefore, we pad the CS82 tiles randomly over the simulated map area. In each tile, the positions and the amplitudes of ellipticities of the galaxies are preserved, but with their orientations being randomized. Because there is no exact redshift information for each galaxy, we assign redshifts to the galaxies following the redshift distribution of equation (1). The galaxy weights and the mask information are also preserved in each tile.

(ii) For a galaxy, its reduced shear $g$ is calculated by interpolating the signals from the pixel positions on simulated maps to the galaxy position. The interpolation is also done in the redshift dimension. Regarding the randomized ellipticity obtained in (i) as its intrinsic ellipticity, we then can construct the mock observed ellipticity for the galaxy by equation (5).

(iii) For each tile of the mock data, we perform the convergence reconstruction with the same procedure for the observed data described in Section 3.2 except we do not correct for the
<figure><img src="image_6.png" /><figcaption>**Figure 3.** CS82 mock simulation results. Upper panels show the peak count distribution in logarithmic scale with equal bins (left), and the peak count distribution in linear scale with unequal bins (right). The three sets of symbols with different colours correspond to the three sets of independent mocks. Within the same colour, five data sets corresponding to five different noise realizations are shown. The blue '*' and the error bars are for the average values and the rms over the 15 mocks. The solid line is for our model predictions. The lower panels show the relative differences between the average values of the 15 mock sets and our model predictions. The error bars are for the rms of the relative differences between each mock and the model predictions.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Cosmological constraints on ($\Omega_m, \sigma_8$) derived from the CS82 mock peak counts with equal bins (left-hand panel) and with unequal bins (right-hand panel), respectively. The blue dot is for the mock simulation input.</figcaption></figure>

multiplicative error in the mock data because our lensing signals are from simulations.

(iv) We perform the peak identifications and peak number counting in the same way as for the observational data described also in Section 3.2.

In our analyses, noise peaks from shape noise are accounted for. Therefore to obtain a good estimate of the average numbers of peaks, we randomly rotate the galaxies five times leading to five realizations of the intrinsic ellipticities for source galaxies. Therefore, we totally generate 3 × 5 mocks for CS82.

In Fig. 3, we show the peak number distributions for the mock data. The upper-left panel is for the results (in logarithmic scale) of equal bins with the bin width $\Delta\nu = 0.5$ in the range of [3., 5.5] and the upper-right panel (in linear scale) is for the case of unequal bins with $\nu = [3, 3.1], (3.1, 3.25], (3.25, 3.5], (3.5, 4], (4, 6]$, respectively. The three sets of symbols with different colours correspond to the mocks from three independent sets of simulated

maps. The five data points within each colour are the results from different shape noise realizations. The blue '*' and error bars are for the average values and rms over the 15 mocks. The solid line is for our model predictions with the shape noise level $\sigma_0$ taken to be the average value over all the tiles. The lower panels show the corresponding relative differences between the average values of the 15 mocks and our model predictions. It is seen that for both binning cases, the averaged mock results agree with our model predictions very well. The relative differences are $\le 10$ per cent, and most often $\le 5$ per cent.

In Fig. 4, we show the derived constraints on ($\Omega_m, \sigma_8$). Here, we use the results averaged over the 15 mocks as the ‘observed’ data, and the covariance matrix is estimated by constructing bootstrap samples each containing 173 tiles as the CS82 data from all the tiles in the 15 mocks. The contours are for $1\sigma$ and $2\sigma$ confidence levels. The left- and right-hand panels are for the cases of equal bins and unequal bins, respectively. The blue symbol indicates the underlying parameters of the mocks. We see that noting the strong
<figure><img src="image_8.png" /><figcaption>**Figure 5.** CS82 observational results. Left-hand panel: the peak count distribution in logarithmic scale with equal bins. Right-hand panel: the peak count distribution in linear scale with unequal bins. The corresponding solid line is the theoretical prediction with the best-fitting cosmological parameters obtained from MCMC fitting. The error bars are the square root of the diagonal terms of the covariance matrix.</figcaption></figure>

degeneracy of the two parameters, the mock constraints recover the input cosmological parameters excellently. The results from the two binning cases agree well.

## 4.3 Results from CS82 observational data

We now proceed to show the results from the CS82 observational data. In Fig. 5, the left- and right-hand panels show the results of peak counts for the equal and unequal bins, respectively. The error bars are from the bootstrap sampling using the CS82 observational tiles. The solid line in each panel is from our theoretical model with the best-fitting cosmological parameters (see Fig. 6). In Fig. 6, we show the cosmological constraints from the CS82 observed weak lensing peak counts with the red and blue contours from equal and unequal bins, respectively. The constraints from the two binning cases are consistent with each other very well. The best-fitting results and the marginalized 1D mean for $\Omega_m$ and $\sigma_8$ are shown in Table 2.

It is seen that similar to weak lensing correlation analyses (e.g. Kilbinger et al. 2013; Fu et al. 2014; Kitching et al. 2014), the constraints on the two parameters from weak lensing peak abundances alone are strongly degenerate. Considering the relation defined by

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Cosmological-parameter constraints for ($\Omega_m$, $\sigma_8$) derived from CS82 observational peak counts with equal bins (red) and with unequal bins (blue), respectively.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Table 2.** Constraints on $\Omega_m$ and $\sigma_8$ from CS82 weak lensing peak abundances. The best-fitting values and the marginalized 1D mean are shown. The errors are 68 per cent confidence intervals.</figcaption></figure>

<table><caption>**Table 3.** Constraints on $\Sigma_8$ and $\alpha$. The errors are 68 per cent confidence intervals. The 2-pt and 3-pt values are the results in Fu et al. (2014) derived from COSEBis, a second-order E-/B-mode measure and the diagonal third-order aperture-mass moment, respectively.</caption><thead><tr><th>Parameter</th><th>Best fit<br>(equal bin)</th><th>(unequal bin)</th><th>1D mean<br>(equal bin)</th><th>(unequal bin)</th></tr></thead><tbody><tr><td>Ω<sub>m</sub></td><td>0.22</td><td>0.27</td><td>0.37<sup>+0.30</sup><sub>-0.24</sub></td><td>0.38<sup>+0.27</sup><sub>-0.24</sub></td></tr><tr><td>σ<sub>8</sub></td><td>0.91</td><td>0.83</td><td>0.83<sup>+0.28</sup><sub>-0.28</sub></td><td>0.81<sup>+0.26</sup><sub>-0.26</sub></td></tr></tbody></table>

<table><thead><tr><th>Parameter</th><th>Equal bin</th><th>Unequal bin</th><th>2-pt</th><th>3-pt</th></tr></thead><tbody><tr><td>Σ<sub>8</sub></td><td>0.82 ± 0.04</td><td>0.82 ± 0.03</td><td>0.79 ± 0.06</td><td>0.73<sup>+0.09</sup><sub>-0.19</sub></td></tr><tr><td>α</td><td>0.43 ± 0.02</td><td>0.42 ± 0.02</td><td>0.70 ± 0.02</td><td>0.58 ± 0.02</td></tr></tbody></table>

$\Sigma_8 = \sigma_8(\Omega_m/0.27)^\alpha$ and using the same estimation method as in Fu et al. (2014), we obtain $\Sigma_8 = 0.82 \pm 0.04$ and $\alpha = 0.43 \pm 0.02$ and $\Sigma_8 = 0.82 \pm 0.03$ and $\alpha = 0.42 \pm 0.02$ for the two binning cases, respectively. In Table 3, we list the constraints on $\Sigma_8$ and $\alpha$ from our peak analyses and the results from 2-pt (COSEBis, a second-order E-/B-mode measure) and 3-pt (3d, diagonal third-order aperture-mass moment) analyses of Fu et al. (2014) using CFHTLenS data. It is noted that for both CS82 and CFHTLenS, the survey areas used in the analyses are ~120 deg$^2$. The mean redshift for CS82 data used in our analyses is $z \sim 0.83$, and the mean redshift for CFHTLenS used in Fu et al. (2014) is $z \sim 0.75$. It is seen that the $\alpha$ value from our peak analyses is significantly smaller than that from 2-pt and 3-pt analyses. This shows a great potential of weak lensing peak analyses in cosmological studies. In Fig. 7, we show the corresponding constraints demonstrating the differences visually.

We note that the above 2-pt COSEBis and 3-pt 3d analyses from Fu et al. (2014) are in $\Lambda$CDM model with five free parameters ($\Omega_m$, $\sigma_8$, $\Omega_b$, $n_s$, $h$). Their results shown in Table 3 and Fig. 7 here are the results marginalized over ($\Omega_b$, $n_s$, $h$). On the other hand, in our peak analyses, we fix the other three parameters and vary only ($\Omega_m$, $\sigma_8$). To see if the constraints, particularly the degeneracy direction between ($\Omega_m$, $\sigma_8$), can be affected significantly by allowing more
<figure><img src="image_11.png" /><figcaption>**Figure 7.** The comparison for the constraints on ($\Omega_m$, $\sigma_8$) between our peak analyses and the results from 2-pt and 3-pt analyses of Fu et al. (2014).</figcaption></figure>

free cosmological parameters, we study the dependence of the peak abundances on ($\Omega_m$, $\sigma_8$, $\Omega_b$, $n_s$, $h$) by calculating the derivatives with respect to these parameters using the model of F10. It is found that the dependences of peak abundances on $\Omega_m$ and $\sigma_8$ are much stronger than the dependences on the other three parameters. Therefore, we do not expect that the inclusion of ($\Omega_b$, $n_s$, $h$) as free parameters can change our results on the constraints of ($\Omega_m$, $\sigma_8$) from peak analyses considerably. We further perform a not-so-rigorous Fisher analysis using the derivatives of the peak abundances with respect to the five parameters and errors corresponding to CS82 data. Assuming Gaussian priors for $\Omega_b$, $n_s$ and $h$ with $\sigma_{\Omega_b} = 0.05$, $\sigma_{n_s} = 0.1$ and $\sigma_h = 0.3$, we find that the constraints on ($\Omega_m$, $\sigma_8$) marginalized over ($\Omega_b$, $n_s$, $h$) are about the same as the results with only ($\Omega_m$, $\sigma_8$) as free parameters. This shows again that our results should not be affected considerably if we include ($\Omega_b$, $n_s$, $h$) in our peak analyses.

It is noted that the degeneracy of ($\Omega_m$, $\sigma_8$) from our peak analyses is comparable to the constraints from cluster studies. From Sunyaev-Zel'dovich effect (SZ) cluster abundance analyses, $\alpha \sim 0.3$ has been obtained (e.g. Hasselfield et al. 2013; Reichardt et al. 2013; Planck Collaboration: Ade et al. 2014b). For X-ray cluster studies, $\alpha \sim 0.5$ (e.g. Vikhlinin et al. 2009; Böhringer, Chon & Collins 2014). Using SDSS 'maxBCG' cluster catalogue, Rozo et al. (2010) derived $\alpha \sim 0.41$. Our constraint is $\alpha \sim 0.42$, which is in good agreement with the cluster studies noting the variations between different analyses largely due to the different observable-mass relation. This is expected because high weak lensing peaks have close associations with clusters of galaxies along lines of sight.

In the very recent studies of LPH2015, they obtain $\Sigma_8 = \sigma_8(\Omega_m/0.27)^{0.60} = 0.84^{+0.03}_{-0.04}$ using the peak analyses alone combining the results from two smoothing scales. The $\alpha$ value is larger than ours of $\alpha = 0.42$ and $\Sigma_8$ is somewhat smaller. There are a number of differences between their analyses and ours. They use the peaks spanning a large range of $\kappa$ value from negative to positive. In our cosmological studies, we only consider high peaks with $\nu \ge 3$, which corresponds to $\kappa \ge 0.066$ for $\theta_G = 1.5$ arcmin (corresponding to $\sim 1$ arcmin in LPH2015 because of the different definition of $\theta_G$). We expect that our results should resemble more those of cluster studies as explained in the previous paragraph. For low peaks, besides the impact of noise, they are related to the projection effects of large-scale structures, which might be

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Constraints from WMAP9 (green) and Planck (red) overplotted on ours from peak counts with unequal bins. The contours are for $1\sigma$ and $2\sigma$ confidence levels.</figcaption></figure>

largely described by the underlying power spectrum. However, in LPH2015, they find that although weak lensing peaks show stronger covariance with the power spectrum than cluster counts, the overall covariance between peaks and the power spectrum is rather weak. Therefore, even low peaks and the power spectrum should contain non-overlapping cosmological information. The error covariance matrix estimation is also different in the two studies. We use the bootstrap approach by using the data themselves. In LPH2015, they use a fiducial simulation to calculate the covariance by randomly rotating and shifting the simulation box of size $240h^{-1}\text{Mpc}$ in ray tracing. To fully explore the cause of the differences from different analyses can be a worthwhile task in the future.

In Fig. 8, we show the constraints in comparison with the results from WMAP9 (Hinshaw et al. 2013, green) and Planck (Planck Collaboration: Ade et al. 2014a, red). It can be seen that our constraints are in good agreement with both.

The above analyses adopt the fiducial redshift distribution given in equation (1) with $a = 0.531$, $b = 7.810$ and $c = 0.517$ for CS82 galaxies. Derived by matching to COSMOS galaxies, this redshift distribution can have significant uncertainties. To understand the impact of the uncertainties on our peak analyses, we follow Hand et al. (2015) to vary the parameters in the redshift distribution and analyse how the peak abundances change using the theoretical model of F10. Considering peaks in the range of $\nu = [3, 6]$, we find that by shifting the peak position of the redshift distribution by $\Delta z = \pm 0.1$, the peak abundances change by $\sim \pm 3$ per cent to $\sim \pm 20$ per cent from low to high peaks. Varying the $b$ parameter by $+30$ per cent ($-30$ per cent) leads to $\sim -2$ per cent ($+5$ per cent) to $\sim -10$ per cent ($+20$ per cent) changes in the peak abundances. A $\pm 30$ per cent change in $c$ parameter causes $\le 5$ per cent changes in peak abundances. The largest impact is from the uncertainties in $a$ parameter. Changing $a$ by $+30$ per cent ($-30$ per cent) leads to $\sim +8$ per cent ($-5$ per cent) to $\sim +30$ per cent ($-20$ per cent) variations in the peak abundances from low to high peaks. With these varied redshift distributions, we estimate the best fit $\sigma_8$ from observed CS82 peak abundances (unequal bins) by keeping $\Omega_m$ to be approximately the best-fitting value from our fiducial analyses. The results are shown in Fig. 9. It is seen that for all the above changes, the best fit $\sigma_8$ are within the $1\sigma$ range of our fiducial constraints. Therefore, we do not expect a considerable impact on our results from the uncertainties
<figure><img src="image_13.png" /><figcaption>**Figure 9.** The impact of the uncertainties in the galaxy redshift distribution on cosmological-parameter constraints derived from CS82 weak lensing peak abundances.</figcaption></figure>

in the galaxy redshift distribution. On the other hand, for future surveys with dramatically improved statistics, the uncertainties in the redshift distribution can be an important source of systematic errors for weak lensing peak statistics, and the required accuracy for redshift measurements needs to be carefully studied.

## 4.4 Other constraints

Besides cosmological-parameter constraints, weak lensing peak statistics may also possibly provide constraints on the density profile of dark matter haloes because massive dark matter haloes are the sources of true weak lensing peaks with high S/N ratios (e.g. Yang et al. 2013). In our theoretical model of F10, the dependence of the peak abundances on the density profile of dark matter haloes is explicit. This allows us to perform constraints on the structural parameters of dark matter haloes simultaneously with cosmological parameters. Similar ideas are also recently proposed in other studies (e.g. Cardone et al. 2014; Mainini & Romano 2014). To show this

<figure><img src="image_14.png" /><figcaption>**Figure 10.** The four-parameter constraints from the mocked (left) and observed (right) unequal binned weak lensing peak counts. The contours are for $1\sigma$ and $2\sigma$ confidence levels</figcaption></figure>

<table><caption>**Table 4.** Marginalized 1D mean from the four-parameter fitting using unequally binned peak counts. The errors are 68 per cent confidence intervals.</caption><thead><tr><th>Parameter</th><th>Mock</th><th>CS82 observation</th></tr></thead><tbody><tr><td>Ω<sub>m</sub></td><td>0.38<sup>+0.30</sup><sub>-0.26</sub></td><td>0.34<sup>+0.30</sup><sub>-0.23</sub></td></tr><tr><td>σ<sub>8</sub></td><td>0.80<sup>+0.28</sup><sub>-0.27</sub></td><td>0.84<sup>+0.28</sup><sub>-0.27</sub></td></tr><tr><td>A</td><td>9.7<sup>+5.4</sup><sub>-4.8</sub></td><td>7.5<sup>+5.2</sup><sub>-4.6</sub></td></tr><tr><td>β</td><td>-0.15<sup>+1.41</sup><sub>-1.45</sub></td><td>-0.15<sup>+1.37</sup><sub>-1.31</sub></td></tr></tbody></table>

feasibility, we assume a power-law mass–concentration relation for NFW haloes in our model calculations with (e.g. Duffy et al. 2008)

$$c_{\text{vir}} = \frac{A}{(1+z)^{0.7}} \left( \frac{M_{\text{vir}}}{10^{14} \text{h}^{-1} \text{M}_{\odot}} \right)^{\beta}, \quad (26)$$

where A and β are regarded as free parameters. The redshift dependence $(1+z)^{0.7}$ is taken to be consistent with recent simulation results (e.g. Bhattacharya et al. 2013). We then perform four-parameter fitting from the observed weak lensing peak counts. The flat priors are [0, 20] for A, [-2, 2] for β, [0.05, 0.95] for $\Omega_m$, and [0.2, 1.6] for $\sigma_8$. The results are shown in Fig. 10 where the unequal binned peak counts are used. The left-hand panels are for the results from mock data described in Section 4.2, and the right-hand panels are from CS82 observational data. We see that although the obtained constraints are mainly for $(\Omega_m, \sigma_8)$, the constraints on the plane of $(A, \beta)$ are apparent even for the current generation of weak lensing surveys. The marginalized 1D mean for both mock and the CS82 observational analyses are shown in Table 4. The results for A and β are in broad agreement with the results from simulated haloes although the error ranges are large. We also note that in this four-parameter fitting, the constraint contours on $(\Omega_m, \sigma_8)$ plane are enlarged somewhat in comparison with the results of the two-parameter fitting case, showing the influence of the uncertainties of the halo structural parameters on the cosmological constraints.

For future surveys with much larger survey areas and the improved depth, statistically we expect significant enhancements of weak lensing peak analyses, which in turn will provide us valuable

<figure><img src="image_15.png" /></figure>
cosmological information complementary to cosmic shear correlation analyses.

## 5 SUMMARY AND DISCUSSIONS

With CS82 weak lensing observations, we study the weak lensing peak abundances for $v$ in the range [3, 6], and derive the first cosmological constraints from peak analyses. We summarize the results as follows.

(i) For flat $\Lambda$CDM, the cosmological constraints on ($\Omega_m, \sigma_8$) from peak analyses are fully consistent with the constraints obtained from cosmic shear correlation studies. On the other hand, the degeneracy direction of the two parameters is flatter than those from the correlation analyses. Quantitatively, with $\Sigma_8 = \sigma_8(\Omega_m/0.27)^\alpha$, we obtain $\alpha = 0.42 \pm 0.02$ in comparison with $\alpha = 0.70 \pm 0.02$ from COSEBis 2-pt correlations and $\alpha = 0.58 \pm 0.02$ from diagonal three-order aperture-mass correlation studies (Fu et al. 2014). This shows a promising potential of weak lensing peak analyses complementary to correlation studies. It is noted that to explore the improvements on the cosmological-parameter constraints from the combined analyses, the full covariance between the peak abundances and the correlation functions should be investigated carefully.

(ii) Our derived cosmological constraints from peak analyses are also consistent with both WMAP9 and Planck results.

(iii) We perform constraints on ($A, \beta$), the power-law form of the mass-concentration relation of dark matter haloes, simultaneously with the cosmological parameters ($\Omega_m, \sigma_8$). For the CS82 survey with relatively large statistical errors, the current constraints are mainly on ($\Omega_m, \sigma_8$). However, the constraints on ($A, \beta$) are already apparent. This shows the capability to constrain the structural parameters of massive structures together with cosmological parameters from weak lensing peak statistics. With much improved data from future surveys, performing simultaneous constraints on the structural and cosmological parameters can potentially allow us to extract important astrophysical effects on the structural evolution of the mass distribution of haloes (e.g. Yang et al. 2013). Meanwhile, it can also avoid the possible bias on cosmological-parameter constraints resulting from the pre-assumption about the halo structures in predicting the cosmological dependence of weak lensing peak abundances.

In this paper, we adopt the theoretical model of F10 for predicting the peak abundances. The model takes into account the effects of shape noise in the calculation, and has been tested extensively with ray-tracing simulations in our previous studies (F10; LWPF2014). For CS82, the shape noise is the dominant source of error. The applicability of the model is further shown with our mock analyses presented in Section 4.2 here. Comparing to the approach fully relying on large simulations, theoretical modelling can help us understand better different effects, and can allow us to explore cosmological- and dark matter halo density profile parameter space efficiently.

On the other hand, for future surveys with much reduced statistical errors, our theoretical model needs to be developed to include the effects neglected in the current treatment, such as the projection effects of large-scale structures, the complex mass distribution of haloes, etc. Unlike the shape noise, such effects themselves also contain cosmological information. With the help of simulations, we are currently working towards improving our model for future cosmological applications.

Weak lensing effects are unique in probing the dark side of the Universe. Current generation of completed surveys, such as CS82

and CFHTLenS, have served as important demonstrations to show the feasibility of weak lensing cosmological studies. Ongoing surveys, such as DES (The Dark Energy Survey Collaboration 2005), HSC (Hyper Suprime-Cam Design Review 2009) and KiDS (de Jong et al. 2013), will expand the survey area to a few thousand square degrees. Future ones, such as Euclid (Amendola et al. 2013) and LSST (LSST Science Collaboration 2009), will target at nearly half of the sky of about 20 000 deg$^2$. The statistical capability of weak lensing studies will increase tremendously. To fully realize the power, however, different systematics, both observational and theoretical ones, need to be understood thoroughly. For this, the current surveys also play important roles in revealing different obstacles that need to be overcome and further paving the road to the future.

## ACKNOWLEDGEMENTS

Based on observations obtained with MegaPrime/MegaCam, a joint project of CFHT and CEA/DAPNIA, at the CFHT, which is operated by the National Research Council (NRC) of Canada, the Institut National des Science de l'Univers of the Centre National de la Recherche Scientifique (CNRS) of France and the University of Hawaii. The Brazilian partnership on CFHT is managed by the Laboratório Nacional de Astrofísica (LNA). This work made use of the CHE cluster, managed and funded by ICRA/CBPF/MCTI, with financial support from FINEP and FAPERJ. We thank the support of the Laboratório Interinstitucional de e-Astronomia (LIneA). We thank the CFHTLenS team for their pipeline development and verification upon which much of this surveys pipeline was built.

We thank the referee for the constructive comments and suggestions. We thank the useful discussions with Z. Haiman and J. Liu. XKL, CZP and ZHF acknowledge the support from NSFC of China under the grants 11333001, 11173001 and 11033005, and from the Strategic Priority Research Program 'The Emergence of Cosmological Structures' of the Chinese Academy of Sciences, grant no. XDB09000000. HYS acknowledges the support by a Marie Curie International Incoming Fellowship within the 7th European Community Framework Programme, and NSFC of China under grants 11103011. LPF acknowledges the support from NSFC grants 11103012 and 11333001, and Shanghai Research grant 13JC1404400 of STCSM. LR acknowledges the support from NSFC grant 11303033 and the support from Youth Innovation Promotion Association of CAS. AL acknowledges the support by World Premier International Research Centre Initiative (WPI Initiative), MEXT, Japan. JPK acknowledges support from the ERC advanced grant LIDA and from CNRS. BM acknowledges financial support from the CAPES Foundation grant 12174-13-0. MM is partially supported by CNPq (grant 486586/2013-8) and FAPERJ (grant E-26/110.516/2-2012).

The MCMC calculations are partly done on the Laohu supercomputer at the Centre of Information and Computing at National Astronomical Observatories, Chinese Academy of Sciences, funded by Ministry of Finance of Peoples Republic of China under grant ZDYZ2008-2. The N-body simulations are performed on the Shuguang cluster at Shanghai Normal University, Shanghai, China.

## REFERENCES

* Albrecht A. et al., 2006, preprint (arXiv:astro-ph/0609591)

* Amendola L. et al., 2013, Living Reviews in Relativity, 16, 6

* Bard D. et al., 2013, ApJ, 774, 49

* Bardeen J. M., Bond J. R., Kaiser N., Szalay A. S., 1986, ApJ, 304, 15

* Bartelmann M., 1995, A&A, 303, 643
Bartelmann M., Schneider P., 2001, Phys. Rep., 340, 291

Bhattacharya S., Habib S., Heitmann K., Vikhlinin A., 2013, ApJ, 766, 32

Böhringer H., Chon G., Collins C. A., 2014, A&A, 570, A31

Bond J. R., Efstathiou G., 1987, MNRAS, 226, 655

Cardone V. F., Camera S., Sereno M., Covone G., Maoli R., Scaramella R., 2015, A&A, 574, A141

Crocce M., Pueblas S., Scoccimarro R., 2006, MNRAS, 373, 369

de Jong J. T. A., Kleijn G. V., Kuijken K., Valentijn E., KiDS, Astro-WISE
Consortiums, 2013, Exp. Astron., 35, 25

Dietrich J. P., Hartlap J., 2010, MNRAS, 402, 1049

Duffy A. R., Schaye J., Kay S. T., Dalla Vecchia C., 2008, MNRAS, 390, L64

Erben T. et al., 2013, MNRAS, 433, 2545

Fan Z., Shan H., Liu J., 2010, ApJ, 719, 1408 (F10)

Fox L., Parker I. B., 1968, Chebyshev Polynomials in Numerical Analysis.
Oxford Univ. Press, Oxford

Fu L. et al., 2008, A&A, 479, 9

Fu L. et al., 2014, MNRAS, 441, 2725

Gavazzi R., Soucail G., 2007, A&A, 462, 459

Geller M., Kurtz M. J., Dell'Antonio I. P., Ramella M., Fabricant D. G., 2010, ApJ, 709, 832

Hamana T., Takada M., Yoshida N., 2004, MNRAS, 350, 893

Hamana T., Oguri M., Shirasaki M., Sato M., 2012, MNRAS, 425, 2287

Hand N. et al., 2015, Phys. Rev. D, 91, 062001

Hartlap J., Simon P., Schneider P., 2007, A&A, 464, 399

Hasselfield M. et al., 2013, J. Cosmol. Astropart. Phys., 7, 008

Hennawi J. F., Spergel D. N., 2005, ApJ, 624, 59

Henry J. P., 2000, ApJ, 534, 565

Heymans C. et al., 2012, MNRAS, 427, 146

Hilbert S., Hartlap J., White S. D. M., Schneider P., 2009, A&A, 449, 31

Hilbert S., Marian L., Smith R. E., Desjacques V., 2012, MNRAS, 426, 2870

Hinshaw G. et al., 2013, ApJS, 208, 19

Hyper Supreme-Cam Design Review, 2009. Available at: http://www.naoj.org/Projects/HSC/report.html

Jauzac M., Jullo E., Kneib J. P., 2012, MNRAS, 426, 3369

Jee M. J., Tyson J. A., Schneider M. D., Wittman D., Schmidt S., Hilbert S., 2013, ApJ, 765, 74

Jullo E., Pires S., Jauzac M., Kneib J. P., 2014, MNRAS, 437, 3969

Kaiser N., Squires G., 1993, ApJ, 404, 441

Kaiser N., Squires G., Broadhurst T., 1995, ApJ, 449, 460

Kilbinger M. et al., 2013, MNRAS, 430, 2200

Kitching T. D. et al., 2014, MNRAS, 442, 1326

Kratochvil J. M., Haiman Z., May M., 2010, Phys. Rev. D, 81, 043519

Lewis A., Challinor A., Lasenby A., 2000, ApJ, 538, 473

Lewis A., Bridle S., 2002, Phys. Rev. D, 66, 103511

Li R. et al., 2014, MNRAS, 438, 2864

Lin C. A., Kilbinger M., 2015, A&A, 576, A24

Liu J., Petri A., Haiman Z., Hui L., Kratochvil J. M., May M., 2015, Phys. Rev. D, 91, 063507 (LPH2015)

Liu X., Wang Q., Pan C., Fan Z., 2014, ApJ, 784, 31 (LWPF2014)

LSST Dark Energy Science Collaboration, 2012, preprint (arXiv:1211.0310)

LSST Science Collaboration, 2009, preprint (arXiv:0912.0201)

Mainini R., Romano A., 2014, J. Cosmol. Astropart. Phys., 8, 63

Marian L., Smith R. E., Bernstein G. M., 2009, ApJ, 698, L33

Marian L., Smith R. E., Bernstein G. M., 2010, ApJ, 709, 286

Marian L., Smith R. E., Hilbert S., Schneider P., 2012, MNRAS, 423, 1711

Maturi M., Angrick C., Pace F., Bartelmann M., 2010, A&A, 519, A23

Miller L., Kitching T. D., Heymans C., Heavens A. F., van Waerbeke L., 2007, MNRAS, 382, 315

Miller L. et al., 2013, MNRAS, 429, 2858

Miyazaki S., Hamana T., Ellis R. S., Kashikawa N., Massey R. J., Taylor J., Refregier A., 2007, ApJ, 669, 714

Navarro J., Frenk C., White S. D. M., 1996, ApJ, 462, 563

Navarro J., Frenk C., White S. D. M., 1997, ApJ, 490, 493

Petri A., Haiman Z., Hui L., May M., Kratochvil J. M., 2013, Phys. Rev. D, 88, 123002

Planck Collaboration XVI, 2014a, A&A, 571, A16

Planck Collaboration XX, 2014b, A&A, 571, A20

Press W. H., 2007, Numerical Recipes, 3rd edn. The Art of Scientific Computing. Cambridge Univ. Press, Cambridge

Reichardt C. L. et al., 2013, ApJ, 763, 127

Rozo E. et al., 2010, ApJ, 708, 645

Rykoff E. S. et al., 2014, ApJ, 785, 104

Schneider P., van Waerbeke L., Jain B., Kruse G., 1998, MNRAS, 296, 873

Seitz C., Schneider P., 1995, A&A, 297, 287

Seitz C., Schneider P., 1997, A&A, 318, 687

Semboloni E., Schrabback T., van Waerbeke L., Vafaei S., Hartlap J., Hilbert S., 2011, MNRAS, 410, 143

Shan H. et al., 2012, ApJ, 748, 56

Shan H. et al., 2014, MNRAS, 442, 2534

Sheth R. K., Tormen G., 1999, MNRAS, 308, 119

Squires G., Kaiser N., 1996, ApJ, 473, 65

Tang J. Y., Fan Z. H., 2005, ApJ, 635, 60

The Dark Energy Survey Collaboration, 2005, preprint (arXiv:astro-ph/0510346)

van Waerbeke L., 2000, MNRAS, 313, 524

van Waerbeke L. et al., 2013, MNRAS, 433, 3373

Vikhlinin A. et al., 2009, ApJ, 692, 1060

Weinberg D. H., Mortonson M. J., Eisenstein D. J., Hirata C., Riess A. G., Rozo E., 2013, Phys. Rep., 530, 87

White M., van Waerbeke L., Mackey J., 2002, ApJ, 575, 640

Wittman D., Dell'Antonio I. P., Hughes J. P., Margoniner V. E., Tyson J. A., Cohen J. G., Norman D., 2006, ApJ, 643, 128

Yang X., Kratochvil J. M., Wang S., Lim E. A., Haiman Z., May M., 2011, Phys. Rev. D, 84, 043529

Yang X., Kratochvil J. M., Huffenberger K., Haiman Z., May M., 2013, Phys. Rev. D, 87, 023511

APPENDIX A:

In this appendix, we describe our programming structure for computing theoretically the peak abundances using the model of F10. For clarity, we copy the relevant equations in Section 3.3 here. For a S/N ratio *v*, the surface number density of peaks is given by

$$
n_{\text{peak}}(v)dv = n_{\text{peak}}^{c}(v)dv + n_{\text{peak}}^{n}(v)dv. \quad (A1)
$$

The part for pure noise peaks $n_{\text{peak}}^n(v)$ can be calculated easily if the regions occupied by haloes are known. Therefore the most computationally heavy part is $n_{\text{peak}}^c(v)$ for peaks in halo regions. In F10, it can be calculated by

$$
n_{\text{peak}}^{c}(v) = \int dz \frac{dV(z)}{dz d\Omega} \int_{M_{\text{lim}}} dM n(M, z) f_{p}(v, M, z), \quad (A2)
$$

and

$$
f_p(v, M, z) = \int_0^{\theta_{vir}} d\theta (2\pi\theta) \hat{n}_{\text{peak}}^c(v, \theta, M, z), \quad \text{(A3)}
$$

where

$$
\begin{equation}
\begin{split}
\hat{n}_{\text{peak}}^{c}(v, \theta, M, z) ={}& \exp\left[-\frac{(K^1)^2 + (K^2)^2}{\sigma_1^2}\right] \\
& \times \left[\frac{1}{2\pi\theta_2^2} \frac{1}{(2\pi)^{1/2}}\right] \exp\left[-\frac{1}{2}\left(v - \frac{K}{\sigma_0}\right)^2\right] \\
& \times \int_0^{\infty} dx_N \left\{ \frac{1}{[2\pi(1 - \gamma_N^2)]^{1/2}} \right. \\
& \qquad \left. \times \exp\left[-\frac{[x_N + (K^{11} + K^{22})/\sigma_2 - \gamma_N(v_0 - K/\sigma_0)]^2}{2(1 - \gamma_N^2)}\right] \times F(x_N) \right\}.
\end{split}
\tag{A4}
\end{equation}
$$
<figure><img src="image_16.png" /><figcaption>Figure A1. Sketch for our code calculations.</figcaption></figure>

It is noted that given a pair of (M, z) for a halo, the function $\hat{n}_{\text{peak}}^{c}(v, \theta, M, z)$ in equation (A4) can be computed independently for different $\theta$ and $v$. We therefore employ GPU for this part of calculations, which improves our computational efficiency enormously.

We illustrate our programming structures in Fig. A1. The specifics are as follows.

(i) We first divide the halo mass range and the redshift range into $N_1 \times N_2$ grid points as shown in the bottom layer of Fig. A1. Memories are allocated for mass $M$ ($N_1$ sample points) and redshift

$z$ ($N_2$ sample points) matrix on host memory (CPU memory). We then calculate the corresponding virial radius $r_{\text{vir}}(M, z)$ on CPU according to equation (17).

(ii) Then the results of $r_{\text{vir}}$ are copied to global memory (GPU memory). We allocate global memory for different $\theta$ ($N_3$ sample points) and SNR $v$ ($N_4$ sample points) as shown in the top layer of Fig. A1.

(iii) We take the advantage of GPU shared memory to evaluate the smoothed convergence field $K$ and its first and second derivatives for an NFW halo using Chebyshev interpolation algorithm for integration (Fox & Parker 1968). The shared memory is highly efficient and suitable to use in this situation.

(iv) Launch the main kernel to evaluate $f_p(v, M, z)$ according to equations (A3) and (A4). This involves massive amount of calculations because of the multidimensional integrations and is benefited greatly by using the GPU feature.

(v) At last, we copy the results of $f_p(v, M, z)$ from global memory to host memory, and calculate the final result of $n_{\text{peak}}^{c}(v)$ using the 11th-order Simpson method on CPU (Press 2007).

The use of GPU improves the model calculation efficiency significantly. For example, our own workstation consists of two Xeon E5-2697v2 CPUs and four Nvidia GTX Titan GPU. Each CPU contains 12 cores while each GPU contains 2880 cores. In this case, we achieved a calculation speed which is more than 20 times faster by using GPU than that by using CPU alone.

This paper has been typeset from a TeX/LATeX file prepared by the author.