The self similarity of weak lensing peaks

TL;DR

This paper investigates the statistics of weak lensing convergence peaks, revealing their self-similar clustering behavior and how their abundance and clustering can jointly constrain cosmological parameters like _m and _8.

Contribution

It introduces accurate models for peak abundance and clustering, demonstrating their potential to break cosmological parameter degeneracies.

Findings

Peak abundance depends on cosmology within 3% accuracy.

Peak 2PCFs exhibit self-similar behavior after rescaling.

Combined peak statistics can break _m-_8 degeneracy.

Abstract

We study the statistics of weak lensing convergence peaks, such as their abundance and two-point correlation function (2PCF), for a wide range of cosmological parameters $\Omega_m$ and $\sigma_8$ within the standard $\Lambda$CDM paradigm, focusing on intermediate-height peaks with signal-to-noise ratio (SNR) of $1.5$ to $3.5$. We find that the cosmology dependence of the peak abundance can be described by a one-parameter fitting formula that is accurate to within $\sim3\%$. The peak 2PCFs are shown to feature a self-similar behaviour: if the peak separation is rescaled by the mean inter-peak distance, catalogues with different minimum peak SNR values have identical clustering, which suggests that the peak abundance and clustering are closely interconnected. A simple fitting model for the rescaled 2PCF is given, which together with the peak abundance model above can predict peak 2PCFs with an accuracy better than $\sim5\%$. The abundance and 2PCFs for intermediate peaks have very different dependencies on $\Omega_m$ and $\sigma_8$, implying that their combination can be used to break the degeneracy between these two parameters.

Equations

$$\text{centering}\kappa =\frac{1}{2}{\nabla }^2\Psi ,\text{@add@centering}$$

$$\Psi (\mathbf{\theta })=\frac{D_{ls}}{D_l{D}_s}\frac{2}{c^2}{\int }_0^{z_s}\Phi (D_l\mathbf{\theta },z)dz.$$

$${\nabla }^2\Phi =4\pi G{a}^2\delta {\rho }_m,$$

$$\kappa (\mathbf{\theta })={\int }_0^{z_s}W(z)\delta {\rho }_m(D_l(z)\mathbf{\theta },z)dz,$$

$$W(z)=\frac{3H_0^2{\Omega }_m}{2c}\frac{1+z}{H(z)}\chi (z){\int }_z^{z_s}\frac{dn}{d{z}_s}d{z}_s\frac{\chi (z_s)-\chi (z)}{\chi (z_s)},$$

$$\nu =\frac{\kappa -\mu }{\sigma },$$

$$\xi_{\rm{LS}}(\theta)=1+\bigg{(}\frac{N_{R}}{N_{D}}\bigg{)}^{2}\frac{DD}{RR}-\bigg{(}\frac{N_{R}}{N_{D}}\bigg{)}\frac{DR}{RR}\,,$$

$$\sigma=m\Sigma_{8}(\alpha)+c\equiv m\sigma_{8}\bigg{(}\frac{\Omega_{m}}{0.3}\bigg{)}^{\alpha}+c\;.$$

$$\log n(>\nu )=B_n\nu +A_n,$$

$$\log n(>\nu )=\nu B_n(A_n)+A_n,$$

$$B_n(A_n)=-0.33A_n+\mathrm{0.28.}$$

$$A_n=-0.30{\sigma }_8({\Omega }_m/0.3{)}^{-0.637}+2.20.$$

$$\xi=\xi_{0}\bigg{(}\frac{\theta}{\theta_{p}}\bigg{)}^{\beta}\,,$$

$${\xi }_0={\xi }_{0,a}{x}^2+{\xi }_{0,b}x+{\xi }_{0,c},$$

$${\xi }_{0,a}=0.253,{\xi }_{0,b}=-0.605,{\xi }_{0,c}=0.514,$$

$${\sigma }_{\mathrm{pix}}^2=\frac{{\sigma }_{\mathrm{int}}^2}{2{\theta }_{\mathrm{pix}}^2{n}_{\mathrm{gal}}},$$

$$\mathrm{d}{P}_{ij}(r)={\bar{n}}^2\left[1+\xi (r)\right]\mathrm{d}{V}_i\mathrm{d}{V}_j,$$

$$\xi_{\rm PH}(\theta)=\bigg{(}\frac{N_{R}}{N_{D}}\bigg{)}^{2}\frac{DD(\theta)}{RR(\theta)}-1,$$

$${\xi }_{\mathrm{DP}}(\theta )=2\frac{N_R}{N_D}\frac{DD(\theta )}{DR(\theta )}-1,$$

$$\xi_{\rm{LS}}(\theta)=1+\bigg{(}\frac{N_{R}}{N_{D}}\bigg{)}^{2}\frac{DD(\theta)}{RR(\theta)}-\bigg{(}\frac{N_{R}}{N_{D}}\bigg{)}\frac{DR(\theta)}{RR(\theta)}.$$

$$⟨{\xi }_{\mathrm{LS}}(\theta ){⟩}_1=⟨1+{\left(\frac{N_R}{N_D}\right)}^2\frac{DD}{RR}-\left(\frac{N_R}{N_D}\right)\frac{DR}{RR}⟩,$$

$$⟨{\xi }_{\mathrm{LS}}(\theta ){⟩}_2=1+{\left(\frac{⟨N_R⟩}{⟨N_D⟩}\right)}^2\frac{⟨DD⟩}{⟨RR⟩}-\left(\frac{⟨N_R⟩}{⟨N_D⟩}\right)\frac{⟨DR⟩}{⟨RR⟩},$$

$$⟨{\xi }_{\mathrm{LS}}(\theta ){⟩}_1\ne ⟨{\xi }_{\mathrm{LS}}(\theta ){⟩}_2.$$

Full text

††pubyear: 2019††pagerange: The self similarity of weak lensing peaks–D

1 Introduction

Gravitational lensing, the deflection of photon trajectory due to the presence of massive objects between the source and the observer, is sensitive to the total matter distribution in the Universe (Bartelmann & Schneider, 2001). The total matter content can be expressed as the sum of baryons ($\Omega_{b}$) and dark matter ($\Omega_{\rm{DM}}$). Combined, these two components make up $\Omega_{m}\approx 30\%$ of the Universe’s energy budget with an abundance ratio of $\Omega_{\rm{DM}}/\Omega_{b}\approx 5$ (Planck Collaboration et al., 2018). While dark matter is the dominant matter component, it is also the more challenging one to detect since it is not directly observable. Gravitational lensing offers a promising probe which allows us to examine the distribution of dark matter in the Universe.

In this paper we are interested in weak lensing (WL), which is the regime of gravitational lensing where the amplitude of light deflections is small. These deflections, caused by inhomogeneities of the total matter distribution, lead to distortions of the images of distant sources, such as background galaxies or the cosmic microwave background. While the distortions are small, they can be extracted with a careful statistical analysis of their correlations (e.g., Bacon et al., 2000; Kaiser et al., 2000; Van Waerbeke et al., 2000; Wittman et al., 2000). Nowadays, WL is widely used as a probe of the large-scale structures (LSS) of the Universe (Albrecht et al., 2006; LSST Dark Energy Science Collaboration, 2012; Amendola et al., 2013; Weinberg et al., 2013).

While WL is only the broad name of a physical phenomenon, a plethora of statistics can be extracted from observations and cosmological simulations. For example, the cosmic shear (the distortion of the shapes of a lensed image) and convergence (the magnification of the magnitude of the lensed source), are usually quantified using their two-point statistics, though higher-order statistics contain extra information about the nonlinear LSS evolution. Analyses performed on the low-redshift universe using shear-shear correlations have led to strong cosmological parameter constraints (e.g., Schneider et al., 2002; Semboloni et al., 2006; Hoekstra et al., 2006; Fu et al., 2008; Heymans et al., 2012; Kilbinger et al., 2013; Hildebrandt et al., 2017). Lensing signatures are also present in the cosmic microwave background for which the lensing potential power spectra can be used as a cosmological test (e.g., Planck Collaboration et al., 2018). WL has also been detected around cosmological objects such as galaxy clusters, which can be used to measure the cluster mass (e.g., Gruen et al., 2014; Umetsu et al., 2014; Applegate et al., 2014; Von der Linden et al., 2014; Hoekstra et al., 2015; Tudorica et al., 2017), and voids, for which void WL profiles can be obtained (e.g., Melchior et al., 2014; Clampitt & Jain, 2015; Sánchez et al., 2017).

Weak lensing can be used to test different cosmological models. For example, the convergence power spectrum, shear-shear correlations as well as WL by voids have been shown to be promising probes for constraining the dark energy equation of state or modified gravity models (e.g. Schmidt, 2008; Tsujikawa & Tatekawa, 2008; Huterer, 2010; Cai et al., 2015; Barreira et al., 2015; Barreira et al., 2017; Baker et al., 2018; Cautun et al., 2018; Van Uitert et al., 2018; Paillas et al., 2019). While WL offers an independant measurement of the absolute cluster mass, this cannot be done for most clusters; one use of WL is to calibrate some observable-cluster-mass scaling relation, which is then used to infer cluster masses and hence the abundance of massive haloes, which is a powerful cosmological probe (e.g. Mantz et al., 2015; Umetsu et al., 2016; Pizzuti et al., 2017; Bocquet et al., 2018). In addition, statistics such as WL Minkowski functionals and WL bispectrum have been used to constrain the sum of neutrino masses (e.g. Marques et al., 2018; Coulton et al., 2018).

The most commonly used WL statistics are the shear correlation function and convergence power spectrum. These two-point statistics alone, however, cannot account for the non-Gaussian features introduced by the nonlinear evolution of structures in the Universe, and other statistics can provide additional and complementary information. In this paper we study one such additional probe, WL peaks, which are the maxima of the WL convergence field. The WL peak abundance is a good example of a statistic that contains complimentary information to two-point statistics (Jain & Van Waerbeke, 2000; Pen et al., 2003; Dietrich & Hartlap, 2010; Shirasaki et al., 2018), and can be used to constrain cosmological parameters within $\Lambda$CDM (Shan et al., 2012; Van Waerbeke et al., 2013; Shan et al., 2014; Liu et al., 2015b), to test alternative cosmological models such as modified gravity (Cardone et al., 2013; Liu et al., 2016b; Higuchi & Shirasaki, 2016; Shirasaki et al., 2017; Peel et al., 2018), dark energy (Giocoli et al., 2018), and to measure the neutrino mass (Li et al., 2018). WL peaks can also be extracted from CMB lensing to provide cosmological constraints (Liu et al., 2016a). Various models have been developed to accurately describe high signal-to-noise-ratio (SNR) WL peaks (e.g., Hamana et al., 2004; Hennawi & Spergel, 2005; Maturi et al., 2005; Fan et al., 2010; Marian et al., 2012; Hamana et al., 2012; Liu & Haiman, 2016; Shan et al., 2018; Wei et al., 2018).

