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

HAL
open science

# On the projected mass distribution around galaxy clusters: a Lagrangian theory of harmonic power spectra

Sandrine Codis, Raphael Gavazzi, Christophe Pichon, Céline Gouin

► To cite this version:

Sandrine Codis, Raphael Gavazzi, Christophe Pichon, Céline Gouin. On the projected mass distribution around galaxy clusters: a Lagrangian theory of harmonic power spectra. Astronomy & Astrophysics - A&A, 2017, 605, pp.A80. 10.1051/0004-6361/201630091. hal-01645464

HAL Id: hal-01645464

https://hal.science/hal-01645464v1

Submitted on 14 Aug 2025

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>

Distributed under a Creative Commons Attribution 4.0 International License
# On the projected mass distribution around galaxy clusters
## A Lagrangian theory of harmonic power spectra

Sandrine Codis<fnref n="1" />, Raphaël Gavazzi<fnref n="2" />, Christophe Pichon<fnref n="2" />,<fnref n="3" />, and Céline Gouin<fnref n="2" />

Received 18 November 2016 / Accepted 30 May 2017

## ABSTRACT

**Aims.** Gravitational lensing allows us to quantify the angular distribution of the convergence field around clusters of galaxies to constrain their connectivity to the cosmic web. We describe the corresponding theory in Lagrangian space in which analytical results can be obtained by identifying clusters to peaks in the initial field.

**Methods.** We derived the three-point Gaussian statistics of a two-dimensional (2D) field and its first and second derivatives. The formalism allowed us to study the statistics of the field in a shell around a central peak, in particular its multipolar decomposition.

**Results.** The peak condition is shown to significantly remove power from the dipolar contribution and to modify the monopole and quadrupole. As expected, higher order multipoles are not significantly modified by the constraint. Analytical predictions are successfully checked against measurements in Gaussian random fields. The effect of substructures and radial weighting is shown to be small and does not change the qualitative picture. The non-linear evolution is shown to induce a non-linear bias of all multipoles proportional to the cluster mass.

**Conclusions.** We predict the Gaussian and weakly non-Gaussian statistics of multipolar moments of a 2D field around a peak as a proxy for the azimuthal distribution of the convergence field around a cluster of galaxies. A quantitative estimate of this multipolar decomposition of the convergence field around clusters in numerical simulations of structure formation and in observations will be presented in two forthcoming papers.

**Key words.** large-scale structure of Universe – galaxies: clusters: general – gravitational lensing: weak – methods: analytical – methods: statistical

## 1. Introduction

Galaxies are not islands uniformly distributed in the Universe. In recent decades and with the increasing precision of both observations and simulations, galaxies have been shown to reside in a complex network made of large filaments surrounded by walls and voids and intersecting at the overdense nodes of this so-called cosmic web (Klypin & Shandarin 1993; Bond et al. 1996). From the pioneering works of Zeldovich in the 1970s to the peak-patch picture of Bond & Myers (1996), the anisotropic nature of the gravitational collapse has been used to explain the birth and growth of the cosmic web. The origin of filaments and nodes lies in the asymmetries of the initial Gaussian random field describing the primordial universe and amplified by gravitational collapse. The above-mentioned works pointed out the importance of non-local tidal effects in weaving the cosmic web. The high-density peaks define the nodes of the evolving cosmic web and completely determine the filamentary pattern in between. In particular, one can appreciate the crucial role played by the study of constrained random fields in understanding the geometry of the large-scale matter distribution.

Galaxy clusters sitting at these nodes are continuously fed by their connected filaments (e.g. Aubert et al. 2004, and reference therein; see also Pogosyan et al., in prep., for a study of the connectivity of the cosmic web). The key role played by this anisotropic environment in galaxy formation is increasingly

underlined. For instance, it has been observed that the properties of galaxies – morphology, colours, luminosities, spins among others – are correlated to their large-scale environment (see Oemler 1974; Guzzo et al. 1997; Tempel & Libeskind 2013; Kovač et al. 2014, among many others).

Numerical simulations allow us to study the details of this large-scale structure of the Universe together with its impact on the formation and evolution of galaxies. Using N-body simulations, Hahn et al. (2007), Gay et al. (2010), Metuki et al. (2015) found that the properties of dark matter halos, such as their morphology, luminosity, colour, and spin parameter, depend on their environment as traced by the local density, velocity, and tidal field. In addition to scalar quantities, it also appears that their shape and spin are correlated to the directions of the surrounding filaments and walls both in dark matter (see for instance Aubert et al. 2004; Bailin & Steinmetz 2005; Brunino et al. 2007; Aragón-Calvo et al. 2007; Sousbie et al. 2008; Paz et al. 2008; Codis et al. 2012; Aragon-Calvo & Yang 2014) and hydrodynamical simulations (Navarro et al. 2004; Hahn et al. 2010; Dubois et al. 2014).

Analytical works provide important insights to understand the results of those simulations in the quasi-linear regime. As already pointed out, the theory of constrained random fields is an important tool that allows analytical calculations in the linear or weakly non-linear regime, which is effective at large scales or early times in the Universe. Virialised halos are the highly

<footnote n="1">¹ Canadian Institute for Theoretical Astrophysics, University of Toronto, 60 St. George Street, Toronto, ON M5S 3H8, Canada
e-mail: codis@cita.utoronto.ca</footnote>

<footnote n="2">² Institut d'Astrophysique de Paris, UMR 7095 CNRS & Université Pierre et Marie Curie, 98 bis Bd Arago, 75014 Paris, France</footnote>

<footnote n="3">³ Korea Institute of Advanced Studies (KIAS) 85 Hoegiro, Dongdaemun-gu, 02455 Seoul, Republic of Korea</footnote>
non-linear result of gravitational dynamics. They tend to form in the high-density peaks of the density field by gravitational instability and as such represent a biased tracer of the density field (Kaiser 1984; Bardeen et al. 1986). Peak statistics has focused a lot of attention in the recent years as it provides a unique way to analytically study the statistics of halos from their spatial distribution to their mass function (Paranjape & Sheth 2012) or their spin (Codis et al. 2015), at least for rare enough objects (Ludlow & Porciani 2011).

Despite clear evidence from numerical simulations, the detection of filaments and cold flows is still a debated but crucial issue as filamentary flows are often depicted as the solution to the missing baryons problem (Persic & Salucci 1992; Fukugita et al. 1998; Davé et al. 2001; Shull et al. 2012). In particular, gravitational lensing has emerged as a potential powerful probe of the filamentary cosmic web despite being challenging because of the systematics and the weakness of the signal (Dietrich et al. 2005; Mead et al. 2010; Martinet et al. 2016).

Gravitational lensing is related to the projected density integrated along the line of sight from distant source to the observer. The so-called convergence $\kappa$ is proportional to the projection of the density contrast $\delta$, and, as such, it inherits its statistical properties. In particular, projection tends to wash the non-Gaussianities of the $\delta$ field out. One would therefore try and enhance the importance of the filamentary structure by looking at the statistical properties of the convergence field at the vicinity of the rarest, most singular events, which are the clusters at the nodes of the web. In this work, we quantify the amount of symmetry of the matter distribution around clusters of galaxies by means of the aperture multipolar moments of the convergence field (Schneider & Bartelmann 1997) and their power spectrum. In particular, this tool should allow us to detect the signature of filaments feeding galaxy clusters in weak lensing surveys. This paper aims to investigate the theory of this observable in the Gaussian regime while a companion paper (Gouin et al. 2017) explores the fully non-linear regime by analyzing clusters of galaxies within cosmological N-body simulations.

This works complements in two dimensions (2D) the three-dimensional (3D) harmonic analysis of infall at the virial radius presented in Aubert & Pichon (2007). The paper is organised as follows. Section 2 describes the mathematical formalism from the general definition of multipolar moments to the statistical description of peaks in Gaussian random fields (GRF) and their impact on the statistics of the multipolar moments. Section 3 then compares the predictions to measurements in Gaussian random fields. Section 4 studies the effect of substructures and Sect. 5 adds a generic radial weight function. We describe the weakly non-linear evolution of the multipolar moment in Sect. 6. Finally, we give preliminary conclusions of this work in Sect. 7 and propose possible follow-up developments. A statistical characterization of the geometry of peaks for 2D Gaussian random fields is given in Appendix A.

## 2. Formalism

### 2.1. Aperture multipolar moments

The focus of this paper lies in the azimuthal mass distribution at various scales around massive galaxy clusters. For a thin gravitational lens plane, the convergence $\kappa$ at a given position $\mathbf{r}$ in the sky corresponds to the projected excess surface density expressed in units of the so-called critical density $\Sigma_{\text{crit}}$