In contrast to high WL peaks, there have been relatively few studies on the abundance of low and intermediate peaks (see, e.g., Yang et al., 2011; Lin & Kilbinger, 2015; Shirasaki, 2017), which nevertheless contains rich cosmological information (Dietrich & Hartlap, 2010; Kratochvil et al., 2010; Yang et al., 2011), and even fewer on the spatial correlation of such peaks (e.g., Marian et al., 2013; Shan et al., 2014). Upcoming wide and deep field galaxy surveys such as euclid (Refregier et al., 2010) and lsst (LSST Science Collaboration et al., 2009) will produce large high-resolution WL maps, with significant improvements compared to the current generation of WL observations. Understanding how WL peak statistics behave will be important if we want to maximise the cosmological information that can be gained from the new surveys. In particular, the higher source number density of these surveys will lead to a reliable determination of peak abundance and clustering down to low SNR values, so it is important to have accurate models to describe the statistics of low- and intermediate-height peaks.

In this work we study properties of WL peak statistics in $\Lambda$CDM, by modelling the peak abundance, peak two-point correlation functions (2PCFs) and the convergence rms fluctuation (convergence map standard deviation). Most importantly, we identify a universal self-similar behaviour in the peak 2PCF, which holds for a large range of peak heights and different cosmologies. The self-similarity is observed when expressing the 2PCF in terms of the angular separation divided by the mean peak separation, with the resulting rescaled 2PCFs lying on top of each other. We propose a general model that describes the abundance and clustering of WL peaks and that allows us to access cosmological information contained on non-linear scales.

The structure of the paper is outlined as follows: in Section 2 we briefly summarise the relevant theory for WL, describe the numerical simulations used to construct the WL maps, and introduce the statistics we use to study WL peaks. Next, in Section 3 we present the WL peak abundance, WL peak 2PCF, and identify a self similarity in the peak 2PCF for a given fiducial cosmology. Then, in Section 4 we give general fitting functions that describe the convergence rms fluctuation, peak abundance, 2PCF and its self similarity in $\Lambda$CDM for a large range of $\Omega_{m}$ and $\sigma_{8}$ values. We then show in Section 5 that our model can accurately reproduce the original peak 2PCF. Finally, in Section 6 we show that the self similarity of the 2PCF is robust to the inclusion of galaxy shape noise.

2 Theory, simulations and analysis pipeline

In this section we briefly summarise the essential elements of WL theory, and describe the simulations and methodology used in this work to study WL peak statistics.

2.1 Theory

The deflection of photon trajectories due to the mass of objects in the lens plane can be quantified using the lensing convergence $\kappa$,

[TABLE]

where $\Psi$ is the 2D lensing potential given by:

[TABLE]

In the above equation, $\bm{\theta}$ is the observed angular position of the lensed image (in the Born approximation); $D_{s}$,$D_{l}$ and $D_{ls}$ are the angular diameter distances between the observer and source, observer and lens, and lens and source; $z_{s}$ is the source redshift; and $c$ is the speed of light. The symbol $\Phi$ denotes the 3D gravitational potential given by the Poisson equation,

[TABLE]

where $\delta\rho_{m}=\rho_{m}-\overline{\rho}_{m}$ with $\rho_{m}$ and $\overline{\rho}_{m}$ respectively the local and the background matter density; $a$ is the scale factor; and $G$ is the gravitational constant. Eq. (1) therefore behaves as a projected version of the Poisson equation weighted by a geometric factor determined by the distances between the lens, source and observer.

The above equations describe the WL effect of a single lens. However, as the light travels between the source and the observer, it experiences gravitational lensing from the entirety of the mass distribution along its path. The single lens equation can be generalised to multiple lenses as

[TABLE]

where $W(z)$ is the lensing kernel that accounts for the redshift distribution of the lenses. This kernel is given by

[TABLE]

where $H(z)$ is the Hubble parameter and $H_{0}$ is its present-day value; $\Omega_{m}$ is the fractional total matter density at present-day; $\chi$ is the comoving distance; and $\frac{dn}{dz_{s}}$ is the redshift distribution of sources.

As can be seen from Eq. (4), the WL convergence corresponds to the projected mass density contrast weighted by a geometric factor, and thus positive and negative $\kappa$ values correspond to overdense and underdense lines of sight. For self consistency across different convergence maps we define the SNR, $\nu$, as

[TABLE]

where $\mu$ is the mean value of the convergence field of a given map, and $\sigma$ is its rms fluctuation. We have subtracted the map mean $\mu$ in the definition of SNR because our maps are relatively small and can have non-zero means due to sample variance that vary from map to map, which can affect the consistency of the SNR definition. Note that the subtraction of $\mu$ does not affect $\sigma$, and it is not needed in the case of $\kappa$ reconstructed from the (directly observable) cosmic shear field. An example of a $\kappa$ map generated from numerical simulations through ray tracing is shown in Fig. 1.

For most of the paper, the lensing quantities are measured in maps obtained from cosmological simulations without added noise, except for Section 6 where galaxy shape noise is included and all quantities are measured from the noise-added maps (for more details see Section 6).

2.2 Numerical simulations

In order to study weak lensing peak statistics, in this paper we use a large suite of WL convergence maps constructed from two sets of N-body simulations. The first are the publicly-available all-sky convergence maps of Takahashi et al. (2017) (hereafter T17). These maps have a source redshift of $z_{s}=1$ and have been generated using the ray tracing algorithm described in Hamana et al. (2015) (see also Shirasaki et al., 2015), with a HEALpix resolution of $N_{\rm{side}}=16384$. To avoid probing the same structures along the line-of-sight, T17 constructed the light cone by stacking a hierarchy of cubic simulation boxes, with comoving sizes $L,2L,3L,\cdots,14L$, where $L=450h^{-1}$Mpc. The simulations had a particle number of $2048^{3}$, where the particle mass depends on the box size, and ranges from $8.2\times 10^{8}$ to $2.3\times 10^{12}M_{\odot}$ (see Table 1 of T17 for more details). Each of the simulation boxes was duplicated 8 times and nested around the observer, such that nests of larger boxes contain nests of smaller boxes at their centers. Ray tracing was then performed on the nested simulation boxes by taking the projected mass distribution in spherical shells of $150~{}h^{-1}\mathrm{Mpc}$ in thickness centred on the observer (see T17 for illustration). The cosmological parameters adopted for the T17 simulations are $\Omega_{m}=0.279$, $\sigma_{8}=0.820$ and $h=0.7$, where $h=H_{0}/100$ km s*-1* Mpc*-1*. Throughout this paper we have split the T17 all-sky map into 184 separate $10\times 10~{}\mathrm{deg}^{2}$ maps with $2048^{2}$ pixels per map, for which we can use the flat sky approximation to simplify our analysis. A detailed description of the method we use to split the all-sky map into smaller squares is given in Appendix A

The second set of WL maps we use are taken from Zorrilla Matilla et al. (2016) (hereafter Z16; see also Gupta et al. 2018) and consist of maps for 96 different cosmologies. It was built with the simulation pipeline described in Petri (2016). For each cosmology, the maps were obtained from an N-body simulation of a periodic box with length $L=240~{}h^{-1}\mathrm{Mpc}$ and $512^{3}$ simulation particles with a particle mass of $\sim 10^{10}h^{-1}M_{\odot}$ (the exact value depends on the actual cosmology). Ray tracing was then performed by using a source redshift of $z_{s}=1$ and by stacking particles into lens planes with a thickness of $80~{}h^{-1}\mathrm{Mpc}$ between the source and the observer. The lens planes were generated by taking a slice along a coordinate axis of the original simulation box and applying a random shift and rotation. This process was repeated to generate $512$ $3.5\times 3.5~{}\mathrm{deg}^{2}$ maps per cosmology with $1024\times 1024$ pixels per map. Note that each of the 512 maps were obtained from the same periodic simulation by varying the orientation of the line-of-sight direction. For a more detailed description we refer the reader to Z16.

In total we have two sets of maps, one with $184$ $10\times 10~{}\mathrm{deg}^{2}$ maps for a fixed cosmology and the other with $512$ $3.5\times 3.5~{}\mathrm{deg}^{2}$ maps for 96 cosmologies with different values of $\Omega_{m}$ and $\sigma_{8}$. Larger maps are ideal for 2PCF studies as the 2PCF cannot be reliably calculated at large separations where pair measurements are affected by the finite size of the map. However, the differences in the two simulation data sets used here bring some benefits for our analysis. First, given the simulations use different ray tracing codes and box tiling methods, if we are able to identify certain features of the WL peak statistics in both simulations, this can be a check that the said features are not an unphysical consequence of the procedure used to generate the convergence maps. Second, the different simulation maps have different angular sizes and resolutions, which can help highlight any potential systemics in our analysis due to the box size or the pixel resolution.

2.3 Weak lensing peaks

WL peaks in this paper are defined as the maxima of the convergence field, which trace local over-densities in a range of environments. To extract the WL peaks, we first smooth the convergence map with a Gaussian filter. The convergence field has power on all scales, so the number and spatial distribution of WL peaks depends on the smoothing scale, with a larger smoothing washing out low contrast peaks and merging neighbouring peaks. We mainly study peaks identified with a Gaussian filter with smoothing length $\theta_{s}=1~{}\mathrm{arcmin}$, a range of smoothing scales have been studied in Liu et al. (2015a), showing that this smoothing scale is ideal for WL peak studies. In some cases (which will be explicitly mentioned), we vary $\theta_{s}$ to understand how the results depend on smoothing scale. Next we identify WL peaks by finding all pixels in the maps whose values are larger than those of their 8 neighbours, and peaks within $3\theta_{s}$ of the map boundary are removed to avoid edge effects where the Gaussian filter is truncated. The height of a peak is given by the $\nu$ value of the smoothed convergence field at the peak position. For a given convergence map, we can generate multiple peak catalogues by imposing a $\nu_{\mathrm{cut}}$ threshold and keeping only peaks with $\nu\geq\nu_{\mathrm{cut}}$.

Fig. 1 illustrates the distribution of peaks (shown as green dots) for three different SNR thresholds, $\nu_{\mathrm{cut}}=1,2$ and $3$. To highlight the distribution of peaks on both small and large scales, we show peaks identified with a Gaussian smoothing scale, $\theta_{s}=2~{}\mathrm{arcmin}$; using a smaller $\theta_{s}$ value would result in many more peaks and make the graph less legible.

Fig. 1 shows that peak catalogues with different $\nu_{\mathrm{cut}}$ values trace different features of the convergence field. The catalogue with $\nu_{\mathrm{cut}}=1$ traces the over-dense regions of the convergence map, whilst avoiding the darker under-dense regions. In particular, many peaks seem to be arranged in a somewhat filamentary pattern. By increasing $\nu_{\mathrm{cut}}$ to 2, we find that the resulting catalogue has a significantly lower number of peaks and the peaks are now more clustered. Most of these peaks are found in highly over-dense regions, with some small filamentary patterns still remaining. Finally, there are few peaks with $\nu_{\mathrm{cut}}=3$, but they show a high degree of spatial clustering and are located in the very over-dense regions of the map.

The description of Fig. 1 above highlights two important features in the behaviour of WL peaks: the number of WL peaks and their clustering, which are respectively quantified by two commonly used statistics, the peak abundance and the peak two point correlation function (2PCF). The former is well studied and has been considered for many cosmological applications (e.g., Liu et al., 2015a; Liu & Haiman, 2016; Liu et al., 2016b; Shirasaki et al., 2017; Shirasaki, 2017; Shan et al., 2018; Li et al., 2018; Wei et al., 2018), whereas weak lensing 2PCFs are usually measured as shear-shear correlations (Fu et al., 2008; Heymans et al., 2012; Kilbinger et al., 2013), with very few studies directly focused on the peak 2PCFs (Marian et al., 2013; Shan et al., 2014).

The two point correlation function measures the probability of finding two points (or in our case, WL peaks) at a given separation ($\theta$ for angular separations on the sky). It can also be interpreted as a measure of the excessive clustering of a distribution of points relative to the clustering of randomly distributed points. To estimate the 2PCFs, we use the Landy-Szalay estimator (Landy & Szalay, 1993), which is a robust way of measuring 2PCFs, especially for small maps and low tracer number densities. Using this estimator requires a catalogue of randomly distributed points, whose role is to account for boundary effects and serves as a proxy for the volume (area in 2D) of the sample. The Landy-Szalay estimator is evaluated as

[TABLE]

where $N_{D}$ and $N_{R}$ are the numbers of data and random points and $DD$, $DR$ and $RR$ are the number of data-data, data-random and random-random pairs in bins $\theta\pm\delta\theta$. We calculate the 2PCFs by taking the average over many small maps (see description in section 2.2). Since the maps are small, taking the average of the $\xi$ values measured for each map leads to biased results and we discuss this subtlety in detail in appendix B. To obtain unbiased results, we calculate the average of the $DD$, $DR$ and $RR$ pair counts over all maps, and then we insert the average pair counts into Eq. (7).

3 Weak lensing statistics

As mentioned above, we are mainly interested in the one- and two-point statistics of WL peaks. In order to gain some first insight into the properties of these quantities, we use the large, $100~{}\mathrm{deg}^{2}$, maps from the T17 simulations for the results shown in this section. In the next section we shall use the small maps from the Z16 simulations to quantify the dependence of peak statistics on cosmology.

3.1 Peak abundance

We start by studying the mean abundance of WL peaks, which is expressed in terms of the cumulative peak abundance, $n(>\nu)$. This represents the number density in deg*-2* of all peaks whose SNRs are higher than $\nu$. The peak abundance is illustrated in Fig. 2, where the results are averaged over the 184 T17 maps. The blue (upper) and orange (lower) curves correspond to a Gaussian smoothing kernel $\theta_{s}=1$ and $2~{}\mathrm{arcmin}$, respectively. According to Eq. (6), the smoothing scale $\theta_{s}$ enters the definition of $\nu$ in two ways, by affecting the pixel values of $\kappa$ and the overall rms $\kappa$ fluctuation, $\sigma$. For the T17 maps we find $\sigma=0.013$ and $0.010$ respectively for $\theta_{s}=1$ and $2~{}\mathrm{arcmin}$, and in the next section we will see that $\sigma$ has a clear cosmology dependence as well.

The qualitative behaviour shown in Fig. 2 is as expected. There are very few peaks with high $\nu$ values since these correspond to massive dark matter structures, which are rare. As $\nu$ decreases, the peak abundance, $n(>\nu)$, increases quickly until $\nu\sim 0$ since lower $\nu$ values correspond to lower mass and thus more abundant dark matter structures. However, for $\nu\lesssim 0$111Note that WL peaks can have $\nu<0$, or equivalently $\kappa<0$. These are local maxima in underdense regions of the convergence map. we see that $n(>\nu)$ flattens, showing that there are few peaks with $\nu<0$. It highlights that there are few structures in underdense regions that are massive enough to lead to a local maximum, especially when smoothing over $1$ and $2~{}\mathrm{arcmin}$. Increasing the smoothing scale $\theta_{s}$ leads to a lower peak abundances at fixed $\nu$, since smoothing over a larger region tends to eliminate some peaks.

We have checked that the WL peak abundance shown in Fig. 2 can be well fitted by the function $n({>}\nu)=-a\left[\tanh(b\nu)-1\right]$ for the entire $\nu$ range. However, for reasons that will become clear in Section 3.3, in this paper we are interested in the range $\nu\in[1.5,3.5]$ (shown as the grey shaded region in Fig. 2), where $n(>\nu)$ can be modelled as a power law (see Section 4.2). Note that Fig. 2 also shows the uncertainties in the $n(>\nu)$ measurement, which are the standard errors of the mean of the 184 T17 maps; however, these error bars are not visible as they are roughly of the same size as the line width.

3.2 Peak two point correlation function and $\nu_{\rm{cut}}$ dependence

For WL peaks it has been suggested that the 2PCF can be well modelled by a power law (Shan et al., 2014). In this section, we will confirm this power-law description using the T17 convergence maps, and show that it works well for peak catalogues with a wide range of $\nu_{\mathrm{cut}}$ thresholds.

The 2PCF dependence on $\nu_{\mathrm{cut}}$ is of particular interest, because by decreasing $\nu_{\mathrm{cut}}$ we are including lower peaks into the analysis, which is equivalent to incorporating smaller dark matter structures into the clustering statistics. In the current standard cosmological paradigm, large-scale structures (LSS) evolve hierarchically, with larger objects forming from higher initial density peaks. This means that by varying $\nu_{\mathrm{cut}}$ we probe the different regimes of nonlinear LSS formation and thus potentially provides more powerful cosmological tests. As an example, in certain modified gravity models, smaller structures experience a stronger boost in their nonlinear growth (e.g., Clifton et al., 2012, and references therein), and we expect this to leave potentially detectable signatures in the peak 2PCFs at different $\nu_{\rm cut}$ values. In addition, as we have seen above, lowering the $\nu_{\mathrm{cut}}$ threshold increases the number of peaks included in the catalogue, and this can help increase the statistical constraining power. We will see shortly that there is a self-similarity in the peak 2PCF, which means that having peak catalogues for multiple $\nu_{\mathrm{cut}}$ values does not require separate modelling for each catalogue; this can potentially strongly improve the constraining power by WL peaks.

The left panels of Fig. 3 show the mean 2PCFs of the T17 maps for a range of $\nu_{\mathrm{cut}}$ values and for two smoothing scales, $\theta_{s}=1$ and $2~{}\mathrm{arcmin}$. The error bars, which are the standard errors of the 184 maps, are shown as shaded regions around the curves, but they are very small and barely visible.

A quick inspection of Fig. 3 by eye confirms that the 2PCFs are well described by power laws. We can see that as $\nu_{\rm{cut}}$ increases the amplitude of the 2PCF, $\xi(\theta)$, also increases. This is intuitive to understand: the high WL peaks correspond to more massive structures which tend to cluster more strongly. Moreover, the 2PCF amplitude is higher for peak catalogues obtained using a larger smoothing length, $\theta_{s}$. This is because a higher $\theta_{s}$ value suppresses peaks originating from low mass dark matter structures that cluster less. The gradient of the 2PCFs from maps with a fixed $\theta_{s}$ increases slightly with $\nu_{\mathrm{cut}}$, but this effect is weak for both smoothing scales shown, and the dominant effect of varying $\nu_{\mathrm{cut}}$ is in the amplitude of the 2PCF.

Note that smoothing can lead to a merging of peaks separated by distances comparable to the smoothing scale $\theta_{s}$, and therefore eliminates some peaks which are close to each other. This leads to a drop off of the 2PCF from the power law on scales $\lesssim\theta_{s}$, which is why we show a different $\theta$ range in the two left panels of Fig. 3. Additionally, for a given WL map size, the 2PCFs cannot be reliably measured at large separations as there are too few peak pairs, which is why in Fig. 3 we adopted a conservative $\theta_{\rm max}$ which is 1/10 the map size. Therefore, the smoothing scale and the map size set a limited range in $\theta$ within which we can measure the 2PCF. More explicitly, while for the T17 maps we use $\theta_{s}\in[1,3]$ arcmin in this study, for the smaller Z16 maps we only use $\theta_{s}=1$ arcmin to avoid having a too narrow $\theta$ range. A larger $\theta_{s}$ is necessary for maps where galaxy shape noises (GSNs) are included, to suppress the biasing effects caused by the latter (e.g., Davies et al., 2018); we will come back to this point in Section 6 below.

3.3 2PCF rescaling and self-similarity

We now move on to one of the most important results of this work: the self similarity of the peak 2PCFs. This has been first studied (very briefly) in Davies et al. (2018) in the context of explaining the self-similar behaviour of the abundances for voids identified from WL peaks with varying $\nu_{\mathrm{cut}}$. As we show later, the 2PCF self similarity is a very useful property that merits the more detailed investigation presented here.