$$ \kappa(\mathbf{r}) = \frac{1}{\Sigma_{\text{crit}}} \int dz (\rho(\mathbf{r}, z) - \bar{\rho}), \quad (1) $$

with the convention that the line of sight corresponds to the z-axis and the plane of the sky $\mathbf{r}$ vector can be defined by polar coordinates $(r, \varphi)$. The critical density involves distance ratios between a fiducial source at an angular diameter distance $D_s$, the distance to the lensing mass $D_l$, and the distance between the lens and the source $D_{ls}$

$$ \Sigma_{\text{crit}} = \frac{c^2}{4\pi G} \frac{D_s}{D_l D_{ls}}. \quad (2) $$

On cosmological scales, the thin lens approximation is generally not valid and the integrated deflections experienced by light rays as they travel from the source to the observer requires numerical treatment, but for most cosmological applications the integration of the deflections along the unperturbed light rays (so-called Born approximation; see e.g. Bartelmann & Schneider 2001) yields a linear integral relation between the convergence $\kappa$ and the density contrast $\delta$. For a known time-varying<fnref n="1" /> 3D power spectrum $P_\delta(k, \chi)$, and for a given source plane redshift $z_s$, one can thus write the convergence power spectrum $P_\kappa(\ell, z_s)$ by means of the Limber approximation (Blandford et al. 1991; Miralda-Escudé 1991; Kaiser 1992; Bartelmann & Schneider 2001; Simon 2007)

$$ P_{\kappa}(\ell, z_s) = \frac{9}{4} \Omega_m^2 \left(\frac{H_0}{c}\right)^4 \int_0^{\chi_s} d\chi \frac{(\chi_s - \chi)^2}{\chi_s^2} \frac{P_\delta(\ell/\chi, \chi)}{a^2(\chi)}. \quad (3) $$

Following early works by Schneider & Bartelmann (1997), we define the aperture multipolar moments of the convergence (projected surface mass density) field $\kappa$ as

$$ Q_m = \int_0^\infty dr \, r^{1+m} w_m(r) \int_0^{2\pi} d\varphi \, e^{im\varphi} \kappa(r, \varphi), \quad (4) $$

with a radial weight function $w_m(r)$ commonly defined on a compact support. Those multipoles aim to quantify possible asymmetries in the mass distribution as probed by gravitational lensing.

The covariance between multipolar moments can straightforwardly be written as

$$ \langle Q_n Q_m^* \rangle = 2\pi i^{n-m} \int kdk U_n(k) U_m(k) P(k), \quad (5) $$

where $U_n(\ell)$ is the Hankel transform of the radial weight function

$$ U_m(\ell) = \int r dr \, r'^m w_m(r) J_m(\ell r), \quad (6) $$

$J_m(x)$ are the first kind Bessel functions and $P(k)$ is the power spectrum of the 2D random field $\kappa$.

In a suite of papers (including Gouin et al. 2017, and Gavazzi et al., in prep.), we propose using the full statistics of these multipolar moments around clusters of galaxies. The covariance of the aperture multipolar moments in specific locations of space, such as the vicinity of clusters, becomes