The quest for a self-similar behaviour in the peak 2PCFs is motivated by the following observations: the 2PCF amplitude is lower for peak catalogues with lower $\nu_{\mathrm{cut}}$; meanwhile, these catalogues have more peaks and hence a smaller mean peak separation, $\theta_{p}$. By expressing the 2PCF in terms of $\theta/\theta_{p}$ the various curves could potentially be brought closer together. The question is whether after this rescaling the 2PCF curves for different $\nu_{\mathrm{cut}}$ thresholds can be made to overlap, in which case their modelling can be significantly simplified.

The right panels in Fig. 3 show the rescaled 2PCF, that is the 2PCF expressed as a function of $\theta/\theta_{p}$ instead of $\theta$. To obtain this result, we calculated the mean peak separation as $\theta_{p}=(N/A)^{-1/2}$, where $N$ is the number of peaks in a catalogue and $A$ is the area of the map. The $\theta_{p}$ value for a peak catalogue can be directly inferred from the peak abundance, $n(>\nu)$, as $\theta_{p}=n(\nu>\nu_{\mathrm{cut}})^{-1/2}$. We find that the rescaled 2PCFs lie on top of each other and thus it indicates that the peak 2PCF is self similar. This shows that the one-point statistic of WL peaks, $n(>\nu)$, can be tied to the amplitude of the 2PCF to achieve the mentioned self similarity. The self-similar behaviour is mainly limited to the range $\nu_{\mathrm{cut}}\in[1.5,3.5]$, with the rescaled 2PCFs starting to peel off from the average relation for $\nu_{\rm{cut}}<1.5$ and $\nu_{\rm{cut}}>3.5$. At this stage it is unclear whether the breakdown of self-similarity at $\nu_{\mathrm{cut}}>3.5$ is physical or due to the small number of high SNR peaks in our maps (which, as discussed in Appendix B, could bias the estimation of the two-point correlation function); this will be investigated in more detail in the future. Bearing this issue in mind, in this work we limit our investigation to the modelling of WL peak statistics for $1.5<\nu_{\mathrm{cut}}<3.5$ only. Note that this happens to be the same range within which the peak abundance can be well described by a power law (see Section 3.1).

The self-similar behaviour holds for both smoothing scales shown in Fig. 3, however, the rescaled 2PCFs for the larger smoothing length ($\theta_{s}=2~{}\mathrm{arcmin}$; bottom right panel of Fig. 3) are shifted to lower $\theta/\theta_{p}$ values than the results for $\theta_{s}=1~{}\mathrm{arcmin}$. It suggests that 2PCFs are self-similar when keeping the smoothing scale constant, and that the self-similarity behaviour does not hold when comparing 2PCFs obtained for peak catalogues with different smoothing scales.

The panels in the right column of Fig. 3 also show that the rescaled 2PCFs are well fitted by a power law, as shown by the black dashed curves with gradients $-0.94$ ($\theta_{s}=1~{}\mathrm{arcmin}$) and $-1.02$ ($\theta_{s}=2~{}\mathrm{arcmin}$) .

4 Cosmology dependence and universal fitting functions

In this section we study the dependence of peak abundance and peak 2PCF on the $\Omega_{m}$ and $\sigma_{8}$ cosmological parameters by analysing these statistics for the set of 96 different cosmologies used for the Z16 maps. The ($\Omega_{m}$,$\sigma_{8}$) parameter space of the Z16 maps is indicated by the points in Fig. 4. The parameter space is densely sampled around $\Omega_{m}=0.26$ and $\sigma_{8}=0.8$, which corresponds to the fiducial cosmology, and only sparsely sampled for models with very different parameter values. In particular, when describing the cosmology dependence of various peak statistics, we will limit our fitting procedure to the ($\Omega_{m}$, $\sigma_{8}$) pairs shown as orange points in the figure. This removes extreme and unrealistic cosmological parameters from our analysis. For comparison, the parameters used for the T17 maps are indicated by the black triangle in Fig. 4.

The two cosmological parameters, $\Omega_{m}$ and $\sigma_{8}$, are degenerate because they can impact the size of the matter fluctuations in similar ways, and the direction of degeneracy depends on the physical quantity which is being studied. In order to better assess the potential and limitations of using WL peak statistics to constrain these cosmological parameters, it is important to know the degeneracy direction for the physical quantities of interest. Following the usual approach, we define the parameter combination, $\Sigma_{8}(\alpha)\equiv\sigma_{8}\left(\Omega_{m}/0.3\right)^{\alpha}$, where $\alpha$ characterises the degeneracy for a given statistic ($\alpha$ is allowed to vary for different statistics since they usually do not have exactly identical degeneracy directions). Note that for studying the cosmology dependence we use only the Z16, and not the T17, convergence maps, and the latter is used as a consistency check of the fitted models.

The fittings carried out in this section are mainly to exemplify the cosmology dependence of the self-similar feature present in the 2PCF in Fig. 3, which applies only to theoretical (simulated) lensing maps with no noises and with a specific smoothing length. Before this approach can be used for observational constraints, further development will be required, notably the inclusion of galaxy shape noise. We discuss briefly the impact the latter has on the self similarity of the 2PCF in Section 6. In order to study the rescaled 2PCF for a range of cosmologies using more realistic noise-added maps, larger simulations for this range of models are required, which we leave to future work.

4.1 Convergence rms fluctuation

We describe peaks in terms of the convergence SNR value at their position. To calculate this, we use the root-mean-square (rms) fluctuations, $\sigma$, of the convergence field (see e.g. Eq. 6). In principle, $\sigma$ is used merely as a normalisation factor and it is not entirely unreasonable to use the same value to define $\nu$ across all cosmologies. However, the standard deviation (or rms fluctuation) of the corresponding WL convergence map, $\sigma$ is a quantity with a clear physical meaning, and hence it is natural to use its correct value for a given cosmology. Therefore, we need a general description of $\sigma$ as a function of input cosmological parameters, $\sigma=\sigma(\Omega_{m},\sigma_{8})$. Having this function is also of interest on its own, since it is useful to know how the rms fluctuation of the WL convergence field depends on the cosmological model.

The dependence of the convergence rms fluctuation, $\sigma$, on $\sigma_{8}$ and $\Omega_{m}$ is illustrated in Fig. 5, where we show the results for two smoothing lengths, $\theta_{s}=1$ (left panel) and $2~{}\mathrm{arcmin}$ (right panel). In both cases we varied $\alpha$ such that $\sigma$ is well described by a linear function of $\Sigma_{8}(\alpha)$, that is:

[TABLE]

We achieved this by performing a $\chi^{2}$ minimisation procedure with three free parameters: $m$, $c$ and $\alpha$. The best fitting parameter values are: $m=2.08\times 10^{-2}$, $c=-3.39\times 10^{-3}$ and $\alpha=6.66\times 10^{-1}$ for $\theta_{s}=1~{}\mathrm{arcmin}$ and $m=1.59\times 10^{-2}$, $c=-2.37\times 10^{-3}$ and $\alpha=6.56\times 10^{-1}$ for $\theta_{s}=2~{}\mathrm{arcmin}$. For these fits, we used the standard errors in the determination of $\sigma$. These errors are shown in Fig. 5 as vertical error bars, however in most cases the error bars are smaller than the size of the data points plotted and are not visible. The bottom panels in Fig. 5 show the residuals, i.e., the deviations of the $\sigma$ values for each cosmology from the best-fitting lines. The fit residuals are small when compared to the absolute values of $\sigma$ and show no systematic trend with $\Sigma_{8}$, indicating that Eq. (8) provides an excellent description for both smoothing lengths. In particular, we notice that the best-fit value of $\alpha$ is similar for $\theta_{s}=1$ and $2~{}\mathrm{arcmin}$, and thus weakly dependent on the smoothing scale,. It is also reassuring to note that T17 is in very good agreement with the best-fitting lines, which are obtained from the Z16 convergence maps only, even though the two sets of maps are constructed from very different simulations.

Even though it is not the main line pursued in this work, we note that the measured standard deviation of the (reconstructed) WL maps is a useful quantity (see, e.g., Van Waerbeke et al., 2013) and therefore a simple fitting formula for $\sigma$ in terms of $\Omega_{m}$ and $\sigma_{8}$ will be useful both theoretically and observationally. However, because GSN introduces a major systematic uncertainty in real WL maps, it is necessary to study the Eq. (8) fitting formula using maps in which realistic GSN is included; this is beyond the scope of this study, because the Z16 maps used in the analysis above have relatively small sizes. In this paper, instead, the primary use of Eq. (8) will be to define $\nu$ for a convergence map with given $\Omega_{m}$ and $\sigma_{8}$ values.

As shown in Figure 5, within the large range of cosmological parameters covered in this study, $\sigma$ varies strongly, by up to a factor of $5$-$7$. By defining $\nu$ relative to the $\sigma$ of the corresponding model, we are able to define the SNR in different models relative to their own clustering amplitude. The alternative way to is use a constant $\sigma$ definition, such as the value for a fiducial model or the rms of the typical noise map. However, in our case this would mean comparing the clustering of map pixels with smaller $\kappa$ values in one model to the clustering of pixels with large $\kappa$ values in another model, when using the same $\nu_{\rm{cut}}$, and so a cosmology dependant $\nu_{\rm{cut}}$ range would have to be applied which is a far more complicated approach. We have explicitly checked by defining the SNR $\nu$ in Eq. (6) using the $\sigma$ of the fiducial model where $\sigma_{8}=0.8$ and $\Omega_{m}=0.3$, and found that the self-similar behaviour still holds though slightly worse than shown here. Later in Section 6, when dealing with noisy maps, we shall use the $\sigma$ measured from the smoothed noisy maps to define the SNR.

4.2 WL peak abundance

The peak abundances for the 96 models in the Z16 simulations are shown in Fig. 6. The colour-bar shows the cosmological parameters for a given curve with the form $\Sigma_{8}=\sigma_{8}(\Omega_{m}/0.3)^{-0.637}$. The spread of the amplitudes of the peak abundances across the 96 cosmologies is up to a factor of two. Most of the curves appear straight, indicating that the peak abundance is well described by a power law. The only exceptions are the black curves with very small values of $\Sigma_{8}$: these curves will be removed in the analysis as they correspond to the extreme cosmological parameter values indicated by the blue points in Fig. 4. As $\Sigma_{8}$ increases, the slope of the peak abundance decreases, while its amplitude increases at larger $\nu$ ($\nu\gtrsim 2.5$) and decreases for small $\nu$ ($\nu\lesssim 2.5$). This ‘rotation’ of the curves about $\nu\approx 2.5$ as $\Sigma_{8}$ changes implies a correlation between the slope and amplitude of the peak abundances as the cosmological parameters vary, and we shall see shortly that this fact can be utilised to reduce the number of fitting parameters in our peak abundance model.

In the range $1.5<\nu_{\rm{cut}}<3.5$, the peak abundance as a function of $\nu$ is well described by the power-law,

[TABLE]

in which the fitting parameters $A_{n}$ and $B_{n}$ depend on the input cosmology.

In order to model $A_{n}=A_{n}(\Omega_{m},\sigma_{8})$ and $B_{n}=B_{n}(\Omega_{m},\sigma_{8})$, we fit Eq. (9) to the peak abundance for each of the cosmological models indicated by the orange points in Fig. 4. The fitting results confirm that $A_{n}$ and $B_{n}$ – which respectively characterise the amplitude and slope of the peak abundance – are strongly correlated, so that $B_{n}$ can be replaced with a function of $A_{n}$, $B_{n}(A_{n})$, and the peak abundance can now be described using a one-parameter power law of the form

[TABLE]

with

[TABLE]

This indicates that a universal model for WL peak abundance that works for a wide range of cosmological models can be obtained if one can fit the cosmology dependence of the single parameter $A_{n}$. The result is shown in the left panel of Fig. 7, where the $A_{n}$ values measured from the 96 Z16 cosmologies are plotted against $\Sigma_{8}(\alpha)$ with $\alpha=-0.637$. The value of $\alpha$ corresponds to the one for which $A_{n}$ is well fitted by a linear function of $\Sigma_{8}(\alpha)$. The latter is shown as the solid line in the left panel of Fig. 7, and is given by

[TABLE]

For fitting the above equation the errors in $A_{n}$ are given by the uncertainties of fitting Eq. (11) to the peak abundance. These errors are small and are not visible in Fig. 7 since they are smaller than the symbol size.

Eqs. (10)-(12) can be used to predict the peak abundance for any input $\Omega_{m}$ and $\sigma_{8}$ values. For this, we first use Eq. (12) to calculate $A_{n}$ for given ($\Omega_{m}$, $\sigma_{8}$), and then infer $n(>\nu)$ using Eq. (11). The accuracy of this prediction is shown by circle symbols in the right panel of Fig. 7. We quantify the success of the method in terms of the mean percentage residual, which is defined as the absolute value of the fractional difference between the measured peak abundance and the predicted one, averaged over all bins in $\nu\in[1.5,3.5]$. We find a mean percentage residual of $1-5$ percent, indicating that the model performs well. To understand what is the major factor affecting the model accuracy, we also calculate the mean percentage residual between the measured $n(>\nu)$ and the direct one-parameter power-law fit to it, which is shown as triangles in the right panel of Fig. 7. There is only a very slight difference between the triangles and the circles, with the former being generally lower. This indicates that our prediction of the peak abundance is very similar in accuracy to the original fit, which is further supported by the fact that the blue line in the left panel fits the various $A_{n}$ data points rather well. In summary, our model is able to predict the peak abundance to within ${\sim}3$ percent accuracy for most of the chosen cosmological models.

If we use the original form of the power law, Eq. (9), to model the peak abundance, with the two parameters, $A_{n}$ and $B_{n}$, both left to vary freely in the fitting, we get fits and predictions that match the raw data at the sub percent level. However, in attempt to minimise the number of parameters in our model we chose the one parameter power law, Eq. (10), at the cost of roughly a $2\%$ loss in accuracy.

4.3 Peak two-point correlation functions

We now move on to check whether peak 2PCFs display the self-similar behaviour described in Section 3.3 for a wide range of cosmological models. The result is shown in Fig. 8, where we plot the rescaled 2PCFs for all the pairs of $(\Omega_{m},\sigma_{8})$ values of the Z16 maps (96 models in total). The re-scaled 2PCFs are normalised to the center of the panels to exemplify the self-similar behaviour. Fig. 8 illustrates that the self-similarity of the 2PCFs is indeed robust against the change of cosmological parameters $\sigma_{8}$ and $\Omega_{m}$. The parameter space in Figs. 4 and 8 is much larger than what is allowed by current constraints. Thus a model describing this self similarity will not only have the potential to provide additional cosmological constraints, but can also be applied to scenarios where predictions for a large parameter space is required, such as generating training sets for machine learning algorithms.

Fig. 8 also shows that the rescaled 2PCFs are well described by power laws which have very similar slopes across all cosmological models. To be more quantitative, we have fitted the following power-law function,

[TABLE]

to the curves in each of the panels, with each data point weighted by its standard error (see Appendix C for details about the error calculation). The best-fitting slope, $\beta$, as a function of cosmological parameters is shown in Fig. 9, where the horizontal axis shows $\sigma_{8}$ while the colour of the points indicates the $\Omega_{m}$ value. The scatter of the points in Fig. 9 does not follow any clear trends, and combinations of $\sigma_{8}$ and $\Omega_{m}$ in the form of $\sigma_{8}(\Omega_{m}/0.3)^{\alpha}$, where $\alpha$ is allowed to vary, does not to lead to improvements in the correlation between $\beta$ and the input cosmology. In Fig. 9 we find that the mean value across the entire sample is $\beta\approx-1.1$. As such, for simplicity, we take $\beta=-1.1$ as the power-law slope of the rescaled 2PCF (we have checked that the results are not particularly sensitive to the value of $\beta$).

On the other hand, we find that the amplitude of the rescaled peak 2PCFs, $\xi_{0}=\xi_{0}(\Omega_{m},\sigma_{8})$ in Eq. (13), shows a systematic dependence on the cosmological parameters. Therefore, in order to have a complete description we also need to model $\xi_{0}(\Omega_{m},\sigma_{8})$. The fitting result for $\xi_{0}$ for the selected cosmological models (orange points in Fig. 4) is displayed in the left panel in Fig. 10, where it is plotted against $\sigma_{8}(\Omega_{m}/0.3)^{0.501}$ with the index $\alpha=0.501$ characterising the degeneracy direction between $\Omega_{m}$ and $\sigma_{8}$ for the rescaled 2PCF amplitude $\xi_{0}$. The value $0.501$ is tuned such that the data points on the left panel of Fig. 10 are fitted using a smooth quadratic curve with the lowest $\chi^{2}$. This is shown as the blue solid line in the left panel of Fig. 10, which takes the form

[TABLE]

where

[TABLE]

and $x=\sigma_{8}(\Omega_{m}/0.3)^{0.501}$. The lower sub-panel of the left panel shows the residual between the measured amplitude $\xi_{0}$ and its fitted values. The residuals show no systematic trends with varying $\sigma_{8}(\Omega_{m}/0.3)^{0.501}$, indicating that the fitting function works equally well for all cosmologies.

In the middle panel of Fig. 10 we have randomly selected one of the cosmologies from the Z16 maps, and compared the rescaled 2PCFs at several $\nu_{\mathrm{cut}}$ values between $1.5$ and $3.5$ (coloured lines), the power-law fit to these rescaled 2PCFs (black solid line, which we call the ‘fitted’ curve), and the predicted rescaled 2PCF for this particular cosmology (the grey straight line, which we call the ‘prediced’ curve). The latter was obtained by calculating $\xi_{0}$ using Eq. (14), and then inferring the 2PCF from Eq. (13). This matches the original fitted power law very closely, indicating that the model described by Eqs. (13) and (14) works very well.

We next quantify the accuracy of our prediction for the rescaled 2PCF. For a given cosmological model, such as the one shown in the middle panel of Fig. 10, we calculate the residuals, i.e., the fractional differences of the ‘fitted’ and ‘reconstructed’ curves with respect to the rescaled measured 2PCFs. This is done for each of the five $\nu_{\mathrm{cut}}$ values shown in Fig. 10, and we define the mean residual as the average over all $\theta/\theta_{p}$ bins and all $\nu_{\mathrm{cut}}$ values. The mean residuals for the fitted and predicted curves are respectively shown by a large triangle and a large circle in the right panel of Fig. 10. We have repeated this procedure for all the cosmological models and have plotted their residuals in the right panel, with the associated $\sigma_{8}$ values shown in the x-axis and $\Omega_{m}$ values shown by the colour bar to the right. We find that the model prediction is almost as accurate as the direct fitting, and is able to match the rescaled 2PCFs at about 5% accuracy level. The large symbols in the right panel correspond to the model shown in the middle panel, to give a visual illustration about how well the 2PCF model in Eqs. (13, 14) works for an ‘average’ cosmology for which the mean residual is $4.8\%$.

5 A pipeline for 2PCF Reconstruction

We can combine the models developed in the previous section for the convergence rms fluctuation, peak abundance and rescaled peak 2PCF, to develop an integrated pipeline that allows us to predict the (un-rescaled) peak 2PCFs, $\xi(\theta)$, as a function of $\nu_{\mathrm{cut}}$. The procedure is schematically illustrated in Fig. 11 and outlined as follows:

For chosen $\Omega_{m}$ and $\sigma_{8}$ values, one can use the models to predict the peak abundance (Section 4.2) and the rescaled 2PCF (Section 4.3). 2. 2.

These two statistics are combined, using $\theta_{p}=1/\sqrt{n(>\nu_{\rm{cut}})}$, to give the 2PCF, $\xi(\theta)$, for peak catalogues with $\nu_{\rm{cut}}\in[1.5,3.5]$. 3. 3.

If needed, the above-predicted peak abundance and 2PCFs can then be expressed in terms of $\kappa$ by using the $\sigma(\Sigma_{8})$ fit in Section 4.1.