$$ \langle Q_n Q_m^* | \text{clusters} \rangle = \int_0^\infty r dr \, r' dr' \int_0^{2\pi} d\varphi d\varphi' \, r'^n w_n(r) r'^m w_m(r') \\ \times e^{i(n\varphi - m\varphi')} \langle \kappa(r, \varphi) \kappa(r', \varphi') | \text{clusters} \rangle, \quad (7) $$

<footnote n="1">¹ Where time variation is captured by an explicit dependence on comoving distance $\chi$.</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Left-hand panel: this paper aims at describing the angular distribution of a 2D Gaussian field $\kappa$ around a peak in $r_z$. We therefore consider two points on the annulus at a distance $r$ from the central peak. Their respective Cartesian coordinates are $\mathbf{r}_x = r(\cos \theta, \sin \theta)$ and $\mathbf{r}_y = r(\cos \theta + \psi, \sin \theta + \psi)$. Right-hand panel: example is shown of such a 2D Gaussian random field with a power-law power spectrum with spectral index $n_s = 0$. Peaks of height $\nu = 3 \pm 0.5$ are highlighted with black dots. We hereby investigate the polar distribution of the field around such peaks.</figcaption></figure>

where $\langle \kappa(r, \varphi)\kappa(r', \varphi')|_{\text{clusters}} \rangle$ is a constrained two-point correlation function as we impose a cluster at the origin of the polar coordinate system.

In order to develop a physical intuition of the effect of this cluster constraint on the statistics of the multipolar moments, we propose in this paper to study analytically this observable for a Gaussian random field in which clusters are identified as high peaks. To simplify the problem, we drop the radial weight function and focus on Gaussian random fields smoothed with a Gaussian kernel on a given scale $R$. In what follows, we investigate the angular distribution of a Gaussian random field around a peak. We therefore need to study the joint statistics of the field in three locations of space: the location of the peak and two arbitrary points on the circle at a distance $r$ away from the central peak. In addition, according to the peak theory originally developed in Bardeen et al. (1986), we need to consider the field and its first and second derivatives at the location of the peak. In Sect. 2.2, we first present the result for the joint PDF of those random variables before computing the resulting multipolar decomposition around a central peak in Sect. 2.6.

## 2.2. Three-point statistics of the field and its derivatives

For a given 2D Gaussian field $\kappa$ (for example, the projected density contrast), we define the moments

$$ \sigma_0^2 = \langle \kappa^2 \rangle, \sigma_1^2 = \langle (\nabla \kappa)^2 \rangle, \sigma_2^2 = \langle (\Delta \kappa)^2 \rangle. \quad (8) $$

From these moments, we use two characteristic lengths $R_0 = \sigma_0/\sigma_1$ and $R_\star = \sigma_1/\sigma_2$, as well as the spectral parameter

$$ \gamma = \frac{\sigma_1^2}{\sigma_0 \sigma_2}. \quad (9) $$

We now define the following normalised random variables

$$ x = \frac{1}{\sigma_0 \kappa}, \quad x_i = \frac{1}{\sigma_1 \nabla_i \kappa}, \quad x_{ij} = \frac{1}{\sigma_2 \nabla_i \nabla_j \kappa}, \quad (10) $$

which have unit variance by construction.

In what follows, $\mathcal{P}(X)$ denotes the one-point probability density (PDF) and $\mathcal{P}(\mathbf{X}, Y, Z)$ designates the joint PDF for the normalised field and its derivatives, $X = \{x\}$, $Y = \{y\}$, and $Z = \{z, z_i, z_{ij}\}$, at three prescribed comoving locations ($\mathbf{r}_x, \mathbf{r}_y$ and $\mathbf{r}_z$) separated by a distance $r = |\mathbf{r}_x - \mathbf{r}_z| = |\mathbf{r}_y - \mathbf{r}_z|$ so that we are considering the density field in two locations, $\mathbf{r}_x$ and $\mathbf{r}_y$ on the same infinitely thin shell around the central peak in $\mathbf{r}_z$; see also the left-hand panel of Fig. 1. The right-hand panel of Fig. 1 shows a Gaussian random field and the position of its peaks. This paper aims to investigate the angular matter distribution around those peaks.

For a Gaussian field (in particular cosmic fields at early times or large scales), the joint PDF is a multivariate normal distribution

$$ \mathcal{N}(\mathbf{X}, \mathbf{Y}, \mathbf{Z}) = \frac{\exp\left[-\frac{1}{2} \begin{pmatrix} \mathbf{X} \\ \mathbf{Y} \\ \mathbf{Z} \end{pmatrix}^\mathrm{T} \cdot \mathbf{C}^{-1} \cdot \begin{pmatrix} \mathbf{X} \\ \mathbf{Y} \\ \mathbf{Z} \end{pmatrix}\right]}{\det|\mathbf{C}|^{1/2} (2\pi)^{(6+3d+d^2)/4}}, \quad (11) $$

where $d$ is the dimension – $d = 2$ here – and $\mathbf{C}$ is the covariance matrix that depends on the separation vectors only because of homogeneity

$$ \mathbf{C} = \begin{pmatrix} 1 & \langle xy \rangle & \mathbf{C}_{XZ} \\ \langle xy \rangle & 1 & \mathbf{C}_{YZ} \\ \mathbf{C}_{XZ}^\mathrm{T} & \mathbf{C}_{YZ}^\mathrm{T} & \mathbf{C}_{ZZ} \end{pmatrix}, \quad (12) $$

with

$$ \mathbf{C}_{XZ} = (\langle xz \rangle, \langle xz_1 \rangle, \langle xz_2 \rangle, \langle xz_{11} \rangle, \langle xz_{12} \rangle, \langle xz_{22} \rangle), \quad (13) $$

$$ \mathbf{C}_{YZ} = (\langle yz \rangle, \langle yz_1 \rangle, \langle yz_2 \rangle, \langle yz_{11} \rangle, \langle yz_{12} \rangle, \langle yz_{22} \rangle), \quad (14) $$

$$ \mathbf{C}_{ZZ} = \begin{pmatrix} 1 & 0 & 0 & -\gamma/2 & 0 & -\gamma/2 \\ 0 & 1/2 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1/2 & 0 & 0 & 0 \\ -\gamma/2 & 0 & 0 & 3/8 & 0 & 1/8 \\ 0 & 0 & 0 & 0 & 1/8 & 0 \\ -\gamma/2 & 0 & 0 & 1/8 & 0 & 3/8 \end{pmatrix}. \quad (15) $$
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Left-hand panel: number density of minima, saddle points, and peaks compared to the approximation of rare events in units of $R_*^{-2}$. Right-hand panel: relative error on the number density of peaks of height $v$ when using the rare event approximation instead of the exact result is shown. Different colours correspond to different spectral indices as labelled.</figcaption></figure>

For instance, for a 2D power-law power spectrum with spectral index $n_s$ smoothed with a Gaussian filter ($r$ is now the separation in units of the Gaussian smoothing length)

$$ \langle xz \rangle = {}_1F_1\left(\frac{n_s}{2} + 1; 1; -\frac{r^2}{4}\right) \equiv \xi(r), \quad (16) $$

$$ \langle xy \rangle = \xi(|\mathbf{r}_x - \mathbf{r}_y| = 2r \sin(\psi/2)), \quad (17) $$

$$ \langle x\nabla z \rangle = \frac{\sqrt{n_s + 2}}{2\sqrt{2}} {}_1F_1\left(\frac{n_s}{2} + 2; 2; -\frac{r^2}{4}\right) \mathbf{r}, \quad (18) $$

$$ \langle xz_{11} \rangle = -\frac{\gamma}{2} \left[ 2 \cos^2(\theta + \psi) {}_1F_1\left(\frac{n_s}{2} + 2; 1; -\frac{r^2}{4}\right) - \cos(2(\theta + \psi)) {}_1F_1\left(\frac{n_s}{2} + 2; 2; -\frac{r^2}{4}\right) \right], \quad (19) $$

$$ \langle xz_{12} \rangle = -\frac{r^2 \gamma(n_s + 4)}{32} \sin(2(\theta + \psi)) {}_1F_1\left(\frac{n_s}{2} + 3; 3; -\frac{r^2}{4}\right), \quad (20) $$

$$ \langle xz_{22} \rangle = -\frac{\gamma}{2} \left[ 2 \sin^2(\theta + \psi) {}_1F_1\left(\frac{n_s}{2} + 2; 1; -\frac{r^2}{4}\right) + \cos(2(\theta + \psi)) {}_1F_1\left(\frac{n_s}{2} + 2; 2; -\frac{r^2}{4}\right) \right]. \quad (21) $$

Here ${}_1F_1(a; b; z)$ is the confluent hypergeometric function, $\xi$ is the two-point correlation function of the field, and the spectral parameter reads $\gamma = \sqrt{(n_s + 2)/(n_s + 4)}$. The correlation matrix $C_{YZ}$ is obviously the same as $C_{XZ}$ once $\psi$ has been set to zero.

## 2.3. Central peak condition

Equation (11) is sufficient to compute the expectation of any quantity involving the fields and its derivatives up to second order in three different locations. This is the case if one wants to implement a peak condition at the $\mathbf{r}_z$ location. Indeed, following Longuet-Higgins (1957), Adler (1981), Bardeen et al. (1986), this peak constraint reads $|\det z_{ij}|\delta_D(z_i)\Theta_H(-\lambda_i)$, where $\delta_D(z_i) \equiv \delta_D(z_1)\delta_D(z_2)$ is a product of Dirac delta functions, which imposes the gradient to be zero, $\Theta_H(-\lambda_i) \equiv \Theta_H(-\lambda_1)\Theta_H(-\lambda_2)$, a Heaviside function forcing the curvatures (equivalently the eigenvalues of the Hessian matrix $\lambda_i$) to be negative. The factor $|\det z_{ij}| = |z_{11}z_{22} - z_{12}^2|$ encodes the volume associated with each peak, in other words the Jacobian, which allows us to go from a smoothed field distribution to the discrete distribution of peaks.

The rareness of the peak $v$ can also be imposed by adding a factor $\delta_D(z - v)$. We therefore denote $n_{pk}(\mathbf{Z})$ the localised density of peaks

$$ n_{pk}(\mathbf{Z}) = \frac{1}{R_*^2} |\det z_{ij}| \delta_D(z_i) \Theta_H(-\lambda_i) \delta_D(z - v). \quad (22) $$

The most difficult part in the peak constraint is often to impose the sign of the curvatures and the positivity of the Jacobian, which can prevent us from getting analytical results as is the case for 3D differential peak counts (Gay et al. 2012), or peak-peak correlation functions (as described in Baldauf et al. 2016, in one dimension and Regos & Szalay 1995, in three dimensions), which can only be solved numerically. A standard approximation to keep results analytical is to drop this sign constraint and remove the absolute values of the determinant factor for high contrasts as one expects sufficiently rare critical points to be essentially peaks. If this approximation is very accurate for one-point statistics, it may not be the case for ($N > 1$) point statistics. For instance, peak-peak correlation functions on small scales are not very well reproduced by this approximation even for large contrasts because the contribution from the other critical points actually dominates at small distance; there is at least one saddle point between two peaks. However, in the context of this work, we impose the peak constraint in one location only and therefore the rare peak approximation is expected to be accurate for $v \gtrsim 2$. As an illustration, Fig. 2 shows the Gaussian mean number density of minima, saddle points, and peaks (Longuet-Higgins 1957; Adler 1981, and later generalised to weakly non-Gaussian fields by Pogosyan et al. 2011)

$$ \begin{aligned} \bar{n}_{\text{pk/min}}(\nu) = & \frac{\gamma^2 \sqrt{2} \exp\left(-\frac{\nu^2}{2}\right)}{16\pi^{3/2} R_*^2} (\nu^2 - 1) \left(1 \pm \text{erf}\left(\frac{\gamma\nu}{\sqrt{2(1-\gamma^2)}}\right)\right) \\ & + \frac{\sqrt{2} \exp\left(-\frac{3\nu^2}{6-4\gamma^2}\right)}{16\sqrt{3-2\gamma^2}\pi^{3/2} R_*^2} \\ & \times \left(1 \pm \text{erf}\left(\frac{\gamma\nu}{\sqrt{2(1-\gamma^2)(3-2\gamma^2)}}\right)\right) \\ & \pm \frac{\sqrt{1-\gamma^2}}{8\pi^2} \gamma\nu \exp\left(-\frac{\nu^2}{2-2\gamma^2}\right), \end{aligned} $$

$$ \bar{n}_{\text{sad}}(\nu) = \frac{\sqrt{2} \exp\left(-\frac{3\nu^2}{6-4\gamma^2}\right)}{8\sqrt{3-2\gamma^2}\pi^{3/2} R_*^2}, $$
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Top panels: expected two-point correlation function $\langle \kappa\kappa'|pk \rangle$ in units of $\sigma_0^2$ for a 2D power spectrum with spectral index $n_s = 0$ and a central peak of height $v_\star = \sqrt{7/3}\gamma$ and eigenvalues $\lambda_{1\star} = (-\sqrt{7/3} + \sqrt{1/3})/2$ and $\lambda_{2\star} = (-\sqrt{7/3} - \sqrt{1/3})/2$ in $\mathbf{r}_z$. Different values of $\theta$, the angle between the major axis of the ellipse (i.e. smaller curvature) and the first point on the annulus between 0 and $\pi/2$ are shown from left to right. $\psi$ is the angle between $\mathbf{r}_x$ and $\mathbf{r}_y$ and $r$ is the separation to the central peak (in units of the smoothing length) so that the shown value corresponds to the correlation function between this point, $\mathbf{r}_y$ and that on the positive x-axis at the same radius, $\mathbf{r}_x$. The values we chose here correspond to the most likely height and curvatures of a peak (and do not depend on the spectral index). Bottom panels: same as top panels for $v_r = 3$. The corresponding most likely curvatures of the peak are $\lambda_{1r} = -0.94$ and $\lambda_{2r} = -1.6$.</figcaption></figure>

and compares the latter to the high-$v$ approximation (related to the genus) $\chi(v) = \langle \det z_{ij} \delta_D(z_i) \delta_D(z_j) \rangle / R_\star^2$ which can be easily computed

$$ \chi(v) = \frac{\gamma^2}{4\sqrt{2\pi^{3/2}}R_{\star}^2} \exp\left(-\frac{v^2}{2}\right) (v^2 - 1). \quad (23) $$

The relative error between the number density of peaks and its high-$v$ approximation is shown on the right-hand panel of Fig. 2.

## 2.4. Density correlations on the circle surrounding a central peak with given geometry

The expected product of projected density $\kappa$ in two locations of space $\mathbf{r}_x$ and $\mathbf{r}_y$ such that $\mathbf{r}_x - \mathbf{r}_z = r(\cos\theta, \sin\theta)$ and $\mathbf{r}_y - \mathbf{r}_z = r(\cos(\theta+\psi), \sin(\theta+\psi))$ and given a peak in $\mathbf{r}_z$ of height $v$ and curvatures $0 > \lambda_1 > \lambda_2$ along the first and second coordinates can be analytically computed. For instance for a power-law power spectrum with $n_s = 0$ (and $\gamma = 1/\sqrt{2}$), we get

$$ \frac{\langle \kappa(r, \theta) \kappa(r, \theta + \psi) | pk \rangle}{\sigma_0^2} = \xi(|\mathbf{r}_x - \mathbf{r}_y|) + \exp\left(-\frac{r^2}{2}\right) [l_0 + l_2 r^2 + l_4 r^4], \quad (24) $$

where

$$ l_0 = (v^2 - 1), \quad (25) $$

$$ l_2 = [v^2 + \sqrt{2\nu}I_1(1 - 2e \cos \psi \cos(2\theta + \psi)) - \cos \psi]/2, \quad (26) $$

$$ l_4 = [v^2 - 2 \cos^2 \psi + 2 \sqrt{2\nu}I_1(1 - 2e \cos \psi \cos(2\theta + \psi)) + 2I_1^2(1 - 2e \cos(2\theta))(1 - 2e \cos(2\theta + 2\psi))] / 16, \quad (27) $$

$l_1 = \lambda_1 + \lambda_2$ is the trace of the density Hessian at the location of the peak, $e = (\lambda_2 - \lambda_1)/(2I_1)$ is the ellipticity of the peak and the unconstrained correlation function is

$$ \xi(|\mathbf{r}_x - \mathbf{r}_y|) = \exp\left(-\frac{r^2}{2}(1 - \cos \psi)\right). \quad (28) $$

To start with, Fig. 4 shows the zero lag contribution to the annulus correlation function (i.e. when the two points are at the same location on the annulus, $\psi = 0$) in the frame of the central peak. As expected the amplitude of fluctuations around the peak have an ellipsoidal shape that is more elongated along the smallest curvature $\lambda_{1\star}$. This zero-lag annulus correlation is dominated by the square of the mean density profile at small separations and by the fluctuations at larger separations.

Figure 3 then shows the full constrained correlation function on the annulus. We use different orientations of the pair $(\mathbf{r}_x, \mathbf{r}_y)$
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Zero lag annulus correlation function $\langle \kappa(r, \theta)^2 | \text{pk}(v) \rangle$ in units of $\sigma_0^2$ for a central peak with height $v_* = \sqrt{7/3}\gamma$ and eigenvalues $\lambda_{1*} = (-\sqrt{7/3} + \sqrt{1/3})/2$ and $\lambda_{2*} = (-\sqrt{7/3} - \sqrt{1/3})/2$ in $r_s$ for a Gaussian random field with power spectrum $P(k) \propto k^0$ smoothed with a Gaussian filter (similar to Fig. 3 when the two points on the annulus are the same). The separation $r$ is given in units of the smoothing length.</figcaption></figure>

with regard to the axis of smaller curvature (corresponding to $\lambda_1$) of the central peak. The orientation of $r_x$ is described by the angle $\theta$, which is taken to be 0, $\pi/4$ and $\pi/2$ from the left-hand to the right-hand panel. On each plot, the angle between $r_x$ and $r_y$, namely $\psi$, vary between 0 and $2\pi$ and the separation to the central peak is described by the value $r$. We show the result for two different peak heights: the most likely value $v_* = \gamma\sqrt{7/3}$ (top panels) and a rarer case $v_r = 3$ (bottom panels), which is more relevant to our study. In each case respectively, we fix the peak curvatures to their most likely values $\lambda_{1*} = (-\sqrt{7/3} + \sqrt{1/3})/2$, $\lambda_{2*} = (-\sqrt{7/3} - \sqrt{1/3})/2$ and $\lambda_{1r} = -0.94$, $\lambda_{2r} = -1.6$, respectively (see Appendix A for a description of the most likely geometry of a peak). As expected, the product of density is larger when the separation vectors are close one to the other and aligned with the major axis of the peak. For the case of a rare peak (bottom panels), the prominence of the peak is obviously larger (increased magnitude and spatial extend of the peak). Conversely, the common peak, $v_*$ occupies a smaller volume and is surrounded by two closer voids and peaks. In what follows, we do not fix the shape of the peak and therefore we marginalise over $\lambda_1$ and $\lambda_2$.

## 2.5. Density correlations around a peak of specified height $v$

If one wants to marginalise over the shape of the peak, which means integrating over the eigenvalues $\lambda_1$ and $\lambda_2$ in the range $\lambda_2 < \lambda_1 < 0$, then the expected product of projected densities on the annulus (the annulus two-point correlation function) becomes

$$ \frac{\langle \kappa(r, \theta) \kappa(r, \theta + \psi) | \text{pk}(v) \rangle}{\sigma_0^2} = \frac{\langle xy \det(z_{ij}) \delta_D(z - \nu) \delta_D(z_i) \Theta_H(-\lambda_i) \rangle}{\langle \det(z_{ij}) \delta_D(z_i) \Theta_H(-\lambda_i) \rangle}, $$

where we marginalise over all variables except $v$, which is fixed. Unfortunately, this expression cannot be analytically computed. For sufficiently rare peaks (high $v$), we drop the constraint on the

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Annulus correlation function $\langle \kappa(r, \theta) \kappa(r, \theta + \psi) | \text{pk}(v) \rangle$ in units of $\sigma_0^2$ for a central peak with height $v = 3$ in a Gaussian random field with power spectrum $P(k) \propto k^0$ smoothed with a Gaussian filter. The separation $r$ is given in units of the smoothing length. The angular anisotropy of the annulus correlation function is quantified using a multipolar decomposition in Sect. 2.6.</figcaption></figure>

sign of the eigenvalues (high critical points are most of the time peaks) and an explicit expression for $\langle \kappa(r, \theta) \kappa(r, \theta + \psi) | \text{pk}(v) \rangle$ can be obtained

$$ \frac{\langle \kappa(r, \theta) \kappa(r, \theta + \psi) | \text{pk}(v) \rangle}{\sigma_0^2} = \xi(|\mathbf{r}_x - \mathbf{r}_y|) + \frac{2 f_{21}^2}{\nu^2 - 1} + 4 f_{11} f_{21} + \frac{\nu^4 - 6\nu^2 + 3}{\nu^2 - 1} f_{11}^2 - \frac{n_s + 2}{4} r^2 f_{22}^2 \cos \psi - \frac{2 \cos(2\psi)}{\nu^2 - 1} (f_{21} - f_{11})^2, \quad (29) $$

where $f_{ij}$ and the unconstrained correlation function $\xi(|\mathbf{r}_x - \mathbf{r}_y|)$ are functions of the following Kummer confluent hypergeometric functions

$$ f_{ij} = {}_1F_1\left(\frac{n_s}{2} + i; j; -\frac{r^2}{4}\right), \quad (30) $$

$$ \xi(|\mathbf{r}_x - \mathbf{r}_y|) = {}_1F_1\left(\frac{n_s}{2} + i; j; -\frac{r^2}{2}(1 - \cos \psi)\right). \quad (31) $$

As an illustration, for a power spectrum $P(k) \propto k^0$, it becomes

$$ \frac{\langle \kappa(r, \theta) \kappa(r, \theta + \psi) | \text{pk}(v) \rangle}{\sigma_0^2} = \xi(|\mathbf{r}_x - \mathbf{r}_y|) + \exp\left(-\frac{r^2}{2}\right) \times \frac{8(\nu^2 - 1)^2 - 8\nu^2 r^2 + r^4 - 4(\nu^2 - 1)r^2 \cos \psi - r^4 \cos 2\psi}{8(\nu^2 - 1)}, \quad (32) $$

where $\xi(|\mathbf{r}_x - \mathbf{r}_y|) = \exp(-\frac{r^2}{2}(1 - \cos \psi))$ is the unconstrained correlation function on the annulus. The apparent singularity at $v = \pm 1$ is due to our high $v$ approximation, which breaks down in this regime as many $v = 1$ critical points are not peaks but saddle points. Figure 5 illustrates the behaviour of $\langle \kappa(r, \theta) \kappa(r, \theta + \psi) | \text{pk}(v) \rangle$ for a central peak with height $v = 3$. Similar to the case where the peak geometry is imposed, here the annulus correlation function is larger when the separation vectors are close one to the other and aligned with the major axis of the peak. The isocontours are close to spherical for small separations but become very anisotropic and elongated along the axis $\psi = 0$ (when the two points overlap) at larger separations.
<figure><img src="image_8.png" /><figcaption>Fig. 6. Multipoles $|Q_m|^2$ for a central peak with height $v = 2$ to 4 as labelled in a Gaussian random field with power spectrum $P(k) \propto k^0$ smoothed with a Gaussian filter and on the annulus at a distance $r = 0.1$ R (top left-hand panel), R (top right), and 2 R (bottom left). The dashed line corresponds to the random case in which we do not impose a central peak. The bottom right-hand panel shows the ratio of those multipoles to the random case for $r = R$.</figcaption></figure>

## 2.6. Multipoles around a peak of specified height $v$

Once the two-point correlation function around a peak – $\langle \kappa(r, \theta)\kappa(r, \theta + \psi)|pk(v)\rangle$ – is known, one can compute the corresponding multipolar moments that we define here as

$$ \langle |Q_m|^2 |pk \rangle (r, v) = \int_0^{2\pi} \frac{d\psi}{2\pi\sigma_0^2} \langle \kappa(r, \theta)\kappa(r, \theta + \psi)|pk \rangle e^{im\psi}. \quad (33) $$

The result is again analytical. As expected, only the first three multipoles are modified by the peak condition, while the rest are unchanged, i.e.,

$$ \langle |Q_m|^2 |pk \rangle = \langle |Q_m|^2 \rangle \quad \text{for all } m \ge 3. \quad (34) $$

For instance, for $P(k) \propto k^0$ power spectra, those multipoles read

$$ \langle |Q_0|^2 |pk \rangle = \langle |Q_0|^2 \rangle + \frac{r^4 - 8\nu^2 r^2 + 8(\nu^2 - 1)^2}{8(\nu^2 - 1)} \exp\left(-\frac{r^2}{2}\right), \quad (35) $$

$$ \langle |Q_1|^2 |pk \rangle = \langle |Q_1|^2 \rangle - \frac{1}{4} r^2 \exp\left(-\frac{r^2}{2}\right), \quad (36) $$

$$ \langle |Q_2|^2 |pk \rangle = \langle |Q_2|^2 \rangle - \frac{1}{16} \frac{r^4}{\nu^2 - 1} \exp\left(-\frac{r^2}{2}\right), \quad (37) $$

$$ \langle |Q_m|^2 |pk \rangle \stackrel{m \ge 3}{=} \langle |Q_m|^2 \rangle \equiv \exp\left(-\frac{r^2}{2}\right) I_m\left(\frac{r^2}{2}\right), \quad (38) $$

where $I_m$ are the modified Bessel functions of the first kind. In particular the correction to the monopole (resp. dipole, quadrupole) is maximal for $r = 0$ (resp. $\sqrt{2}$, 2). It can easily be checked that the condition of zero gradient only affects

the dipole, while the constraint on the peak height changes the monopole and the Hessian modifies both the monopole and quadrupole.

Figure 6 shows the amplitude of the multipoles for various peak heights and separations. There is a significant drop of power in the dipole while the change in the monopole and quadrupole is much less pronounced. The dependence on the peak height is rather small. Those predictions are checked against GRF realisations in Sect. 3.

## 2.7. Dependence on the slope of the power spectrum

In this work, we have shown results for a power-law power spectrum ($n_s = 0$) but the qualitative conclusions can be shown to be almost independent from the spectral index. To illustrate this property, we computed the multipoles for different slopes of the power spectrum from -1.5 (close to the effective spectral index of the convergence field at cluster scale) to 1 as shown in Fig. 7. The correction to the monopole and dipole are quasi-linearly suppressed when $n_s$ increases while the quadruple is constant for a wide range of slopes $n_s \lesssim 1$ and shows only a decrease at very low spectral indices. Overall, this shows that the qualitative picture described in this paper does not depend significantly on the slope of the power spectrum. Investigating the effect of its running (i.e the scale-invariance of the slope) is left for future works as no analytical results can be obtained in this case. The study of a more realistic $\Lambda$CDM power spectrum in the non-linear regime will be presented elsewhere.
<figure><img src="image_9.png" /><figcaption>Fig. 7. Same as Fig. 6 for central peaks of height $v = 3$ and separation $r = 1$ (in units of the smoothing length) as a function of the spectral index $n_s$.</figcaption></figure>

### 3. Comparison with direct measurements in GRF

Let us generate 10 maps of a $2048^2$ GRF with power spectrum $P(k) \propto k^0$. Each map is then smoothed with a Gaussian kernel on $R = 8$ pixels. A portion of such a map is shown in the right-hand panel of Fig. 1.

Peaks are then found using the code `map2ext` (Colombi et al. 2000; Pogosyan et al. 2011): for every pixel a segment of quadratic surface is fit in the tangent plane based on the field values at the pixel of origin and its neighbours. The position of the extremum of this quadratic surface, its height, and its Hessian are computed. The extremum is counted into the tally of the type determined by its Hessian (two negative eigenvalues for peaks) if its position falls within the original pixel. Several additional checks are performed to preclude registering extrema in the neighbouring pixels and to minimise missing extrema due to jumps in the fit parameters as region shifts to the next pixel. This procedure performs with better than 1% accuracy when the map is smoothed with a Gaussian filter whose full width at half maximum exceeds 6 pixels.

The field is then interpolated at 100 equally spaced points on the circle located at $r = R$ around each peak and Fourier transformed. Only the square modulus of the Fourier coefficients are stored. For comparison, a similar procedure is followed to estimate the multipolar decomposition around the same number of random points in the field.

The resulting multipolar decomposition measured in GRF is shown in Fig. 8 for various peak heights and separations. Those measurements are in very good agreement with the theoretical predictions described in Sect. 2.6. The high-$v$ approximation used to derive the prediction is therefore shown to be very accurate in the regime $v \ge 2.5$. Below this threshold, some departures – in particular in the quadrupole – are seen and would require a numerical integration of the equation with the correct peak curvature constraints.

### 4. Effect of substructures

In practice, measurements in simulations and observations of the angular distribution of the convergence field around clusters naturally involve two separate scales: the relatively large scale of the cluster and the smaller scale of the convergence field, or the dark matter density field in a N-body simulation, around it. Even if those scales are not identical, they are necessarily highly correlated and the effect described in this paper should persist. To

study the effect of substructures, we redo the analysis but introduce two different smoothing lengths, one $R_1$ for the field $z$ at the location of the peak and one $R_2$ at the location of the anulus. In this section only, we denote $R = R_2/R_1 \le 1$ as the corresponding (dimensionless) ratio. The same formalism as described above applies but all the coefficients of the covariance matrix are changed. Let us first redefine the random variables as

$$z = \frac{1}{\sigma_0} \kappa(\mathbf{r}_z), \quad z_i = \frac{1}{\sigma_1} \nabla_i \kappa(\mathbf{r}_z), \quad z_{ij} = \frac{1}{\sigma_2} \nabla_i \nabla_j \kappa(\mathbf{r}_z), \quad (39)$$

$$x = \frac{1}{\sigma_0} \kappa(\mathbf{r}_x), \quad y = \frac{1}{\sigma_0} \kappa(\mathbf{r}_y), \quad (40)$$

where the factors $\sigma_i$ are the respective variances of the field, gradient, and Laplacian smoothed on scale $R_1$. With this definition, one can easily recompute the coefficients of the covariance matrix. For instance,

$$\langle xz \rangle = \langle yz \rangle = \beta^{-n_s-2} {}_1F_1 \left( \frac{n_s}{2} + 1; 1; -\frac{r^2}{4\beta^2} \right), \quad (41)$$

$$\langle xy \rangle = R^{-\frac{n_s+2}{2}} {}_1F_1 \left( \frac{n_s}{2} + 1; 1; -\frac{r^2}{2R^2} (1 - \cos \psi) \right), \quad (42)$$

$$\langle xz_{22} \rangle = -\frac{\gamma}{2} \beta^{-n_s-4} \left[ 2 \sin^2(\theta + \psi) {}_1F_1 \left( \frac{n_s}{2} + 2; 1; -\frac{r^2}{4\beta^2} \right) + \cos(2(\theta + \psi)) {}_1F_1 \left( \frac{n_s}{2} + 2; 2; -\frac{r^2}{4\beta^2} \right) \right], \quad (43)$$

where the separation $r$ is again a dimensionless quantity (expressed in units of $R_1$) and $\beta$ is the dimensionless quadratic mean of the two smoothing lengths

$$\beta = \sqrt{(1 + R^2)/2}, \quad (44)$$

which ranges from $1/\sqrt{2}$, when $R_2$ goes to zero, to 1, when the two smoothing scales are equal $R_2 = R_1$.

An analytical solution for the mean amplitude of the multipoles of the field around a central peak can again be computed. For the same example of a power-law power spectrum $P(k) \propto k^0$, those multipoles read

$$\langle |Q_0|^2 |pk \rangle = \langle |Q_0|^2 \rangle + \frac{a_0^{pk}}{8\beta^{12}(\nu^2 - 1)} \exp\left(-\frac{r^2}{4\beta^2}\right), \quad (45)$$

$$\langle |Q_1|^2 |pk \rangle = \langle |Q_1|^2 \rangle - \frac{r^2}{4\beta^8} \exp\left(-\frac{r^2}{4\beta^2}\right), \quad (46)$$

$$\langle |Q_2|^2 |pk \rangle = \langle |Q_2|^2 \rangle - \frac{r^4}{16(\nu^2 - 1)\beta^{12}} \exp\left(-\frac{r^2}{4\beta^2}\right), \quad (47)$$

$$\langle |Q_m|^2 |pk \rangle \stackrel{m \ge 3}{=} \langle |Q_m|^2 \rangle \equiv \frac{1}{R^2} \exp\left(-\frac{r^2}{2R^2}\right) I_m\left(\frac{r^2}{2R^2}\right), \quad (48)$$

where $I_m$ are the modified Bessel functions of the first kind and

$$a_0^{pk} = r^4 - 8\beta^2 [1 + \beta^2(\nu^2 - 1)]r^2 + 8\beta^4 [2 + 4\beta^2(\nu^2 - 1) + \beta^4(\nu^4 - 6\nu^2 + 3)].$$

The limit $R = \beta = 1$ trivially reduces to the former Eqs. (35)–(38).

Those small-scale multipoles are shown in Fig. 9 for $R = 1/100$ to $R = 1$. The multiscale approach described in this section does not modify the $m > 2$ multipoles. As expected, the correction due to the peak decreases when $R$ goes to 0 as the scales decorrelate. In addition, we expect that non-linearities and corrections beyond the Hessian change the power of higher order multipoles.
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Left-hand panel: same as the bottom right-hand panel of Fig. 6 for measurements in 10 realisations of a 2048² 2D GRF smoothed with a Gaussian filter on 8 pixels. The height of the peaks are binned as labelled and the separation considered here is r = R = 8 pixels. Right-hand panel: same as left-hand panel when we vary the separation r instead of the peak height, which is set to v = 3 here. We overplotted the theoretical predictions with dashed lines that are almost indistinguishable from the measurements.</figcaption></figure>

## 5. Beyond the thin shell approximation

The effect of the radial weight function in Eq. (7) can be studied by relaxing the assumption that $r_x$ and $r_y$ are on a same infinitely thin shell around the central peak in $r_z$. Let us therefore consider the general setting for which $r_x$ is at a distance r from the central peak and $r_y$ at a distance r'. In this case, the constrained two-point correlation function reads

$$ \frac{\langle \kappa(r, \theta)\kappa(r, \theta + \psi) |pk(v)\rangle}{\sigma_0^2} = \xi(|r_x - r_y|) + \frac{2f_{21}f'_{21}}{v^2 - 1} + 2(f'_{11}f_{21} + f_{11}f'_{21}) + \frac{v^4 - 6v^2 + 3}{v^2 - 1} f_{11}f'_{11} - \frac{n_s+2}{4}r^2\cos\psi f_{22}f'_{22} - \frac{2\cos(2\psi)}{v^2 - 1}(f_{21} - f_{11})(f'_{21} - f'_{11}), \quad (49) $$

where f, f', and the unconstrained correlation function $\xi(|r_x - r_y|)$ are functions of the following Kummer confluent hypergeometric functions

$$ f_{ij} = {}_1F_1\left(\frac{n_s}{2} + i; j; -\frac{r^2}{4}\right), \quad (50) $$

$$ f'_{ij} = {}_1F_1\left(\frac{n_s}{2} + i; j; -\frac{r'^2}{4}\right), \quad (51) $$

$$ \xi(|r_x - r_y|) = {}_1F_1\left(\frac{n_s}{2} + i; j; -\frac{r^2 + r'^2 - 2rr' \cos \psi}{4}\right). \quad (52) $$

It can easily be checked that Eq. (49) trivially reduces to Eq. (29) when $r' = r$. As an illustration, for a power-law power spectrum $P(k) \propto k^0$, the correction to the unconstrained correlation function $\xi(|r_x - r_y|)$ reads

$$ \frac{\langle \kappa(r, \theta)\kappa(r', \theta + \psi) |pk\rangle}{\sigma_0^2} - \xi(|r_x - r_y|) = \exp\left(-\frac{r^2 + r'^2}{4}\right) \delta\xi, \quad (53) $$

with

$$ \delta\xi = \frac{8(v^2 - 1)^2 - 4v^2(r^2 + r'^2) + r^2r'^2(1 - \cos 2\psi) - 4(v^2 - 1)rr'\cos\psi}{8(v^2 - 1)} $$

<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Same as the bottom right-hand panel of Fig. 6 when the field is smoothed at two different scales whose ratio R goes from 1/100 to 1.</figcaption></figure>

From Eq. (49), one can now easily compute the statistics of the multipoles including the radial weight function $w_n(r)$ that appears in Eq. (7). We find again that only the amplitude of the first three multipoles are affected by the peak constraint

$$ \langle |Q_0|^2 |pk \rangle = \langle |Q_0|^2 \rangle + \int_{S_0} \frac{r^4 - 8v^2r^2 + 8(v^2 - 1)^2}{8(v^2 - 1)} \exp\left(-\frac{r^2}{2}\right), \quad (54) $$

$$ \langle |Q_1|^2 |pk \rangle = \langle |Q_1|^2 \rangle - \int_{S_1} \frac{r^2}{4} \exp\left(-\frac{r^2}{2}\right), \quad (55) $$

$$ \langle |Q_2|^2 |pk \rangle = \langle |Q_2|^2 \rangle - \int_{S_2} \frac{r^4}{16v^2 - 1} \exp\left(-\frac{r^2}{2}\right), \quad (56) $$

$$ \langle |Q_m|^2 |pk \rangle \stackrel{m \ge 3}{=} \langle |Q_m|^2 \rangle, \quad (57) $$

where $\int_{S_m}$ stands for the following 2D radial integral

$$ \int_{S_m} f(r, r') = (2\pi)^2 \sigma_0^2 \int r dr r' dr' r^n w_n(r) r'^m w_m(r'). \quad (58) $$

Figure 10 shows the resulting multipoles for a radial weight function defined following Schneider & Bartelmann (1997) as

$$ R_{\max}^{1+m} w_m(r) = \frac{1}{x^{1+m} + \alpha^{1+m}} - \frac{1}{1 + \alpha^{1+m}} + \frac{(1+m)(x-1)}{(1 + \alpha^{1+m})^2} \quad (59) $$

over the range $x = r/R_{\max} \in [\alpha, 1]$ and zero elsewhere (which was found to be optimal for an isothermal mass distribution).
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Same as the bottom right-hand panel of Fig. 6 when we apply a radial weight function defined by Eq. (59), where $R_{\max}$ equals the smoothing length and $\alpha = 0.5$, which means that the minimum radius considered is half the smoothing length.</figcaption></figure>

The qualitative picture does not change; the most affected multipole is the dipole whose power is significantly reduced by the peak constraint, the monopole and quadrupole are slightly affected in a $v$-dependent way, and all other coefficients are unaffected.

# 6. A non-linear theory of harmonic power spectra

In this section, we study the weakly non-linear evolution of the multipolar moments. We therefore no longer assume that the PDF is Gaussian $P(x, y, z) = G(x, y, z)$. Instead, we expand the PDF around a Gaussian by means of the so-called Gram-Charlier expansion (Cramér 1946; Pogosyan et al. 2009a). For simplicity, we restrict ourselves to the case in which we only impose the height of the cluster but not the rest of the peak condition (no zero gradient or constraint on the eigenvalues of the Hessian). We show that this effect dominates the high multipoles.

## 6.1. The Gram Charlier expansion

The Gaussian PDF has zero means and covariance matrix

$$C = \begin{pmatrix} 1 & b & a \\ b & 1 & a \\ a & a & 1 \end{pmatrix}, \qquad (60)$$

with $a = \xi(r)$ and $b = \xi(2r \sin(\psi/2))$.

We first diagonalise this matrix and use a new set of variables $(u, v, z)$, where

$$w = \frac{y - az}{\sqrt{1 - a^2}} \qquad (61)$$

$$u = \frac{x(1 - a^2) + y(a^2 - b) + za(b - 1)}{\sqrt{(1 - b)(1 - 2a^2 + b)(1 - a^2)}} \qquad (62)$$

so that

$$P(u, w, z) = N(u)N(w)N(z) \qquad (63)$$

with $N$ a normal distribution of zero mean and unit variance.

Following Gay et al. (2012), Codis et al. (2013), we then use a Gram-Charlier expansion of the PDF

$$\mathcal{P}(u, w, z) = G(u, w, z) \left[ 1 + \sum_{i+j+k=3}^{\infty} \frac{H_i(u)H_j(w)H_k(z)}{i! j! k!} \langle u^i w^j z^k \rangle_{\text{GC}} \right],$$

where $H_i$ represent probabilistic Hermite polynomials and the Gram-Charlier coefficients are given by

$$\langle u^i w^j z^k \rangle_{\text{GC}} = \langle H_i(u) H_j(w) H_k(z) \rangle. \qquad (64)$$

Once the joint PDF is known, we can compute the annulus two-point correlation function $\langle xy | z = v \rangle$ as

$$\langle xy | z = v \rangle = \frac{\int dx dy \mathcal{P}(u(x, y, z), w(y, z), v) xy}{\mathcal{P}(z = v)}, \qquad (65)$$

which can be rewritten

$$\langle xy | z = v \rangle = \frac{\int du dw \mathcal{P}(u, w, v) x(u, w, v) y(u, w, v)}{\mathcal{P}(z = v)}, \qquad (66)$$

where $xy$ is a polynomial of $u$ and $w$

$$xy = a^2 z^2 + uza \sqrt{\frac{(1-b)(1-2a^2+b)}{1-a^2}} + wza \frac{1+b-2a^2}{\sqrt{1-a^2}} + uw \sqrt{(1-b)(1-2azx+a^2z^2+a^2)} + (b-a^2)w^2, \qquad (67)$$

which is the sum of four terms proportional respectively to $H_0(u)H_0(w)$, $H_1(u)H_0(w)$, $H_0(u)H_1(w)$, $H_1(u)H_1(w)$ and $H_0(u)(H_2(w) + H_0(w))$ as $H_0(x) = 1$, $H_1(x) = x$ and $H_2(x) = x^2 - 1$, respectively. Using the property of orthogonality of Hermite polynomials, it is then easy to compute Eq. (66) so that eventually

$$\langle xy | z = v \rangle = b + a^2 (v^2 - 1) + \Delta_{\text{NL}}, \qquad (68)$$

where the non-linear contribution reads

$$\Delta_{\text{NL}} = \frac{\sum_{k=1}^{\infty} \frac{H_k(v)}{k!} \langle H_k(z) [2av(x-az) + xy - b - 2azx + a^2z^2 + a^2] \rangle}{1 + \sum_{k=3}^{\infty} \frac{H_k(v)}{k!} \langle z^k \rangle_{\text{GC}}},$$

which, at first order in $\sigma_0$, the amplitude of fluctuations, is given by<fnref n="2" />

$$\Delta_{\text{NL}}^{(1)} = a^2(2v - v^3) \langle z^3 \rangle + a(v^3 - 3v) \langle xz^2 \rangle + v \langle xyz \rangle.$$

In terms of multipoles, it means that for $m > 0$, we get a non-linear bias given by

$$\frac{\langle |Q_m|^2 | v \rangle}{\langle |Q_m|^2 \rangle} = 1 + v \frac{\langle xyz \rangle_m}{b_m} + O(\sigma_0^2), \qquad (69)$$

where the subscript $m$ refers to the associated multipole of order $m$. The multipoles near a high density cluster are therefore biased compared to random locations; this bias is proportional to the height $v$ with a proportionality coefficient related to the ratio between the isosceles three-point function $\langle xyz \rangle$ and the two-point correlation function of its base $b = \langle xy \rangle = \xi(2r \sin(\psi/2))$. The all-order expression is also easily obtained once it is realised that the only terms that depend on the angle $\psi$ are $b$ and the cumulants involving the product $xy$

$$\frac{\langle |Q_m|^2 | v \rangle}{\langle |Q_m|^2 \rangle} = 1 + \frac{\sum_{k=1}^{\infty} \frac{H_k(v)}{k!} \langle H_k(z)(xy-b) \rangle_m}{b_m \left( 1 + \sum_{k=3}^{\infty} \frac{H_k(v)}{k!} \langle z^k \rangle_{\text{GC}} \right)}, \qquad (70)$$

<footnote n="2">² An easy way to get this expression is to only keep Gram-Charlier coefficients $\langle u^i w^j z^k \rangle_{\text{GC}}$ for which $i + j + k = 3$, which were shown to be equivalent to cumulants and correspond exactly to the first-order correction that is proportional to $\sigma_0$ (Gay et al. 2012).</footnote>
<figure><img src="image_14.png" /><figcaption>Fig. 11. Slice through the fields $n_{step} = 0, 1, 2$, and 3 obtained by a Zeldovich displacement of an initial GRF. We measure the multipolar moments around all points of height $v > 2$ in those maps.</figcaption></figure>

The monopole is also easy to compute

$$ \frac{\langle |Q_0|^2 |v| \rangle}{\langle |Q_0|^2 \rangle} = 1 + \frac{a^2(v^2-1)}{b_0} + \delta q_0^{\text{NL}}, \qquad (71) $$

where the first order non-linear correction reads

$$ \delta q_0^{\text{NL}} = +v \frac{a^2(2-v^2)\langle z^3 \rangle + a(v^2-3)\langle xz^2 \rangle + \langle xyz \rangle_0}{b_0} + O(\sigma_0^2). $$

## 6.2. Comparison with simulations

To test the $m > 0$ prediction, we have generated various GRF and displaced the density field following a Zeldovich displacement with different time steps denoted $n_{step} = 0, 1, 2$, and 3 (from Gaussian to more evolved fields). In practice, we compute the displacement field as the gradient of the gravitational potential by FFT and we multiply by a constant times $n_{step}$. We then move the mass in each pixel according to this displacement and distribute it to the eight closest pixels. Those fields are illustrated in Fig. 11. We measure the multipoles around field point of height $v > 2$ together with the mean height of those peaks (resp. $\bar{v} = 2.36, 2.52, 2.63, 2.86$) and the multipolar decomposition of the bispectrum $\langle xyz \rangle_m$. The result is shown in Fig. 12 and shows a fair agreement of the $m > 0$ multipoles with the prediction given in Eq. (69).

<figure><img src="image_15.png" /><figcaption>Fig. 12. Ratio between the multipolar moments at a distance $r = 8$ pixels from a field point of height $v > 2$ compared to random locations, for a GRF displaced following the Zeldovich approximation for different time steps between 0 and 20 as labelled and smoothed over 3 pixels. The solid line is the measurements and the dashed line is the prediction for $m > 0$ given by Eq. (69).</figcaption></figure>

## 7. Conclusions

We have computed the statistics of the multipolar moments around a peak for a generic 2D Gaussian field as a proxy for
the azimuthal distribution of matter around clusters seen by weak gravitational lensing experiments. For rare enough peaks ($v \gtrsim 2.5$), all results are completely analytical. It is shown that only the monopole, dipole, and quadrupole are affected by the central peak while higher order multipoles are essentially left unchanged by the peak constraint. Overall, the dominant effect we find is a significant drop in the dipole coefficient as expected from the zero gradient condition. Substructures in the Gaussian field and the addition of a radial weighting function do not change this qualitative picture.

This feature in the dipole can also be detected in numerical simulations of structure formation as will be shown in a forthcoming paper (Gouin et al. 2017). We anticipate that higher order corrections will also emerge from the non-linear evolution of the density field in the vicinity of peaks beyond the Gaussian picture described here, and also from a possible departure from the peak model itself which, as we showed in this paper, boils down to modifying the power in the monopole, dipole, and quadrupole only. As an illustration, we computed the non-linear bias of the multipolar moments due to the height of the cluster. This bias is proportional to the height $v$ and to the variance of the field $\sigma$ by means of the rescaled bispectrum. This approach based on the statistics of multipolar moments in the convergence field around clusters will soon be applied to data (Gavazzi et al., in prep.).

Extensions of this analytical work in the future might include i) an investigation of the accuracy of the large $v$ approximation and a precise numerical integration in the regime of intermediate contrasts where this approximation breaks down; and ii) a study of the effect of the scale-dependence of the power spectrum.

**Acknowledgements.** This work is partially supported by the grants ARN-13-BS05-0005 of the French Agence Nationale de la Recherche. This work has made use of the Horizon cluster on which the GRF maps were generated, hosted by the Institut d'Astrophysique de Paris. We warmly thank D. Pogosyan for insightful discussions, his careful reading of the manuscript, and for providing us with his code map2ext to detect extrema in 2D maps. We also thank S. Rouberol for running the Horizon cluster for us and D. Munro for freely distributing his Yorick programming language and opengl interface (available at yorick.sourceforge.net).

## References

* Adler, R. J. 1981, The Geometry of Random Fields (Chichester: Wiley)

* Aragón-Calvo, M. A., van de Weygaert, R., Jones, B. J. T., & van der Hulst, J. M. 2007, ApJ, 655, L5

* Aragon-Calvo, M. A., & Yang, L. F. 2014, MNRAS, 440, L46

* Aubert, D., & Pichon, C. 2007, MNRAS, 374, 877

* Aubert, D., Pichon, C., & Colombi, S. 2004, MNRAS, 352, 376

* Bailin, J., & Steinmetz, M. 2005, ApJ, 627, 647

* Baldauf, T., Codis, S., Desjacques, V., & Pichon, C. 2016, MNRAS, 456, 3985

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

* Bartelmann, M., & Schneider, P. 2001, Phys. Rep., 340, 291

* Blandford, R. D., Saust, A. B., Brainerd, T. G., & Villumsen, J. V. 1991, MNRAS, 251, 600

* Bond, J. R., & Myers, S. T. 1996, ApJS, 103, 1

* Bond, J. R., Kofman, L., & Pogosyan, D. 1996, Nature, 380, 603

* Brunino, R., Trujillo, I., Pearce, F. R., & Thomas, P. A. 2007, MNRAS, 375, 184

* Codis, S., Pichon, C., Devriendt, J., et al. 2012, MNRAS, 427, 3320

* Codis, S., Pichon, C., Pogosyan, D., Bernardeau, F., & Matsubara, T. 2013, MNRAS, 435, 531

* Codis, S., Pichon, C., & Pogosyan, D. 2015, MNRAS, 452, 3369

* Colombi, S., Pogosyan, D., & Souradeep, T. 2000, Phys. Rev. Lett., 85, 5515

* Cramér, H. 1946, Mathematical Methods of Statistics (Princeton Univ. Press)

* Davé, R., Cen, R., Ostriker, J. P., et al. 2001, ApJ, 552, 473

* Dietrich, J. P., Schneider, P., Clowe, D., Romano-Díaz, E., & Kerp, J. 2005, A&A, 440, 453

* Doroshkevich, A. G. 1970, Astrophysics, 6, 320

* Dubois, Y., Pichon, C., Welker, C., et al. 2014, MNRAS, 444, 1453

* Fukugita, M., Hogan, C. J., & Peebles, P. J. E. 1998, ApJ, 503, 518

* Gay, C., Pichon, C., Le Borgne, D., et al. 2010, MNRAS, 404, 1801

* Gay, C., Pichon, C., & Pogosyan, D. 2012, Phys. Rev. D, 85, 023011

* Gouin, C., Gavazzi, R., Codis, S., et al. 2017, A&A, 605, A27

* Guzzo, L., Strauss, M. A., Fisher, K. B., Giovanelli, R., & Haynes, M. P. 1997, ApJ, 489, 37

* Hahn, O., Porciani, C., Carollo, C. M., & Dekel, A. 2007, MNRAS, 375, 489

* Hahn, O., Teyssier, R., & Carollo, C. M. 2010, MNRAS, 405, 274

* Kaiser, N. 1984, ApJ, 284, L9

* Kaiser, N. 1992, ApJ, 388, 272

* Klypin, A., & Shandarin, S. F. 1993, ApJ, 413, 48

* Kovač, K., Lilly, S. J., Knobel, C., et al. 2014, MNRAS, 438, 717

* Longuet-Higgins, M. S. 1957, Phil. Trans. R. Soc. London Ser. A, 249, 321

* Ludlow, A. D., & Porciani, C. 2011, MNRAS, 413, 1961

* Martinet, N., Clowe, D., Durret, F., et al. 2016, A&A, 590, A69

* Mead, J. M. G., King, L. J., & McCarthy, I. G. 2010, MNRAS, 401, 2257

* Metuki, O., Libeskind, N. I., Hoffman, Y., Crain, R. A., & Theuns, T. 2015, MNRAS, 446, 1458

* Miralda-Escudé, J. 1991, ApJ, 380, 1

* Navarro, J. F., Abadi, M. G., & Steinmetz, M. 2004, ApJ, 613, L41

* Oemler, Jr., A. 1974, ApJ, 194, 1

* Paranjape, A., & Sheth, R. K. 2012, MNRAS, 426, 2789

* Paz, D. J., Stasyszyn, F., & Padilla, N. D. 2008, MNRAS, 389, 1127

* Persic, M., & Salucci, P. 1992, MNRAS, 258, 14

* Pogosyan, D., Gay, C., & Pichon, C. 2009a, Phys. Rev. D, 80, 081301

* Pogosyan, D., Pichon, C., Gay, C., et al. 2009b, MNRAS, 396, 635

* Pogosyan, D., Pichon, C., & Gay, C. 2011, Phys. Rev. D, 84, 083510

* Regos, E., & Szalay, A. S. 1995, MNRAS, 272, 447

* Schneider, P., & Bartelmann, M. 1997, MNRAS, 286, 696

* Shull, J. M., Smith, B. D., & Danforth, C. W. 2012, ApJ, 759, 23

* Simon, P. 2007, A&A, 473, 711

* Sousbie, T., Pichon, C., Colombi, S., & Pogosyan, D. 2008, MNRAS, 383, 1655

* Tempel, E., & Libeskind, N. I. 2013, ApJ, 775, L42
## Appendix A: Typical peak geometry

The typical geometry of a Gaussian peak in two dimensions can easily be computed. Starting from the Gaussian joint PDF of the field value $v$ and local curvatures $\lambda_1 > \lambda_2$ (Doroshkevich 1970; Pogosyan et al. 2009b)

$$ \mathcal{P}(v, \lambda_1, \lambda_2) = \frac{2\sqrt{J_2}}{\pi\sqrt{1-\gamma^2}} \exp\left(-\frac{1}{2}\left(\frac{v+\gamma I_1}{\sqrt{1-\gamma^2}}\right)^2 - \frac{1}{2}I_1^2 - J_2\right), \quad (A.1) $$

where $I_1 = \lambda_1 + \lambda_2$ and $J_2 = (\lambda_1 - \lambda_2)^2$, one can show that the PDF for a peak to have height $v$ and geometry $0 > \lambda_1 > \lambda_2$ reads (Bardeen et al. 1986; Codis et al. 2015)

$$ \mathcal{P}(v, \lambda_1, \lambda_2|pk) = \frac{8\sqrt{3}(\lambda_1 - \lambda_2)\lambda_1\lambda_2}{\pi\sqrt{1-\gamma^2}} \times \exp\left(-\frac{1}{2}\left(\frac{v+\gamma(\lambda_1+\lambda_2)}{\sqrt{1-\gamma^2}}\right)^2 - \frac{1}{2}(\lambda_1+\lambda_2)^2 - (\lambda_1-\lambda_2)^2\right). \quad (A.2) $$

It has to be emphasised that here we do impose exactly the peak constraint given by Eq. (22). The most likely value of the peak height and curvatures is therefore given by $v_{\star} = \sqrt{7/3}\gamma$, $\lambda_{1,2\star} = (-\sqrt{7/3} \pm \sqrt{1/3})/2$, which corresponds to an ellipticity $e_{\star} = 1/(2\sqrt{7})$.

If $v$ is fixed, for example to a rare value $v_r = 3$, the maximum of the PDF given by Eq. (A.2) is changed to $\lambda_{1r} = -0.94$ and $\lambda_{2r} = -1.6$ so that the ellipticity of the peak is given by $e_r = 0.13$, independently from the spectral parameter $\gamma$. The evolution of the most likely peak curvatures as a function of height is shown in Fig. A.1. In particular, this illustrates the well-known result that high peaks are increasingly spherical.

<figure><img src="image_16.png" /><figcaption>Fig. A.1. Most likely peak curvatures as a function of the peak height $v$ (top panel) and the corresponding most likely peak ellipticity (bottom panel). These results do not depend on the spectral index.</figcaption></figure>