This pipeline offers a simple apparatus to make predictions of the one- and two-point statistics for intermediate ($\nu$ $\in$ [1.5,3.5]) WL peaks, which can be used (on its own or together with other cosmological probes) to constrain the parameters $(\sigma_{8},\Omega_{m})$ using observational data. It will be interesting to see if these new statistics are complimentary to other probes, such as the shear-shear correlation, when constraining ($\Omega_{m}$,$\sigma_{8}$), but this will be left for future follow up works. In the next section we will discuss further aspects which need to be checked before applying this method.

As a proof of concept, we show an example of this 2PCF reconstruction pipeline in the left panel of Fig. 12. The solid curves show the 2PCFs measured from the simulation data for an arbitrarily selected cosmology, with shaded regions showing the (under) estimated standard error (see Appendix C for more detail). The dashed lines show the predictions by our 2PCF reconstruction pipeline. We find a reasonably good agreement between the simulated and reconstructed 2PCFs, with the latter mostly lying within or just outside the (under)estimated errors bars. The second panel in Fig. 12 shows the mean percentage difference between the reconstructed and measured 2PCFs, averaged over the 5 plotted 2PCFs and all $\theta$ bins with $\nu_{\rm{cut}}\in[1.5,3.5]$ separated by a $\Delta\nu_{\rm{cut}}=0.5$ increment. The model that has been selected to exemplify the reconstruction is indicated by the large symbol in the left panel of Fig. 12, which is an "average" one in terms of the performance the reconstruction (there are many models for which the reconstruction works better). We can see that for all of the selected cosmologies, our model is able to predict the 2PCF to within a roughly $6\%$ uncertainty on average. Relative to the estimated errors bars, the quality of our reconstruction is reasonably good.

We find that generally the amplitude of the 2PCFs is overestimated for the larger $\nu_{\rm{cut}}$ catalogues. This could be a fundamental aspect of the 2PCF evolution, however due to small map sizes and low peak number densities (at approximately 7 deg*-2*) it is likely that 2PCFs with $\nu>3.5$ are biased. The true amplitude of the 2PCFs with larger $\nu_{\rm{cut}}$ could be measured more accurately with larger weak lensing maps, which we leave to further study.

6 The impact of galaxy shape noise

Up to here, we have discussed the WL peak abundance and 2PCF in a theoretical context with the aim of having a model that allows us to accurately describe and predict these statistics in an idealised situation. Whilst this theoretical model can have useful applications in, e.g., mock WL peak catalogue generation, to be more useful for cosmological constraints, we need to investigate the self similarity of the 2PCF in more realistic situations. One of the things we have not included in our analysis so far is galaxy shape noise (GSN).

GSN is a source of uncertainty in WL observations, where the measured ellipticity of galaxies is dominated by their random orientation, and only weakly correlated due to gravitational lensing on scales much larger than the galaxy-galaxy separation. Observations of cosmic shear, and therefore cosmic convergence, is contaminated by this noise. One usually uses large smoothing lengths to suppress this noise in order to recover statistics more reliably. However, large smoothing lengths could either dampen the amplitude of the measured statistics, which is evident from the decrease of the WL peak abundance with increasing smoothing scales in Fig. 2, or increase the noise in the measurements, which can be seen to a small extent in the 2PCFs for different smoothing lengths in Fig. 3. Therefore, a trade-off has to be struck between using a large enough smoothing length in order to suppress the galaxy shape noise, and not over smoothing so that interesting statistics are not suppressed more than they need to be. With convergence maps from N-body simulations, we can test the difference in the peak abundance and 2PCF for cases with and without GSN for a range of smoothing lengths.

For this section, we include GSN in the T17 convergence maps that match LSST specifications by adding to each pixel within a map random values drawn from a Gaussian distribution with a standard deviation $\sigma_{\rm{pix}}$ given by

[TABLE]

where $\sigma_{\rm{int}}$ is the dispersion of the intrinsic source galaxy ellipticity, $\theta_{\rm{pix}}$ is the angular width of the pixel to which noise is added and $n_{\rm{gal}}$ is the number density of source galaxies. To match LSST specification we use $\sigma_{\rm{int}}=0.4$ and $n_{\rm{gal}}=40$ arcmin*-2* (LSST Science Collaboration et al., 2009)

After GSN is added to the pixels, we smooth the maps, identify peaks in the noise-added smoothed WL maps using Eq. (6) with $\sigma$ also directly measured from the noisy maps, recalculate the peak abundance and peak 2PCFs, and compare these statistics to the case with no GSN, with the same smoothing.

The impact of GSN on the WL peak abundance is shown in Fig. 13, where the solid and dashed lines respectively correspond to peaks identified in WL maps with and without GSN. Here we study four smoothing scales, $\theta_{s}=1$, $2$, $2.5$ and $3$ arcmin. In each instance, $\nu$ is defined relative to the WL map in which the peaks are identified, so for the GSN added case, $\sigma$ in Eq. (6) includes contributions to the rms fluctuations from both GSN and the underlying convergence signal, while for the no GSN case $\nu$ is defined by taking $\sigma$ as the rms convergence fluctuation.

For all smoothing lengths, by adding GSN, the peak abundance increases at low $\nu$ and decreases at high $\nu$, with a crossover between $\nu=1.5$ and $2.5$ depending on the smoothing scale. GSN has the largest impact on the peak abundance for the smallest smoothing length, while for larger $\theta_{s}$ the agreement between the peak abundances in the GSN and no GSN cases is better, although substantial difference remains even in the case of $\theta_{s}=3$ arcmin. This means that the fitting formulae, Eqs. (10, 11), which describe the cosmology dependence of peak abundance, need to be recalibrated by using peaks extracted from GSN-added maps. Due to the small size of the Z16 WL maps, this will be left as future work when larger simulations of different cosmologies are available.

Note that in Fig. 13 the peaks are defined using Eq. (6), where $\sigma$ is the total rms convergence that includes contributions from the physical rms convergence and from the rms of noise. This explains the crossover mentioned above: because $\sigma$ is increased, for the high peaks their $\nu$ values actually decrease, and the number of such high peaks does not increase quickly enough to maintain $n(<\nu)$ at large $\nu$, which causes the latter to drop compared with the no GSN case. We have explicitly checked (not shown here) that, if one defines $\nu$ in Eq. (6) by using the same $\sigma$ for the GSN and no GSN cases, then the peak abundance is consistently higher in the former case, due to artificial peaks created by noise.

In order to closely inspect the impact of using different smoothing lengths on the self similarity of the rescaled 2PCFs, we have tried four different $\theta_{s}$ values, respectively $\theta_{s}=1$, $2$, $2.5$ and $3$ arcmin. The results are shown in Fig. 14, where the peaks are all identified from, and the $\sigma$ used to define the SNR $\nu$ are all measured by using, the smoothed noisy maps. Interestingly, we find that the rescaled 2PCFs are still on top of each other for all four smoothing lengths. With $\theta_{s}=1~{}\mathrm{arcmin}$, the agreement between the rescaled 2PCFs is weaker, where only the curves with $2\leq\nu_{\rm{cut}}\leq 3.5$ appear to be self similar. For $2~{}\mathrm{arcmin}{}$ smoothing the 2PCFs appear to be self similar in the entire $1.5\leq\nu_{\rm{cut}}\leq 3.5$ range, and shows that the self-similarity of 2PCFs is robust against GSN. With $2.5~{}\mathrm{arcmin}{}$ smoothing the overall self similarity appears to be tighter, however the $\nu_{\rm{cut}}=1.5$ appears to be outside the self similar range. Finally, for $3~{}\mathrm{arcmin}{}$ smoothing we see that the self similarity of the 2PCFs holds up to $\nu_{\rm{cut}}=3$, after which the rescaled 2PCFs drop off in amplitude. It is possible that this drop in amplitude is caused by the small map size ($10\times 10~{}\mathrm{deg}^{2}$) and low number density of tracers ($\approx$ 0.5 deg*-2*), rather than a breakdown of the self similarity. As $\theta_{s}$ increase it also appears that the overall gradient of the rescaled 2PCFs decreases.

Having verified that the 2PCFs remain self similar in the presence of GSN, next we want to see how including the latter affects the power law of the rescaled self-similar peak 2PCFs. In each subpanel of Fig. 14, we have overplotted, as the grey dashed lines, the best-fit power-law functions for the rescaled 2PCFs of the peaks extracted from the T17 maps smoothed using the same $\theta_{s}$ values but without adding GSN (the grey dashed lines in the top two panels of Fig. 14 are the same as the black dashed lines in Fig. 3). The two cases are in good agreement for all four smoothing scales, which shows that the impact of GSN on the rescaled 2PCF is minor. This is a nice property, since it indicates that GSN will not significantly contaminate the underlying cosmology dependence of the rescaled 2PCF if the same observation applies to other cosmologies. However, due to the limited map size from Z16 we leave this investigation to future study.

In short, we conclude that the prevalence of the self similarity in the 2PCFs for peaks extracted from GSN-added WL maps shows that this feature is robust to this observational systematic, and therefore has the potential to be used in cosmological constraints.

7 Discussion and Conclusions

In this paper, we have investigated the one- and two-point statistics for intermediate peaks, with SNR values $\nu\in[1.5,3.5]$, from weak lensing convergence maps. These WL peaks contain useful information about the LSS formation, and the analyses of them are expected to place complementary constraints on the cosmological model. However, unlike high peaks, the intermediate WL peaks are not individually associated to the most massive dark matter structures, making the modelling of their statistical properties more challenging. To overcome this difficulty, we rely on WL convergence maps constructed from a large number of N-body simulations with varying cosmological parameters and technical specifications, to attempt to find patterns of the peak statistics and their cosmology dependence. Our main findings are summarised as follows:

• •

The rms fluctuation of WL convergence, $\sigma$, has a linear dependence on a particular combination of $\Omega_{m}$ and $\sigma_{8}$ via $\sigma_{8}\left(\Omega_{m}/0.3\right)^{\alpha}$, with the parameter $\alpha$ weakly dependent on the smoothing length of the convergence map, cf. Fig. 5. This linear dependence is given in Eq. (8), and highlights a universal behaviour within $\Lambda$CDM which may be exploited to make cosmological constraints.

• •

A universal one-parameter power law function is found, which can describe the WL peak count for $\nu\in[1.5,3.5]$ with an accuracy of within $\approx 1$-$5\%$, for a large range of $\Omega_{m}$ and $\sigma_{8}$ values, cf. Fig. 2 and Eq. (10). The accuracy of the power-law description of the peak abundance can reach the sub-percent level if two free parameters are used in the power-law function.

• •

A self-similar behaviour of the WL peak 2PCF has been found by rescaling the angular separation, $\theta$, between a pair of peaks by the mean inter-peak separation, $\theta_{p}$. While the amplitude of the original 2PCF increases with $\nu_{\mathrm{cut}}$, the rescaled 2PCFs for $\nu_{\mathrm{cut}}\in[1.5,3.5]$ lie on top of each other cf. Fig. 3.

• •

This self-similar behaviour holds for a very wide range of $(\Omega_{m},\sigma_{8})$ values, and we find a simple quadratic dependence of the amplitude of the rescaled 2PCFs on $\sigma_{8}\left(\Omega_{m}/0.3\right)^{\alpha}$, while the slope of the rescaled 2PCFs have negligible dependence on $\Omega_{m}$ and $\sigma_{8}$, cf. Figs. 8, 9 and 10. A fitted model to predict the peak 2PCF for any chosen $\Omega_{m}$ and $\sigma_{8}$ is given in Eq. (13).

• •

A pipeline is developed which combines the above three fitted models, for the convergence rms fluctuation, WL peak abundance and rescaled peak 2PCF respectively, to predict the raw peak 2PCF $\xi(\theta;\nu_{\mathrm{cut}})$ for $\nu_{\mathrm{cut}}\in[1.5,3.5]$ and any given $\Omega_{m}$ and $\sigma_{8}$ with good accuracy, cf. Fig. 12.

• •

We found that the self similarity of the peak 2PCF holds in the presence of galaxy shape noise and larger smoothing lengths, cf. Fig. 14.

The most important application of the results presented in this work is in constraining the $\Omega_{m}$ and $\sigma_{8}$ cosmological parameters. As demonstrated above, the pipeline integrating the models for WL peak abundance and self-similar rescaled 2PCFs is able to reconstruct the raw, unrescaled, peak 2PCFs for various $\nu_{\mathrm{cut}}$ values with a typical accuracy of better than $6\%$. Furthermore, we have seen that the WL peak abundance and 2PCFs depend on very different combinations of $\Omega_{m}$ and $\sigma_{8}$, one with $\sigma_{8}\left(\Omega_{m}/0.3\right)^{-0.638}$ and the other $\sigma_{8}\left(\Omega_{m}/0.3\right)^{0.501}$. This indicates that a simultaneous use of these statistics already holds the potential of breaking the degeneracy between $\Omega_{m}$ and $\sigma_{8}$ before including other cosmological probes. Marian et al. (2013) found that the 2PCFs of high WL peaks only provide weakly complimentary constraints on ($\Omega_{m}$,$\sigma_{8}$) when combined with the peak abundance. In this work we investigate the 2PCFs of WL peaks with intermediate heights and above, as well as combining the 2PCFs from multiple peak catalogues in the form of a rescaled 2PCF described by a single power law. The powerlaw describing the rescaled 2PCF may be more sensitive to cosmology than the 2PCF of high peaks.

We note that the degeneracy direction of the peak abundance of intermediate height peaks, which are studied in this paper, are very different to that of low and high peaks, which has also been observed in Liu et al. (2015a) and explained in Yang et al. (2011). Therefore, using the counts of intermediate height peaks may be complimentary to using the full peak abundance and could aid in breaking the $\Omega_{m}$ and $\sigma_{8}$ degeneracy.

Another potential application of our results is the generation of mock WL peak catalogues. For a given input cosmological model, the pipeline can be used to predict the WL peak counts and 2PCFs as described above. Random realisations of peaks can then be generated with the peaks arranged such that they have the desired number density and spatial clustering. One technique to do this is point process (see, e.g., Öztireli & Gross, 2012, for some recent progress and applications). This is a Monte Carlo approach where a candidate point (e.g., a WL peak) is placed in a field, which is accepted if its inclusion into the field pushes the measured 2PCF closer to the input one, and rejected otherwise. Point process is a well-developed and widely-used technique to generate point catalogues. In the WL peak case, the situation is slightly more complicated, because the generated catalogue should have peaks of different SNR (or $\nu$ values), which simultaneously have the desired 2PCFs at different $\nu_{\mathrm{cut}}$ values. We expect that the good agreement between the rescaled peak 2PCFs will prove useful in dealing with this issue, though a detailed investigation into this interesting question will be left for a future work. The fast generation of mock WL peak catalogues can be used for evaluating covariance matrices and studying other cosmological quantities, such as voids identified from WL peaks (Davies et al., 2018).

The proof-of-concept study in this work has also left various possible further extensions of the analyses presented here. One of the most important considerations for future WL surveys and their cosmological applications is the effect of galaxy shape noise. Using the all-sky maps from T17, we have shown that (i) the inclusion of GSN necessitates a larger smoothing length than used in the bulk of this paper, $\theta_{s}=2$-$3$ arcmin, to suppress its impact on the extracted cosmological statistics, and (ii) with a suitable smoothing, the self-similarity of the peak 2PCFs still holds for the cosmology used in the T17 simulations. While we expect these conclusions to apply for other cosmological models, Fig. 14 shows that the use of GSN and larger $\theta_{s}$ does indeed affect the slope of the rescaled peak 2PCF. Therefore, in the presence of GSN our fitted models need to be re-analysed before it can be directly useful for cosmological tests.

Unfortunately, the 96 Z16 maps with varying cosmological parameters have a relatively small size, at $3.5\times 3.5$ deg2. Including GSN in these maps and increasing the smoothing length will reduce both the number of peaks in the maps and the dynamical range over which the 2PCFs can be reliably studied. This consideration makes a compelling case that larger convergence maps, constructed from N-body simulations with larger boxes and varying cosmologies, are a natural next step, to re-calibrate our peak models so that they can be readily applied for upcoming WL surveys. Again, we leave these to a future, more comprehensive, study.

The planned larger simulations will have other applications as well. For example, they will allow us to study low/intermediate WL peaks and the high peaks, as well as other statistics such as the WL shear power spectrum, simultaneously. It will also be possible to look at source galaxies with a certain redshift distribution compared to the currently idealised case with a single source redshift, $z_{s}=1$. Larger WL maps will also allow us to more accurately estimate the errors on the 2PCFs, with large-scale modes properly included. Further more, in future studies we will try methods of extracting WL peaks that are more similar to approaches taken in observations, such as starting with the shear field and adding GSN to this before we then transform to the convergence field.

Finally, it will also be interesting to analyse the rescaled WL peak 2PCFs in cosmological models beyond $\Lambda$CDM. We can envisage two possible scenarios here: the first is that the rescaled 2PCFs may not be self similar, which would offer a potentially strong constraint on these models. Alternatively, the detailed properties of the self-similarity in the 2PCFs may change, in the form of a different amplitude or slope, which can also be used to test models with observational data. Therefore, it will be important to consider models which are expected to alter the large-scale clustering of matter. These include the various dark energy models which may couple to dark matter or have different equation-of-state $w$ parameters. The neutrino mass is another interesting possibility, as massive neutrinos tend to dampen structure formation, which leaves signatures in the WL peak abundance and 2PCF. Modified gravity models can also be potentially tested since they generally introduce fifth forces on cosmological scales, which modify the clustering of matter or even the geodesics of photons. The studies of these topics will require new simulations and dedicated effort, and will be deffered for the future.

Acknowledgements

We thank Zuhui Fan, Daniel Gruen, Zoltan Haiman, Catherine Heymans, Chieh-An Lin, Jia Liu, Xiangkun Liu, Peder Norberg, Masato Shirasaki, for useful discussions and comments on this draft.

CTD is funded by a Science and Technology Facilities Council (STFC) PhD studentship through grant ST/R504725/1. MC acknowledges support by the European Research Council through an ERC Advanced Investigator Grant, DMIDAS [GA-786910]. BL is supported by an ERC Starting Grant, ERC-StG-PUNCA-716532. MC and BL are additionally supported by the STFC Consolidated Grants [Nos. ST/I00162X/1, ST/P000541/1]

This work used the DiRAC facility, hosted by Durham University, managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K00042X/1, ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University and STFC operations grant ST/R000832/1. DiRAC is part of the National e-Infrastructure.

We thank the Columbia lensing group (http://columbialensing.org) for making their suite of simulated lensing maps available, and the US National Science Foundation (NSF) for supporting the creation of those maps through grant AST-1210877 and XSEDE allocation AST-140041.

Appendix A Partitioning an all-sky map into smaller non-overlapping maps

The T17 maps are all-sky maps with a HEALPix data structure, in which pixels are stored on the surface of a sphere. To simplify our analysis, we used the flat-sky approximation and thus we needed to partition the all-sky map into smaller, and preferably, non-overlapping maps.

To achieve this, we capitalise on the HEALPix data structure and first define a set of coarser HEALPix pixels with a resolution of $N_{\rm{side}}=4$, which corresponds to a pixel area of roughly $215~{}\mathrm{deg}^{2}$. We then assign each of the (higher resolution) data pixels to the coarser pixel that they are enclosed by. This is shown by the illustration on the left in Fig. 15 (using Mollweide projection), where each coloured patch shows a course pixel. Next, for each sub region defined by the coarse pixels, we define a (flat) plane tangential to the centre of the coarse pixels and project the data pixels onto that plane. We then extract a square of $10\times 10~{}\mathrm{deg}^{2}$ (centred on the centre of the plane) from each plane giving us 184 $10\times 10~{}\mathrm{deg}^{2}$ flat maps. The HEALPix pixels that are projected onto the flat maps are converted into regular square pixels, where we interpolate between HEALPix pixels for square pixels that overlap with multiple HEALPix pixels. The benefit of this approach is that there is no overlap between any two maps as illustrated in the right panel of 15.

We note that a HEALPix resolution of $N_{\rm{side}}=4$ actually gives 192 pixels, however due to the irregular shapes of HEALPix pixels (which arises from the requirement that all pixels have the same area), we find that 8 pixels have to be discarded since they cannot enclose squares of size $10\times 10~{}\mathrm{deg}^{2}$.

Appendix B Biased 2PCF estimation for small maps

Estimation of 2PCFs is straightforward in idealised situations. The 2PCF, $\xi(r)$, characterises the excess probability of finding a pair of tracers in two volume elements, ${\rm d}V_{i}$ and ${\rm d}V_{j}$, that are separated by a distance $r$:

[TABLE]

where $\bar{n}$ represents the expected tracer number density. In N-body simulations with periodic boundary conditions, as an example, $\bar{n}$ is the known mean number density and so the excess probability ${\rm d}P_{ij}$ can be evaluated by counting the number of pairs that are separated by a distance $r-\Delta r$ to $r+\Delta r$ and comparing that against $\bar{n}^{2}$. In realistic situations, $\bar{n}$ is not always known – this can for example be due to the geometry, mask, fibre collision and redshift failure in a galaxy redshift survey, or the small map size with boundaries in our WL peak catalogues. The uncertainty in the expected number of tracers in a given volume can cause biased 2PCF estimations. It is known that, for examples, the Peebles & Hauser (1974) estimator

[TABLE]

and the Davis & Peebles (1983) estimator

[TABLE]

have errors that depend to the first order on the uncertainty of the expected tracer number density. On the other hand, the Hamilton (1993) and the Landy & Szalay (1993) estimators have errors which are second order in this uncertainty and are more commonly used. In the above $N_{D}$ and $N_{R}$ are the numbers of data and random points, and $DD$, $DR$ and $RR$ are the numbers of data-data, data-random and random-random pairs in bins $\theta\pm\delta\theta$ respectively.

The Landy-Szalay estimator is given by,

[TABLE]

When analysing $n$ maps, there are $n$ different $N_{D}$, $DD$ and $DR$ values, that is one per map ($N_{R}$ and $RR$ can be taken as constants since the same random catalogue can be used for each map). We checked that our results are stable to a change in the number of randoms used.

Given the expression of Eq. (20), there are two possible ways to calculate the mean 2PCF, $\langle\xi\rangle$, where $\langle\cdot\rangle$ denotes the mean value over the $n$ maps, given respectively by

[TABLE]

and

[TABLE]

where we have dropped the $\theta$ dependence of $DD,DR$ and $RR$ to lighten the notations. For Eq. (21), we calculate $n$ 2PCFs from the $n$ maps and take the mean value. For Eq. (22), we first calculate the mean over all maps of $N_{D}$, $DD$ and $DR$, and then use these mean values to calculate the mean 2PCF. In general, these two approaches do not give identical results, that is

[TABLE]

Naively, it seems natural to calculate the mean 2PCF using Eq. (21) – after all, if we only had a single map, we would use this formula (excluding the outer $\langle\cdot\rangle$) to estimate the 2PCF. However, we found that this approach actually leads to biased estimates when the number of tracers varies considerably between the different maps. This is particularly the case when the number of peaks in a catalogue is low, and the effect is particularly strong for peaks with high $\nu$ values, for which the number density is low, and for small maps, such as the $3.5\times 3.5~{}\mathrm{deg}^{2}$ ones.

To see this point, we compare these two approaches as follows. First, we split each of the 184 T17 maps, whose size is $10\times 10~{}\mathrm{deg}^{2}$, into four $5\times 5~{}\mathrm{deg}^{2}$ maps and nine $3.3\times 3.3~{}\mathrm{deg}^{2}$ maps, which give us in total 184, 736 and 1656 maps of the three sizes respectively. Then, using respectively Eqs. (21) and (22), we calculate the mean 2PCF for the three different maps sizes, and the results are shown in Fig. 16.

Fig. 16 clearly shows that, as the T17 maps are split into progressively smaller sections, the mean 2PCF calculated using Eq. (21) drops in amplitude, whereas using Eq. (22) leads to a constant amplitude. The difference between the two approaches is small for the $10\times 10~{}\mathrm{deg}^{2}$ maps and only becomes significant for the smaller maps. This implies that the bias from Eq. (21) depends on the map size, or more exactly the number of tracers used for the 2PCF estimation. We have performed similar tests for 3D galaxy 2PCFs and found a similar bias effect when using small box sizes. Finally, the mean 2PCFs from Eq. (22) for the different maps sizes do not line up exactly, which is due to some pairs being lost at the sub-map boundaries as large maps are split up into smaller maps. We checked for this and found that the inclusion of cross sub-map pairs restores the original 2PCFs.

Physically, the reason why Eq. (21) leads to biased 2PCF estimations is that the number of WL peaks per map is small and this translates into a large uncertainty in the mean tracer number density when estimated individually for each map. Even though this uncertainty enters the 2PCF estimation only at second order for the Landy-Szalay estimator, it can still strongly affect the latter. In contrast, Eq. (22) essentially treats the $n$ maps as a single (combined) one, for which the uncertainty in the expected mean peak number is small.

The biased 2PCF estimation using the Landy-Szalay estimator caused by the small tracer number is important for this study, since the WL maps from the Z16 simulations have a map size of $3.5\times 3.5~{}\mathrm{deg}^{2}$, which is in the regime where the biasing effect is strong. As a result, in this paper we calculate the mean 2PCF using Eq. (22).

Appendix C Error estimates

For each of the Z16 cosmologies we used the $N=$ 512 $3.5\times 3.5$ deg2 maps to evaluate the mean 2PCF $\langle\xi\rangle$. We estimate the standard error for $\langle\xi\rangle$ using the jackknife method, by calculating $N-1$ mean values from sequentially removing individual maps from the sample, and taking the standard deviation of the $N-1$ means from the 512 maps. However, we found the error to be significantly smaller than expected, of roughly $<1\%$ of $\langle\xi\rangle$ itself. On the other hand, when repeating the same practice on 512 $3.5\times 3.5$ deg2 maps extracted from the T17 all-sky map, we found the standard errors were larger and more reasonable, of roughly $2$-$3\%$ of the mean 2PCF.

This discrepancy in the magnitudes of the standard errors in the two different suites of maps is likely caused by the way in which the multiple convergence maps were generated. In the Z16 case, the 512 maps were generated from multiple lines of sight by shifting, reorienting and tiling a single simulation box of size $240h^{-1}$Mpc, which means that the scatter in the different maps is likely to only contain an error component representing the line-of-sight variation. In contrast, the T17 maps were all-sky maps generated using much larger boxes with minimal repetition of structures along the lines of sight, which means that these maps better sample the variation due to large-scale modes. The additional source of variance in the T17 maps can explain the increase in their measured standard error.

In order to have a more realistic estimate of the standard error associated to the Z16 maps, we extract 512 $3.5\times 3.5~{}\mathrm{deg}^{2}$ maps from the T17 all-sky map and use jackknife to find the error of the mean 2PCFs, $\langle\xi\rangle$. For illustration purposes, the resulting $\langle\xi\rangle$ and their errors for a few values of $\nu_{\mathrm{cut}}$ are shown in Fig. 17 as respectively lines and shaded regions.

We then take this relative error as our estimate of the standard error for the mean 2PCFs from the Z16 maps, as a way to (approximately) include the contributions to the error from large-scale modes.

The above estimate of the error associated to the Z16 maps is likely to be an underestimate since the estimated error corresponds to the case when each of the 512 Z16 maps would have been obtained from a different N-body simulation. However, this is not the case since all the Z16 maps were obtained from the same simulation. Thus, the errors used in this paper serve only as a way to gain rough indications of the quality of our models for the WL peak statistics, which we present as a proof of concept. In a future work, we plan to run a suite of large simulations similar to those used by T17, for different cosmological models, to further study the self-similar properties of the rescaled peak 2PCFs.

Appendix D The independence of self-similarity on the $\nu$ definition

In Eq. (6) we choose to define the SNR, $\nu$, in terms of a cosmology-dependent rms convergence, $\sigma$, which is analytically parameterised through a simple dependence on $(\Omega_{m},\sigma_{8})$, as exemplified in Fig. 5. Besides having a readily-predictable $\sigma$, this approach has the added benefit of allowing us to more naturally define the amplitude of WL peaks for a given cosmology relative to its own convergence rms, bearing in mind that the wide coverage of cosmological parameters means that the $\sigma$ values can vary by a factor of a few across the Z16 maps; cf. Fig. 5.

One can argue that given an observational WL map, the value of $\sigma$ receives contributions from both the physical convergence rms and the GSN, and that the actual value of $\sigma$ as measured from such noisy maps is a natural choice that can be used to define $\nu$. Such is the logic followed in Section 6 where we analysed the rescaled peak 2PCFs in the GSN-added maps. Alternatively, one may argue that in real observations we do not necessarily know the true cosmology, but we do understand the survey specifications well enough to know the expected noise level. This leads to another natural way to define $\nu$, namely by using $\sigma\equiv\sigma_{\rm GSN}$. Given this flexibility in $\nu$ definition, we would like to check that the self-similar behaviour of the 2PCFs for the resulting peak catalogues is not affected by it. This is done in Fig. 18, which is similar to Fig. 8 but for a subset of cosmological models (the ones represented by the orange symbols in Fig. 4), and where $\nu$ is defined as $\nu=(\kappa-\mu)/\sigma_{\rm{GSN}}$ with $\sigma_{\rm{GSN}}=0.013$, which corresponds to the rms of a GSN only map smoothed with $\theta_{s}=1$ arcmin, is used for all models.

We find that the cosmology-dependent description of WL peak amplitude results in a (marginal) improvement of the self similarity of the 2PCFs for all cosmologies (shown in Fig. 8) compared to using a $\sigma$ definition that is constant across all cosmologies, which is shown in Fig. 18. It can be seen that the self similarity of the 2PCFs worsens notably for some of the panels, which correspond to models with more extreme $(\Omega_{m},\sigma_{8})$ values. This is not surprising because, as mentioned above, the models studied here vary wildly in their $\sigma$ values, and by using a constant $\sigma_{\rm GSN}$ to define $\nu$ one is essentially selecting very different peak populations in them – in the more extreme models the peaks that end up being selected do not possess the self-similarity (remember that this property is only present for a limited range of peak heights). Hence, by using the cosmology dependent form of $\nu$ in Eq. (6), the 2PCF is self similar for a larger range of cosmologies. However, if one focuses on the more realistic $(\Omega_{m},\sigma_{8})$ parameters, then using a constant $\sigma_{\rm GSN}$ to define $\nu$ should not affect the potential of the rescaled peak 2PCF as a cosmological probe.